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

T型微通道內氣泡生成大小影響因素研究

2017-01-03 06:27:53呂明明劉志剛管寧王樹眾
山東科學 2016年5期
關鍵詞:界面

呂明明,劉志剛*,管寧,王樹眾

(1. 山東省科學院流動與強化傳熱重點實驗室,山東省科學院能源研究所,山東 濟南 250014;2. 西安交通大學能源與動力工程學院,陜西 西安 710049)

?

【能源與動力】

T型微通道內氣泡生成大小影響因素研究

呂明明1,劉志剛1*,管寧1,王樹眾2

(1. 山東省科學院流動與強化傳熱重點實驗室,山東省科學院能源研究所,山東 濟南 250014;2. 西安交通大學能源與動力工程學院,陜西 西安 710049)

采用流體體積函數(shù)(VOF)方法,對T型微通道中氣泡形成過程進行數(shù)值模擬研究,根據(jù)氣泡形成機理,分析了氣液流速、流體性質和微通道尺寸等因素對生成氣泡大小的影響。研究結果表明,T型微通道內生成氣泡長度隨氣體份額的增加呈指數(shù)增加趨勢,而在相同氣體份額下氣液流速對氣泡長度影響不大;比較而言,液體粘度和表面張力對生成氣泡大小的影響較小,當液相表面張力從0.072 N·m-1降低到0.01 N·m-1時,T型微通道內生成氣泡的長度減小了18%,主要是因為在阻塞階段,最大頸部寬度和塌陷時間減小了;氣泡長度隨微通道直徑的增加而增大,而氣泡的無量綱長度基本不受微通道直徑的影響。

微通道;數(shù)值模擬;氣泡;卡斷機理;表面張力

泡沫流體應用于油田開發(fā)在國內外已有40多年的歷史,如泡沫驅油就是一項重要的強化采油(Enhanced Oil Recovery, EOR)技術[1-4]。氣液在表面活性劑作用下形成的低流度泡沫可以實現(xiàn)對高低滲透層原油的等流度驅替,避免水驅或氣驅應用中常見的指進現(xiàn)象,從而提高原油采收率[5-6]。

泡沫滲流特性和驅油效果受氣液流速、溶液性質和孔隙結構等參數(shù)的影響,研究者對此進行了大量的研究,得出了泡沫流動壓降或驅油效率與各因素之間的趨勢關系,然而不同研究者得出的結論有一致性也存在矛盾[7-9]。泡沫結構(即泡沫體系中氣泡的尺寸大小和分布)[10]是影響泡沫滲流特性和驅油效果的一個重要參數(shù),尺寸大的氣泡會阻塞大孔隙,降低氣體的滲透率。所以,研究注入速率、氣液比、溶液性質和多孔介質結構等因素對生成氣泡大小的影響具有重要意義。泡沫在多孔介質中主要通過3種作用機理生成,即氣泡卡斷、分叉和液膜滯后[11-12]。其中卡斷作用的過程是氣液接觸,在喉道處截斷形成不連續(xù)的氣泡,該作用容易形成強泡,在較高的注入速度下,是泡沫產(chǎn)生的主要機理。

本文利用Ansys 14.0中的流體計算軟件FLUENT來模擬研究氣泡在T型微通道內的卡斷形成過程,采用多相流流體體積函數(shù)法(volume of fluid,VOF)[13]捕捉氣液界面。通過對T型微通道中氣泡形成過程中的形變和受力分析,研究氣液流速、流體性質和微通道尺寸等因素對氣泡生成大小的影響機制。

1 計算模型

1.1 物理模型

