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

激光輻照下黑索今固相與氣相反應(yīng)模型對(duì)比

2016-11-03 05:35:23張家雷劉倉(cāng)理譚福利王偉平
固體火箭技術(shù) 2016年5期
關(guān)鍵詞:區(qū)域模型

張家雷,劉倉(cāng)理,譚福利,王偉平

(1.中國(guó)工程物理研究院 流體物理研究所,綿陽(yáng) 621900;2.中國(guó)工程物理研究院,綿陽(yáng) 621900)

?

激光輻照下黑索今固相與氣相反應(yīng)模型對(duì)比

張家雷1,劉倉(cāng)理2,譚福利1,王偉平1

(1.中國(guó)工程物理研究院 流體物理研究所,綿陽(yáng)621900;2.中國(guó)工程物理研究院,綿陽(yáng)621900)

為了研究激光輻照下含能材料的點(diǎn)火和燃燒機(jī)理,建立了典型含能材料黑索今(RDX)的固相與氣相反應(yīng)模型。氣相模型考慮了多組分氣體化學(xué)反應(yīng),采用有限差分方法求解控制方程,模擬了黑索今激光點(diǎn)火全過(guò)程。固相模型只考慮固相區(qū)域的傳熱和化學(xué)反應(yīng),采用有限元方法求解控制方程,獲得了黑索今溫度變化和點(diǎn)火時(shí)刻。最后,將2種模型計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)進(jìn)行了比較。結(jié)果表明,2種模型的激光點(diǎn)火時(shí)刻,計(jì)算值均與試驗(yàn)值吻合;氣相模型更能全面描述含能材料的點(diǎn)火過(guò)程,但計(jì)算過(guò)程相對(duì)復(fù)雜;固相模型不能描述點(diǎn)火過(guò)程中氣相產(chǎn)物分解和反應(yīng)過(guò)程,但計(jì)算過(guò)程相對(duì)簡(jiǎn)單,能模擬的入射熱流范圍更大,預(yù)測(cè)的點(diǎn)火時(shí)刻精度也滿足工程需要。

激光;點(diǎn)火;有限差分法;有限元法

0 引言

含能材料是固體推進(jìn)劑的重要組成部分,其燃燒機(jī)理研究對(duì)于火箭推進(jìn)劑的工藝設(shè)計(jì)和安定性評(píng)估都具有重要意義。相比火焰加熱和電加熱等傳統(tǒng)加載手段,激光能夠提供更加精確可測(cè)的輻照熱流,在材料燃燒機(jī)理和診斷方面具有獨(dú)特優(yōu)勢(shì)。目前,國(guó)外學(xué)者已經(jīng)廣泛應(yīng)用連續(xù)激光來(lái)研究含能材料的燃燒機(jī)理[1],即通過(guò)激光輻照來(lái)引燃含能材料,再用微型熱電偶、質(zhì)譜儀和光譜儀來(lái)測(cè)試火焰溫度和反應(yīng)氣體組分濃度等數(shù)據(jù)。這種激光點(diǎn)火技術(shù)為含能材料燃燒機(jī)理研究提供了新的試驗(yàn)方法。與試驗(yàn)相對(duì)應(yīng),含能材料反應(yīng)動(dòng)力學(xué)理論研究也在快速發(fā)展。特別是詳細(xì)化學(xué)反應(yīng)動(dòng)力學(xué)模型逐步替代了傳統(tǒng)的宏觀反應(yīng)動(dòng)力學(xué)模型,使得計(jì)算模型能夠描述含能材料氣相產(chǎn)物之間復(fù)雜的反應(yīng)動(dòng)力學(xué)過(guò)程。Liau等[2]較早基于詳細(xì)化學(xué)反應(yīng)動(dòng)力學(xué),建立了黑索今(RDX)激光點(diǎn)火模型,模擬了激光輻照下黑索今的加熱、熔化、蒸發(fā)、氣體反應(yīng)和點(diǎn)火的全過(guò)程。Meredith[3]建立HMX氣相反應(yīng)模型,模擬了HMX激光點(diǎn)火過(guò)程和快烤燃過(guò)程。Davidsonm[4]模擬了RDX、HMX和GAP的點(diǎn)火燃燒過(guò)程,含能材料分為固相、液相和氣相3部分,細(xì)致分析了穩(wěn)定燃燒波的結(jié)構(gòu)。Beckstead[5]總結(jié)了近年來(lái)固體推進(jìn)劑點(diǎn)火和燃燒理論工作進(jìn)展,特別是含能材料激光點(diǎn)火數(shù)值模擬進(jìn)展。國(guó)內(nèi)激光與含能材料相互作用研究,主要集中在脈沖激光對(duì)含能材料的起爆試驗(yàn)技術(shù)研究,有關(guān)激光輻照下含能材料燃燒行為研究較少。田占東等[6]建立了黑索今激光點(diǎn)火的一維氣相模型,模擬了典型工況下黑索今的氣相點(diǎn)火過(guò)程,點(diǎn)火時(shí)間計(jì)算結(jié)果與國(guó)外文獻(xiàn)數(shù)據(jù)一致。

