999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

深海水下采油樹流場分析及流阻計算

2016-01-12 15:08:49張鵬舉吳曉冬李博吳恭興
計算機輔助工程 2015年6期
關鍵詞:有限元

張鵬舉++吳曉冬++李博++吳恭興

摘要: 建立500 m水深的采油樹的流場簡化模型,對采油樹流場在給定的2種工況下進行分析,評定流體流動和分布特點并計算流阻.分析結果表明:采油樹本體在2種工況下的壓力和速度變化較均勻,只是在變直徑處或直角彎管處會出現湍流現象.

關鍵詞: 采油樹; 流場; 流阻; 湍流; 有限元; FLUENT

中圖分類號: TE53;TB115.1文獻標志碼: B

0引言

21世紀是海洋的世紀,隨著海洋的開發與利用,海洋石油成為熱點話題.我國南海海域寬達2×106 km2,是世界上四大海洋油氣聚集中心之一.我國海洋石油儲量約為2.30×1010~3.00×1010 t,天然氣儲量3.38×1014 m3.[12]目前,國際海洋石油工程界普遍認可的深水定義是水深300 m.[3]我國水深300 m以上的海域達1.537×106 km2[4],但大部分尚未被勘探,因此,深海石油的勘探將成為我國未來海洋資源開發的主要方向.

作為海洋油氣開發不可或缺的設備,水下采油樹在水下生產系統中起著至關重要的作用.[58]但是,由于我國深水石油的開采起步晚,技術不全面,深水采油樹的開發受到限制.長期以來,水下采油樹的關鍵技術被美國FMC和Cameron以及挪威Aker Kvaerner Subsea 三大廠商所壟斷,這3家企業占有世界采油樹市場90%以上的份額.由于水下采油樹的結構復雜,對材料性能和密封技術的要求很高,控制系統和閥等單元部件容易出現問題.目前,國內海洋石油的裝備主要靠進口,因此對水下采油樹相關技術進行國產化研究,對突破國外技術封鎖有重大意義.[911]

1CAE分析簡介

1.1分析目的和要求

水下采油樹可分為立式和臥式2種,兩者的主要區別在于油管懸掛器和閘閥閥組的安裝位置不同.[12]本文研究的水下采油樹為臥式采油樹,其特點是油管掛在采油樹本體內,采油樹的閥組位于油管掛的側面,生產油管和油管掛的安裝要后于采油樹自身.

采用通用CFD軟件FLUENT 6.3對采油樹流場進行CAE分析,評定流體流動和分布特點,并計算流阻.CAE分析和結果要滿足合適的API 6A的溫度要求:上游18~121 ℃,下游46~121 ℃.

1.2采油樹工況

采油樹以節流閥作為整棵樹的上下游分界,以流進節流閥的一端為上游,流出節流閥的一端為下游.上、下游油管流道由130.175 0 mm(518英寸)變為103.187 5 mm(4116英寸).

在采油樹實際工作過程中,流道內一般為氣液混合相(CO2含量較高,基本沒有H2S),其設計使用壽命為20 a,設計水深為500 m.在流體特別是氣體流過節流閥時,由于壓差的存在,將產生絕熱節現象,即焦耳湯普生效應,節流閥及節流閥下游的設備會受此低溫效應的影響,因此上游主要設備(指與流體接觸的設備)額定設計溫度為18~121 ℃,下游設備額定設計溫度為46~121 ℃.

采油樹的額定壓力為69.0 MPa,測試壓力為103.5 MPa.實際上,節流閥的節流作用及流道孔徑由大變小導致下游壓力變小,也就是說,在上游負荷達到額定工作壓力的時候,下游設備并沒有達到額定工作壓力.

1.3分析內容和假設

根據給定的工作環境,分析各組件關鍵部位在各種工況下流場的變化規律.

以當前的計算機計算速度和內存,直接求解NS方程非常困難,故采用雷諾平均方程.

1.4計算方法

可壓縮流體流動和傳熱的控制方程可用雷諾平均NS方程表示,把流動變量放入連續性方程和動量方程,并且取一段時間的平均值,即

ρt+xi(ρui)=0 (1)

t(ρui)+xj(ρuiuj)=-pxi+

xjμuixj+ujxi-23δijulxl+xj(-ρu′iu′j)(2)