本文選用T型微通道作為物理模型,研究氣泡在微通道中的生成大小,如圖1所示。T型微通道由氣液混合區(qū)和氣泡流動區(qū)兩部分組成。氣相和液相分別從上下兩個入口進入混合區(qū)域,形成氣泡后進入流動區(qū)域。T型微通道入口和出口的直徑相同,用d來表示,氣液混合區(qū)域長度為6d,氣泡流動區(qū)域長度為50d。本文中,結合泡沫驅油實際應用地層孔隙的大小以及滲流速度,微通道直徑d選取0.05 mm、 0.1 mm和0.2 mm,氣液流速范圍為0.000 6~0.02 m·s-1。流體物性采用標準狀況下的物性,不考慮能量傳遞。

圖1 T型微通道二維物理模型Fig.1 Schematic of T-Junction microchannel

1.2 數(shù)學模型

1.2.1 控制方程

1.2.1.1 連續(xù)方程

氣液兩相在流動過程中滿足質量守恒定律,在微尺度條件下,忽略重力,連續(xù)性方程可簡化為:

(1)

1.2.1.2 動量方程

(2)

式中,u=(u,v)為流體速度,m·s-1;ρ為流體密度,kg·m-3;μ為流體動力粘性系數(shù),Pa·s;F為由集中在相界面上的表面張力轉化為的體積力項,稱為表面體積力。在對N-S方程求解過程中,如何處理表面張力是一個關鍵問題。Brackbill等[14]假設界面有一定厚度,將界面上的表面張力轉化為體積力(即表面體積力),解決了表面張力的離散問題,由此建立的模型稱為連續(xù)表面張力模型(CSF)。本文采用CSF模擬表面張力,通過添加附加分布函數(shù)將表面張力拓展到整個計算區(qū)域。上式(2)中的表面體積力F由下式計算:

(3)

(4)

1.2.1.3 相函數(shù)輸運方程

研究氣液兩相流動過程的關鍵是如何有效地捕捉兩相界面。本文選取VOF方法捕捉氣液相界面,在VOF方法中定義流體體積函數(shù),其輸運方程如下:

(5)

在VOF方法中,氣液相界面的捕捉是通過計算各網(wǎng)格單元內的兩相體積分數(shù)αl和αg來實現(xiàn)的。αl=1(αg=0),表示網(wǎng)格單元全部被液相占據(jù);αl=0(αg=1),表示網(wǎng)格單元全部被氣相占據(jù);氣液相界面存在于0<αl<1的網(wǎng)格單元內。在氣液兩相混合的網(wǎng)格單元中,公式(1)、(2)中的密度和粘度皆為混合屬性,通過下式計算:

ρ=α1ρ1+(1-αg)ρg,

(6)

μ=α1μ1+(1-αg)μg。

(7)

氣液相界面采用分段線性界面算法(PLIC)重構。Youngs[15]提出的PLIC型重構技術是一種完全的重構技術,在單個網(wǎng)格內用直線段近似界面,它通過界面附近網(wǎng)格內流體的細微輸運,幾何地精確給定界面的位置變化。

1.2.2 主要假設及邊界條件

該數(shù)值模擬主要假設為:

(1)忽略重力對孔隙內流體流動的影響;

(2)兩相均不可壓縮。

模型邊界條件為:

入口1:速度入口velocity inlet,液相,速度0.000 6~0.02 m·s-1;

入口2:速度入口velocity inlet,氣相,N2,速度0.005~0.02 m·s-1;

出口:充分發(fā)展出口outflow。

1.2.3 計算設置

在動量方程的求解中,壓力差值方案選用PRESTO!(pressure staggering option)算法,壓力-速度耦合采用PISO(pressure-implicit with splitting of operators)算法,動量方程采用二階迎風格式。Courant數(shù)、時間步長和亞松弛迭代因子等參數(shù)在計算過程中根據(jù)計算結果的穩(wěn)定性以及收斂性選取。

2 模型驗證

