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

小尺度受限空間內瓦斯湍流爆燃大渦模擬

2016-08-22 02:45:10溫小萍余明高鄧浩鑫陳俊杰王發輝劉志超
化工學報 2016年5期
關鍵詞:模態實驗

溫小萍,余明高,鄧浩鑫,陳俊杰,王發輝,劉志超

(1河南理工大學機械與動力工程學院,河南 焦作 454003;2河南理工大學安全科學與工程學院,河南 焦作 454003)

?

小尺度受限空間內瓦斯湍流爆燃大渦模擬

溫小萍1,余明高2,鄧浩鑫1,陳俊杰1,王發輝1,劉志超1

(1河南理工大學機械與動力工程學院,河南 焦作 454003;2河南理工大學安全科學與工程學院,河南 焦作 454003)

構建了150 mm × 150 mm × 500 mm小尺度受限空間三維模型,基于火焰表面密度模型和Charlette湍流燃燒模型,對兩側連續障礙物條件下瓦斯爆燃火焰與湍流耦合過程進行了大渦模擬(LES)。模擬結果均與實驗結果進行了比較。結果表明:大渦模擬可以很好預測瓦斯爆燃過程中的火焰結構、火焰鋒面位置、火焰傳播速度及超壓,驗證了大渦模擬及湍流燃燒模型對于瓦斯爆燃的適用性。此外,通過Karlovitz數定量描述了瓦斯爆燃火焰與湍流之間的相互作用及其變化規律,并對不同時刻的火焰模態進行了判別,在兩側連續障礙物條件下瓦斯湍流爆燃火焰先后經歷波紋小火焰和薄反應區兩種模態。

爆燃;湍流;模型;數值模擬;超壓;火焰

DOI:10.11949/j.issn.0438-1157.20151219

Lee等[5]通過實驗研究指出,大尺度障礙物形成的湍流是爆燃火焰加速機制的最重要因素,并將這種在大尺度湍流激勵下的爆燃過程稱為湍流爆燃。國內外不少學者研究過瓦斯爆燃火焰加速機理,并進行了相關實驗[6-10]或數值模擬[11-18]的研究,但由于現有實驗方法和測試手段的限制,均未能很好捕捉詳實的火焰傳播過程和湍流流場。特別在數值模擬方面,對于瓦斯爆燃過程的研究大都采用傳統的雷諾平均模擬(RANS)方法,其主要缺點在于它只能提供湍流的平均信息,這對于研究瓦斯爆燃的復雜湍流過程遠遠不夠,因此無法有效預測瓦斯爆燃火焰-湍流耦合的行為細節。

大渦模擬(LES)將耗散尺度的脈動進行過濾,只對大尺度脈動進行求解,與雷諾平均模擬相比計算量更大,但其空間分辨率更高,因此可以有效捕捉瓦斯爆燃過程中的湍流特征。目前國內外關于瓦斯爆燃的有效大渦模擬仍然較少,未能揭示瓦斯爆燃火焰與湍流的耦合特性,更缺少這方面的實驗佐證。本文對兩側連續障礙物條件下瓦斯爆燃火焰與湍流耦合過程進行大渦模擬,將數值模擬得到的火焰形狀、火焰鋒面位置、火焰傳播速度及超壓等數據與筆者實驗結果[19]進行比較,進而揭示瓦斯爆燃過程中火焰-湍流耦合機理。

1 數學模型

1.1大渦模擬控制方程

實驗[19]中涉及小尺度受限空間內的瓦斯爆燃氣體流速不超過100 m·s-1,即Ma<0.3,故能量方程中的壓力梯度項可以忽略。但由于反應區的溫度變化較大,爆燃氣流應視為可壓縮理想氣體。一般情況下,可壓縮湍流的大渦模擬控制方程可以通過N-S及其他方程的過濾導出,然而直接過濾可壓縮N-S方程將得到十分復雜的可解尺度湍流的運動方程,因此本文采用Favre過濾方法,得到封閉的可壓縮湍流大渦模擬方程