激光輻照下含能材料點(diǎn)火模型主要包括氣相模型和固相模型2種。氣相模型基于詳細(xì)反應(yīng)動(dòng)力學(xué),能夠處理多種氣體組分的化學(xué)反應(yīng),并給出溫度、速度和氣體濃度等關(guān)鍵參數(shù)。固相模型基于宏觀反應(yīng)動(dòng)力學(xué),假設(shè)點(diǎn)火發(fā)生在固相區(qū)域,主要計(jì)算固相的溫升和反應(yīng)過(guò)程。這兩類模型在文獻(xiàn)中都具有廣泛應(yīng)用,且各自具有優(yōu)缺點(diǎn)。氣相模型能夠描述的詳細(xì)地點(diǎn)火過(guò)程,但計(jì)算過(guò)程非常復(fù)雜,目前還只能針對(duì)典型的含能材料(RDX和HMX等)進(jìn)行一維數(shù)值模擬。固相模型只考慮固相點(diǎn)火,不能反映氣相反應(yīng),但計(jì)算相對(duì)簡(jiǎn)單。迄今為止,國(guó)內(nèi)外文獻(xiàn)對(duì)這2種反應(yīng)模型的適用條件并沒(méi)有明確規(guī)定,對(duì)計(jì)算結(jié)果也沒(méi)有細(xì)致比較。因此,有必要對(duì)比研究2種反應(yīng)模型的預(yù)測(cè)結(jié)果,并確定其適用判據(jù)。

本文建立了激光輻照下黑索今的固相和氣相反應(yīng)模型。氣相模型考慮了氣體組分之間的詳細(xì)化學(xué)反應(yīng),采用一維有限差分法求解。固相模型采用簡(jiǎn)化的宏觀反應(yīng)模型,采用有限元方法求解。基于這2種反應(yīng)模型,模擬了激光輻照下典型含能材料黑索今的點(diǎn)火燃燒過(guò)程,并將二者結(jié)果與試驗(yàn)值進(jìn)行了比較。

1 激光輻照下黑索今反應(yīng)模型

1.1氣相反應(yīng)模型

激光輻照下黑索今點(diǎn)火過(guò)程大致如下:首先,激光輻照下固體黑索今的表面區(qū)域溫度迅速升高,最先到達(dá)熔點(diǎn),轉(zhuǎn)變?yōu)橐簯B(tài)黑索今;接著,液態(tài)黑索今溫度繼續(xù)升高,一部分蒸發(fā)為氣相黑索今,一部分熱分解為氣體;最后,氣相黑索今吸收激光能量,溫度進(jìn)一步升高,同時(shí)氣體產(chǎn)物組分之間迅速反應(yīng),并最終形成可自持的燃燒。因此,黑索今的穩(wěn)態(tài)燃燒波包含了固相、固液混合相以及液態(tài)相,如圖1所示。其中,泡沫層是液相和氣相的混合區(qū)域,泡沫層與固相之間為熔化線;氣相區(qū)域發(fā)生劇烈化學(xué)反應(yīng),形成可持續(xù)的穩(wěn)定燃燒波。

圖1 黑索今燃燒波結(jié)構(gòu)示意圖

氣相模型以液體和氣體分界面作為坐標(biāo)原點(diǎn),分別建立凝聚相(包含了固相和液相)和氣相的控制方程組。

