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

瓦斯爆炸壓力與波前瞬態(tài)流速演化特征及其定量關(guān)系

2015-04-14 06:56:32林柏泉洪溢都朱傳杰江丙友孫豫敏
爆炸與沖擊 2015年1期
關(guān)鍵詞:實驗

林柏泉,洪溢都,朱傳杰,江丙友,劉 謙,孫豫敏

(1.中國礦業(yè)大學(xué)安全工程學(xué)院,江蘇 徐州 221116; 2.煤炭資源與安全開采國家重點實驗室,江蘇 徐州,221008; 3.安徽理工大學(xué),安徽 淮南 232001; 4.龍巖學(xué)院,福建 龍巖 364012)

?

瓦斯爆炸壓力與波前瞬態(tài)流速演化特征及其定量關(guān)系

林柏泉1,2,洪溢都1,2,朱傳杰1,2,江丙友3,劉 謙4,孫豫敏1,2

(1.中國礦業(yè)大學(xué)安全工程學(xué)院,江蘇 徐州 221116; 2.煤炭資源與安全開采國家重點實驗室,江蘇 徐州,221008; 3.安徽理工大學(xué),安徽 淮南 232001; 4.龍巖學(xué)院,福建 龍巖 364012)

為了建立爆炸波前的瞬態(tài)流速和超壓的定量關(guān)系,采用數(shù)值模擬的方法分別研究了開口型方管內(nèi)瓦斯爆炸超壓和瞬態(tài)流速傳播特征。研究結(jié)果表明:開口型方管內(nèi),波前瞬態(tài)流速和超壓的波形曲線的峰值個數(shù)不一樣,而且超壓峰值總是早于波前瞬態(tài)流速峰值出現(xiàn)。大部分情況下,方管橫截面邊長越大,其超壓峰值相對較小,并且超壓峰值沿傳播方向呈現(xiàn)降低趨勢。波前瞬態(tài)流速峰值沿傳播方向呈不斷增長趨勢,而且方管橫截面邊長越大,其峰值也相對較小。長徑比(方管長度與橫截面邊長之比)小于125時,超壓峰值與波前瞬態(tài)流速峰值的定量關(guān)系始終呈現(xiàn)線性反比關(guān)系;大于或等于125時,超壓峰值和波前瞬態(tài)流速峰值的定量關(guān)系呈分段關(guān)系。研究結(jié)果可為爆炸沖擊波揚塵的研究提供基礎(chǔ)數(shù)據(jù)。

爆炸力學(xué);瞬態(tài)流速;超壓;瓦斯爆炸

可燃性煤塵爆炸是煤礦災(zāi)害的重要形式之一,但是單一的煤塵爆炸在煤礦中很少見,它往往由瓦斯爆炸引發(fā),作為一種次生災(zāi)害發(fā)生。爆炸沖擊波能否揚起煤塵,是發(fā)生煤與瓦斯爆炸的關(guān)鍵。B.Fletcher[1]通過研究發(fā)現(xiàn)沖擊波前的高速氣流是粉塵揚起的關(guān)鍵因素,目前這種觀點已經(jīng)得到學(xué)界的認(rèn)可。因此,研究沖擊波揚塵問題關(guān)鍵在于是否能夠準(zhǔn)確了解爆炸沖擊波的波前流場特征。

由于當(dāng)前實驗條件有限,波前瞬態(tài)流速的分布特征仍然難以通過實驗方法獲得,目前能得到的最多的是用紋影或陰影法拍到的流場特征[2-4]。但這種方法歸根結(jié)底是一種定性方法,只能配合定量測試方法來定性描述爆炸火焰或沖擊波的傳播特征[5-7]。考慮到爆炸波的超壓是相對容易在實驗室測得的參數(shù),如果能建立爆炸波前瞬態(tài)流速與超壓的定量耦合關(guān)系,對認(rèn)識不同爆炸強度的揚塵能力具有重要意義。但目前很少有學(xué)者對爆炸的波前流速進行單獨研究。楊書召等[8]給出了沖擊波和波前流速關(guān)系的理論公式,但是該公式只適用于惰性沖擊波,也難以準(zhǔn)確估算沖擊波對氣體介質(zhì)所做的功。

基于上述研究現(xiàn)狀以及課題組此前所做的研究工作[9-10],本文中研究開口型系統(tǒng)內(nèi)的爆燃波前瞬態(tài)流速與爆炸超壓的變化規(guī)律,并嘗試建立兩者之間的關(guān)系,為預(yù)測開口型系統(tǒng)內(nèi)的沖擊波前流場特征提供一種新的方法,從而為評估不同強度下爆炸波沖擊波揚塵特征提供依據(jù)。