式中,ρ為密度;p為壓力;T為溫度;t為時間;i、j為坐標方向;ui是i方向速度;xi為直角坐標參量;σij為黏性應力張量;h為焓;λ為熱導率;c為反應進程變量;Tω.為化學反應熱;cω.為歸一化的化學反應速率;D為擴散系數;上標“—”及“~”分別表示物理空間過濾量和Favre過濾量;亞網格應力、亞網格物質流、亞網格熱流和濾波后的反應速率采用亞網格湍流模型和亞網格燃燒模型進行封閉,其中亞網格Prandtl數PrSGS和亞網格Schmidt數ScSGS均取0.7[11]。

1.2亞網格燃燒模型

過濾后的反應進程變量分子擴散及燃燒速率可寫為

式中,ρu為未燃混合氣體密度;Sl為層流火焰速度。將式(8)、式(11)代入式(4)中可得到

為了考慮火焰鋒面與亞網格湍流之間的相互作用,Charlette等[20]給出了湍流燃燒模型,即式(12)中亞網格火焰褶皺系數ΞΔ表達式為

式中,u'Δ為亞網格湍流脈動速度;Δ為過濾尺度;亞網格湍流Reynolds數ReΔ=u'ΔΔ/ν;層流火焰厚度δf通過δfSl/ν=4估算得到,ν為運動黏度。

2 數值模擬

圖1 實驗裝置及物理模型(單位:mm)Fig.1 Experimental setup and physical model (unit: mm)

為了提高管道出口邊界條件的準確性,減小出口回流對管道內部壓力的影響,物理模型在管道出口往x、y及z方向均延伸至1000 mm,即增加了一個尺寸為1000 mm × 1000 mm × 500 mm的擴展區域。這個擴展區域與原管道之間相通,在模擬過程中允許壓力波在擴展區域繼續傳播,因此可以更加真實地反映管道出口條件(擴展區出口設為壓力出口)。計算網格采用六面體非均勻網格,即在障礙物和點火源的附近選擇較細的網格,而在遠離障礙物和點火源的區域(特別是擴展區域)采用較粗的網格,其三維計算網格如圖2所示。

圖2 三維計算網格(單位:m)Fig.2 3D computational grids (unit: m)

在大渦模擬中,亞網格尺度的物理擴散是隨著網格的加密而減小的,因此,對于大渦模擬而言,并不存在通常數值計算中所謂“網格無關性”的概念,較好的做法是采用十分精細的網格和很小的時間步長,以確保其數值擴散遠小于雷諾平均方法。本文在障礙物和點火源附近區域的網格大小為δx = δy = δz = 1.25 mm,遠離障礙物和點火源的區域的網格逐漸增大。計算單元總數約為180萬,其中管道區的單元數約為150萬,約占計算單元總數的83.3%。數值求解是利用CFD求解器ANSYS FLUENT實現的,整個計算區域采用有限容積法離散控制方程,空間離散采用二階中心差分格式,選取二階隱式時間推進法來提高計算精度,利用SIMPLE算法對速度和壓力耦合方程組進行解耦。質量方程、動量方程、能量方程、反應進程變量方程殘差分別小于1×10-6、2.5×10-5、1×10-6及1×10-5。

3 結果與分析

3.1火焰結構