網(wǎng)格的劃分對于模擬計算至關重要,網(wǎng)格質量的好壞直接影響計算的準確性,因此,對于網(wǎng)格數(shù)的確定需要反復計算調試和比較。在相同的條件下對網(wǎng)格做加密處理,如對數(shù)值計算基本沒有影響,則判定網(wǎng)格數(shù)的劃分是合理的。利用Gambit軟件進行物理建模,劃分網(wǎng)格。本文T型微通道二維模型采用四邊形網(wǎng)格,圖2所示為直徑d= 0.1 mm 的T型微通道內,相同工況下不同網(wǎng)格數(shù)時形成的氣泡輪廓圖,深色表示液相,淺色表示氣相。為了模擬表面活性劑溶液起泡,氣液表面張力設為0.032 N·m-1,氣體選為N2,由于添加表面活性劑對液相粘度影響不大,所以液相粘度設為水的粘度。由圖可以看出在考察的網(wǎng)格密度下,流動情況相似,主要的不同在于氣液界面。模型網(wǎng)格越密,氣液界面越清晰。鑒于計算量和界面的清晰度,研究中選用的網(wǎng)格數(shù)為7 168。

d=0.1 mm, UG =UL =0.02 m·s-1a 網(wǎng)格數(shù)3584; b 網(wǎng)格數(shù)7168; c 網(wǎng)格數(shù)14336圖2 網(wǎng)格獨立性驗證Fig.2 Verification of grid independence

Laplace壓力是指曲面所受的內外壓力差,是由氣/液界面的表面張力引起的,計算公式為:

(8)

圖3所示為d=0.05 mm,0.1 mm和0.2 mm 3種不同直徑二維微通道中氣液界面所受的Laplace壓力,直線所示為由式(8)計算的理論值,散點為由數(shù)值計算得到的結果。數(shù)值計算情況下,對應每一個尺寸的微通道,隨機選取了兩個氣泡,每個氣泡前后界面的Laplace壓力平均值對應圖3中的一個數(shù)據(jù)點。從圖中可以看出,數(shù)值計算結果和理論計算值基本吻合,由于數(shù)值計算不能完全準確地模擬氣液界面形狀,所以存在一定誤差。

圖3 數(shù)值模擬和理論計算的Laplace壓力對比Fig.3 Laplace pressure comparison between numerical simulation and theoretical computation results

3 結果與討論

3.1 氣泡形成過程

氣泡在T型微通道內形成過程如圖4(1)所示,分為三個階段[16-17]:(1)阻塞階段(a~c),隨著流體流入,氣體量增多,氣體前端占據(jù)水平通道入口,阻塞液體流入;(2)塌陷階段(c~g),由于氣體堵塞水平通道,氣泡后端的液體向氣體通道內流動,氣泡被逐漸拉長,在后端出現(xiàn)頸部,本文定義頸部寬度Wn為拐點A到氣液界面的垂直距離,如圖4(2)所示,隨著氣泡繼續(xù)向下游運動,頸部逐漸變細;(3)卡斷階段(h),氣泡在頸部斷裂,一部分氣體收縮到下部氣體通道內,另一部分氣體則收縮到右側水平通道內,在水平微通道內形成兩端界面呈弧形的氣泡。

(1) d=0.1 mm,UG=UL=0.01 m·s-1, σ=0.032 N·m-1;(2) 頸部寬度Wn圖4 T型微通道中氣泡形成過程及頸部寬度定義Fig.4 Bubble formation process in T-Junction microchannel and the definition of neck width

3.2 氣液流速的影響

圖5給出了不同氣體份額下T型微通道內生成氣泡的大小。固定氣相流速UG=0.01 m·s-1,隨著液相流速的降低(即氣體份額的增加),T型微通道內生成氣泡的長度增大。

d=0.1 mm,UG=0.01 m·s-1, σ=0.032 N·m-1圖5 不同氣體份額下生成氣泡的大小Fig.5 Bubble lengths at different gas fractions

