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

微傾管中低含液率氣液分層流臨界攜液流速預(yù)測模型

2019-04-09 09:12:16蒲雪雷王武杰閆敏敏王亮亮
天然氣工業(yè) 2019年12期
關(guān)鍵詞:界面模型

潘 杰 蒲雪雷 王武杰 閆敏敏 王亮亮

1.西安石油大學(xué)石油工程學(xué)院 2.上海理工大學(xué)新能源科學(xué)與工程研究所 3.中國石油長慶油田公司第四采氣廠

1 研究背景

氣井采出的天然氣往往伴隨著游離水、水蒸氣和液態(tài)烴等成分,被稱為濕天然氣[1]。濕天然氣在管道輸送過程中隨著管道沿線壓力、溫度等參數(shù)的變化,會在管道內(nèi)部析出液體[2]。析出的液體在重力的作用下逐漸在管道低洼處沉積形成積液。管道積液會影響天然氣管道的輸送效率、腐蝕管道甚至?xí)氯艿酪鹗鹿蔥3]。因此,有必要對濕天然氣管道的積液特性開展深入研究。

在濕天然氣管道輸送過程中,增加氣體流速能夠?qū)⒐艿赖屯萏幍姆e液攜帶出去。一般將能夠攜帶管道積液的最小氣相表觀流速稱為臨界攜液流速[4]。圖1給出了通過增加氣體表觀流速攜帶管道積液的過程。當(dāng)氣相表觀流速為零時(shí),管道中析出的液體沉積在管道底部形成積液(圖1-a);當(dāng)氣相表觀流速小于臨界攜液氣體流速時(shí),管道積液在氣相剪切力的作用下有沿著管壁向上運(yùn)動的趨勢,同時(shí)積液在重力的作用下不均勻地分布在管道內(nèi)壁上,管道低洼處液相層較厚,越靠近管道上部液相層越薄(圖1-b);當(dāng)氣相表觀流速達(dá)到臨界攜液流速時(shí),管道中的氣液兩相流動過程處于臨界狀態(tài),液相的凈流量和平均流速均為零[5],管道積液在剪切力和重力的共同作用下均勻分布在管壁上(圖1-c)。當(dāng)氣相的表觀流速大于氣相臨界攜液流速時(shí),天然氣中析出的液體受到的剪切力大于重力回流,會沿著管道向上運(yùn)動,此時(shí)濕天然氣中析出的液體可以被氣體攜帶出去。因此,準(zhǔn)確預(yù)測臨界攜液流速對濕天然氣管道輸送具有重要意義[6-8]。

圖1 積液示意圖

早在1949年,Lockhart和Martinelli[9]就針對管內(nèi)氣液兩相流提出了壓降計(jì)算關(guān)聯(lián)式,并初步建立了流型過度準(zhǔn)則。1959年,Hoogendoorn[10]對水平光滑管中氣液混合物的流動特性進(jìn)行了實(shí)驗(yàn)研究,并給出了預(yù)測流動壓降的經(jīng)驗(yàn)關(guān)聯(lián)式。1973年,Beggs和Brill[11]通過引入含液率和壓力梯度關(guān)聯(lián)式研究了傾斜管中氣液兩相流的壓降特性。之后,國內(nèi)外學(xué)者先后建立了不同的界面形狀模型(圖2)來研究氣液兩相分層的流動特性[12-17]。1976年,Taitel和Dukler[12]最早提出了針對氣液兩相分層流的界面形狀模型——FLAT模型,并用于研究氣液兩相流型轉(zhuǎn)變機(jī)理。1989年,Hart等[13]提出了新的界面形狀模型——ARS(Apparent Rough Surface)模型,并在此基礎(chǔ)上建立了水平氣液管流摩擦壓降和含液率關(guān)聯(lián)式。1997年,Chen等[14]提出了雙圓環(huán)(Double-Circle)模型,并通過該模型對低含液率下水平管中氣液兩相分層波狀流的含液率和壓降進(jìn)行了理論研究。Grolman和Fortuin[15]在1997年提出了MARS(Modified Apparent Rough Surface)模型,并用來預(yù)測微傾管中氣液兩相流含液率和壓降。2014年,Birvalski[16]通過MARS模型對微傾管中氣液兩相分層流臨界攜液流速、臨界含液率和臨界壓降進(jìn)行了預(yù)測。Banafi和Talaie[17]在2014年通過實(shí)驗(yàn)研究提出了一個(gè)新的氣液兩相分層流界面形狀模型,并用于預(yù)測含液率和壓降特性。2015年,徐英等[18]通過ARS模型預(yù)測了水平管內(nèi)的氣液兩相流動壓降特性。