1 數(shù)值模型及驗證

1.1 數(shù)值模型

數(shù)值模擬軟件AutoReaGas能夠較好地模擬氣體爆炸,模擬的可靠性較高[11-13]。湍流作為氣體燃燒爆炸的重要因素,采用k-ε模型。燃燒反應(yīng)過程簡化成基元反應(yīng),而燃燒速率Rc表示為[14]

(1)

式中:Ct是量綱一常數(shù),ρ是密度,Γ是質(zhì)量擴散系數(shù),Rmin是可燃物、氧氣和燃燒產(chǎn)物中的最小值。 湍流燃燒速度St表示為[15]

(2)

式中:ut湍流強度,Lt湍流的特征長度,Sl是層流燃燒速度,ν是運動黏度。數(shù)值模擬中的其余設(shè)置參數(shù)參考W.P.M.Mercx等[16]的研究設(shè)定。初始階段的層流燃燒速率按照準(zhǔn)層流模型處理。Fs是另外一個重要的修正系數(shù),主要是為了修正壓力、溫度和火焰前沿褶皺對層流燃燒速度的影響。Sl,eff是湍流火焰速度,其與Fs、火焰半徑r和理論層流火焰速度Sl的關(guān)系表示為[17]

(3)

數(shù)值計算的初始參數(shù)按參考文獻[18]進行設(shè)置。

1.2 實驗驗證

有學(xué)者曾利用實驗礦井進行了大量實驗,發(fā)現(xiàn)只要數(shù)值模擬結(jié)果與實際實驗誤差值在±47%以內(nèi),數(shù)值模擬的結(jié)果就能夠滿足工程現(xiàn)場的需要[19]。C.J.Lea等[20]也通過實驗證實了AutoReaGas軟件模擬的可靠性。本文中通過相關(guān)實驗來驗證實驗數(shù)據(jù)與數(shù)值結(jié)果,從而驗證網(wǎng)格劃分的合理性以及模型的合理性。

1.2.1 實驗設(shè)備

實驗設(shè)備如圖1所示。實驗管道為方形不銹鋼管,管長L=5 m,橫截面邊長a=8 cm, 如圖1(a)所示。數(shù)值模擬所建模型管道與實驗管道一樣。管道左端為封閉點火端,管道右端開放。圖1(b)是活塞式壓力計,主要用于標(biāo)定傳感器,盡量避免測量誤差。壓力測點自管道左端0.5 m開始布置,每隔0.5 m安置一個,總共9個,火焰測點安裝在管道左端0.25 m處,起觸發(fā)開關(guān)用,如圖1(c)所示。點火能量是2 J。實驗瓦斯氣體體積分?jǐn)?shù)為9.5%。在對比參數(shù)的選取上,采用實驗方法易于獲得的爆炸超壓值[21]。

圖1 實驗設(shè)備Fig.1 Schematic of experimental equipment

1.2.2 實驗結(jié)果與數(shù)值模擬結(jié)果對比

圖2 管道實驗數(shù)據(jù)和數(shù)值模擬結(jié)果的對比Fig.2 Comparison of numerical and experimental results

數(shù)值模擬實驗中選取了2種網(wǎng)格進行對比選擇。一種網(wǎng)格尺寸為2 cm×2 cm×2 cm,另一種網(wǎng)格尺寸為4 cm×4 cm×4 cm,模擬結(jié)果與實驗結(jié)果對比如表1中所示。應(yīng)該說明的是,表格中對比的數(shù)據(jù)是超壓峰值pmax,測點與點火端的距離為l,數(shù)值模擬結(jié)果與管道實驗數(shù)據(jù)之間的偏差程度為E。從中可以發(fā)現(xiàn)2 cm×2 cm×2 cm的網(wǎng)格劃分所得的模擬結(jié)果與實際更吻合。因此用2 cm×2 cm×2 cm的網(wǎng)格進行更深入的模擬研究。圖2為數(shù)值模擬結(jié)果與管道實驗數(shù)據(jù)的對比,可以看出數(shù)值模擬結(jié)果與實驗結(jié)果吻合較好。表2為數(shù)值模擬結(jié)果與管道實驗數(shù)據(jù)以及兩者的偏差,其中最大的偏差是8.35%,明顯小于47%。因此,我們認(rèn)為在這種數(shù)值模型以及這種網(wǎng)格劃分方法下的模擬結(jié)果具有較大的可靠性。