圖3為采用大渦模擬和Charlette湍流燃燒模型時不同時刻火焰結構與實驗結果的比較,其中模擬結果中的火焰結構采用反應進程變量c=0.5等值面表示[14]。可以清楚地看到,Charlette湍流燃燒模型不僅很好地預測了火焰結構演變過程(包括火焰結構和火焰鋒面達到各個位置的時間),而且準確捕捉到了瓦斯湍流爆燃火焰的結構特征。在爆燃初期(如t=20 ms),火焰結構主要為半球形,火焰傳播速度較慢,致使火焰鋒面到達第1個障礙物的時間約為25 ms,約占用火焰在整個管道內傳播時間的59%。至t=32 ms時,火焰鋒面從第1對障礙物之間的間隙噴出,并逐漸向兩側卷曲,說明火焰明顯受到兩側障礙物后的湍流影響。火焰繞過第1對障礙物之后,繼續向第2對和第3對障礙物傳播。至t=42 ms時,火焰鋒面不斷發生褶皺,出現了明顯的火焰湍流現象,這種現象是由于受到多個障礙物持續作用后,火焰被不斷褶皺、變形,導致未燃氣體與已燃氣體能夠快速相互摻混,由此出現明顯的火焰湍流特征。在整個過程中,大渦模擬和實驗結果非常相似。

圖3 大渦模擬的不同時刻火焰結構與實驗結果的比較Fig.3 Comparison between LES flame structures at different instants and experimental data

3.2火焰鋒面位置與火焰傳播速度

圖4、圖5分別表示利用Charlette燃燒模型模擬的不同時刻火焰鋒面位置及火焰鋒面傳播速度與實驗結果的對比圖。其中,火焰鋒面位置是指每張火焰圖片中火焰鋒面與點火點的最大距離,火焰速度則是通過前后連續兩個時間點的火焰鋒面位置差與時間間隔之比計算而得。從圖4可以看出,模擬與實驗在爆燃初期的火焰鋒面位置存有偏差,實驗數值略高于模擬結果,這是由于實驗中高速攝像的角度偏差造成的。但是,在瓦斯湍流爆燃的整個過程中,模擬預測的火焰鋒面位置與實驗吻合較好。圖5則反映了模擬的火焰傳播速度與實驗也非常接近,尤其是很好再現了火焰繞過障礙物時的加速、減速特征。稍有區別的是數值模擬中火焰加速程度略低于實驗數據,這是由于數值模擬的計算網格精度有限,使得大渦模擬不能捕捉到極小尺度的湍流渦團對火焰加速的激勵作用。顯然,這種影響對模擬的精確度影響較小。

圖4 大渦模擬的不同時刻火焰鋒面位置與實驗結果的比較Fig.4 Comparison between LES time histories of flame front position and experimental data

圖5 大渦模擬的不同時刻火焰鋒面位置與實驗結果的比較Fig.5 Comparison between LES simulated time histories offlame front position and experimental data

3.3爆燃超壓

圖6為大渦模擬得出的瓦斯爆燃超壓度與實驗結果的對比。從實驗結果來看,爆燃超壓在15 ms ≤t ≤ 25 ms時間段出現了一個小幅度峰值,該峰值是由于管道頂端薄膜破裂形成的,而數值模擬并未考慮薄膜的影響,因此模擬超壓并未出現這一小幅度峰值。還可以看到,數值模擬預測的最大超壓為10.1 kPa,與實驗最大超壓11.6 kPa的誤差為12.9%。這里值得注意的是,數值模擬和實驗的最大超壓的時間非常吻合,均出現在t=41 ms,再對照圖3中的火焰結構可以看出,此時火焰鋒面已經越過最后一對障礙物并迅速往兩側發展,火焰在障礙物后卷曲變形呈現蘑菇狀,附近的未燃氣體則以更快的速度摻混其中并燃燒,導致超壓的急速上升。在41 ms ≤t ≤ 42 ms時間段內,雖然第3對障礙物下游的火焰繼續膨脹,但管道底部的火焰在抵達壁面后開始出現局部燃盡熄滅現象,這可能使管道內火焰面密度減小,并導致超壓開始下降。此后,由于火焰鋒面已經越過容器出口,管道超壓迅速回落。

圖6 模擬的不同時刻爆燃超壓與實驗結果的比較Fig.6 Comparison between LES time histories of deflagration overpressure and experimental data

4 爆燃湍流與火焰模態分析