在凝聚相區(qū)域,主要的物理化學(xué)過(guò)程是激光輻照下,含能材料快速升溫熔化,液態(tài)含能材料表面發(fā)生蒸發(fā)形成氣態(tài)物質(zhì),同時(shí)液態(tài)物質(zhì)內(nèi)部發(fā)生分解反應(yīng),產(chǎn)生氣體。固相的熱交換形式為熱傳導(dǎo),而液相的熱交換形式為熱傳導(dǎo)和對(duì)流換熱。因此,凝聚相區(qū)域的控制方程組包括了能量方程和組分連續(xù)方程,即

(1)

(2)

能量方程左邊包含了熱對(duì)流和熱傳導(dǎo),右邊為熱源項(xiàng),即化學(xué)反應(yīng)熱源以及激光體吸收熱源。

在氣相區(qū)域,主要物理化學(xué)過(guò)程主要是氣相產(chǎn)物流動(dòng)換熱,產(chǎn)物之間發(fā)生化學(xué)反應(yīng)。控制方程組包括了質(zhì)量守恒、動(dòng)量守恒、能量守恒和組分連續(xù)方程,分別為

(3)

(4)

(5)

(6)

式中P為混合物的壓力;μ為粘性系數(shù);Vk為組分k在混合物中的擴(kuò)散系數(shù);Sg為氣相吸收的激光能量(其余變量意義與凝聚相相同)。

激光點(diǎn)火過(guò)程中,黑索今氣相區(qū)域發(fā)生了多組分氣體反應(yīng),反應(yīng)過(guò)程較復(fù)雜。本文采用詳細(xì)的反應(yīng)動(dòng)力學(xué)氣體模型,即45組分232反應(yīng)道模型,化學(xué)反應(yīng)數(shù)據(jù)來(lái)源于文獻(xiàn)[4]。黑索今固相、液相和氣相熱物性參數(shù)采用文獻(xiàn)[5]數(shù)據(jù)。同時(shí),為了使方程封閉,還需要?dú)怏w狀態(tài)方程、邊界條件和激光吸收等。

1.2固相反應(yīng)模型

固相模型忽略了氣相化學(xué)反應(yīng)過(guò)程,只考慮凝聚相的傳熱和化學(xué)反應(yīng)。為了簡(jiǎn)化分析,可將物理問(wèn)題簡(jiǎn)化為軸對(duì)稱模型。在柱坐標(biāo)系統(tǒng)中,激光輻照下黑索今的能量平衡方程為

(7)

式中S=S1+S2,S1代表黑索今自熱反應(yīng)源項(xiàng),S2代表激光體吸收源項(xiàng)。

S1可用Arrhenius定律來(lái)表述為

(8)

式中ρ為黑索今密度;Q為反應(yīng)熱;A為指前因子;E為活化能;R為普適氣體常數(shù)。

S2代表激光體吸收源項(xiàng),即

(9)

式中f為材料表面的激光反射率;β為材料對(duì)激光的吸收系數(shù);I0為激光功率密度。

固相模型所用的黑索今熱物理性能參數(shù)如表1所示。考慮波長(zhǎng)為10.6 μm的CO2連續(xù)激光,黑索今的激光能量吸收系數(shù)1.2×105,吸收率為0.86。

表1 黑索今熱物理性能參數(shù)

2 求解方法

氣相模型采用有限差分法求解控制方程。在氣相區(qū)域,對(duì)流換熱效應(yīng)顯著;而在凝聚相區(qū)域,熱傳導(dǎo)是主要的換熱方式。因此,氣相計(jì)算區(qū)域比固相計(jì)算區(qū)域大。氣相模型計(jì)算過(guò)程比較復(fù)雜,詳細(xì)計(jì)算方法主要參考了文獻(xiàn)[3],相關(guān)計(jì)算程序是參考了文獻(xiàn)[6]工作,并對(duì)計(jì)算流程進(jìn)行了優(yōu)化,改進(jìn)了溫度場(chǎng)顯示方式。固相點(diǎn)火模型采用有限元方法求解,計(jì)算程序?yàn)樽跃幍腇ortran程序,能夠模擬激光輻照下金屬和復(fù)合材料的傳熱-燒蝕過(guò)程[7],本文在原程序基礎(chǔ)上,擴(kuò)充了含能材料化學(xué)反應(yīng)計(jì)算功能。

3 結(jié)果分析與討論

3.1點(diǎn)火過(guò)程溫度場(chǎng)結(jié)果

