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

充氣鉆井環(huán)空彈狀流穩(wěn)態(tài)流動模型與實例驗證

2020-04-24 09:24:26徐新紐曹光福
科學技術與工程 2020年5期
關鍵詞:模型

徐新紐, 阮 彪, 曹光福, 黃 鴻, 楊 虎

(1.中國石油新疆油田分公司,克拉瑪依 834000;2.自然資源部深部地質鉆探技術重點實驗室, 克拉瑪依 834000;3.中國石油大學(北京)克拉瑪依校區(qū),克拉瑪依 834000)

目前,在氣液兩相流體力學中常用的模型有流動型態(tài)、均相流動、分相流動和漂移流動等模型[1]。國外許多學者[2-6]對氣液兩相管流的流動型態(tài)判別及不同流型下不同管路的流動參數、含氣率、流動壓降等方面進行實驗和半經驗性研究,建立了許多數學模型[7]。其中,對于鉆井管柱和環(huán)空內的泡狀流、分散泡狀流,研究的學者和成果較多,部分計算模型已廣泛應用于石油工業(yè)。但是,由于各種模型的假設條件與實際流動的差異,流動模型的計算結果誤差不一[4,7-8]。

充氣鉆井井筒內存在多種流型,對于彈狀流(段塞流)、環(huán)狀流、攪拌流等流型的判別和轉換關系,國內外尚未統(tǒng)一標準[9-11]。現(xiàn)有穩(wěn)態(tài)多相流動所開展的實驗均以低黏牛頓流體(清水)為實驗介質,而鉆井條件下高黏、非牛頓流體鉆井液與其性質大為不同[12-13]。因此,井下復雜環(huán)境的流態(tài)演化與工業(yè)界通用實驗圖版有較大差異[14-16]。

充氣鉆井液作為欠平衡鉆井流體,其氣液比和實現(xiàn)的井底當量密度范圍較大,在鉆井環(huán)空上部易出現(xiàn)彈狀流形態(tài)[17]。國內外多種軟件均采用常規(guī)穩(wěn)定泡狀流力學模型,井筒參數的計算誤差較大[14-16]。因此,有必要調研國內外氣液兩相流流型研究成果,結合充氣鉆井流體流動實際情況,建立較準確的充氣鉆井環(huán)空彈狀流穩(wěn)態(tài)力學模型,用于井筒流體參數的預測和設計。

1 鉆井氣液兩相管流形態(tài)描述

兩相流體同時在管道內流動可能產生的不同相交界面簡稱兩相流型,不僅與每一相是層流或湍流有關,而且與兩相交界面的變化和組合有關。實驗表明[3-6],垂直向上的兩相管流經常出現(xiàn)5種流型,即泡狀流、彈狀流、攪拌流、環(huán)狀流、霧狀流,如圖1所示。

圖1 垂直氣液兩相管流流型Fig.1 Gas-liquid two-phase flow patterns in vertical pipe

彈狀流是當泡狀流中氣相流量增大到一定值,將發(fā)生氣泡聚合,甚至聚合成接近管徑的彈狀氣泡(泰勒泡),即泡狀流過渡為彈狀流。大塊彈狀氣泡與含有彌散小氣泡的液塊在管內間隔流動,彈狀氣泡外圍液相以液膜狀向下流動。攪拌流僅出現(xiàn)在通徑較大的管流中,液相呈不定型向上或向下振蕩運動,呈攪拌狀態(tài)。然而對于鉆井常規(guī)尺寸環(huán)空不會出現(xiàn)攪拌流,可能會由彈狀流向環(huán)狀流直接平穩(wěn)過渡。環(huán)狀流是指液相沿管壁呈膜狀流動,氣相在管道芯部流動。實際管流中呈現(xiàn)環(huán)狀流的工況參數范圍很窄。