表1 不同網(wǎng)格劃分方法下的數(shù)值模擬結(jié)果的對比

注:表中的偏差以管道實驗結(jié)果為基準(zhǔn),其順序與數(shù)值模擬結(jié)果的呈對應(yīng)關(guān)系。

表2 數(shù)值模擬結(jié)果與實驗數(shù)據(jù)的對比

2 數(shù)值計算結(jié)果與分析

2.1 超壓和波前瞬態(tài)流速的變化規(guī)律

圖3 超壓和波前瞬態(tài)流速隨時間的變化規(guī)律Fig.3 Overpressure and flow speed versus time

圖3為開口型管道內(nèi)的超壓p和波前瞬態(tài)流速v隨時間t的演化特征曲線。由于管道中火焰始終處于緩燃階段,前驅(qū)沖擊波尚未完全形成,其超壓值很小,所以超壓波形曲線只體現(xiàn)火焰鋒面的影響,而波前瞬態(tài)流速演化波形曲線卻表現(xiàn)出了前驅(qū)沖擊波和火焰鋒面兩者的影響。因此,前驅(qū)沖擊波和火焰鋒面的影響下,波前瞬態(tài)流速在26.9、43.3 ms時依次出現(xiàn)2個正向峰值,而由于波后回流現(xiàn)象的影響,在32.6、53.2 ms出現(xiàn)2個反向峰值;然后由于封閉端的反射作用,波前瞬態(tài)流速再一次加速,在83.3 ms時出現(xiàn)第3個正向峰值;最后由于氣體動能不斷耗散,波前瞬態(tài)流速不斷減速直至趨于零。比較而言,超壓的演化波形曲線稍微簡單一點,在33.1 ms時才出現(xiàn)首個正向峰值,而后不斷降壓,在53.2 ms時出現(xiàn)首個反向峰值,接著再次升壓,在78.7 ms時出現(xiàn)第2個正向峰值,而后不斷降壓,直至趨于一個定值。沖擊波傳播過程中,在到達某處時會首先對該處的氣體進行壓縮,使密度、壓力升高,而后由于密度增加、體積減小,流體質(zhì)點被沖擊波推動,這時才產(chǎn)生運動速度[22]。所以超壓峰值總是在波前瞬態(tài)流速峰值之前出現(xiàn),這在圖中清晰可見。

2.2 超壓峰值與波前瞬態(tài)流速峰值沿傳播方向上的分布規(guī)律

圖4~5為開口型管道內(nèi)超壓峰值pmax和波前瞬態(tài)流速峰值vmax沿傳播方向的演化曲線。管道橫截面邊長a為8、20和30 cm等3種,管道長度L有5、10、15和20 m等4種。從圖4中可以看到,L和a都能夠明顯影響超壓峰值的變化趨勢。當(dāng)L相同時,a越小,超壓值越大。而當(dāng)a為20和30 cm,超壓峰值在點火端要比管道末端大,沿傳播方向上呈降低趨勢。而當(dāng)a為8 cm時,其變化趨勢與L的值相關(guān)。5和10 m管道中超壓峰值沿傳播方向呈降低趨勢,最小值在管道末端;但是在15和20 m管道中超壓峰值卻呈現(xiàn)先降低后增長的變化趨勢,最小值在管道中間。

圖4 管道尺寸不同時超壓峰值沿傳播方向的分布Fig.4 Overpressures in different tubes

圖5 不同管道尺寸中流速峰值沿傳播方向的分布Fig.5 Flow velocities in different tubes

從圖5中可以看出,a值對波前瞬態(tài)流速峰值的影響與對超壓峰值的影響類似。a越小,流速峰值在不同管道的相同測點位置上越大。但是L的值對兩者的影響卻有明顯的不同,當(dāng)a相同時,波前瞬態(tài)流速峰值沿傳播方向上只呈明顯的增大趨勢,而不是超壓峰值的多樣變化。因此,L和a的綜合影響使得波前瞬態(tài)流速峰值在不同情況下有明顯區(qū)別。當(dāng)L為20 m、a為8 cm時,波前瞬態(tài)流速最大可達到852 m/s;而當(dāng)L為5 m、a為30 cm時,波前最大流速最大才278 m/s,兩者相差數(shù)倍。

2.3 超壓峰值與波前瞬態(tài)流速峰值的定量關(guān)系