此方程適用于變密度的氣液混合流動,其中-ρu′iu′j稱為雷諾壓力,可利用Boussinesq假設把雷諾壓力與平均速度梯度連續起來,即-ρu′iu′j=μiuiuj+ujxi-23ρk+μtuixiδij(3)Boussinesq假設用在SA模型,kε模型和kω模型中,采用標準的kε湍流方程描述流體在管道內的流動狀態,該方程主要基于湍流動能k和擴散率ε,即

t(ρk)+xi(ρkixi)=

xjμ+μiσkkxj+Gk+Gb-ρε-YM+Sk(4)

t(ρε)+xi(ρεμi)=xjμ+μiσεεxj+

C1εεk(Gk+C3εGk)-C2ερε2k+Sε (5)

式中:ρ為流體密度,kg/m3;μ為流體動力黏度;C1ε,C2ε和C3ε為經驗常數,取值分別為1.44,1.92和0.99;σk和σε分別為k和ε的普朗特數,取值分別為1.0和1.3.[13]

此模型包含低雷諾數影響和可壓縮影響,適用于混合邊界層計算和受壁面限制的流動計算.

1.5采油樹整體分析

采油樹主要組成部分見圖1和2.H4連接器用于采油樹本體和井口的連接.油管掛坐于采油樹本體,并與采油樹本體的生產主閥和環空主閥相通.油管掛垂直中心孔以中心堵頭密封.油管掛往上為內置樹帽(仍位于采油樹本體內),再往上是垃圾帽,防止落入泥沙等.采油樹本體內的環空主閥連接環空翼閥閥組,環空翼閥閥組連接環空跨接管.環空跨接管和采油樹本體的生產主閥都接入生產翼閥閥組.通過連接管1,生產翼閥閥組和雙孔連接器流入口連接,流入口經過變徑〔130.175 0 mm(518英寸)變為103.187 5 mm(4116英寸)〕進入RPM模塊內的節流閥,然后又流回雙孔連接器的出口,經過生產隔離閥,通過連接管2將直徑101.6 mm(4英寸)管道連接器連接起來.endprint

圖 1采油樹總裝圖

Fig.1Christmas tree assembly diagram

圖 2采油樹總裝圖剖面

Fig.2Profile of christmas tree of assembly diagram

在采油樹內,流體的流通路線見圖3,通過H4連接器將采油樹本體和井口鎖緊固定,通過井下的油管將流體引入油管掛,再通過采油樹上的各類閥、閥組、連接管和連接器最終通過直徑101.6 mm(4英寸)管道連接器連接到外接管線.

圖 3流體流通路線

Fig.3Fluid circulation route

2CAE流場分析

2.1流體域模型

采油樹流體域模型見圖4.抽取流道壁并進行適當簡化,如忽略螺紋、小孔、小插拴和倒角等.圖 4采油樹流體域模型

Fig.4Fluid domain model of christmas tree

2.2網格劃分

網格劃分見圖5,共劃分2 381 521個四面體網格.由Gambit網格質量檢查功能知:最差網格歪斜度為0.784 991,小于標準值0.97,所以所劃分的網格質量較好,所有網格在整個流體計算域內分布均勻、整齊,所以,用此網格進行計算能夠滿足工程精度要求,具體的網格質量檢查見圖5c.

a)采油樹網格劃分

b1)節流閥b2)油管掛

b)節流閥和油管掛本體局部網格放大

c)采油樹網格質量檢查

圖 5采油樹網格劃分和質量檢查

Fig.5Christmas tree mesh and quality check

2.3邊界條件設定

計算域入口采用質量流量入口,出口采用壓力出口,其他位置設置為固壁條件.認為管道中流動介質為完全甲烷氣體.壓力差分格式采用標準離散差分格式,動量方程、動能方程和湍動能耗散率均采用2階迎風差分格式,通過SIMPLE算法耦合求解速度和壓力方程.

2.4計算工況

環境溫度:取海底深度為-340 m以下,中部深度約-170 m,平均環境溫度見表1.

表 1環境溫度

Tab.1Temperature of environment℃項目空氣海水海面中部海底最高36.031.423.519.5最低15.922.017.09.5

根據已知資料得到采油樹的最大日產量為6×105 m3,因此得到甲烷氣體在標準大氣壓下15.6 ℃時的體積流量Q=6×10524×3 600=6.944 m3/s (6)質量流量M=ρ標Q (7)式中:M為質量流量,kg/s;ρ標為甲烷氣體在標準大氣壓15.6 ℃時的密度,ρ標=0.677 6 kg/m3.具體計算過程見式(3).計算工況見表2.

表 2計算工況

Tab.2Calculation conditions工況質量流量

