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

障礙物對預混火焰特性影響的大渦數(shù)值模擬*

2017-04-10 13:20:35王公忠張建華李登科陳先鋒
爆炸與沖擊 2017年1期
關鍵詞:實驗模型

王公忠,張建華,李登科,陳先鋒

(1.河南工程學院安全工程學院,河南鄭州451191;2.武漢理工大學安全科學與工程系,湖北武漢430070)

障礙物對預混火焰特性影響的大渦數(shù)值模擬*

王公忠1,2,張建華2,李登科2,陳先鋒2

(1.河南工程學院安全工程學院,河南鄭州451191;2.武漢理工大學安全科學與工程系,湖北武漢430070)

障礙物在預混氣體火焰?zhèn)鞑ミ^程中對其結構及傳播特性造成較大影響,對火焰的加速和爆燃轉(zhuǎn)爆轟過程(deflagration-to-detonation transition,DDT)起到直接的促進作用。通過障礙物條件下可視管道中甲烷/空氣預混火焰?zhèn)鞑嶒灒东@其火焰微觀結構變化。采用三維物理模型,采用壁面自適應局部渦黏模型(wall-adapting local eddy-viscosity,WALE)的大渦模擬(large eddy simulation,LES),并用火焰增厚化學反應模型(thickened flame model,TFM)對實驗過程進行重現(xiàn)。分析開口管道中預混火焰翻越障礙物后的復雜流場變化,并分析層流向湍流轉(zhuǎn)變過程的特點。揭示了在障礙物影響下預混火焰擾動失穩(wěn)現(xiàn)象的直接原因,是由障礙物引發(fā)的3個氣流渦團同時作用而形成Kelvin-Helmholtz不穩(wěn)定及Rayleigh-Taylor不穩(wěn)定現(xiàn)象耦合作用所導致。

爆炸力學;火焰增厚模型;大渦模型;預混火焰;障礙物

在火災或爆炸災害現(xiàn)場常有不同形式的障礙物,這些障礙物不但會影響火焰?zhèn)鞑サ耐緩竭€會對火焰的傳播速度、溫度及環(huán)境壓力造成一定的影響。由G.Ciccarelli等[1]、Chen Xianfeng等[2]的研究表明,DDT過程的主要促成因素是粗糙壁面的影響和爆燃波的誘導作用,而DDT過程則是氣體爆炸中危害最大的形式。因此研究障礙物對可燃氣體爆燃影響,對降低災害造成的破壞和損失有著重要的意義。

甲烷/空氣預混火焰在障礙物控制條件下在開口管道的燃燒實驗可以初步了解到由于障礙物的作用,使預混火焰在管道中傳播過程變得更加復雜,特別是管道內(nèi)的一些特性參數(shù)例如管內(nèi)流場等特性參量無法在實驗中觀測到。隨著流體力學數(shù)值計算精度和準確性的大幅提高,數(shù)值模擬已成為普遍的研究手段。S.Sklavounos等[3]、B.J.Arntzen[4]基于CFD軟件,對可燃氣云的擴散模型及氣體爆炸過程進行了模擬研究,得到與實驗相符合的數(shù)據(jù),驗證了數(shù)值模擬方法的可行性。而對障礙物對火焰影響的數(shù)值模擬研究,大多利用雷諾平均模型(Reynolds average Navier-Stokes,RANS),其中k-ε模型應用最為廣泛。但在應用RANS方法計算時,所得的火焰形狀與實驗中的火焰形狀不能很好吻合。利用LES模型能夠更清晰地反映出火焰及流場的大渦結構,特別是在受限空間中模擬火焰的傳播,可以得到較好的火焰精細結構和火焰內(nèi)外部流場特性。V.D.Sarli等[5]、A.R.Masri等[6]、S.R.Gubba等[7]應用LES模型進行了大量的三維模擬,結果與實驗十分相似,并且得出不同亞網(wǎng)格選擇對模擬結果影響較大的結論。另外,障礙物影響對預混火焰DDT過程有著明顯的促進作用,在對障礙物誘發(fā)的DDT過程的數(shù)值模擬研究也做了大量的工作。D.A.Kessler等[8]、V.N.Gamezo等[9]、T.Ogawa等[10]利用數(shù)值模擬研究在障礙物作用下DDT的形成過程及危害影響。預混火焰發(fā)生的火焰微觀結構變形及火焰失穩(wěn)的流場特性對預混火焰后期的加速及DDT過程的影響較大。本文中通過實驗發(fā)現(xiàn)并捕捉低壓環(huán)境下預混火焰通過障礙物的過程,并利用LES方法結合TFM模型對實驗過程進行重現(xiàn),以揭示預混火焰在障礙物影響下發(fā)生的復雜流場變化。