圖7為火焰鋒面分別經過3對障礙物不同時刻的速度矢量分布。當火焰經過第1對障礙物時(t=29 ms),3對障礙物之后均形成了一對左右對稱的大尺度旋渦,左側旋渦逆時針旋轉,而右側旋渦順時針旋轉。這是由于燃燒區的高溫氣流不斷膨脹,迫使下游未燃氣體向管道出口方向流動,而這些未燃氣體經過障礙物時,由于主流流速與障礙物后氣體流速存在較大的速度差,從而在障礙物后形成旋渦。這些旋渦的尺度大小與障礙物尺度、障礙物位置、流體流速、流體黏度有關。當火焰繼續經過第2對和第3對障礙物時,隨著氣流流速的增大,旋渦尺度和速度大小都在增加。

為了揭示爆燃過程中火焰與湍流的耦合關系,本文引入Karlovitz數(Ka)來表征湍流與燃燒相互作用的程度,從而定性分析瓦斯爆燃火焰在連續障礙物條件下的火焰模態。Ka定義為火焰面的時間尺度τf與流動中最小渦(Kolomogrov渦)的時間尺度τη之比

圖7 不同時刻速度矢量分布Fig. 7 Velocity vector maps when flame passes sequentially over three obstacles

對于預混湍流燃燒的大渦模擬,Pitsch等[22]提出Ka按式(15)近似計算

圖8表示3個不同時刻(29、37及41 ms)火焰面上Ka的空間分布,火焰面取反應進程變量c = 0.5的等值面。在這3個時刻,火焰鋒面正在分別繞過3個障礙物。由圖8可以看出,由于位于點火源附近的湍流強度很小,其火焰面比較光滑,Ka始終處于較低水平。但隨著火焰鋒面先后繞過3個障礙物時,位于障礙物附近的火焰鋒面上Ka出現明顯增長,其中35 ms時刻位于第3個障礙物附近火焰面上的Ka較高,最大值可達到10~11。其原因是在障礙物湍流激勵作用下,障礙物附近流體的應變率和亞網格黏度較高,使其亞網格湍流脈動處于較高水平造成的,根據式(15)可知,較高的湍流脈動所對應的Ka也就越大。因此,可進一步利用整個流場的Ka最大值(在障礙物附近)來反映火焰與湍流之間的作用程度。

圖8 火焰面(c=0.5)上Ka在不同時刻的空間分布Fig.8 Spatial distribution of Ka over flame surface (c=0.5) at different time instants

Pitsch等[22]將大渦模擬過濾尺度和Ka引入之后,對燃燒模態區間作了進一步改進,改進后燃燒模態區間的劃分如圖9所示。圖中的方形點表示瓦斯爆燃火焰在中間連續障礙物條件下不同時刻的Ka-Δ/δl數據點,對應時刻分別為2、14、34、37、41及42 ms,其中Ka均取對應時刻整個流場的最大值。由圖9可以看出,在瓦斯爆燃初期Ka略低于1,還處于波紋火焰面模式;隨著火焰向下游傳播,Ka不斷升高,火焰開始向薄反應區轉變;當t = 42 ms時(此時火焰鋒面抵達管道出口),Ka達到9.8,但仍然屬于薄反應區。由此可知,在連續障礙物條件下,瓦斯爆燃火焰先后經歷了波紋小火焰和薄反應區兩種火焰模態,但對于火焰-湍流耦合作用的時間段(34 ms < t < 42 ms)而言,則主要處于薄反應區火焰模態。

圖9 大渦模擬結果在湍流火焰模態區間上的分布Fig.9 Distribution of LES data on turbulent flame regimes

5 結 論

(1)利用大渦模擬及Charlette湍流燃燒模型對兩側障礙物條件下的瓦斯湍流爆燃過程進行數值模擬,并將模擬獲得的火焰結構、火焰鋒面位置、火焰速度及超壓等參數特征與實驗進行比較,驗證了大渦模擬及Charlette湍流燃燒模型對于瓦斯爆燃的適用性。