入口/(kg/s)出口

壓力/MPa入口

溫度/°C出口

溫度/°C環境

溫度/°C工況14.706 3.81032019工況24.70612.11032019

2.5設計壓力

采油樹額定工作壓力為69.0 MPa,測試壓力為103.5 MPa.

2.6分析計算

1)選擇FLUENT求解器.選擇三維單精度求解器進行分析.

2)網格的相關操作.讀入Gambit軟件生成的網格文件;檢查網格,確認最小網格體積小于零,否則需重新劃分網格;設置計算區域大小、模型和實際幾何尺寸的單位換算.

3)選擇計算模型.定義求解器,指定計算模型:選擇能量方程,選擇流態為湍流,湍流模式為kε兩方程模型,設置模型參數.

4)設置操作環境.

5)定義流體的物理性質.將流入采油樹中的流體看作純甲烷氣體,將材料的物理性質從數據庫中調出,并給定2種工況下的密度.

6)設置邊界條件.設置流體計算區域,設置質量流量入口邊界條件和壓力出口邊界條件,其他邊界條件設置為固壁,設置入口溫度和出口溫度(見表2).

7)求解方法的設置.設置求解參數,進行初始化;打開殘差圖,設置精度為0.000 01,在迭代計算時可動態顯示殘差,保存計算文件;設置迭代500次,進行迭代計算.

8)計算結果顯示.對給定的2種工況分別進行計算,完成后提取壓力云圖和速度云圖,見圖6~9,計算結果匯總見表3.

a)壓力云圖

b)節流閥和計量閥處壓力云圖局部放大

圖 6工況1采油樹壓力云圖

Fig.6Pressure contours of christmas tree under working condition 1

由圖6可知:最大壓力Pmax=24.1 MPa,最小壓力Pmin=3.2 MPa,壓差ΔP=20.9 MPa;由于節流閥和計量閥的作用,整個采油樹的壓力出現3段明顯的變化.由圖7可知:最大速度umax=739 m/s,最小速度umin=0;流體在經過節流閥過程中速度增大,最大速度出現在節流閥和計量閥處,其他位置速度都很低.由圖8可知:最大壓力Pmax=24.4 MPa,最小壓力Pmin=7.95 MPa,壓差ΔP=16.45 MPa;由于節流閥的節流作用,壓力發生急劇變化.由圖9可知:最大速度umax=494 m/s,最小速度umin=0;由于節流閥的節流作用,流體在經過節流閥時速度急劇增加.

a)速度云圖endprint

b)速度云圖局部放大

圖 7工況1采油樹速度云圖

Fig.7Velocity contours of christmas tree under working condition 1

a)壓力云圖

b)壓力云圖局部放大

圖 8工況2采油樹壓力云圖

Fig.8Pressure contours of christmas tree under working condition 2

a)速度云圖

b)速度云圖局部放大

圖 9工況2采油樹速度云圖

Fig.9Velocity contours of christmas tree under working condition 2

表 3計算結果

Tab.3Calculation results工況最大壓力/

MPa最小壓力/

MPa最大速度/

(m/s)最小速度/

(m/s)工況13.846 73.790 770.722 50工況212.138 712.093 121.594 50

3流阻計算

計算甲烷氣體在采油樹內流動過程中密度隨壓力和溫度的變化.

由表3可知,甲烷氣體在流動過程中壓力變化不大,溫度變化也不大,所以密度基本不變.只需算在2種工況下的甲烷氣體密度,即可算出壓力損失.

真實氣體狀態方程為P=R(T+273.15)V-b-aV2 (8)式中:P為壓力,MPa;V為摩爾體積,m3/kmol;R為通用氣體常數或摩爾氣體常量,各種氣體的通用氣體常數相同,R=8.314 J/(mol·°C);T為攝氏溫度;a為反應分子相互作用強度的常數,與物質的性質有關,一般由試驗獲得,對于甲烷氣體,a=0.228 589 m3/kmol;b為與分子體積有關的常數,與物質的性質有關,一般由試驗獲得,對于甲烷氣體b=0.042 71.

甲烷氣體在壓力12.1和3.8 MPa,溫度103 ℃時的摩爾體積可以根據式(8)進行計算,因此密度可以根據式(9)獲得.ρ甲烷=MmolV (9)式中:Mmol為甲烷氣體在相應壓力和溫度下的摩爾質量.計算12.1和3.8 MPa時甲烷氣體的密度ρ12.1=66.147 kg/m3,ρ3.8=20.793 kg/m3.