1 實 驗

采用高速紋影拍攝技術和高頻壓力傳感器對開口管道中甲烷/空氣預混火焰在障礙物影響下傳播過程中的火焰結構及傳播特性進行捕捉,如圖1所示。開口管道是為預混火焰?zhèn)鞑ヌ峁┑蛪涵h(huán)境的受限空間,設計為尺寸是8cm×8cm×50cm的長方體不銹鋼管道,尾端中心設有正方形泄壓口。任少鋒[11]認為泄壓口比率為30%可作為有效泄爆設計的重要參考值,因此為了有效泄爆和便于實驗的操作,這里設計泄壓口的尺寸為4cm×4cm,泄壓口比率為25%。實驗中障礙物設置于距管道首端30cm處,其尺寸為(即阻塞比為50%)。甲烷/空氣預混氣的體積比=1∶0.42,初始壓力為p0=1.0×105Pa,初始溫度T0=273K。

圖1 管道中心軸切面結構圖Fig.1 Pipe axis structure in cross-section

2 物理模型

物理模型的建立主要考慮甲烷/空氣的層流預混火焰厚度,通過對模型局部做了0.5、1.0、2.0和4.0mm的網(wǎng)格尺寸獨立性分析,模型選用網(wǎng)格尺寸為2mm×2mm×2mm的submap格式的六面體網(wǎng)格布置障礙物設計放置在x=30cm處。三維物理模型的單元格數(shù)量為726 628。

3 數(shù)學模型

3.1 大渦模型控制方程

在預混氣體燃燒中大渦模擬的基本原理是對其控制方程通過濾波方程在流場中區(qū)別出的大尺度渦團和Kolmogorov尺度渦團,然后把相關的流場變量區(qū)分成大尺度量和小尺度量。對大尺度量通過直接數(shù)值模擬方法計算,而對小尺度量采用亞網(wǎng)格模型(sub-grid scale model,SGS model)進行模型假設來模擬。大渦模型的控制方程為[12]:

連續(xù)性方程為:

動量方程為:

式中:LES濾波的參量標注為橫線上標,質(zhì)量權重濾波的參量標注為波浪線上標;ρ為密度,ui、uj為速度分量,t為時間,應力張量σij由分子黏度μ決定,可表示為,τij為亞網(wǎng)格尺度應力,將其定義為,下標k表示為Kolmogorov尺度。

能量方程為:

式中:hs為顯焓,λ為熱導率。亞網(wǎng)格熱焓通量可通過梯度假設近似為:

式中:μSGS為亞網(wǎng)格黏度,Pr為亞網(wǎng)格的普朗特數(shù),cp為定壓比熱,T為溫度。

亞網(wǎng)格模型選用WALE(wall-adapting local eddy-viscosity,WALE),WALE模型增強了大渦計算中的壁面處理,而且在計算層流流動狀態(tài)時,湍流黏度的返回值為零,比較適合本文中所研究的工況。

在WALE模型中將渦流黏度定義為:

式中:Ls為亞網(wǎng)格尺寸和應變率張量的混合長度量,可以定義為表示網(wǎng)格尺寸,Sij可定義為,為von Kármán常數(shù),d為流體質(zhì)點離最近壁面的距離,Cw為WALE常數(shù)。

3.2 增厚火焰模型

甲烷/空氣層流預混火焰厚度約為0.26mm,對網(wǎng)格尺寸要求苛刻。為降低計算成本,選用增厚火焰模型(TFM),它適用于較粗的網(wǎng)格來模擬預混火焰的傳播,且與LES模型配合計算結果較好[13]。

增厚火焰模型是通過層流火焰在不改變層流火焰?zhèn)鞑ニ俣鹊那闆r下被人為加厚,通過增加物質(zhì)擴散率和降低反應速率來實現(xiàn)。用ul表示層流預混火焰?zhèn)鞑ニ俣龋瑄l與成比例,D為分子擴散系數(shù),R為平均反應速率。層流火焰厚度與D/ul成比例。在ANSYS Fluent軟件中,增厚因子其中火焰內(nèi)指定的網(wǎng)格單元數(shù)N=5,層流火焰厚度σ=0.26;Δ為為網(wǎng)格尺寸,由定義,V為單元體積,d為空間維度。層流火焰厚度σ可為常數(shù)、也可自定義函數(shù),或由D/ul計算得到。其中熱擴散系數(shù)D,由決定,k為未燃混合氣體導熱系數(shù)。

所有組分擴散系數(shù),包括導熱系數(shù),均乘以增厚因子F;所有反應速率均除以F。然而,遠離火焰的區(qū)域會由于擴散率增加而出現(xiàn)混合及熱傳導錯誤,因此,火焰增厚區(qū)域應限制在火焰前鋒附近的狹窄區(qū)域。火焰增厚區(qū)域還要乘以系數(shù)Ω:

3.3 化學反應機理

偏重研究火焰在受限空間內(nèi)的傳播行為和流場特性,為了減少計算成本,并保證計算準確性,因此化學反應機理選用簡化后的甲烷/空氣兩步反應機理。其反應機理步驟如下:

3.4 Spark點火模型

為了更精準地模擬實驗中的電擊點火過程,模擬中使用的是Spark點火模型,它是根據(jù)網(wǎng)格的尺寸,假設一個完美的球形火焰面,球形火焰面的半徑r隨著時間軸的推移而逐漸增大。控制方程如下:

式中:ρu、ρb分別為未燃區(qū)和已燃區(qū)氣體密度,St為湍流火焰速度。為了控制火星的大小,給出了r的限制條件:,Δ為網(wǎng)格尺度,lt為湍流尺度,初始火星半徑r0設定為0.5mm。

4 結果分析

實驗通過高速攝影對預混火焰在開口管道中傳播和翻越障礙物過程的紋影圖像進行拍攝。圖2所示為預混火焰翻越障礙物過程紋影圖片。紋影圖片記錄了從點火后t=32~42ms火焰鋒面的傳播過程。從圖可以看出t=32~38ms的過程中預混火焰處于層流氣泡膨脹狀態(tài)。在t=36ms的紋影圖像中,球形火焰鋒面由于障礙物的阻礙作用球心開始向上傾斜,下部火焰鋒面開始拉伸。在t=38ms時火焰鋒面到達塊狀障礙物的橫軸坐標x=30cm處,但是火焰鋒面并未接觸障礙物壁面,而是因為誘導氣流的作用下,越過障礙物從障礙物上部的狹窄通道穿過。此時,火焰的鋒面彎折成約45°的銳角,尖端仍是直徑為1mm的半球火焰鋒面。當t=39ms時,火焰上端鋒面越過障礙物,而且從圖中可以發(fā)現(xiàn)在其下表面層流火焰厚度明顯增大。接著在t=40ms時,紅色標記的火焰增厚處其火焰鋒面開始出現(xiàn)失穩(wěn),并伴隨有細小的渦團出現(xiàn),同時在標記處出現(xiàn)擠壓褶皺現(xiàn)象,這里稱之為“誘導褶皺”。其后渦團隨時間擴大,逐漸發(fā)展成湍流狀態(tài)。