圖6所示為不同氣相流速下,氣體份額εG(εG≈UG/(UL+UG))對T型微通道內生成氣泡長度的影響。由圖可以看出氣泡長度隨氣體份額的增加呈指數(shù)增加趨勢。相同氣體份額下,氣液流速對生成氣泡長度影響不大,氣泡長度隨流速增大呈減小趨勢。

d=0.1 mm, σ=0.032 N·m-1圖6 不同流速下生成的氣泡長度Fig.6 Bubble lengths at different flow velocities

從氣泡的生成過程分析氣液流速對生成氣泡大小的影響。定義塌陷速率Vc,Vc=ΔWn/Δt為塌陷階段頸部寬度的變化速率,塌陷時間tc為頸部寬度由最大值變?yōu)?(即斷裂時)的時間。圖7所示為不同流速下氣體份額對氣泡生成過程中氣泡后端頸部塌陷速率和塌陷時間的影響。由圖可以看出,固定氣相流速下,隨著氣體份額的增大(即液相流速的減小),氣泡后端頸部塌陷速率減低,塌陷時間增大;相同氣體份額情況下,隨著流速的增加,氣泡后端頸部塌陷速率增大,塌陷時間減小。氣泡后端頸部塌陷速率越低,塌陷時間越大,則氣相阻塞水平通道向前推進的距離越大,從而生成的氣泡越大。當氣體份額增大,即液相流速減小,液相內由于擠壓形成的壓力小,則阻礙氣相前緣向Y軸推進的速度慢,所以氣相形成的最大頸部寬度較大。并且由于液相內部形成的壓力小,則在塌陷階段擠壓頸部斷裂的速度慢。因而,固定氣相流速下,隨著氣體份額的增大(即液相流速的減小),氣泡后端頸部塌陷速率減低,塌陷時間增大,生成的氣泡增大。

(1) 塌陷速率 (2) 塌陷時間圖7 流速對氣泡形成過程中塌陷速率和塌陷時間的影響Fig.7 Impact of flow velocity on collapse rate and collapse time in bubble formation process

3.3 液體粘度的影響

為了增強泡沫的穩(wěn)定性,實際應用中,常在起泡劑溶液中添加泡沫穩(wěn)定劑,以增加泡沫液相的粘度。圖8所示為液相粘度對泡沫生成氣泡長度的影響,液相粘度μ范圍為0.001~0.01 Pa·s。氣液入口流速為0.02 m·s-1,氣液界面張力為0.032 N·m-1,液相接觸角為0°。由圖可以看出,液相粘度對生成氣泡大小影響較小,氣泡長度隨液相粘度增大呈減小趨勢。T型微通道內氣泡在形成過程中主要受3種力,即由于氣體阻塞引起的液體擠壓力、粘性剪切應力和表面張力,液相粘度主要影響粘性剪切應力。而在氣泡的形成過程中,粘性剪切應力比液體的擠壓力和表面張力小1~2個數(shù)量級,所以流體粘度對生成氣泡大小的影響相對較小。粘性剪切應力主要在氣泡生成過程中的阻塞階段起作用,隨著液相粘度增大,氣相很快占據(jù)水平通道阻塞液相流入,所以氣相在微通道Y方向的擴張停止,形成的最大頸部寬度Wnmax相對減小,所以塌陷階段的塌陷時間減小,最終生成的氣泡長度減小。

UG=UL=0.02 m·s-1, σ=0.032 N·m-1圖8 液相粘度對生成氣泡長度的影響Fig.8 Impact of liquid viscosity on bubble length

3.4 表面張力的影響

泡沫驅油過程中通過向液相內加入表面活性劑來降低溶液的表面張力,從而使氣液易于形成泡沫。界面張力在泡沫形成和微通道流動中起著重要的作用,所以本節(jié)研究氣液界面張力對T型微通道內形成氣泡大小的影響。由于泡沫驅油中加入的表面活性劑降低了溶液表面張力,所以主要考察溶液表面張力小于水的表面張力(σ=0.072 N·m-1)的情況。氣液入口流速為0.02 m·s-1,液相粘度0.001 Pa·s,液相接觸角為0°。圖9為表面張力對生成氣泡長度的影響。隨著表面張力的增大,氣泡長度呈線性增大的趨勢。當表面張力從0.072 N·m-1降低到0.01 N·m-1時,T型微通道內氣泡的長度減小了18%。