(2)通過Karlovitz數定量描述了瓦斯爆燃火焰與湍流之間的相互作用及其變化規律,對中間連續障礙物條件下不同時刻的火焰模態進行了判別。在兩側連續障礙物條件下,瓦斯湍流爆燃火焰先后經歷了波紋小火焰和薄反應區兩種模態。

References

[1] 周心權. 煤礦井下瓦斯爆炸的基本特性[J]. 中國煤炭, 2002, 28(9):8-11. DOI: 10.3969/j.issn.1006-530X.200209002. ZHOU X Q. Basic characteristics of gas explosion in coal mines [J]. China Coal, 2002, 28(9): 8-11.DOI: 10.3969/j.issn.1006-530X. 200209002.

[2] ZHU C J, LU Z G, LIN B Q, et al. Effect of variation in gas distribution on explosion propagation characteristics in coal mines [J]. Mining Science and Technology, 2010, 4(20): 516-519.

[3] JOHANSEN C T, GICCARELLI G. Visualization of the unburned gas flow field ahead of an accelerating flame in an obstructed square channel [J]. Combustion and Flame, 2009, 156(2): 405-416.

[4] KIM W K, MOGI T, DOBASHI R. Flame acceleration in unconfined hydrogen/air deflagrations using infrared photography [J]. Journal of Loss Prevention in the Process Industries, 2013, 26(6): 1501-1505.

[5] LEE J H S, KNYSTAUTAS R, FREIMAN A. High speed turbulent deflagrations and transition to detonation in H2/air mixtures [J]. Combustion and Flame, 1984, 56(2): 227-239.

[6] BYCHKOV V, VALIEV D, AKKERMAN V, et al. Gas compression moderates flame acceleration in deflagration-to-detonation transition[J]. Combustion Science and Technology, 2012, 184(7/8): 1066-1079.

[7] IBRAHIM S, HAEGRAVE G, WILLIAMS T. Experimental investigation of flame/solid interactions in turbulent premixed combustion [J]. Experimental Thermal and Fluid Science, 2001, 24:99-106.

[8] CICCARELLI G, JOHANSEN C, PARRAVANI M. The role of shock-flame interactions on flame acceleration in an obstacle laden channel [J]. Combustion and Flame, 2010, 157(11): 2125-2136.

[9] WEN X P, YU M G, JI W T, et al. Methane-air explosion characteristics with different obstacle configurations [J]. International Journal of Mining Science and Technology, 2015, 25(2): 213-218.

[10] 鄭立剛, 呂先舒, 鄭凱, 等. 點火源位置對甲烷-空氣爆燃超壓特征的影響[J]. 化工學報, 2015, 66(7): 2749-2756. DOI:10.11949/j.issn.0438-1157.20141789. ZHENG L G, Lü X S, ZHENG K, et al. Influence of ignition position on overpressure of premixed methane-air deflagration [J]. CIESC Journal, 2015, 66(7): 2749-2756. DOI: 10.11949/j.issn.0438-1157. 20141789.

[11] SARLI V, Benedetto A, RUSSO G, et al. Large eddy simulation and PIV measurements of unsteady premixed flames accelerated by obstacles [J]. Flow Turbulence and Combustion. 2009, 83: 227-250.

[12] SARLI V, Benedetto A. Effects of non-equidiffusion on unsteady propagation of hydrogen-enriched methane/air premixed flames [J]. International Journal of Hydrogen Energy, 2013, 38(18): 7510-7518.

[13] BI M S, DONG C, ZHOU Y. Numerical simulation of premixed methane-air deflagration in large L/D closed pipes [J]. Applied Thermal Engineering, 2012, 40: 337-342.

[14] SARLI V, Benedetto A, RUSSO G. Large eddy simulation of transient premixed flame-vortex interactions in gas explosions [J]. Chemical Engineering Science, 2012, 71: 539-551.