圖2 甲烷/空氣預混火焰?zhèn)鞑ミ^程紋影圖片F(xiàn)ig.2 High-speed schlieren images of premixed methane/air flame propagation

為了觀察三維預混火焰結構,選用T=1 800K的溫度等值面近視作為預混火焰鋒面,得到LES模型計算三維預混火焰結構時間序列圖,如圖3所示。將圖3與圖2進行比較,可以看出LES的計算結果與實驗的紋影圖像是比較吻合的。在t=5~10ms時,認為是火焰形成初期,火焰呈球形火焰成長,當火焰繼續(xù)成長階段,因受管道壁面擠壓,火焰形態(tài)拉長轉(zhuǎn)而向管道兩端傳播。對于預混火焰翻越障礙物后層流向湍流轉(zhuǎn)變的過程,LES模型的計算結果表現(xiàn)出色,可以很明顯地觀察到處在t=37ms時火焰失穩(wěn)并在火焰鋒面前段出現(xiàn)向內(nèi)卷曲,同時可以看到“誘導褶皺”的出現(xiàn)。而在t=42ms時可以清楚的看到,此時的預混火焰一部分已經(jīng)從泄壓口瀉出,另一部分則不斷向內(nèi)卷吸,最終形成1個C字形的火焰結構,火焰翻越障礙物時管道內(nèi)流場結構示意圖可參考圖4。

圖3 LES模型計算三維預混火焰結構時間序列圖Fig.3 Sequence diagram three-dimensional premixed flame structure by LES model

另一方面從圖3中可以觀察到預混火焰的三維結構中z軸的火焰結構特點,而這是實驗的紋影圖片捕捉不到的火焰結構特征。其中在t=30ms時看出在火焰遇到障礙物之前,火焰鋒面銳角化變形并向上偏移,而同時火焰鋒面下端,即面對障礙物的火焰鋒面出現(xiàn)向內(nèi)凹陷的現(xiàn)象。在t=39ms時發(fā)現(xiàn)在火焰越過障礙物后,火焰的失穩(wěn)造成的湍流使火焰頭部產(chǎn)生較大的扭曲變形,湍流強度有明顯增大的趨勢。另外從圖3中發(fā)現(xiàn)在開口管道中,障礙物的作用還會影響z軸方向的火焰鋒面形狀。

圖4 火焰翻越障礙物時管道內(nèi)流場結構示意圖Fig.4 The pipe flow field structure diagram when the flame climb over the obstacle

圖5所示為傳感器所在位置點(圖1中所示)的壓力時程曲線與LES模型計算值的對比,從圖中可以看出LES的計算結果與實驗值相當吻合,其上升趨勢和數(shù)據(jù)基本一致,但其峰值略高于實驗值20.8%。結合圖3中t=37ms時“誘導褶皺”效應也強于實驗值,其原因可以認為是LES模型運用WALE亞網(wǎng)格濾波,對壁面效應有所增強,使其對預混火焰翻越障礙物的過程中障礙物通道空間內(nèi)的壁面受沖擊壓力后的反作用壓力較真實值偏高,產(chǎn)生一定的放大效應。由此可以認為預混火焰?zhèn)鞑毫ES數(shù)值模擬與實驗結果較為吻合。

圖5 特征點處壓力時程曲線實驗與LES模型計算值對比Fig.5 Comparison of pressure histories at characteristic point between experiment and LES

截取y軸上截面對開口管道內(nèi)的流場進行分析。圖6展示的是在t=35,37,39ms時刻開口管道內(nèi)的氣體速度矢量以及流場情況。其中圖6(a)所示為管道內(nèi)的速度矢量圖,為了更好觀察管道內(nèi)流場變化,將圖6(a)中的速度矢量通過計算得到圖6(b)。