液滴夾帶現(xiàn)象在氣液兩相流中廣泛存在[19]。當(dāng)發(fā)生液滴夾帶時(shí),部分液體以液滴的形式分布在氣相中,會對氣液兩相流動特性產(chǎn)生較大的影響[20]。然而,上述研究均未考慮液滴夾帶的影響。

筆者針對微傾管中低含液率氣液兩相分層流,基于氣液兩相流動量平衡方程和新的氣—液界面形狀閉合關(guān)系式,建立了考慮液滴夾帶的臨界攜液流速預(yù)測模型。

圖2 不同氣—液界面形狀示意圖

2 模型建立

2.1 動量方程

當(dāng)微傾管中氣液兩相分層流處于臨界狀態(tài)時(shí),管內(nèi)液體在剪切應(yīng)力和重力回流的共同作用下,處于動態(tài)平衡狀態(tài),均勻分布在管壁上,液相平均流速為零。氣相沿著管道上部向上運(yùn)動,氣/液界面處的液相在氣相剪切應(yīng)力的作用下向上運(yùn)動,管道下部液相在重力作用下產(chǎn)生重力回流向下運(yùn)動。此時(shí),氣液兩相所受剪切應(yīng)力如圖3所示。

圖3 臨界狀態(tài)剪切應(yīng)力示意圖

氣相動量方程為:

式中(-dp/dx)GFE表示氣相流動壓力梯度,Pa/m;AGFE表示考慮液滴夾帶時(shí)的氣相流道面積,m2;SG表示氣相潤濕周長,m;Si表示氣液界面處潤濕周長,m;ρGFE表示考慮液滴夾帶時(shí)的氣相密度,kg/m3;g表示重力加速度,m/s2。

液相動量方程為:

式中(-dp/dx)LFE表示液相流動壓力梯度,Pa/m;ALFE表示考慮液滴夾帶時(shí)的液相流道面積,m2;SL表示液相潤濕周長,m;ρL表示液相密度,kg/m3。

氣/液界面處的剪切應(yīng)力相等,即:

在氣液兩相流動過程中,氣液兩相的流動壓力梯度相等,即:

聯(lián)立式(1)~(4),得到考慮液滴夾帶時(shí)的氣液兩相動量平衡方程為:

2.2 液滴夾帶率

液滴夾帶率表示分散在氣相中的液體質(zhì)量流量所占總液體質(zhì)量流量的比例,即液相在氣相中分散的液體量所占總液體量的份額。在氣液兩相分層流中,液滴夾帶率雖然較小,但仍對流動特性有較大的影響。

在氣液兩相分層流動過程中,氣相流速較高時(shí),氣液界面波峰處的液體在氣相剪切力的作用下,以液滴的形式進(jìn)入氣相,這個(gè)過程稱為霧化。氣相中的液滴在重力的作用下沉積到液膜中,這個(gè)過程稱為沉積。Pan和Hanratty[21]在2002年基于液滴霧化與沉積的動態(tài)過程建立了液滴夾帶率計(jì)算關(guān)聯(lián)式。采用該公式進(jìn)行液滴夾帶率計(jì)算:

式中FE表示液滴夾帶率,無量綱;FE,max表示最大液滴夾帶率,無量綱,取FE,max=1.0;vG,cr表示開始夾帶時(shí)的臨界氣體速度,m/s;d表示管道直徑,m;ρG表示未產(chǎn)生夾帶時(shí)的氣相密度,kg/m3;σ表示液相表面張力,N/m。

2.3 氣—液界面形狀

2014年,Banafi和Talaie[17]根據(jù)氣液兩相分層流動機(jī)理以及通過絲網(wǎng)電氣系統(tǒng)捕捉到的實(shí)驗(yàn)數(shù)據(jù),提出了一個(gè)新的氣液兩相分層流界面形狀預(yù)測模型。筆者采用了Banafi和Talaie[17]提出的界面形狀模型。當(dāng)考慮液滴夾帶時(shí),其氣—液界面形狀如圖4所示。