通常,井筒兩相流型取決于兩相流體流量、流體性質和井眼幾何參數等。液體欠平衡鉆井和充氣鉆井考慮到井控和地面流體處理設備的安全性,需附加一定的井口回壓減小氣體返速,以避免井口出現(xiàn)高速氣流并對返出管線造成損壞。對于充氣鉆井,氣、液兩相的流量范圍分別為10~50 m3/min和0.189~1.325 m3/min。Lage[5]繪制的環(huán)空兩相流流型分布如圖2所示。由圖2可知,對于多數充氣鉆井(井深1 500~5 500 m),井口環(huán)空氣體在一定井口回壓下,井口附近環(huán)空流型以彈狀流為主,攪拌流或環(huán)狀流出現(xiàn)概率很小,可視為彈狀流處理,而井筒環(huán)空中下部以泡狀流或分散泡狀流為主[1,6]。此觀點與文獻[18-20]的環(huán)空充氣液主要出現(xiàn)泡狀流和彈狀流的實驗結果一致。

圖2 環(huán)空氣液兩相流流型隨壓力變化曲線Fig.2 Dominant gas-liquid two-phase flow patterns’ map for annular geometries

2 充氣鉆井環(huán)空彈狀流力學模型

2.1 假設條件

(1)充氣鉆井流體在環(huán)空中混合有巖屑和地層流體(天然氣、原油或地層水)。模型假設環(huán)空鉆井流體中氣相為注入氣體和地層氣體的混合物,液相為注入液體、地層液體與巖屑的混合物[21-23]。

(2)假設注入氣體與地層產出氣的流動性無差異,兩種氣體環(huán)空返速相同。同樣,假設注入液體和地層產出液環(huán)空返速相同。氣相和液相的密度為相應組分的混合密度,其中,假設巖屑為顆粒流[24]。

(3)考慮到環(huán)空氣液兩相流將產生較高的摩阻壓耗和表觀黏度[25-26],假設井眼凈化理論是在常規(guī)鉆井液流變學的基礎上,將氣相的環(huán)空返速、表觀黏度等參數特殊處理。

(4)模型中所有與位移有關的參數(包括井深、速度、壓力等)均為沿井眼的測量深度進行計算。因此力學模型適用于水平井[24]。

(5)井筒流體穩(wěn)定流動時,流體溫度沿井深為線性分布。

2.2 模型基本參數

由于井筒內氣、液兩相流體的密度和流速存在差異,易出現(xiàn)兩相流體間的滑脫現(xiàn)象[24-26]。因此,本節(jié)定義的基本參數有助于理解和建立環(huán)空彈狀流穩(wěn)態(tài)力學模型。

折算速度為管流的全部過流斷面僅被兩相流中的單相占據時的流動速度[6,14,24]。氣相和液相的折算速度分別為

(1)

(2)

由此,混合兩相流的速度vm可定義為

vm=vSL+vSG

(3)

式中:QL、QGSC分別為液體流量和氣體地面流量,m3/s;T、TSC分別為井筒溫度和地面溫度,K;P、PSC分別為井筒壓力和地面壓力,Pa;Ap為流道斷面面積,m2;Z為氣體壓縮因子。

真實含液率HL(截面含液率或持液率)為兩相流動的過流斷面中,液相面積占過流總斷面的份額。其取值為0~1,受氣、液相流體性質、流型和井筒幾何參數的影響。

(4)

若兩相流中氣、液體的流速相同,無滑脫現(xiàn)象,則無滑脫含液率λL定義為

(5)

實際速度為一相流體隨其他相流體流動的真實流速,此相流體的實際過流斷面小于過流總斷面。

(6)

(7)

對于環(huán)空同時存在鉆井液、原油和地層水時,通常這些液體間的滑脫現(xiàn)象與氣、液體間的滑脫現(xiàn)象相比微乎其微[27]。因此,假設環(huán)空各種液體間無滑脫,且鉆井基液的體積系數

(8)

式(8)中:QG、Qoil、Qw分別為氣體流量、原油流量和地層水流量,m3/s。

同樣,地層產出天然氣和鉆井注入氣體間的體積系數也如此計算。

2.3 環(huán)空彈狀流型判別

泡狀管流中分散氣泡的上升類似于泰勒泡[28]。文獻[5-7,22,24]給出泡狀流中分散氣泡的絕對舉升速度

v

(9)

而泰勒泡舉升速度的表達式為

(10)

式中:v∞為氣泡絕對舉升速度,m/s;vTB為泰勒泡舉升速度,m/s;Dep為環(huán)空等圓周直徑,m;ρG、ρL、ρm分別為氣體、液體和混合流體密度,kg/m3;g為重力加速度;σ為氣泡表面張力,N/m。