可以清楚的看到當t=35ms時火焰鋒面前部有較強的引導氣流繞過障礙物向x軸正向流動,同時對火焰鋒面有法向的擠壓作用,以致火焰鋒面向上移動并發(fā)生銳角化變形。在t=37ms時可以發(fā)現(xiàn)由于引導氣流作用下火焰鋒面通過障礙物通道時,障礙物后方形成回流氣流向x軸負方向翻越障礙物,與引導氣流形成對流,因此在障礙物臺階上方形成小渦團,對火焰鋒面形成擠壓與拉伸,并導致“誘導褶皺”的形成。同時這個小渦團使障礙物上方空間內(nèi)產(chǎn)生較大的速度梯度,從而促使預混火焰鋒面在翻越障礙物時,滿足Kelvin-Helmholtz不穩(wěn)定的條件導致火焰鋒面失穩(wěn)引發(fā)湍流。同時火焰鋒面的向右傳播和火焰氣泡膨脹的作用下,高溫低密度的燃燒產(chǎn)物流以火焰的傳播速度流入低溫高密度的未燃氣體,使火焰鋒面上形成了特定的密度梯度,同時在火焰鋒面與障礙物壁面之間的引導氣流的擠壓作用下,出現(xiàn)如圖2中火焰鋒面的小幅振動,引起Rayleigh-Taylor不穩(wěn)定現(xiàn)象。其流場結構如圖4所示,由于2個同向的渦流之間的擠壓,再加上火焰鋒面的切向移動速度,致使“誘導褶皺”的產(chǎn)生。而在2個渦流作用的同時,使障礙物后方的火焰鋒面上產(chǎn)生足夠大的速度與密度梯度,以致Rayleigh-Taylor不穩(wěn)定和Kelvin-Helmholtz不穩(wěn)定的耦合作用,使預混火焰鋒面失穩(wěn),失穩(wěn)火焰鋒面的擾動隨時間逐漸增大,火焰鋒面與未燃氣體快速卷吸,促使渦流形成,直至渦流破碎形成湍流。

圖6 預混火焰翻越障礙物過程中管道內(nèi)的速度矢量及流場變化圖Fig.6 Velocity vector and the flow field when premixed flame climbed over obstacle

在LES的計算結果中選取管道的截面上x=30cm處的截面進行分析得到圖7。從圖7中發(fā)現(xiàn)預混火焰在z軸方向上即火焰?zhèn)鞑サ姆ㄏ蚍较蛏弦矔霈F(xiàn)擠壓、拉伸褶皺。可以從t=39ms的火焰溫度云圖中發(fā)現(xiàn)火焰形態(tài)在法向方向上出現(xiàn)2道向內(nèi)的折痕。其產(chǎn)生原因則是由于陡坡的上升氣流,即前文中的提到的引導氣流在y軸上的速度矢量對預混火焰下端鋒面造成的擠壓。同時由于火焰處于受限空間中,側面壁面對火焰存在限制,以致于造成火焰在法向方向上出現(xiàn)褶皺變形。此現(xiàn)象對火焰的失穩(wěn)也存在促進作用。

圖7 預混火焰法向截面的溫度云圖及流場結構Fig.7 Temperature contours and the flow field structure of premixed flame normal section

5 結 論

利用大渦模型的數(shù)值模對障礙物影響預混火焰的實驗進行了重現(xiàn),通過實驗所得到的預混火焰微觀結構的紋影圖像及火焰?zhèn)鞑サ膲毫?shù)據(jù)驗證了模擬的可靠性。由大渦模型的數(shù)值模擬結果揭示了預混火焰在障礙物條件下的傳播過程中發(fā)生的火焰結構變化及火焰失穩(wěn)直接原因,并得到以下結論:

(1)通過大渦模擬的流場分析,在實驗現(xiàn)象中發(fā)現(xiàn)“誘導褶皺”是預混火焰翻越障礙物后有規(guī)律性的火焰褶皺,而非火焰鋒面隨機的擾動的結果;