d=0.1 mm,UG=UL=0.02 m·s-1圖9 表面張力對生成氣泡長度的影響Fig.9 Impact of surface tension on bubble length

由公式ΔpS=2σ/d可知,氣泡在形成過程中的阻塞階段,隨著界面張力的增加,界面氣液兩側的壓力差增大,此壓力差促使氣液界面向微通道Y軸方向移動,所以在阻塞階段,達到的最大頸部寬度Wnmax增大。從氣泡生成過程分析,表面張力對塌陷時間tc和塌陷速率Vc的影響如圖10所示,隨著表面張力的增加,塌陷時間延長,塌陷速率呈先增大后減小的趨勢。結合最大頸部寬度的變化趨勢,表面張力的增加使生成的氣泡長度變大。

圖10 表面張力對氣泡生成過程中塌陷速率和塌陷時間的影響Fig.10 Impact of surface tension on collapse rate and collapse time in bubble formation process

3.5 微通道直徑的影響

泡沫在多孔介質中的生成受孔隙大小的影響。本節(jié)研究了T型微通道直徑對氣泡生成大小的影響。氣液入口流速為0.02 m·s-1,氣液界面張力為0.032 N·m-1,液相粘度0.001 Pa·s,液相壁面接觸角為0°。圖11描述了3種不同微通道直徑下,無量綱氣泡長度LG/d與氣體份額εG的關系。由圖可以看出,在微通道直徑0.05 ~0.2 mm范圍內,氣泡的無量綱長度基本相同,不受微通道直徑的影響。氣泡無量綱長度隨氣體份額的增加,呈指數(shù)增大趨勢。在氣體份額0.5~0.9范圍內,對應的氣泡無量綱長度為4~26。

圖11 不同微通道直徑下無量綱氣泡長度隨氣體份額的變化Fig.11 Variation of dimensionless bubble length with gas fraction in the microchannel with different diameters

4 結論

本文利用數(shù)值計算的方法研究了氣液流速、流體性質等因素對T型微通道內氣泡生成大小的影響,由分析可知,氣體份額和微通道尺寸是影響氣泡生成大小的主要因素,液相粘度和表面張力對生成氣泡大小具有一定影響,具體結論如下:

(1)固定氣速下,氣泡長度隨氣體份額的增加呈指數(shù)增加趨勢;相同氣體份額下,氣液流速對生成氣泡長度影響不大,氣泡長度隨流速增大呈減小趨勢。

(2)液相粘度主要影響粘性剪切應力,在氣泡的形成過程中,粘性剪切應力比液體的擠壓力和表面張力小1~2個數(shù)量級,所以液體粘度對生成氣泡大小的影響相對較小,氣泡長度隨液相粘度增大呈減小趨勢。

(3)當液相表面張力從0.072 N·m-1降低到0.01 N·m-1時,生成氣泡的長度減小了18%,這是因為氣泡在形成過程中的阻塞階段,達到的最大頸部寬度和塌陷時間減小了。

(4)氣泡長度隨微通道直徑的增加而增大,而氣泡的無量綱長度基本不受其影響。在微通道直徑0.05 mm~0.2 mm,氣體份額0.5~0.9范圍內,對應的無量綱氣泡長度為4~26。

[1]孫曉,王樹眾,白玉.交聯(lián)氮氣泡沫壓裂液的攜砂性能研究[J].工程熱物理學報,2011,32(1):67-70.

[2]劉祖鵬,李兆敏.CO2驅油泡沫防氣竄技術實驗研究[J].西南石油大學學報(自然科學版),2015,37(5):117-122.