實際流體在運動過程中必然會產生能量損失,在伯努利方程中將單位流體的機械能損失用hw表示,稱為水頭損失.根據產生水頭損失的流動邊界條件和機理不同,水頭損失可分為沿程水頭損失和局部水頭損失.

流體在進行均勻流動時,由于黏性作用,過流斷面上的流速分布不均勻,相鄰兩流層間存在相對運動,從而使流體流層之間及流體與邊界之間存在切應力,形成流動阻力.這種在均勻流段上產生的流動阻力稱為沿程阻力,由于沿程阻力做功而引起的水頭損失稱為沿程水頭損失,其計算公式為hf=fLQm/db (10)式中:f為沿程阻力因數,與管道的種類和材質有關,一般由試驗獲得,此處根據相關資料取5 650;L為管道長度,m;Q為體積流量,m3/s,根據式(6)進行計算;m為流量指數,與管道的類型有關,一般由試驗獲得,此處根據相關資料取1.85;d為管徑,m;b為管徑指數,與管道的類型有關,一般由試驗獲得,此處根據相關資料取5.04.[14]

流體在流動的局部區域,如流體流經管道的彎頭、突擴、突縮和閘門等處,由于固體邊界的急劇改變而引起速度分布的變化,甚至使主流脫離邊界,形成漩渦區從而產生的阻力稱為局部阻力.由于局部阻力做功引起的水頭損失稱為局部水頭損失,其計算公式為 hj=ξu22g (11)式中:ξ為局部阻力因數,雷諾數對ξ的影響很小,可以忽略不計,一般情況下ξ只取決于流道的邊界形式變化,可以根據經驗選取;u為損失后的速度,m/s;g為重力加速度,m/s2.

由甲烷氣體在2種工況下的密度,可以根據質量流量和管徑計算流速, M=ρuS (12)式中:M為質量流量;ρ為甲烷氣體在相應壓力和溫度下的密度,kg/m3;u為流速,m/s;S為管道橫截面積,m2.

以直徑101.6 mm(4英寸)彎管1為例說明流阻計算過程.流體介質完全為甲烷氣體,壓力為3.8 MPa,溫度為103.15 ℃.

將彎管1從油氣進口至油氣出口劃分為小段,依次編號為34,35,36和37,見圖10.34段為直管,只需計算沿程損失.體積流量Q=0.226 3 m3/s,直徑d34=0.102 m,管長L34=0.390 m,代入公式得沿程損失hf,34=0.040 m.37段為直管,同理Q=0.226 3 m3/s,d37=0.102 m,管長L37=0.652 m,代入公式得沿程損失hf,37=0.067 m.35和36段都是彎管,局部損失為主要部分,而沿程阻力相對局部阻力而言可以忽略.質量流量M=4.706 kg/s,直徑d35=d36=0.101 m,彎管截面積S35=S36=0.008 2 m2.在壓力為3.8 MPa,溫度為103.15 ℃時,速度u=27.435 m/s.局部阻力因數根據參考資料分別取ξ35=0.30和ξ36=0.18.代入公式得局部損失hj,35=11.509 m,hj,36=6.905 m.

圖 10直徑101.6 mm(4英寸)彎管1

Fig.10Pipe 1 of 101.6 mm(4 m) diameter

根據上例,可以得到直徑101.6 mm(4英寸)彎管1的總水頭損失,匯總見表4.

表 4壓力為3.8 MPa,溫度為103.15 ℃時彎管1總損失

Tab.4Total pressure loss of pipe 1 under 3.8 MPa pressure and 103.15 ℃ temperature劃分段34353637總損失直徑d/mm102.0102.0截面積S/m20.008 20.008 2長度L/m0.390.39流量Q/(m3/s)0.226 30.226 3速度v/m/s27.43527.435沿程阻力因數5 6505 650局部阻力因數0.300.18沿程水頭損失/m0.0400.0670.107局部水頭損失/m11.5096.90518.414endprint

同理,采油樹整體的壓力損失計算結果見表5和6,采油樹內流體流動過程見圖4.其中壓力損失ΔP=ρgh (12)式中:h為總水頭損失,單位m.由表6可知,當管道中流動的介質完全是甲烷氣體時候,最大的壓力損失為1.596 3 MPa.

表 5工況1采油樹總流阻計算結果