Taitel等[28]認為若分散氣泡的舉升速度大于泰勒泡的舉升速度,則分散氣泡合并成泰勒泡,此時泡狀流將轉變?yōu)閺棤盍鳌H魞烧吲e升速度大小相反,泰勒泡將穿過許多小氣泡,相對運動產生的席卷作用使泰勒泡分散成小氣泡,不能合并。因此,結合式(9)和式(10)為

(11)

文獻[5,20,22,27]通過實驗驗證非牛頓流體的含氣率應采用α=0.20,則推導出泡狀流與彈狀流間轉換界限的判別式(12),方程曲線如圖2所示。

vSL=3.167vSG-0.745v

(12)

采用式(12)針對注入氣液比較高(10∶1)的充氣鉆井進行模擬計算,當井深超過1 800 m時環(huán)空含氣率僅為5%,氣液比僅為0.05∶1。只有在井深約300 m的環(huán)空含氣率才超過50%,氣液比超過1.4∶1,才會出現(xiàn)彈狀流。

2.4 環(huán)空彈狀流物理模型

彈狀流表現(xiàn)為氣體和液體的交替流動[21,28],氣體為近似彈頭形狀的大氣泡(泰勒泡)。泰勒泡幾乎占據了整個流道截面,并均勻地上行流動。而彈狀流中的液體也包括兩部分:泰勒泡之間的液塞和管壁與泰勒泡間的下行液膜。液塞內包含許多小的球形氣泡,近似為泡狀流。通常,泰勒泡的氣泡帽長度LC很小,可以忽略不計,液膜厚度變化不大,計算時可采用整個液膜的平均厚度。

圖3所示為彈狀流體中生長完全泰勒泡的流動參數和物理模型。假設環(huán)空彈狀流的水動力學參數的推導過程是基于環(huán)空兩相流體為一維穩(wěn)定流動,流動參數公式推導如下:

圖3 彈狀管流中泰勒泡的流動物理模型Fig.3 Taylor bubble flow physical model in Slug flow unit

由于環(huán)空彈狀流(段塞流)表現(xiàn)為間歇性流動,任一流動斷面上的流動參數隨時間不斷變化,呈不穩(wěn)定流動狀態(tài)。通常,模型建立時假設一個彈狀流體單元由一個泰勒泡及其液膜和一段相鄰的液塞組成。因此,泰勒泡和液塞的絕對平移流速

(13)

式(13)中:um為液相和氣相的平均流速,m/s。

在一固定的坐標系中,以絕對平移速度流動的彈狀流體單元流動參數并不隨時間變化,僅在流道空間上變化。因此,假設環(huán)空向上流動為正方向,則液塞與泰勒泡之間液相的質量平衡方程為

(uT-uLLS)HLLS=(uT+|uLTB|)HLTB

(14)

在一個固定截面A—A(圖3)的流道內,一個流體單元勻速流動時,一個泰勒泡和一段液塞流過該截面所需時間分別為

(15)

(16)

于是泰勒泡區(qū)域和液塞區(qū)域液體體積分別為

(17)

(18)

由于一完整彈狀流體單元流過截面A—A所需總時間為ΔtSU=ΔtTB+ΔtLS,則彈狀流單元內液體總體積可由下式計算:

VLSU=uSLApΔtSU

(19)

式(17)和式(18)中的流道截面積均為持液率的函數,將式(15)和式(16)代入式(19),求得

(20)

圖3及公式中:uGTB為泰勒泡中氣體的流速,m/s;uLTB為泰勒泡中液體的流速,m/s;uGLS為液塞中氣體的流速,m/s;uLLS為液塞中液體的流速,m/s;HLTB為泰勒泡段兩相流的持液率;HLLS為液塞段兩相流的持液率;ALTB為泰勒泡段液流截面積,m2;ALLS為液塞段液相截面積,m2;LC為泰勒泡帽的長度,m;LTB為泰勒泡的長度,m;LLS為液塞段的長度,m;LSU為彈狀流單元的長度,m;(下標:TB表征泰勒泡區(qū)域;LS表征液塞區(qū)域;SU表征整體區(qū)域)。

考慮到液塞中的氣液兩相流動形態(tài)類似于穩(wěn)定的泡狀流,液塞內的氣體流速為