[3]呂明明,王樹眾.二氧化碳泡沫穩(wěn)定性及聚合物對其泡沫性能的影響[J],化工學報,2014,65(6):2219-2224.

[4]LI S Y, LI Z M , LI B F. Experimental study on foamed gel flow in porous media[J]. Journal of Porous Media, 2015, 18(5):519-536.

[5]TALEBIAN S H, MASOUDI R, TAN I M, et al. Foam assisted CO2-EOR: A review of concept, challenges, and future prospects [J]. Journal of Petroleum Science and Engineering, 2014, 120(8):202-215.

[6]FARAJZADEH R, ANDRIANOW A, BRUINING H, et al. Comparative study of CO2and N2foams in porous media at low and high pressure and temperatures[J]. Industrial and Engineering Chemistry Research, 2009, 48(9): 4542-4552.

[7]ASGHARI K, KHALIL F. Operation parameters on CO2-foam process[J]. Petroleum Science and Technology, 2005, 23(2):189-198.

[8]SHAN D, ROSSEN W R. Optimal injection strategies for foa m IOR[J]. SPE Journal, 2004, 9(2):132-150.

[9]KIM J S, DONG Y, ROSSEN W R. Steady-state flow behavior of CO2foam[J]. SPE Journal, 2005,10(4): 405-415.

[10]ETTINGER R A, RADKE C J. Influence of texture on steady foam flow in Berea sandstone[J]. SPE Reservoir Engineering, 1992, 7(1): 83-90

[11]KOVSCEK A R, TANG G Q, RADKE C J. Verification of roof snap off as a foam-generation mechanism in porous media at steady state[J]. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 2007, 302(1/2/3):251-260

[12]GAUTEPLASSA J, CHAUDHARYB K, KOVSCEK A R, et al. Pore-level foam generation and flow for mobility control in fractured systems[J]. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 2015, 468:184-192.

[13]HIRT C W, NICHOLS B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J]. Journal of Computational Physic, 1981, 39(1):201-225.

[14]BRACKBILL J U, KOTHE D B, ZEMACH C. A continuum method for modeling surface tension[J]. Journal of Computational Physics, 1992, 100(2):335-354.

[15]YOUNGS D L. Time-dependent multi-material flow with large fluid distortion[M]∥Numerical methods for fluid dynamics. New York: Academic Press, 1982: 273-285.

[16]DAI L, CAI W F, XIN F. Numerical study on bubble formation of a gas-liquid flow in a T-junction microchannel[J]. Chemical Engineering & Technology, 2009, 32(12): 1984-1991.

[17]GARSTECKI P, FUERSTMAN M J, STONE H A, et al. Formation of droplets and bubbles in a microfluidic T-junction-scaling and mechanism of break-up[J]. Lab on A Chip, 2006, 6(3): 437-446.

Influential factors of bubble length in a T-junction microchannel

Lü Ming-ming1, LIU Zhi-gang1*, GUAN Ning1, WANG Shu-zhong2

(1.Key Lab for Flow & Enhanced Heat Transfer of Shandong Academy of Sciences, Energy Research Institute,Shandong Academy of Sciences, Jinan 250014, China; 2. School of Energy and Power Engineering,Xi′an Jiaotong University, Xi′an 710049, China)

∶We performed numerical simulation for bubble formation process in a T-junction microchannel with volume of fluid method (VOF). We also analyzed the impact of gas/liquid fluid velocity, fluid properties and microchannel diameter on bubble length based on bubble formation mechanism. Results demonstrate that the bubble length exponentially increases with the increase of gas fraction in T-junction microchannel, but gas/liquid fluid velocity has little effect on bubble length for fixed gas fraction. Liquid viscosity and surface tension comparatively have less effect on bubble length. When surface tension of liquid phase reduces from 0.072 N·m-1to 0.01 N·m-1, the bubble length in T-junction microchannel decreases by 18%. This is because that maximum neck width and collapse time decrease in expansion stage of bubble formation process. The bubble length increases with the increase of microchannel diameter, but dimensionless length of the bubble is less influenced by microchannel diameter.