[15] JOHANSEN C, CICCARELLI G. Modeling the initial flameacceleration in an obstructed channel using large eddy simulation [J]. Journal of Loss Prevention in the Process Industries, 2013, 26(4):571-585.

[16] 朱傳杰, 林柏泉, 江丙友. 受限空間內爆燃波瞬態流速與超壓的耦合關系[J]. 燃燒科學與技術, 2012, 18(4): 326-330. ZHU C J, LIN B Q, JIANG B Y. Coupled relationship between gas velocity and peak overpressure of deflagration wave in confined space [J]. Journal of Combustion Science and Technology (China),2012, 18(4): 326-330.

[17] WEN X P, YU M G, LIU Z C, et al. Effects of cross-wise obstacle position on methane-air deflagration characteristics [J]. Journal of Loss Prevention in the Process Industries, 2013, 26(6): 1335-1340.

[18] XIAO H H, MAKAROV D, SUN J H, et al. Experimental and numerical investigation of premixed flame propagation with distorted tulip shape in a closed duct [J]. Combustion and Flame, 2012, 159(4):1523-1538.

[19] 溫小萍, 武建軍, 解茂昭. 瓦斯爆炸火焰結構與壓力波的耦合規律[J]. 化工學報, 2013, 64(10): 3871-3877. DOI: 10.3969/j.issn. 0438-1157.2013.10.052. WEN X P, WU J J, XIE M Z. Coupled laws between flame structure and pressure wave of gas explosion [J]. CIESC Journal, 2013, 64(10):3871-3877. DOI: 10.3969/j.issn.0438-1157.2013.10.052.

[20] CHARLETTE F, MENEVEAU C, VEYNANTE D. A power-law flame wrinkling model for LES of premixed turbulent combustion(Ⅰ): Non-dynamic formulation and initial tests [J]. Combustion and Flame, 2002, 131(1): 159-180.

[21] IBRAHIM S S, GUBBA S R, MALALASEKERA W. A parametric study on large eddy simulations of turbulent premixed flames [C]// The 8th Asia-Pacific Conference on Combustion. 2010: 87-94.

[22] PITSCH H, DUCHAMP L. Large-eddy simulation of premixed turbulent combustion using a level-set approach [J]. Proceedings of the Combustion Institute, 2002, 29(2): 2001-2008.

Large eddy simulation of gas turbulent deflagration in small-scale confined space

WEN Xiaoping1, YU Minggao2, DENG Haoxin1, CHEN Junjie1, WANG Fahui1, LIU Zhichao1
(1School of Mechanical and Power Engineering, Henan Polytechnic University, Jiaozuo 454003, Henan, China;
2School of Safety Science and Engineering, Henan Polytechnic University, Jiaozuo 454003, Henan, China)

A 3D model of small-scale confined space with an inner size of 150 mm × 150 mm × 500 mm was set up. Based on the flame surface density model and the turbulent combustion model by Charlette et al., a large wddy simulation (LES) had been carried out on the process of gas deflagration flame-turbulence interaction with continuous obstacles at two sides of the chamber. All numerical results have been compared to experimental data. It showed that the LES is capable to predict the flame structure, position, speed, and overpressure in the process of gas deflagration, and the applicability of the LES and turbulent combustion model on gas deflagration was verified. In addition, the interaction between gas deflagration and turbulence and the relationship were described quantitatively by the Karlovitz number, and the transient flame regimes also identified. Under condition of continuous obstacles at double sides of the chamber, the gas turbulent deflagration flame experienced in the subsequent states of corrugated flamelets zone and thin reaction zone.

deflagration; turbulent flow; model; numerical simulation; overpressure; flame

引 言