基于氣相模型模擬了黑索今激光點(diǎn)火全過(guò)程(圖2和圖3)。圖2給出了不同功率密度激光輻照下黑索今的點(diǎn)火過(guò)程。外部環(huán)境壓力為1個(gè)標(biāo)準(zhǔn)大氣壓,激光功率密度分別為200、400、600、800 W/cm2,初始溫度為300 K,零點(diǎn)坐標(biāo)為氣相和凝聚相界面,負(fù)值區(qū)域[-0.01 cm,0]代表凝聚相,正值區(qū)域[0,1.0 cm]代表氣相。

(a) 200 W/cm2  (b) 400 W/cm2

(c) 600 W/cm2  (d) 800 W/cm2

由圖2可見,激光輻照下凝聚相黑索今迅速升溫達(dá)到熔點(diǎn),開始發(fā)生分解和蒸發(fā),分解產(chǎn)物和氣體蒸氣進(jìn)行氣相區(qū)域,并開始輕微的放熱反應(yīng);當(dāng)氣相區(qū)域的化學(xué)反應(yīng)釋放大量熱量后,溫度升高到1 500 K,形成初級(jí)火焰;隨著反應(yīng)繼續(xù),化學(xué)反應(yīng)釋放出大量能量,氣相溫度達(dá)3 000 K以上,形成次級(jí)火焰,并最終達(dá)到穩(wěn)定燃燒狀態(tài)。比如,激光功率密度為600 W/cm2工況,當(dāng)輻照時(shí)間4.5 ms時(shí),氣相區(qū)域達(dá)到1 500 K,形成初級(jí)火焰;當(dāng)輻照時(shí)間4.7 ms時(shí),次級(jí)火焰開始準(zhǔn)備;當(dāng)輻照時(shí)間4.8 ms時(shí),次級(jí)火焰形成,氣相區(qū)域最高溫度超過(guò)3 000 K;當(dāng)輻照時(shí)間4.8 ms以后,氣相發(fā)生點(diǎn)火,火焰中心向表面移動(dòng),然后再遠(yuǎn)離表面移動(dòng),并達(dá)到穩(wěn)定燃燒狀態(tài)(對(duì)應(yīng)5.5 ms)。由此可見,氣相模型預(yù)測(cè)的黑索今激光點(diǎn)火過(guò)程清晰,與點(diǎn)火過(guò)程試驗(yàn)觀測(cè)非常相符[5]。

由圖2還可發(fā)現(xiàn),隨激光輻照功率密度增大,黑索今點(diǎn)火位置距離氣/液界面更近,點(diǎn)火前的氣相反應(yīng)區(qū)域尺寸更小。

火焰溫度分布是表征含能材料燃燒特性的重要數(shù)據(jù)。Lee[1]開展了CO2激光輻照下黑索今的燃燒特性試驗(yàn)研究,測(cè)試了火焰溫度分布和氣相反應(yīng)組分濃度。Parr[8]也開展了類似的試驗(yàn)工作。本文給出了凝聚相表面以上的黑索今火焰溫度分布情況,并與試驗(yàn)值進(jìn)行了比較,如圖3所示。由圖3可見,火焰溫度分布計(jì)算值與試驗(yàn)值吻合較好,試驗(yàn)也存在1 500 K左右的初級(jí)火焰和3 000 K左右的次級(jí)火焰,但高溫階段計(jì)算值比試驗(yàn)值稍高,Beckstead[5]在綜述文章中也提到了這個(gè)問(wèn)題,認(rèn)為這可能來(lái)源于微型熱電偶的測(cè)量誤差。

圖3 黑索今穩(wěn)定燃燒狀態(tài)溫度分布

固相模型計(jì)算的黑索今點(diǎn)火溫度曲線如圖4所示,激光功率密度為200、400、600、800 W/cm2,初始溫度為300 K。由圖4可見,隨激光輻照功率密度提高,黑索今表面溫升率越大,從緩慢加熱快速放熱反應(yīng)的時(shí)間越短。固相模型化學(xué)反應(yīng)相對(duì)簡(jiǎn)單,只是將含能材料反應(yīng)熱加入了控制方程,從宏觀角度模擬了外熱流作用下含能材料從加熱到點(diǎn)火的過(guò)程。

3.2點(diǎn)火時(shí)刻比較