圖6為超壓峰值與波前瞬態(tài)流速峰值的定量關(guān)系。圖中r為長徑比,即L與a的比值。從圖4~5可知,除了在a為8 cm、L為15和20 m時超壓先減小后增大的情況外,當(dāng)a的值相同時,波前瞬態(tài)流速沿傳播方向不斷增大,而超壓沿傳播方向不斷減小。當(dāng)火焰處于緩燃階段時,氣體爆炸釋放的能量較小,超壓和波前瞬態(tài)流速會出現(xiàn)此消彼長的變化;而當(dāng)火焰處于爆燃階段時,氣體爆炸釋放的能量足以維持超壓和波前瞬態(tài)流速的同時增長[23]。所以,當(dāng)r較小時,管道中的火焰始終處于緩燃階段,超壓和波前瞬態(tài)流速的耦合關(guān)系是反比關(guān)系;而當(dāng)r較大時,火焰在經(jīng)歷過一段緩燃階段之后會進入爆燃階段,所以波前瞬態(tài)流速和超壓的耦合關(guān)系早期是反比關(guān)系,但是到后期是正比關(guān)系,這在圖6中有著直接的體現(xiàn)。從圖中可以清晰地看出,當(dāng)r<125時,超壓峰值和波前瞬態(tài)流速峰值都呈單一線性反比關(guān)系;當(dāng)r=125時,兩者之間的定量關(guān)系開始呈分段關(guān)系,并且兩段都呈線性反比關(guān)系;等到r=187.5、250時,則早期依然呈反比關(guān)系,但是后期就變?yōu)檎汝P(guān)系。

圖7為圖6中超壓與波前瞬態(tài)流速耦合得到的斜率和截距。圖形中左軸是截距h,右軸是斜率k,橫坐標(biāo)是r。當(dāng)r<100時,k與r呈反比關(guān)系,如圖中B所示。當(dāng)r>125時,即當(dāng)a=8 cm、L=15和20 m時,k與r在早期仍然呈反比關(guān)系,如圖中C所示,而在后期則轉(zhuǎn)換為正比關(guān)系,如圖中D所示。而對于h而言,在r<100時,h與r呈正比關(guān)系,如圖中A所示,但是在r>125時,h與r卻都是反比關(guān)系,如圖中E和F所示。

圖6 超壓峰值與波前瞬態(tài)流速峰值的定量關(guān)系Fig.6 Relationships between flow velocities and overpressures

圖7 斜率和截距的分布關(guān)系Fig.7 Fitting curves of slope and intercept

3 結(jié) 論

(1)開口型管道內(nèi),波前瞬態(tài)流速和超壓在前驅(qū)沖擊波和火焰鋒面的不同影響下表現(xiàn)出不一樣的演化波形曲線。超壓演化波形曲線只有2個正向峰值和1個反向峰值,而波前瞬態(tài)流速演化波形曲線卻有3個正向峰值和2個反向峰值。此外,超壓峰值總是早于波前瞬態(tài)流速峰值出現(xiàn)。

(2)管道長度L和管道橫截面邊長a對超壓峰值和流速峰值都有明顯影響。a值較大的管道,其超壓峰值比a值較小管道中的超壓峰值小,而且在a值較大管道中超壓峰值的降低趨勢要比a值較小管道中的降低趨勢更為平緩,同時超壓峰值沿傳播方向上呈降低趨勢。但是當(dāng)a=8 cm、L=15和20 m的管道中,超壓峰值沿傳播方向上表現(xiàn)出先下降后上升的變化趨勢。波前瞬態(tài)流速峰值在a值較大的管道中比a值較小的管道中要小,而且其沿傳播方向上呈增長趨勢。

(3)在對超壓峰值和流速峰值的定量關(guān)系分析中發(fā)現(xiàn),當(dāng)長徑比小于125時,超壓峰值和波前瞬態(tài)流速峰值的定量關(guān)系始終呈單一線性反比關(guān)系;當(dāng)長徑比是125時,兩者之間的定量關(guān)系開始呈分段關(guān)系,并且兩段都呈線性反比關(guān)系;等到長徑比為187.5和250時,則定量關(guān)系早期依然呈反比關(guān)系,但是后期就變?yōu)檎汝P(guān)系。

[1] Fletcher B. The interaction of a shock with a dust deposit[J]. Journal of Physics, D: Applied Physics, 1976,9(2):197-202.

[2] Kuznetsov M, Alekseev V, Matsukov I, et al. DDT in a smooth tube filled with a hydrogen-oxygen mixture[J]. Shock Waves, 2005,14(3):205-215.