圖4中的氣—液界面形狀是對FLAT模型中的氣/液界面形狀的改進(jìn)。當(dāng)氣相和液相以極低的速度在管道中接觸時(shí),會形成一個(gè)扁平的界面;隨著氣相表觀流速增大,氣—液界面在氣相的剪切作用下會產(chǎn)生彎曲,且有隨界面速度移動的傾向[17]??拷鼙谔幍囊合嘣跓o滑移條件下速度接近于零,速度水頭被轉(zhuǎn)化為靜壓能,使液相在管壁處形成薄層;除管壁附近的薄膜外,氣液界面基本是水平的。

圖4 氣—液界面形狀示意圖

液滴夾帶會改變氣液兩相的流道面積。當(dāng)存在液滴夾帶時(shí),液相流道面積為:

式中π表示圓周率,無量綱;hL表示底層液膜高度,m。

式中A表示管道橫截面面積,m2。

管道橫截面面積為:

式中S表示管道橫截面周長,m。

式中HL表示考慮液滴夾帶時(shí)液膜部分對應(yīng)的含液率,無量綱。

管壁處形成薄膜的長度為:

式中R表示管道半徑,m。

由于壁面效應(yīng),管壁附近的液相層以較低的速度流動,靠近壁面的位置,速度水頭被轉(zhuǎn)化為靜壓能,使液相在管壁處形成一層薄膜[17],因此:

式中Vsg表示氣相表觀流速,m/s。

由式(19)可以得到管壁處液膜的高度:

當(dāng)發(fā)生液滴夾帶時(shí),部分液體以液滴的形式分布在氣相中,從而改變氣相密度。考慮液滴夾帶時(shí)的氣相密度為:

Soleimani和Talaei[22]的實(shí)驗(yàn)研究表明,管壁處液膜的厚度是恒定的。管壁處液膜厚度為:

式中dL表示液相水力直徑,m。

液相水力直徑(dL)為:

將式(23)代入式(22)可得到管壁處液膜的厚度為:

2.4 剪切應(yīng)力

2.4.1 氣—壁剪切應(yīng)力

通常,剪切應(yīng)力被認(rèn)為是流體動能(單位體積)與界面摩擦因子的乘積[23],可以表示為:

式中fG表示氣/壁摩擦因子,無量綱。

筆者采用Gregory和Fogarasi[24]提出的方法計(jì)算氣/壁摩擦因子:

式中k/d表示管道相對粗糙度,無量綱;ReG表示氣相雷諾數(shù),無量綱;A5、A6分別表示經(jīng)驗(yàn)系數(shù),無量綱。

氣相雷諾數(shù)為:

式中μG表示氣相的動力黏度,Pa·s;θ表示液相濕壁分?jǐn)?shù),無量綱;si表示氣/液界面處的濕壁分?jǐn)?shù),無量綱。

2.4.2 氣—液界面剪切應(yīng)力

氣—液界面剪切應(yīng)力為:

式中fi表示反映氣—液界面摩擦因子,無量綱。

筆者采用Banafi和Talaie[17]提出的方法計(jì)算氣—液相界面處的摩擦因子(fi):

筆者采用Grolman和Fortuin[15]提出的方法計(jì)算管道相對粗糙度(k/d):

式中f1、f2、Fn表示中間變量,無量綱;μL表示液相的動力黏度,Pa·s。

管道相對粗糙度的詳細(xì)計(jì)算步驟如圖5所示。

圖5 管道粗糙度算法示意圖

2.4.3 液—壁剪切應(yīng)力

采用Biberg提出的方法計(jì)算臨界狀態(tài)下的液/壁剪切應(yīng)力:

式中δL表示一個(gè)量化液相含率的角度,(°);函數(shù)f(δL)表示不具有分析解的積分項(xiàng)[25],當(dāng)HL從0到0.5時(shí),δL從0到π/2。

當(dāng)處于臨界狀態(tài)時(shí)液相的平均真實(shí)流速等于零即vL=0 m/s。因此,在臨界狀態(tài)時(shí)液/壁的剪切應(yīng)力為:

Birvalski[16]給出了傾角1.3°~2.1°時(shí)f(δL)的平均值為0.25,但在實(shí)驗(yàn)中可發(fā)現(xiàn)f(δL)的值隨著管道的傾角增大而不斷增大,并且受液相介質(zhì)影響,實(shí)驗(yàn)測得其變化范圍為0~1.0。筆者依據(jù)本文參考文獻(xiàn)[16]中的實(shí)驗(yàn)數(shù)據(jù),擬合出了f(δL)隨管道傾角β變化的函數(shù)關(guān)系式。

當(dāng)液相介質(zhì)為密度ρW=997 kg/m3、黏度μW=0.001 Pa·s的水時(shí):

當(dāng)液相介質(zhì)為密度ρG/W=1 150 kg/m3、黏度μG/W=0.010 8 Pa·s的60%甘油/水時(shí):

3 模型求解方法

由于所建立的氣液兩相流動動量平衡方程是關(guān)于含液率的隱函數(shù),即

因此,需要對其進(jìn)行迭代求解[26]。圖6給出了特定運(yùn)行條件下函數(shù)F(HL)1/3的值隨含液率(HL)的變化規(guī)律。如圖6所示,當(dāng)氣相表觀速度(vsg)為8.700 m/s時(shí),函數(shù)曲線位于圖中水平線[方程F(HL)1/3=0對應(yīng)的直線]的上方,兩者沒有交點(diǎn),表明函數(shù)沒有零解,此時(shí)氣相表觀流速大于臨界攜液流速;當(dāng)vsg=8.600 m/s時(shí),函數(shù)曲線與圖中水平線有2個(gè)交點(diǎn),此時(shí)函數(shù)有2個(gè)零解,表明氣相表觀流速小于臨界攜液流速;當(dāng)vsg=8.649 m/s時(shí),函數(shù)曲線與圖中水平線只有1個(gè)交點(diǎn),此時(shí)函數(shù)只有1個(gè)零解,此時(shí)氣相表觀流速即為臨界攜液流速。因此,可以通過迭代氣相表觀流速使函數(shù)F(HL)1/3有且只有1個(gè)零解,來得到氣相臨界攜液流速。所建立模型的詳細(xì)求解過程如圖7所示。

圖6 函數(shù)F(HL)1/3的值隨含液率(HL)的變化曲線

4 模型驗(yàn)證

采用本文參考文獻(xiàn)[16]中的空氣—水、空氣—60%甘油/水兩相流實(shí)驗(yàn)數(shù)據(jù)對筆者提出的模型以及FLAT模型、ARS模型、雙圓環(huán)模型、MARS模型進(jìn)行了驗(yàn)證和對比,其計(jì)算結(jié)果如圖8~11所示。附具體實(shí)驗(yàn)參數(shù)如下:管徑(d)為0.050 8 m,傾角分別為1.3°、1.7°、2.1°;空氣的密度(ρG)為1.19 kg/m3,黏度(μG)為1.83×10-5Pa · s;水的密度(ρW)為997 kg/m3,黏度(μW)為0.001 Pa·s,表面張力(σW)為0.072 8 N/m;60%甘油/水的密度(ρG/W)為1 150 kg/m3,黏度(μG/W)為0.010 8 Pa·s,表面張力(σG/W)為0.067 6 N/m。

從圖8~11可以看出,當(dāng)液相分別為水和60%的甘油/水混合物時(shí),模型的臨界攜液流速預(yù)測結(jié)果均與實(shí)驗(yàn)值符合較好,F(xiàn)LAT模型、ARS模型、雙圓環(huán)模型以及MARS模型的臨界攜液流速預(yù)測結(jié)果均小于實(shí)驗(yàn)值;ARS模型和雙圓環(huán)模型的臨界攜液流速預(yù)測結(jié)果最小,與實(shí)驗(yàn)值偏差較大,F(xiàn)LAT模型和MARS模型的臨界攜液流速預(yù)測結(jié)果略高于ARS模型和雙圓環(huán)模型。當(dāng)液相為水時(shí),模型的臨界含液率預(yù)測結(jié)果略高于實(shí)驗(yàn)值,F(xiàn)LAT模型和MARS模型的臨界含液率預(yù)測結(jié)果與實(shí)驗(yàn)值符合較好,ARS模型和雙圓環(huán)模型預(yù)測的臨界含液率均低于實(shí)驗(yàn)值。當(dāng)液相為60%的甘油/水的混合物時(shí),模型的臨界含液率預(yù)測結(jié)果與實(shí)驗(yàn)值符合較好,其他模型的臨界含液率預(yù)測結(jié)果均低于實(shí)驗(yàn)值。綜合評價(jià)結(jié)果表明,模型的預(yù)測結(jié)果與實(shí)驗(yàn)值符合最好,可用于預(yù)測微傾管道中氣液兩相分層流的臨界攜液流速。