(21)

式(21)中:n為氣泡密集效應指數;C0為兩相管流速度剖面系數。

假設液塞的持液率為一常數,許多學者[29-30]提出了一些計算該持液率的方法。考慮到研究對象為鉆井環(huán)空中的彈狀流體流動,建議采用以下計算公式:

(22)

由于泰勒泡周圍的液膜呈自由下滑流動,其厚度δ最終為一固定值。

(23)

式(23)中:μL、μG分別為液相、氣相的黏度,Pa·s。

依據環(huán)空彈狀流(段塞流)的幾何參數,泰勒泡區(qū)域的持液率可表示為

(24)

式(24)中:DIC、DOT分別為井眼或套管內徑和鉆柱外徑,m。

通過大量試驗均證實,大多數向上管流的彈狀流體單元的液塞長度近似為管道直徑的16倍。后來,采用了水力直徑的概念,并且Caetano[31]證明該規(guī)律也適用于環(huán)空流道,于是,液塞長度

LLS=16Dh

(25)

式(25)中:Dh為環(huán)空水力直徑,m。

總之,已知環(huán)空氣體和液體流量、流體性能參數,可由式(22)、式(21)、式(13)和式(25)分別計算出液塞持液率、液塞中氣體流速、環(huán)空彈狀流單元平移速度和液塞長度;然后,依據式(14)、式(20)、式(23)和式(24)分別計算出泰勒泡長度、泰勒泡內液體流速、液塞內液體流速和泰勒泡區(qū)域的持液率。

2.5 環(huán)空彈狀流力學模型

Caetano[31]提出了垂直和水平管路中彈狀流的流動模型。近幾年,其他學者[32-35]建立的模型非常復雜,難以付諸于現(xiàn)場應用。為此,本模型充分考慮到彈狀流中液塞和泰勒泡持液率的變化,將液塞視為泡狀流,將泰勒泡視為環(huán)狀流,并利用經驗式(22)和式(24)分別計算兩部分的持液率,由此可計算液塞中氣泡的速度。因此,本模型認為,與彈狀流中整個泰勒泡長度相比,生長完全的泰勒泡的氣泡帽長度很小,可以忽略不計(圖3)。泰勒泡液膜厚度變化不大,計算時可采用整個液膜的平均厚度。由此可簡化壓降計算過程。

根據機械能守恒定律,環(huán)空穩(wěn)定彈狀流的總壓降由重力壓降、摩阻壓降和加速壓降組成[9]。

(26)

重力壓降為

(27)

摩阻壓降為

(28)

液塞混合區(qū)域的加速壓降為

(29)

(30)

式(30)中:ρR為地層巖石密度,kg/m3;vDR為機械鉆速,m/s;β為氣泡長度相對百分比。

上述方程中的摩阻系數由雷諾數計算得出,雷諾數定義為

(31)

3 模型數值解法

考慮到鉆井流體沿井筒(鉆柱和環(huán)空內)流道流動為一維穩(wěn)定流動。可利用一維數值迭代的解法進行計算。由于是穩(wěn)態(tài)模型,算法中無需考慮時間因素[36]。

將井筒流道進行一維迭代網格劃分(圖4),計算步長設定為1~10 m;邊界條件:井口溫度和井口壓力;迭代計算路徑為:井口→環(huán)空→鉆頭噴嘴→鉆柱內→井口鉆柱(立管內)。該鉆井方案的整個計算路徑均為氣液兩相流。

圖4 模型的計算路徑和迭代網格劃分Fig.4 Iterative calculation path and grid division of the model

(1)確定井口環(huán)空處的邊界條件:PL=0=PS,TL=0=TS,QGL=0=QGSC,QL等,L為井深,PL=0為井深為零時的環(huán)空壓力,PS為井口回壓,TL=0為井深為零時的環(huán)空溫度,TS為環(huán)空出口溫度;QGL=0為井深為零時的環(huán)空氣體流量,QGSC為環(huán)空出口氣體流量,QL為環(huán)空液體流量;確定井身結構和鉆具組合數據以計算井筒內的有關幾何尺寸;確定鉆頭參數:噴嘴總流面積或各噴嘴直徑等;確定水平井井眼軌跡數據:斜井段井斜角對應測深、垂深數據。