Tab.5Total flow resistance calculation results of christmas tree under condition 1構件沿程水頭損失/m局部水頭損失/m總水頭頭損失/m壓力損失/MPa油管掛本體0.056 7.7267.7820.002直徑127.0 mm(5英寸)彎管0.112 36.48136.5930.007雙孔連接器0.6806 370.5166 371.1961.298直徑101.6 mm(4英寸)彎管10.12418.41418.5380.004生產隔離閥0.04300.0430直徑101.6 mm(4英寸)管20.41519.81520.230.004總損失1.4306 452.9526 454.3821.315

表 6工況2采油樹總流阻計算結果

Tab.6Total flow resistance calculation results of christmas tree under condition 2構件沿程水頭損失/m局部水頭損失/m總水頭頭損失/m壓力損失/MPa油管掛本體0.0061.0891.0950.000 7127.0 mm(5英寸)彎管0.0135.1205.1330.003 2雙孔連接器0.0792 450.2372 450.3161.588 4直徑101.6 mm(4英寸)彎管10.0142.7882.8020.001 8生產隔離閥0.00500.0050直徑101.6 mm(4英寸)管20.0483.3993.4470.002 2總損失0.1652 462.6332 462.7981.596 3

4結束語

利用FLUENT對采油樹內流體流動特點和壓力的變化進行模擬,并計算流阻.結果表明:2種工況的壓力和速度變化較均勻,只是在變直徑處或直角彎管處會出現湍流現象.由流阻計算可知壓力損失約1.596 3 MPa.參考文獻:

[1]楊致政, 張瀛, 胡啟國. 近海邊際油田勘探開發中若干關鍵問題探討[J]. 特種油氣藏, 2005, 12(5): 59.

YANG Zhizheng, ZHANG Ying, HU Qiguo. Approach to key points in exploration and development of offshore marginal oilfields[J]. Special Oil & Gas Reservoirs, 2005, 12(5): 59.

[2]郭西水, 張林, 肖良平, 等. 深水完井技術研究[J]. 重慶科技學院學報(自然科學版), 2011, 13(4): 7476.

GUO Xishui, ZHANG Lin, XIAO Liangping, et al. Research on deep water well completion technology[J]. J Chongqing Univ Sci & Technol (Nat Sci), 2011, 13(4): 7476.

[3]馬彥金, 荊穎, 王亮. 對我國未來幾種深海油氣田開發新方案的探究[C]// 左其華, 竇希萍. 第十五屆中國海洋(岸)工程學術討論會論文集(上). 北京: 海洋出版社, 2011: 522524.

[4]方華燦. 對我國深海油田開發工程中幾個問題的淺見[J]. 中國海洋平臺, 2006, 21(4): 18.

FANG Huacan. A few superficial views about some problems of engineering of oil eield development in deep water in china[J]. China Offshore Platform, 2006, 21(4): 18.

[5]任克忍, 王定亞, 周天明, 等. 海洋石油水下裝備現狀及發展趨勢[J]. 石油機械, 2008, 36(9): 151153.

REN Keren, WANG Dingya, ZHOU Tianming, et al. Status and development of offshore oil and underwater equipment[J]. China Petroleum Machinery, 2008, 36(9): 151153.

[6]楊進, 劉書杰, 周建良, 等. 深水石油鉆采工程模擬試驗裝置的研制[J]. 石油機械, 2011, 39(8): 13.

YANG Jin, LIU Shujie, ZHOU Jianliang, et al. development of the simulation experimental device for deepwater petroleum drilling and production engineering[J]. China Petroleum Machinery, 2011, 39(8): 13.

[7]趙宏林, 程華榮, 田紅平, 等. 深水采油樹油管懸掛器下放工具關鍵參數設計[J]. 石油機械, 2014, 42(4): 1619.

ZHAO Honglin, CHENG Huarong, TIAN Hongping, et al. Design of the key parameters for the running tool of subsea christmas tree tubing hanger[J]. China Petroleum Machinery, 2014, 42(4): 1619.endprint

[8]龔銘煊, 劉再生, 段夢蘭, 等. 深海水下采油樹安放過程分析與研究[J]. 石油機械, 2013, 41(3): 5054.

GONG Mingxuan, LIU Zaisheng, DUAN Menglan, et al. Research on the running installment process of deep sea underwater christmas tree[J]. Petroleum Machinery, 2013, 41(3): 5054.

[9]劉文霄, 鄧平, 金弢.海洋水下立式采油樹用油管懸掛器: 中國, 201020530991[P]. 20110330.