圖7 模型求解過程示意圖

5 計(jì)算結(jié)果分析

根據(jù)提出的模型對微傾管道中天然氣—水、天然氣—60%甘油/水分層流臨界攜液流速和臨界含液率進(jìn)行了計(jì)算分析。其中管徑(d)為0.050 8 m,工質(zhì)物性參數(shù)如下:天然氣的密度和黏度采用NISTRefprop軟件包計(jì)算,溫度取25 ℃;水的密度(ρW)為997 kg/m3,黏度(μW)為0.001 Pa·s,表面張力(σW)為0.072 8 N/m;60%甘油/水的密度(ρG/W)為1 150 kg/m3,黏度(μG/W)為0.010 8 Pa·s,表面張力(σG/W)為0.067 6 N/m。

圖8 不同模型的臨界攜液流速計(jì)算結(jié)果對比圖(空氣—水)

圖9 不同模型的臨界含液率計(jì)算結(jié)果對比圖(空氣—水)

圖10 不同模型的臨界攜液流速計(jì)算結(jié)果對比圖

圖11 不同模型的臨界含液率計(jì)算結(jié)果對比圖

5.1 管道傾角的影響

圖12給出了運(yùn)行壓力為0.24 MPa、氣相為天然氣(95%甲烷+5%乙烷)、液相分別為水和60%甘油/水時(shí),不同管道傾角條件下的臨界攜液流速和臨界含液率。從圖12可以看出,隨著管道傾角的增大,臨界攜液流速持續(xù)增大,臨界含液率逐漸減小。這是由于重力導(dǎo)致的液相回流作用隨管道傾角的增大而增強(qiáng),不利于氣體攜液。

圖12 管道傾角的影響圖

5.2 運(yùn)行壓力的影響

圖13給出了管道傾角為2.1°、氣相為天然氣(95%甲烷+5%乙烷)、液相分別為水和60%甘油/水時(shí),不同運(yùn)行壓力條件下的臨界攜液流速和臨界含液率。從圖13可以看出,隨著管道運(yùn)行壓力的增大,臨界攜液流速持續(xù)減小,臨界含液率逐漸增大。主要原因是隨著運(yùn)行壓力的增大,天然氣的密度和黏度逐漸增大,導(dǎo)致其攜液能力增強(qiáng)。

圖13 天然氣管道壓力的影響圖

5.3 液相密度的影響

圖14給出了運(yùn)行壓力為0.24 MPa、管道傾角為2.1°、氣相為天然氣(95%甲烷+5%乙烷)、液相分別為水和60%甘油/水時(shí),不同液相密度條件下的臨界攜液流速和臨界含液率。從圖14可以看出,隨著液相密度的增加,臨界攜液流速持續(xù)增大,臨界含液率逐漸減小。這是由于隨著液相密度的增大,重力導(dǎo)致的液相回流作用逐漸增強(qiáng),不利于氣體攜液。

圖14 液相密度的影響圖

5.4 天然氣組成的影響

圖15給出了運(yùn)行壓力為0.24 MPa、管道傾角為2.1°、氣相為天然氣、液相分別為水和60%甘油/水混合物時(shí),不同天然氣組分條件下的臨界攜液流速和臨界含液率。圖1中組分1為100%甲烷;組分2為95%甲烷+5%乙烷;組分3為90%甲烷+10%乙烷;組分4為90%甲烷+5%乙烷+5%丙烷;組分5為85%甲烷+10%乙烷+5%丙烷;組分6為85%甲烷+5%乙烷+5%丙烷+5%丁烷。從圖15可以看出,隨著甲烷含量的逐漸降低和重組分含量的逐漸增加,臨界攜液流速持續(xù)減小,臨界含液率略有增大。這是由于隨著天然氣組分的改變,天然氣的密度和黏度隨重組分比例的增大而增大,提高了攜液能力。

圖15 天然氣組分的影響圖

6 結(jié)論