點(diǎn)火時(shí)刻對(duì)于含能材料燃燒非常重要。這里將激光輻照下黑索今點(diǎn)火時(shí)刻計(jì)算值和試驗(yàn)值進(jìn)行了比較,如圖5所示。

圖4 固相模型黑索今前表面溫度-時(shí)間曲線

圖5 計(jì)算值和試驗(yàn)點(diǎn)火時(shí)刻對(duì)比

由圖5可見,當(dāng)激光功率密度為102W/cm2量級(jí)時(shí),氣相模型和固相模型給出的點(diǎn)火時(shí)刻預(yù)測(cè)值非常接近,均與試驗(yàn)值比較吻合,特別是在高功率水平范圍內(nèi)誤差更小。在連續(xù)激光輻照下,含能材料固相區(qū)域的加熱、熔化和蒸發(fā)過(guò)程所需時(shí)間較長(zhǎng),而點(diǎn)火前的氣相劇烈反應(yīng)所需時(shí)間較短。因此,固相模型忽略了氣相的熱流和化學(xué)反應(yīng),仍能給出較為準(zhǔn)確的點(diǎn)火時(shí)刻預(yù)測(cè)值。當(dāng)激光功率密度較高時(shí),點(diǎn)火前的氣相劇烈反應(yīng)所需時(shí)間更短,固相模型預(yù)測(cè)的點(diǎn)火時(shí)刻將更加準(zhǔn)確。

綜上所述,氣相模型計(jì)算過(guò)程復(fù)雜,但能更好描述含能材料的點(diǎn)火過(guò)程,適用于需要了解氣體組分反應(yīng)過(guò)程的情況,如含能材料配方優(yōu)化;固相模型計(jì)算過(guò)程簡(jiǎn)單,能給出較為準(zhǔn)確的點(diǎn)火時(shí)刻預(yù)測(cè)值,但不能描述氣體分解和反應(yīng)過(guò)程,適用于預(yù)估特定熱流下含能材料點(diǎn)火時(shí)刻的情況。

4 結(jié)論

(1)氣相模型和固相模型給出的激光點(diǎn)火時(shí)刻預(yù)測(cè)值非常接近,均與試驗(yàn)值比較吻合。

(2)氣相模型能給出黑索今氣相產(chǎn)物的化學(xué)反應(yīng)動(dòng)力學(xué)過(guò)程,描述的點(diǎn)火過(guò)程更加清晰合理。

[1]Lee Young-joo, Tang Ching-jen, Thomas. A study of the chemical and physical processes governing CO2laser-induced pyrolysis and combustion of RDX[J]. Combustion and Flame, 1999, 117(3): 600-628.

[2]Liau Y C, Kim E S, Yang V. A comprehensive analysis of laser-Induced ignition of RDX monopropellant[J]. Combust. Flame, 2001, 126(3): 1680-1698.

[3]Meredith K V. Ignition modeling of HMX in laser-induced and fast-cookoff environments[D]. PhD. Dissertation, Provo UT: Brigham Young University, 2003.

[4]Davidson J E. Combustion modeling of RDX, HMX and GAP with detailed kinetics[D]. PhD. Dissertation, Provo UT: Brigham Young University, 1996.

[5]Beckstead M W, Puduppakkam K, Thakre P. Modeling of combustion and ignition of solid-propellant ingredients[J]. Progress in Energy and Combustion Science, 2007, 33(6): 497-551.

[6]田占東, 張震宇, 盧芳云, 等. RDX激光點(diǎn)火的一維氣相模型[J]. 高壓物理學(xué)報(bào), 2011, 25(2): 138-142.

[7]張家雷, 劉國(guó)棟, 王偉平, 等. 激光對(duì)碳纖維增強(qiáng)復(fù)合材料的熱燒蝕數(shù)值模擬[J]. 強(qiáng)激光與粒子束, 2013, 25(8): 1888-1892.

[8]Parr T P, Hanson-Parr D M. RDX ignition flame structure[C]//The 27th Symposium (International) on Combustion. Pittsburgh: The Combustion Institute, 1998: 2301-2308.

(編輯:劉紅利)

Comparative analysis of solid-phase and gas-phase reaction model for RDX under laser irradiation

ZHANG Jia-lei1,LIU Cang-li2,TAN Fu-li1, WANG Wei-ping1

(1.Institute of Fluid Physics, CAEP, Mianyang621900,China;2.China Academy of Engineering Physics, Mianyang621900,China)