發生在煤礦井下的瓦斯爆炸絕大多數屬于爆燃[1]。當發生瓦斯爆燃時,火焰在傳播過程中往往會遇到各種障礙物,如礦車、采掘機電設備、液壓支架、管道、巷道風門等,火焰前方未燃氣體不斷受到熱膨脹波的壓縮推動作用,使得在障礙物下游流場形成大尺度渦團,其湍流尺度遠大于層流火焰厚度,湍流脈動將誘導緊隨而至的火焰產生明顯褶皺、卷吸,從而使火焰鋒面附近的已燃高溫氣體和未燃氣體快速摻混,燃燒速率和火焰表面積增大,進而使火焰傳播獲得加速,并伴有爆燃超壓的增加,引起更強烈的湍流[2-4]。

date: 2015-07-29.

WEN Xiaoping, wenxiaoping666@163.com

supported by the National Natural Science Foundation of China (51176021).

TD 712

A

0438—1157(2016)05—1837—07

2015-07-29收到初稿,2016-01-25收到修改稿。

聯系人及第一作者:溫小萍(1977—),男,博士,副教授。

國家自然科學基金項目(51176021);河南省教育廳科學技術研究重點項目(14A410007)。

猜你喜歡
模態實驗
記一次有趣的實驗
微型實驗里看“燃燒”
做個怪怪長實驗
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
車輛CAE分析中自由模態和約束模態的應用與對比
國內多模態教學研究回顧與展望
高速顫振模型設計中顫振主要模態的判斷
航空學報(2015年4期)2015-05-07 06:43:35
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 免费看美女自慰的网站| 欧美一区二区三区香蕉视| 91口爆吞精国产对白第三集| 蜜桃视频一区二区| 亚洲天堂网在线观看视频| 国产成人91精品免费网址在线| 97人人模人人爽人人喊小说| 1769国产精品免费视频| 欧美色99| 中文字幕人成乱码熟女免费| 亚洲国产中文精品va在线播放| 色婷婷亚洲十月十月色天| 丝袜国产一区| 国模在线视频一区二区三区| 国产不卡网| 亚洲精品中文字幕无乱码| 中文字幕第1页在线播| 国产成人综合日韩精品无码不卡| 日韩久草视频| 亚洲最黄视频| 福利在线不卡一区| 久久精品视频亚洲| 国产欧美视频在线| 国产精品短篇二区| 亚洲国产成人精品青青草原| 国产精品999在线| 久久99国产乱子伦精品免| 日韩黄色大片免费看| 日韩在线观看网站| 青青久视频| 精品国产香蕉伊思人在线| 97青青青国产在线播放| 国产精品久线在线观看| 鲁鲁鲁爽爽爽在线视频观看| 无码福利日韩神码福利片| 毛片视频网址| 国产91小视频在线观看| 啪啪国产视频| vvvv98国产成人综合青青| 亚洲高清资源| 狠狠综合久久久久综| 亚洲综合网在线观看| 亚洲天堂视频在线播放| 欧美亚洲第一页| 日韩av电影一区二区三区四区| 国产靠逼视频| 亚洲中文字幕无码爆乳| 久久久久久久97| 中文字幕久久波多野结衣| 国产精品99久久久| 亚洲永久色| 亚洲欧洲日本在线| 国产丰满大乳无码免费播放| 在线欧美国产| 欧美日韩精品综合在线一区| 国产不卡国语在线| 中文字幕 91| 亚洲女人在线| 国产jizzjizz视频| 免费一级无码在线网站| 国产在线观看成人91| 黄色免费在线网址| 中文字幕66页| 喷潮白浆直流在线播放| 2021天堂在线亚洲精品专区| 国产精品太粉嫩高中在线观看| 国产成人区在线观看视频| 色婷婷成人| 国产精品白浆在线播放| 国产高潮流白浆视频| 欧美黄网在线| 香蕉伊思人视频| 中文字幕啪啪| 欧美精品成人| 亚洲欧美日本国产综合在线| 中文字幕一区二区视频| 91在线播放国产| 亚洲婷婷丁香| 免费啪啪网址| 亚洲成A人V欧美综合天堂| 综合亚洲色图| 久久99国产视频|