1)針對微傾管道中低含液率氣液兩相分層流,建立了考慮液滴夾帶的臨界攜液氣體流速預(yù)測模型。該模型引入了液滴夾帶率計(jì)算公式和新的氣—液界面形狀閉合關(guān)系式,同時(shí)擬合得到了考慮管道傾角影響的液/壁剪切應(yīng)力修正系數(shù)計(jì)算公式。

2)結(jié)合實(shí)驗(yàn)數(shù)據(jù),將本文模型和FLAT模型、ARS模型、雙圓環(huán)模型、MARS模型進(jìn)行了驗(yàn)證和對比。結(jié)果表明,當(dāng)液相為水時(shí),模型的臨界攜液流速預(yù)測結(jié)果均與實(shí)驗(yàn)值符合較好,臨界含液率預(yù)測結(jié)果略高于實(shí)驗(yàn)值;當(dāng)液相為60%甘油/水時(shí),模型的臨界攜液流速和臨界含液率預(yù)測結(jié)果均與實(shí)驗(yàn)值符合較好。

3)管道傾角、運(yùn)行壓力、液相密度以及天然氣組分均會對臨界攜液流速產(chǎn)生較大影響。隨著管道傾角和液相密度的增大,臨界攜液流速持續(xù)增大,臨界含液率逐漸減小。隨著運(yùn)行壓力和天然氣中重組分含量的增大,臨界攜液流速持續(xù)減小,臨界含液率逐漸增大。

猜你喜歡
界面模型
一半模型
重要模型『一線三等角』
國企黨委前置研究的“四個(gè)界面”
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
基于FANUC PICTURE的虛擬軸坐標(biāo)顯示界面開發(fā)方法研究
空間界面
金秋(2017年4期)2017-06-07 08:22:16
電子顯微打開材料界面世界之門
人機(jī)交互界面發(fā)展趨勢研究
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 国产欧美精品午夜在线播放| 日本欧美视频在线观看| 91免费观看视频| 久草美女视频| 国产AV毛片| 波多野结衣的av一区二区三区| 99久久无色码中文字幕| 999精品在线视频| 影音先锋丝袜制服| 国产原创第一页在线观看| 99久久免费精品特色大片| 国产原创演绎剧情有字幕的| 91在线一9|永久视频在线| 素人激情视频福利| 亚洲精品无码AⅤ片青青在线观看| 大陆精大陆国产国语精品1024| 不卡网亚洲无码| 国产成人做受免费视频| 欧美日韩在线亚洲国产人| 久久一色本道亚洲| 51国产偷自视频区视频手机观看| 91av国产在线| 中文字幕有乳无码| 亚洲国产欧美中日韩成人综合视频| 无码中文AⅤ在线观看| 色偷偷一区二区三区| 亚洲欧美成人在线视频| 久久久久无码国产精品不卡| 伊人色综合久久天天| 国产精品手机视频一区二区| 久久香蕉国产线看观看式| 久久无码高潮喷水| 91年精品国产福利线观看久久| 午夜视频在线观看免费网站| 丁香婷婷久久| 欧美精品不卡| 国产91蝌蚪窝| 婷婷色在线视频| 亚洲第一成年网| 国外欧美一区另类中文字幕| 亚洲AV无码一二区三区在线播放| 一区二区三区四区日韩| 99热这里只有免费国产精品 | 美女无遮挡免费网站| 国产日韩久久久久无码精品| 国产成人91精品免费网址在线| 成人在线观看一区| AV在线天堂进入| 在线国产毛片| 欧美日韩成人在线观看| 国产av无码日韩av无码网站| 妇女自拍偷自拍亚洲精品| 无码在线激情片| 萌白酱国产一区二区| 久久精品国产电影| 免费 国产 无码久久久| 国产呦视频免费视频在线观看| 亚洲天堂精品在线| 国产va免费精品观看| 在线观看视频99| 欧美综合成人| 免费视频在线2021入口| 亚洲日本www| 国产区免费| 中文字幕av一区二区三区欲色| 久久国语对白| 亚洲天堂久久久| 国产欧美在线| 国产精品99在线观看| 一区二区三区成人| 亚洲国产成人精品一二区| 国产91视频免费观看| 精品成人一区二区三区电影| 国产成人a毛片在线| 欧美不卡视频一区发布| 成人午夜天| 日本免费一区视频| 在线免费a视频| 国产69精品久久| 视频一区视频二区中文精品| 国产不卡网| 欧美v在线|