In order to investigate laser-induced ignition and combustion mechanism of energetic material, solid-phase and gas-phase reaction models for RDX were proposed. Taking into account the detailed chemical reaction of gaseous products, control equations were solved by finite different method and the ignition process was simulated with gas-phase reaction models. Only considering heat transfer and reaction in solid RDX, control equations were solved by finite element method and temperature and ignition time were obtained with solid-phase reaction models. Finally the predictions were compared with experiment data. It was found that ignition times predicted by these two models are all in good agreement with experiment. Though the detailed ignition process can be simulated by gas-phase reaction models, the calculation process is very complex. On the other hand, though the decomposing and chemical reaction of gaseous products can not be described by solid-phase reaction models, the calculation process is very simple and the precision of ignition time can meet the needs of engineering design.

laser;ignition;finite different method;finite element method

2015-06-04;

2016-04-11。

中國(guó)工程物理研究院科學(xué)技術(shù)發(fā)展基金(2014B0401055)。

張家雷(1980—),男,博士生,研究方向?yàn)榧す馀c物質(zhì)相互作用。E-mail:zhangjialei21@126.com

V512

A

1006-2793(2016)05-0638-04

10.7673/j.issn.1006-2793.2016.05.007

猜你喜歡
區(qū)域模型
一半模型
永久基本農(nóng)田集中區(qū)域“禁廢”
分割區(qū)域
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
關(guān)于四色猜想
分區(qū)域
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
基于嚴(yán)重區(qū)域的多PCC點(diǎn)暫降頻次估計(jì)
主站蜘蛛池模板: 91久久国产综合精品女同我| 亚洲精品天堂自在久久77| 亚洲人成网线在线播放va| 国产精品白浆在线播放| 中文天堂在线视频| h网址在线观看| 97国产在线播放| 亚洲精品动漫在线观看| 一本无码在线观看| 国产成人精品亚洲日本对白优播| 五月天丁香婷婷综合久久| 中文纯内无码H| 国产成人精品视频一区视频二区| 国产美女无遮挡免费视频| 欧美成人午夜影院| 亚洲欧美在线精品一区二区| 91激情视频| 天堂久久久久久中文字幕| 欧美激情视频二区| 一区二区午夜| 久热精品免费| 亚洲欧美另类日本| 3p叠罗汉国产精品久久| 91综合色区亚洲熟妇p| 国产一区二区视频在线| 国产AV毛片| 一级黄色欧美| 一级毛片在线免费看| 国产尹人香蕉综合在线电影| 永久免费无码日韩视频| 国产成人综合在线观看| 日韩欧美91| 黄色网页在线播放| 欧美不卡视频在线观看| 日韩欧美综合在线制服| 五月天综合网亚洲综合天堂网| 国产成人精品第一区二区| 亚洲国产日韩在线成人蜜芽| 亚洲无码不卡网| 欧日韩在线不卡视频| 国产欧美日韩在线一区| 亚洲视频免| 成人亚洲国产| 99re在线观看视频| 国产小视频在线高清播放 | 国产69精品久久久久孕妇大杂乱 | 无码中文AⅤ在线观看| 亚洲天堂网视频| 日本欧美视频在线观看| 成人亚洲视频| 九九九九热精品视频| 免费aa毛片| 3344在线观看无码| 97se亚洲综合| 欧美一区二区三区欧美日韩亚洲| 欧美日韩国产精品综合| 国产亚洲欧美日本一二三本道| 国产一级毛片网站| 日本www在线视频| 波多野结衣一二三| 思思热精品在线8| 久久一日本道色综合久久| 永久免费AⅤ无码网站在线观看| 精品国产成人av免费| 高清无码手机在线观看 | 欧美精品不卡| 久久99精品久久久久纯品| 色欲色欲久久综合网| 久久国产乱子伦视频无卡顿| 久久久波多野结衣av一区二区| 国产一在线| 91丨九色丨首页在线播放| 日韩精品久久无码中文字幕色欲| 色哟哟国产精品一区二区| 亚洲综合激情另类专区| 青青操国产| 中文字幕日韩丝袜一区| 人妻丰满熟妇αv无码| 久久综合一个色综合网| 亚洲欧美日韩综合二区三区| 精品久久人人爽人人玩人人妻| 免费国产小视频在线观看|