[3] Ilbas M, Crayford A P, Ylmaz I, et al. Laminar-burning velocities of hydrogen-air and hydrogen-methane-air mixtures: An experimental study[J]. International Journal of Hydrogen Energy, 2006,31(12):1768-1779.

[4] Johansen C T, Ciccarelli 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.

[5] Ciccarelli G, Johansen C T, 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.

[6] Ciccarelli G, Johansen C T, 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.

[7] Ciccarelli G, Dorofeev S. Flame acceleration and transition to detonation in ducts[J]. Progress in Energy and Combustion Science, 2008,34(4):499-550.

[8] 楊書召,景國勛,賈智偉.礦井瓦斯爆炸沖擊氣流傷害研究[J].煤炭學(xué)報,2009,34(10):1354-1358. Yang Shu-zhao, Jing Guo-xun, Jia Zhi-wei. Injury study on impact current of gas explosion in coal mine[J]. Journal of China Coal Society, 2009,34(10):1354-1358.

[9] Zhu C J, Lin B Q, Hong Y D, et al. Numerical simulations on relationships between gas velocity and overpressure of gas explosions in ducts[J]. Disaster Advances, 2013,6(S1):217-227.

[10] Lin B Q, Hong Y D, Zhu C J, et al. Effect of length on the relationships between the gas velocity and peak overpressure of gas explosion disasters in closed-end pipes[J]. Disaster Advances, 2013,6(S2):176-184.

[11] Jiang B Y, Lin B Q, Shi S L, et al. A numerical simulation of the influence initial temperature has on the propagation characteristics of, and safe distance from, a gas explosion[J]. International Journal of Mining Science and Technology, 2012,22(3):307-310.

[12] Maremonti M, Russo G, Salzano E, et al. Numerical simulation of gas explosions in linked vessels[J]. Journal of Loss Prevention in the Process Industries, 1999,12(3):189-194.

[13] Pang L, Zhang Q, Wang T, et al. Influence of laneway support spacing on methane/air explosion shock wave[J]. Safety Science, 2012,50(1):83-89.

[14] Janovsky B, Selesovsky P, Horkel J, et al. Vented confined explosions in Stramberk experimental mine and AutoReaGas simulation[J]. Journal of Loss Prevention in the Process Industries, 2006,19(2):280-287.

[15] Bray K N C. Studies of the turbulent burning velocity[J]. Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences, 1990,431(1882):315-335.

[16] Mercx W P M, Van den Berg A C, Hayhurst C J, et al. Developments in vapour cloud explosion blast modeling[J]. Journal of Hazardous Materials, 2000,71(1):301-319.

[17] Bakke J R. Numerical simulation of gas explosions in two-dimensional geometries[J]. Christian Michelsen Institute, 1986:865403-865408.

[18] AutoReaGas user manual version 3.1[Z]. England: Century Dynamics and TNO, 2002.

[19] Zipf R K, Sapko M J, Brune J F. Explosion pressure design criteria for new seals in US coal mines[S]. 2007.

[20] Lea C J, Ledin H S. A review of the state-of-the-art in gas explosion modelling[M]. Health and Safety Laboratory, 2002.

[21] Salzano E, Marra F S, Russo G, et al. Numerical simulation of turbulent gas flames in tubes[J]. Journal of Hazardous Materials, 2002,95(3):233-247.

[22] 吳望一.流體力學(xué)(下冊)[M].北京:北京大學(xué)出版社,2011:414-432.

[23] Baker W E, Cox P A, Kulesz J J, et al. Explosion hazards and evaluation[M]. Access Online via Elsevier, 1983.

(責(zé)任編輯 曾月蓉)

Quantitative relationship between flow speed and overpressure of gas explosion in the open-end square tube

Lin Bai-quan1,2, Hong Yi-du1,2, Zhu Chuan-jie1,2, Jiang Bing-you3, Liu Qian4, Sun Yu-min1,2

(1.FacultyofSafetyEngineering,ChinaUniversityofMiningandTechnology,Xuzhou, 221116Jiangsu,China; 2.StateKeyLaboratoryofCoalResourcesandSafeMining,Xuzhou221008,Jiangsu,China; 3.SchoolofMiningandSafetyEngineering,AnhuiUniversityofScienceandTechnology,Huainan232001,Anhui,China; 4.SchoolofResourceEngineering,LongyanUniversity,Longyan364012,Fujian,China)