[10]曾鳴, 孫亮, 鐘朝廷, 等.海洋鋪管船用張緊器內懸架的設計與分析[J].石油礦場機械, 2010, 39(6): 3235.

ZENG Ming, SUN Liang, ZHONG Zhaoting, et al. Design and analysis in suspension of pipe tensioner machine used in offshore pipelay vessel[J]. Oil Field Equipment, 2010, 39(6): 3235.

[11]朱江, 周文勝.中國近海油氣田區域開發戰略思考[J].中國海上油氣, 2009, 21(6): 380382.

ZHU Jiang, ZHOU Wensheng. Pondering over regional development strategy for oil and gas fields offshore China[J]. China Offshore Oil & Gas, 2009, 21(6): 380382.

[12]曾憲錦. 海上油氣田生產系統[M]. 北京: 石油工業出版社, 1993: 139156.

[13]李萬平. 計算流體力學[M]. 武漢: 華中科技大學出版社, 2004: 5458.

[14]劉玉林. 管道沿程水力損失的試驗研究[J]. 排灌機械, 1988, 6(5): 2326.

LIU Yulin. Experimental study on hydraulic loss along the course of the pipeline[J]. Drainage & Irrigation Machinery, 1988, 6(5): 2326.(編輯武曉英)endprint

猜你喜歡
有限元
基于擴展有限元的疲勞裂紋擴展分析
非線性感應加熱問題的全離散有限元方法
TDDH型停車器制動過程有限元分析
新型有機玻璃在站臺門的應用及有限元分析
上海節能(2020年3期)2020-04-13 13:16:16
基于I-DEAS的履帶起重機主機有限元計算
基于有限元模型對踝模擬扭傷機制的探討
10MN快鍛液壓機有限元分析
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 久久久久亚洲AV成人人电影软件| 欧美日本一区二区三区免费| 久久这里只精品国产99热8| 国产精品毛片在线直播完整版| 国产白浆在线| 91久久国产综合精品女同我| 午夜日b视频| 国产一区三区二区中文在线| 久久综合国产乱子免费| 欧美自慰一级看片免费| 伊人成人在线| 免费观看国产小粉嫩喷水| 国产人妖视频一区在线观看| 国内精品免费| 亚洲欧洲综合| 国产呦精品一区二区三区网站| 99尹人香蕉国产免费天天拍| 亚洲日本中文综合在线| 老色鬼久久亚洲AV综合| 中文字幕有乳无码| 日韩人妻精品一区| 国内精品伊人久久久久7777人| 国产在线观看99| 国产白浆一区二区三区视频在线| 一本二本三本不卡无码| 嫩草在线视频| 欧美亚洲国产视频| 亚洲综合经典在线一区二区| 狠狠亚洲五月天| 国产a在视频线精品视频下载| 蜜臀av性久久久久蜜臀aⅴ麻豆| Jizz国产色系免费| 国产毛片久久国产| 亚洲成人黄色在线| 国产高潮流白浆视频| 国产一级α片| 亚洲日本在线免费观看| 尤物成AV人片在线观看| 欧洲亚洲欧美国产日本高清| 国产最爽的乱婬视频国语对白| 国产情侣一区二区三区| 黄色网在线免费观看| 色综合久久88| 无码人妻免费| 国产免费福利网站| 97超爽成人免费视频在线播放| 久久久久国产一级毛片高清板| 中文字幕日韩丝袜一区| 欧美一区二区啪啪| 三级国产在线观看| 人妻91无码色偷偷色噜噜噜| 色香蕉网站| 四虎成人精品在永久免费| 99热这里都是国产精品| 亚洲日韩精品欧美中文字幕| 97影院午夜在线观看视频| 国产 在线视频无码| 欧美国产精品不卡在线观看| 亚洲AV无码乱码在线观看代蜜桃| 又黄又爽视频好爽视频| 中文无码日韩精品| 无码免费的亚洲视频| 亚洲第一成年网| 超碰91免费人妻| 日韩成人免费网站| 国产精品 欧美激情 在线播放 | 国产最爽的乱婬视频国语对白| 亚洲成人在线免费| 国产无码在线调教| 欧美伦理一区| 亚洲美女高潮久久久久久久| 日本欧美午夜| 国产永久无码观看在线| 日韩av手机在线| 亚洲另类色| 国产国产人在线成免费视频狼人色| 欧美 亚洲 日韩 国产| 欧美成人看片一区二区三区 | 欧美五月婷婷| 国产理论最新国产精品视频| 亚洲人成网址| 青草免费在线观看|