(2)在1#網格單元內視兩相流體參數穩(wěn)定取平均值,然后利用流體單元的邊界條件,判別兩相流流型,并沿積分路徑進行計算壓降,得出單元下端的壓力、溫度(P2,T2)。

(3)令2#單元(前一單元)的P1=P2,T1=T2,計算出2#單元的流體參數(密度、折算速度,持液率、含氣率、黏度、氣體壓縮因子、雷諾數、摩擦系數等),判別兩相流流型,并沿積分路徑進行計算壓降,得出單元下端的壓力、溫度(P2,T2)。

(4)重復程序(3)直至井底環(huán)空。注意:計算前確定環(huán)空幾何尺寸和兩相流流型的變化,采用相應的流動模型進行計算。

(5)利用氣液兩相流鉆頭壓降模型計算出鉆頭噴嘴上游壓力。

(6)重復程序(2)、(3)。注意:計算前確定鉆柱幾何尺寸和兩相流流型的變化,應采用相應的流動模型進行計算。鉆柱內計算路徑為井底向井口方向,井深差值為負值。計算過程中應考慮馬達、測量儀器、限流器及單流閥等的壓降。

(7)沿鉆柱內迭代計算結果:P2=Pin,Pin為鉆井泵壓(注入壓力)。

4 實例驗證

山西晉城大寧煤礦煤層氣DNP02井實現(xiàn)了兩井連通直井注空氣鉆水平羽狀分支井,水平井段為13個分支,總長為8 019 m[1,24]。該井注氣點位于水平段A點附近,其造斜段和直井段均為充氣鉆井液,環(huán)空內為氣液兩相流[37-38]。由于該井A點垂深僅為195.73 m,流型判別認為該井施加少量井口回壓的情況下連通點以上充氣環(huán)空均為彈狀流(圖5)。

表1所示為DNP02井環(huán)空充氣鉆井的基本數據[1,37]。空氣由直井內的7 in(1 in=25.4 mm)套管注入到水平井環(huán)空,注氣點以上環(huán)空為氣液兩相流,而注氣點以下的環(huán)空流體和鉆柱內流體均為單相液流。考慮到兩井連通的環(huán)空注氣方式下氣體注入點以上環(huán)空內的氣液兩相流動為一維穩(wěn)定流動,可利用數值解法進行計算。而對于環(huán)空注氣點以下的環(huán)空流體和鉆柱內流體均為單相液流,可利用單相液流的水力模型進行計算。另外,由于氣體注入管路中均為單相氣流,其水力參數可利用氣體鉆井水力模型進行計算。因此,該實例中井筒水力參數的計算可分為3部分:注氣點以上環(huán)空兩相流水力參數、注氣點以下環(huán)空及鉆柱內單相液流水力參數和注氣管內單相氣流水力參數。

圖5 DNP02井兩井連通環(huán)空充氣鉆井流程Fig.5 Two well-connected annulus aerated drilling process from DNP02 well

表1 DNP02井環(huán)空充氣鉆井的基本數據

注:①為動力黏度單位,1 cp=10-3Pa·s。

將環(huán)空兩相流流道進行一維網格劃分,其迭代計算路徑為:井口→氣體注入點以上環(huán)空→環(huán)空單相液流→鉆頭噴嘴→鉆柱內單相液流→井口鉆柱(立管內)。該實例中部分計算路徑為氣液兩相流,需采用數值迭代計算;而其余部分的單相液流,可利用液流水力模型直接計算,無須迭代計算。

實例計算步長為1~10 m,邊界條件為井口溫度和井口壓力,結合該井基本參數(表1),計算出該井環(huán)空流體參數隨井深的變化情況(表2)。

由于直井為空氣注入井,井深僅195.73 m。通過計算,井口注入壓力與井底壓力差僅為20 kPa,可忽略不計。因此,將地面空氣注入壓力視為水平井段兩井連通點的環(huán)空壓力。

將空氣注入管線上的PQT計顯示的空氣流量和壓力與模型計算的井底壓力結果進行對比(表3),該模型計算結果非常接近實際環(huán)空壓力(符合率達到95%),證實建立的環(huán)空彈狀流力學模型的準確性。

表2 DNP02井充氣鉆井環(huán)空流體參數