The main objective of this study is to establish the quantitative relationship between overpressure and flow speed in the open-end square tube by the numerical simulation. It is found that the numbers of the peak value of overpressure and flow speed at the same measured point are different. The peak overpressure always appears earlier than peak flow speed in time scale. In mast cases, larger side lenth of square tube corresponds to smaller peak overpressure, and the peak flow speed goes down slowly along the propagation direction. The peak overpressure decreases with the increasing of the distance far from ignition end. However, the peak flow speed increases with the stream wise distance. When normalized distance is less than 125, the peak overpressure and peak flow speed always presents an inverse relationship. Otherwise, the relationship is piecewise-linear. The results may provide reference for the study on evaluating the dust lifting ability behind shock wave in the limited spaces.

mechanics of explosion; flow speed; overpressure; gas explosion

10.11883/1001-1455(2015)01-0108-08

2013-04-25;

2013-11-13

國家重點基礎(chǔ)研究發(fā)展計劃(973計劃)項目(2011CB201205);國家自然科學(xué)基金青年科學(xué)基金項目(51204174);中國礦業(yè)大學(xué)人才引進&青年教師啟航計劃(2011RC07);中央高校基本科研業(yè)務(wù)費專項基金項目(2012QNB01)

林柏泉(1960— ),男,博士,教授,lbq21405@126.com。

O389;TD712.7 國標(biāo)學(xué)科代碼: 13035

A

猜你喜歡
實驗
我做了一項小實驗
記住“三個字”,寫好小實驗
我做了一項小實驗
我做了一項小實驗
記一次有趣的實驗
有趣的實驗
小主人報(2022年4期)2022-08-09 08:52:06
微型實驗里看“燃燒”
做個怪怪長實驗
NO與NO2相互轉(zhuǎn)化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
主站蜘蛛池模板: 久久美女精品| 日韩少妇激情一区二区| 久久香蕉国产线看观看精品蕉| 91毛片网| 免费视频在线2021入口| 国产美女丝袜高潮| 亚洲一区二区日韩欧美gif| 美女啪啪无遮挡| 热伊人99re久久精品最新地| 日韩在线影院| 动漫精品啪啪一区二区三区| 午夜少妇精品视频小电影| 欧亚日韩Av| 91在线一9|永久视频在线| 最新日本中文字幕| 国产99精品久久| 国产又粗又爽视频| 亚洲天堂在线免费| 激情午夜婷婷| 国产最新无码专区在线| 国内精品小视频在线| 国产精品一区在线观看你懂的| 日韩毛片在线播放| 亚洲色图另类| 中文字幕丝袜一区二区| 久久精品女人天堂aaa| 精品久久香蕉国产线看观看gif| 日韩免费中文字幕| 伊人色天堂| 四虎影视无码永久免费观看| 国产高清国内精品福利| 天堂网亚洲综合在线| 日韩无码视频网站| 99在线国产| 日韩资源站| 九九久久99精品| 四虎影视国产精品| 国产区精品高清在线观看| 人人91人人澡人人妻人人爽| 亚洲一区无码在线| 国产精品偷伦视频免费观看国产| 亚洲婷婷六月| 精品国产欧美精品v| 亚洲第一成年网| 欧美日韩激情在线| 欧美成人精品高清在线下载| 国产91精品久久| 久久人搡人人玩人妻精品一| 中文字幕亚洲电影| 国产成人综合亚洲欧洲色就色| 激情国产精品一区| 国产福利在线免费观看| 亚洲性网站| 亚洲第一网站男人都懂| 亚洲人成影院在线观看| 国产尤物视频在线| 亚洲男人的天堂久久香蕉网| 亚洲美女高潮久久久久久久| 亚洲日本中文综合在线| 国产激情无码一区二区免费| 精品国产自| 久久成人免费| 亚洲AⅤ无码国产精品| 国产91丝袜在线播放动漫 | 三上悠亚在线精品二区| 亚洲欧美日韩动漫| www.国产福利| 亚洲精品男人天堂| av在线无码浏览| 中文字幕在线看| 波多野结衣无码视频在线观看| 香蕉久久永久视频| 扒开粉嫩的小缝隙喷白浆视频| 久久精品国产精品青草app| 福利在线一区| 久久精品国产精品国产一区| 青青草国产精品久久久久| 亚洲欧美精品日韩欧美| 欧美a网站| 久久这里只精品国产99热8| 国产国产人免费视频成18| 国产精品亚欧美一区二区三区|