(2)氣體預混火焰在遇障礙物后發(fā)生火焰失穩(wěn)的直接原因是由障礙物引發(fā)的3個氣流渦團同時作用而形成Kelvin-Helmholtz不穩(wěn)定及Rayleigh-Taylor不穩(wěn)定現(xiàn)象耦合作用,導致火焰的失穩(wěn)變形;

(3)在障礙物引發(fā)的上升氣流的拉伸、擠壓影響下,預混火焰在傳播的法向方向也會產(chǎn)生褶皺變形。

[1]Ciccarelli G,Johansen C,Kellenberger M.High-speed flames and DDT in very rough-walled channels[J].Combustion and Flame,2013,160(1):204-211.

[2]Chen Xianfeng,Zhang Yin,Zhang Ying.Effect of CH4-Air ratios on gas explosion flame microstructure and propa-gation behaviors[J].Energies,2012,5(10):4132-4146.

[3]Sklavounos S,Rigas F.Validation of turbulence models in heavy gas dispersion over obstacles[J].Journal of hazardous materials,2004,108(1):9-20.

[4]Arntzen B J.Modelling of turbulence and combustion for simulation of gas explosions in complex geometries[J].Journal of Loss Prevention in the Process Industries,1998,18(4/5/6):225-237.

[5]Sarli V D,Benedetto A D,Long E J,et al.Time-resolved particle image velocimetry of dynamic interactions between hydrogen-enriched methane/air premixed flames and toroidal vortex structures[J].International Journal of Hydrogen Energy,2012,37(21):16201-16213.

[6]Masri A R,Ibrahim S S,Cadwallader B J.Measurements and large eddy simulation of propagating premixed flames[J].Experimental Thermal and Fluid Science,2006,30(7):687-702.

[7]Gubba S R,Ibrahim S S,Malalasekera W,et al.Measurements and LES calculations of turbulent premixed flame propagation past repeated obstacles[J].Combustion and Flame,2011,158(12):2465-2481.

[8]Kessler D A,Gamezo V N,Oran E S.Simulations of flame acceleration and deflagration-to-detonation transitions in methane-air systems[J].Combustion and Flame,2010,157(11):2063-2077.

[9]Gamezo V N,Ogawa T,Oran E S.Flame acceleration and DDT in channels with obstacles:Effect of obstacle spacing[J].Combustion and Flame,2008,155(1/2):302-315.

[10]Ogawa T,Gamezo V N,Oran E S.Flame acceleration and transition to detonation in an array of square obstacles[J].Journal of Loss Prevention in the Process Industries,2013,26(2):355-362.

[11]任少峰.可燃性氣體泄爆動力學機理研究[D].武漢:武漢理工大學,2012.

[12]Lesieur M.Large-eddy simulations of turbulence[M].Cambridge:Cambridge University Press,2005.

[13]肖華華.管道中氫/空氣預混火焰?zhèn)鞑恿W實驗與數(shù)值模擬研究[D].合肥:中國科學技術大學,2013.

Large eddy simulation of impacted obstacles’effects on premixed flame’s characteristics

Wang Gongzhong1,2,Zhang Jianhua2,Li Dengke2,Chen Xianfeng2
(1.School of Safety Engineering,Henan Institute of Engineering,Zhengzhou 451191,Henan,China;2.Department of Safety Science and Engineering,Wuhan University of Technology,Wuhan 430070,Hubei,China)