表3 DNP02井環(huán)空壓力模擬結果與實測數據對比

5 結論

(1)鉆井環(huán)空氣液兩相流主要有5種流型。其中,攪拌流僅出現(xiàn)在通徑較大的管流中,對于鉆井常規(guī)尺寸環(huán)空不會出現(xiàn)攪拌流。實際鉆井管流中呈現(xiàn)環(huán)狀流的工況參數范圍很窄。對于多數充氣鉆井,由于現(xiàn)場作業(yè)主動施加環(huán)空回壓,井口附近的攪拌流或環(huán)狀流出現(xiàn)的概率很小。

(2)調研國外學者的研究成果和實驗數據,結合充氣鉆井具體參數,認為只有在井深300 m以上的環(huán)空含氣率才超過50%,以彈狀流為主,而下部環(huán)空以泡狀流或分散泡狀流為主。

(3)本模型充分考慮到彈狀流中液塞和泰勒泡持液率的變化,將液塞視為泡狀流,將泰勒泡視為環(huán)狀流。為了簡化壓降模型,認為生長完全的泰勒泡的氣泡帽長度很小,可忽略不計,且泰勒泡液膜厚度變化不大,計算時可采用整個液膜的平均厚度。

(4)基于穩(wěn)態(tài)多相流理論和機械能守恒定律,詳細推導出的環(huán)空彈狀流穩(wěn)態(tài)流動模型,利用數值解法進行求解、程序化。DNP02井采用直井與水平井兩井連通注氣技術,該實例證實距井口小于200 m的環(huán)空中,氣液兩相流為彈狀流。通過與實測值的對比,證實本模型計算精度完全滿足充氣鉆井。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产真实乱子伦精品视手机观看 | 狼友av永久网站免费观看| 欧美国产在线看| 青青草原国产| 国产精品美女自慰喷水| 91欧美亚洲国产五月天| 97国产精品视频自在拍| 九九热视频精品在线| 久久国语对白| 成人精品在线观看| 久久一本日韩精品中文字幕屁孩| 亚洲婷婷丁香| 欧美在线黄| 黄色国产在线| 精品视频一区在线观看| 欧美成一级| 99在线视频免费观看| 狠狠色狠狠色综合久久第一次| 亚洲成年人网| 国内精品小视频在线| 国产亚洲精品91| 欧美一区日韩一区中文字幕页| 久久久久九九精品影院| 亚洲女同欧美在线| 国产精品不卡永久免费| 91尤物国产尤物福利在线| 久久人体视频| 精品无码人妻一区二区| 国产成人AV综合久久| 99视频在线看| 国产精品精品视频| 国产极品嫩模在线观看91| 91黄色在线观看| 日韩一区二区在线电影| 亚洲第一精品福利| 国产精品久久久久婷婷五月| 国产a v无码专区亚洲av| 看你懂的巨臀中文字幕一区二区 | 免费一级毛片| 免费人成黄页在线观看国产| 国产97视频在线| 亚洲色图欧美激情| 午夜性刺激在线观看免费| 免费看a毛片| 亚洲日韩精品无码专区| 日韩一级二级三级| 91精品综合| 亚洲av成人无码网站在线观看| 无码精品一区二区久久久| 亚洲av中文无码乱人伦在线r| 国产乱人伦偷精品视频AAA| 精品成人一区二区| 在线看免费无码av天堂的| 色噜噜在线观看| 国产91小视频在线观看| 国产极品美女在线观看| 国产午夜福利亚洲第一| 91精品免费高清在线| 99re热精品视频中文字幕不卡| 欧美无专区| 国产欧美另类| 国产第一页第二页| 蜜臀AVWWW国产天堂| 成人精品免费视频| 国产探花在线视频| 视频一区亚洲| 精品无码人妻一区二区| 久久久久久久97| 永久免费av网站可以直接看的| 精品伊人久久久香线蕉| 久久精品中文字幕少妇| 精品一区二区三区水蜜桃| 国产成人91精品| 99在线视频精品| 久草美女视频| 亚洲色大成网站www国产| 91久久精品国产| 国产18在线| 青青极品在线| 久久人人97超碰人人澡爱香蕉 | 9999在线视频| 美女被操黄色视频网站|