∶ microchannel; numerical simulation; bubble; snap-off; surface tension

10.3976/j.issn.1002-4026.2016.05.014

2016-07-25

呂明明(1987—),女,助理研究員,研究方向為微通道內多相流動研究。

*通信作者。Email:zgliu9322@gmail.com

TK124

A

猜你喜歡
界面
聲波在海底界面反射系數(shù)仿真計算分析
微重力下兩相控溫型儲液器內氣液界面仿真分析
國企黨委前置研究的“四個界面”
當代陜西(2020年13期)2020-08-24 08:22:02
基于FANUC PICTURE的虛擬軸坐標顯示界面開發(fā)方法研究
西門子Easy Screen對倒棱機床界面二次開發(fā)
空間界面
金秋(2017年4期)2017-06-07 08:22:16
鐵電隧道結界面效應與界面調控
電子顯微打開材料界面世界之門
人機交互界面發(fā)展趨勢研究
手機界面中圖形符號的發(fā)展趨向
新聞傳播(2015年11期)2015-07-18 11:15:04
主站蜘蛛池模板: 欧美区一区| 一区二区三区四区精品视频 | 欧美激情伊人| 国产亚洲视频在线观看| 高清无码手机在线观看| 欧美一区二区三区不卡免费| 国产成人啪视频一区二区三区| 精品在线免费播放| 亚洲成人播放| 亚洲免费毛片| 日韩欧美综合在线制服| 伊人无码视屏| 狠狠v日韩v欧美v| 三区在线视频| 国产理论一区| 婷婷色婷婷| 欧美国产精品拍自| 多人乱p欧美在线观看| 亚洲人成网站观看在线观看| 69视频国产| 久久青草视频| 久久情精品国产品免费| 亚洲欧美国产五月天综合| 色综合天天娱乐综合网| 亚洲Va中文字幕久久一区 | 男人天堂亚洲天堂| 无码国内精品人妻少妇蜜桃视频| 国产亚洲精久久久久久无码AV| 99人体免费视频| 青青青草国产| 激情在线网| a毛片基地免费大全| 亚洲天堂精品在线| 一级毛片在线免费看| 老汉色老汉首页a亚洲| 亚洲不卡av中文在线| 97视频在线观看免费视频| 亚洲日本一本dvd高清| 国产亚洲欧美另类一区二区| 免费a在线观看播放| 香蕉伊思人视频| 亚洲 欧美 日韩综合一区| 国产在线视频福利资源站| 92精品国产自产在线观看| 在线观看精品国产入口| 欧美高清国产| 日本亚洲欧美在线| 婷婷开心中文字幕| 在线免费观看AV| 国产99热| 午夜啪啪福利| 久操中文在线| 久久一色本道亚洲| 国产理论精品| 日韩在线1| 国产一在线| 日韩a级片视频| 国产精品欧美激情| 亚洲天堂色色人体| 亚洲欧美国产五月天综合| 久久精品波多野结衣| 国产丝袜第一页| 久久大香香蕉国产免费网站| 香蕉视频在线精品| 国产成人亚洲综合A∨在线播放| 尤物亚洲最大AV无码网站| av无码一区二区三区在线| 青青网在线国产| 国产呦视频免费视频在线观看| 午夜激情福利视频| 激情亚洲天堂| 老汉色老汉首页a亚洲| 亚洲乱码视频| 色屁屁一区二区三区视频国产| 精品一区二区久久久久网站| 国产精品污污在线观看网站| 国产精品白浆在线播放| 欧美日韩免费观看| 久久亚洲国产一区二区| 亚洲黄色激情网站| 国产熟睡乱子伦视频网站| 日韩小视频在线观看|