In the process of the premixed gas flame propagation,obstacles have vital influence on the flame's structure and propagation,and will enhance the flame acceleration and the DDT process.Through the methane/air premixed flame propagation experiment with obstacles in the visual pipe,the flame microstructure changes and the propagation characteristics were captured.By means of the three-dimensional physical model,a large eddy simulation mode(LES)with the WALE sub-grid scale models,and the thickened flame model(TFM)were used to repeat experiment process.The complex changes of the flow field were obtained when the premixed flame climbed over the obstacle,and the characteristics of the flow turbulence transition were analyzed.Finally,the direct cause of the premixed flame disturbance instability under the influence of the obstacles was revealed.It was induced by the coupling effect of the Kelvin-Helmholtz instability and the Rayleigh-Taylor instability,which in turn was affected by three vortexes as a result of the obstacle.

mechanics of explosion;thickened flame model;large eddy simulation;premixed flame;obstacles

O383國標學科代碼:13035

A

10.11883/1001-1455(2017)01-0068-09

(責任編輯 王易難)

2015-05-12;

2015-11-18

國家自然科學基金項目(51174153,51374164);國家重點研發(fā)計劃項目(2016YFC0802801);

建筑消防工程技術公安部重點實驗室開放課題項目(KFKT2014ZD03)

王公忠(1978— ),男,博士,副教授;通信作者:陳先鋒,cxf618@whut.edu.cn。

猜你喜歡
實驗模型
一半模型
記一次有趣的實驗
微型實驗里看“燃燒”
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
做個怪怪長實驗
3D打印中的模型分割與打包
NO與NO2相互轉(zhuǎn)化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 狠狠色狠狠色综合久久第一次| 日韩国产亚洲一区二区在线观看| 久久亚洲AⅤ无码精品午夜麻豆| 国产日韩久久久久无码精品| 国产自在线播放| а∨天堂一区中文字幕| 亚洲精品国产精品乱码不卞| 亚洲成人精品| 成人精品午夜福利在线播放 | 国产Av无码精品色午夜| 欧美黑人欧美精品刺激| 久热99这里只有精品视频6| 中文字幕66页| 久久99国产乱子伦精品免| 国产91色在线| 色悠久久综合| 久久久噜噜噜| 欧洲日本亚洲中文字幕| 国产免费久久精品99re不卡| 女人一级毛片| 亚洲精品无码人妻无码| 成人噜噜噜视频在线观看| 精品撒尿视频一区二区三区| 欧美综合在线观看| 亚洲国产AV无码综合原创| 亚洲无码高清免费视频亚洲| 国产成a人片在线播放| 欧美精品成人| 日韩免费毛片视频| 中文天堂在线视频| 成人字幕网视频在线观看| 国产黑丝一区| 国产靠逼视频| 91久草视频| 日韩在线播放欧美字幕| 欧美日本中文| 亚洲视频免| 国产区人妖精品人妖精品视频| 午夜日韩久久影院| 国产特一级毛片| 欧美激情视频二区| 欧美特黄一免在线观看| 成人午夜精品一级毛片| 丝袜国产一区| 国模在线视频一区二区三区| 无码中文字幕精品推荐| 欧美区一区| 波多野结衣久久高清免费| 日本人真淫视频一区二区三区| 国产精品一线天| 婷婷六月色| 57pao国产成视频免费播放| 日本福利视频网站| 女人av社区男人的天堂| 国产三级视频网站| 伊人久久福利中文字幕| 亚洲色图欧美| 天天躁夜夜躁狠狠躁图片| 在线播放国产99re| 亚洲福利视频一区二区| 久久精品国产999大香线焦| 久久 午夜福利 张柏芝| 欧美色视频日本| 久久一本日韩精品中文字幕屁孩| 亚洲精品天堂自在久久77| 国产成a人片在线播放| 免费看久久精品99| 精品夜恋影院亚洲欧洲| 国产午夜精品鲁丝片| 国产精品免费入口视频| 最新午夜男女福利片视频| 欧日韩在线不卡视频| 精品久久久久久久久久久| 成人免费一区二区三区| 免费毛片a| 伊人丁香五月天久久综合 | 亚洲二区视频| 亚洲第一视频免费在线| 99久久精彩视频| 热99精品视频| 毛片一级在线| 欧美人人干|