• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    RM 不穩(wěn)定過程中預(yù)混火焰界面演化及混合區(qū)增長預(yù)測1)

    2020-12-23 01:16:56汪洋董剛
    力學(xué)學(xué)報(bào) 2020年6期
    關(guān)鍵詞:混合區(qū)激波火焰

    汪洋 董剛

    (南京理工大學(xué)瞬態(tài)物理重點(diǎn)實(shí)驗(yàn)室,南京 210094)

    引言

    Richtmyer-Meshkov(RM)不穩(wěn)定是指不同密度且互不浸潤的介質(zhì)間所形成的有一定曲率的界面,在脈沖加速(impulsive acceleration)作用下,其兩側(cè)介質(zhì)發(fā)生混合并導(dǎo)致其混合區(qū)(mixing zone,MZ)增長的行為,其中脈沖加速通常可由激波沖擊來實(shí)現(xiàn).RM不穩(wěn)定的機(jī)制來自斜壓效應(yīng)(baroclinic effect),即,當(dāng)激波與密度界面作用時(shí),界面形成的密度梯度與激波形成的壓力梯度的方向不一致會導(dǎo)致界面上的渦量沉積,在渦量的誘導(dǎo)下界面上的混合區(qū)逐漸增長,后期這種增長還會逐漸發(fā)展為湍流混合.當(dāng)界面本身具有化學(xué)反應(yīng)活性時(shí)(如預(yù)混火焰界面),化學(xué)反應(yīng)會進(jìn)一步增加RM 不穩(wěn)定的復(fù)雜性.作為一種重要而又復(fù)雜的力學(xué)現(xiàn)象,帶化學(xué)反應(yīng)的RM 不穩(wěn)定頻繁地出現(xiàn)在自然現(xiàn)象(如宇宙超新星爆炸[1])和人類生產(chǎn)生活(如航空發(fā)動機(jī)超燃混合[2]、慣性約束聚變[3]、工業(yè)爆炸災(zāi)害[4]等)中.因此,開展帶化學(xué)反應(yīng)的RM 不穩(wěn)定過程的研究,對認(rèn)識、控制和利用這一現(xiàn)象具有非常重要的現(xiàn)實(shí)意義.

    對惰性界面條件下的RM 不穩(wěn)定過程,國內(nèi)外已經(jīng)開展了大量的理論、實(shí)驗(yàn)和數(shù)值模擬研究.界面混合區(qū)增長的理論預(yù)測模型包含了早期的線性模型[5]以及隨后針對激波單次作用下界面RM 不穩(wěn)定非線性發(fā)展過程的漸進(jìn)模型[6]和勢流模型[7-8]等,基于實(shí)驗(yàn)結(jié)果的適于激波多次作用下的界面混合區(qū)增長模型也得到了進(jìn)一步研究[9-11].文獻(xiàn)[12-13]對這些預(yù)測模型進(jìn)行了很好的總結(jié).此外,有學(xué)者也采用實(shí)驗(yàn)測試的方法針對不同入射激波強(qiáng)度和形狀[14-17]、不同初始界面形態(tài)[18-21]條件下,RM 不穩(wěn)定過程的發(fā)展進(jìn)行了細(xì)致的研究.數(shù)值模擬方面,蔣華等[22-23]采用高精度數(shù)值模擬的方法研究了界面初始形態(tài)對混合區(qū)早期發(fā)展過程的影響,并采用理論模型對混合區(qū)的增長進(jìn)行了預(yù)測[23];Schilling 等[24]采用類似的數(shù)值格式模擬了激波及其反射激波作用下的單模形態(tài)界面RM 不穩(wěn)定過程; Thornber 等[25]采用三維數(shù)值模擬,考察了不同初始波長下的界面混合區(qū)的增長行為; Lombardini 等[26]利用三維大渦模擬方法研究了單次激波作用下,激波強(qiáng)度對界面混合區(qū)發(fā)展后期的影響,重點(diǎn)考察了后期界面混合發(fā)展特性;Tritschler 等[27]采用兩種不同的大渦模擬方法研究了激波及其反射波作用下界面混合區(qū)的發(fā)展,對界面發(fā)展后期小尺度湍流混合特性進(jìn)行了詳盡的分析.

    盡管激波與惰性界面相互作用的研究得到了廣泛深入地開展,但反應(yīng)性界面在激波誘導(dǎo)下的發(fā)展過程的研究卻遠(yuǎn)遠(yuǎn)不夠.由于特定形態(tài)的反應(yīng)性界面在實(shí)驗(yàn)控制上比較困難,因此實(shí)驗(yàn)研究幾乎未見報(bào)道,相關(guān)文獻(xiàn)僅限于少量的數(shù)值模擬研究.Khokhlov等[28]采用二維帶化學(xué)反應(yīng)的Navier-Stokes (NS) 方程,數(shù)值模擬了單次激波與單模態(tài)預(yù)混火焰界面的相互作用過程,指出RM 不穩(wěn)定可促進(jìn)燃燒放熱率,并增加火焰的面積;Kilchyk 等[29]采用帶單步化學(xué)反應(yīng)的二維NS 方程,研究了單次激波或膨脹波與單模預(yù)混火焰界面相互作用的過程,指出激波壓縮效應(yīng)和火焰長度對反應(yīng)性RM 不穩(wěn)定有顯著影響.Massa和Jha[30]采用推導(dǎo)得到的二維線化擾動方程,研究了單次激波作用下預(yù)混火焰界面的RM 不穩(wěn)定過程,指出化學(xué)反應(yīng)誘導(dǎo)的壓力梯度所形成的斜壓效應(yīng)可以減少火焰表面的小尺度結(jié)構(gòu)、減弱火焰界面的增長,因而很可能會弱化激波驅(qū)動的小尺度混合過程.Dong 等采用帶單步化學(xué)反應(yīng)的NS 方程對激波及其反射激波與預(yù)混火焰界面的多次作用的二維[31-32]和三維[33]過程開展了系列研究,研究了大尺度流動、化學(xué)反應(yīng)以及小尺度擴(kuò)散等幾個(gè)因素對反應(yīng)性RM 不穩(wěn)定過程的影響,發(fā)現(xiàn)化學(xué)反應(yīng)在一定程度上可以削弱大尺度流動帶來的界面混合效應(yīng),從而證實(shí)了文獻(xiàn)[30]的推測.此外,Attal 和Ramaprabhu[34]還報(bào)道了非預(yù)混單?;鹧娼缑娴腞M 不穩(wěn)定過程,考察了反射激波對界面的影響,其結(jié)果亦表明激波的多次作用明顯增加了燃燒效率,該工作同時(shí)也反映了反應(yīng)性RM不穩(wěn)定研究開始向多樣性發(fā)展.

    上述研究表明,化學(xué)反應(yīng)的存在進(jìn)一步增加了RM 不穩(wěn)定現(xiàn)象的復(fù)雜性.雖然目前已有眾多的理論模型可以預(yù)測激波誘導(dǎo)下惰性界面的混合區(qū)增長過程,但對于反應(yīng)性界面的混合區(qū)在激波誘導(dǎo)下的增長規(guī)律目前并不明確,例如,文獻(xiàn)[31-33]主要關(guān)注的是RM 不穩(wěn)定過程中流動和化學(xué)反應(yīng)的時(shí)間尺度的變化規(guī)律及其對RM 不穩(wěn)定的影響,而化學(xué)反應(yīng)作用下混合區(qū)增長的定量預(yù)測并未得到研究,也未見其它文獻(xiàn)報(bào)道.為此,本文采用高精度數(shù)值模擬的方法,研究平面入射激波及其反射激波與正弦形狀的預(yù)混火焰界面的相互作用過程,重點(diǎn)考察不同化學(xué)反應(yīng)活性條件下界面混合區(qū)的發(fā)展過程,對化學(xué)反應(yīng)條件下的界面混合區(qū)增長速率進(jìn)行預(yù)測,以期為反應(yīng)性RM 不穩(wěn)定的混合區(qū)發(fā)展模型的建立提供參考.

    1 數(shù)值方法

    1.1 控制方程及其計(jì)算方法

    控制方程采用二維帶單步化學(xué)反應(yīng)的Navier-Stokes(NS)方程,表達(dá)如下

    其中

    以上各式中,t代表時(shí)間,x和y代表空間坐標(biāo)長度,ρ代表體系的密度,u和v分別為x和y方向的速度分量,Y代表反應(yīng)物的質(zhì)量分?jǐn)?shù),p為體系壓力,e是單位質(zhì)量的總能,表達(dá)為

    式(8)中,γ 為絕熱指數(shù),Q代表單位質(zhì)量的化學(xué)反應(yīng)放熱.qx和qy分別代表x和y方向的熱流量,表達(dá)為qx=?k?T/?x和qy=?k?T/?x,T是溫度.τij(i=x,y;j=x,y)為黏性應(yīng)力張量分量,表達(dá)為

    其中,δi j為Kronecker 函數(shù).以上各式中出現(xiàn)的輸運(yùn)系數(shù),如導(dǎo)熱系數(shù)k,擴(kuò)散系數(shù)D和運(yùn)動學(xué)黏度υ 可統(tǒng)一表達(dá)為[4]

    本文計(jì)算使用的預(yù)混氣體參照了文獻(xiàn)[35] 中采用的C2H4+3O2+4N2氣體,同時(shí)假定該預(yù)混氣的Lewis數(shù)、Prandtl 數(shù)和Schemidt 數(shù)均為1,則式(10)中的k0,D0和υ0為相同的常數(shù),根據(jù)文獻(xiàn)[4],這些輸運(yùn)系數(shù)取值為k0=D0=υ0=3.2×10?7kg/(s·m·K0.7),n取值為0.7.這些取值能夠保證輸運(yùn)系數(shù)隨溫度的變化關(guān)系很好地符合實(shí)際C2H4+3O2+4N2氣體的這種變化.

    本文采用了基于Arrhenius 形式的單步化學(xué)反應(yīng)描述燃燒現(xiàn)象,因此式(7) 中的反應(yīng)速率˙ω 可以表達(dá)為

    式中β 代表指前因子,Ea為活化能,R代表通用氣體常數(shù).

    對控制方程(1)采用了分裂算法進(jìn)行求解,其中對流通量的空間導(dǎo)數(shù)項(xiàng)?F/?x和?G/?y采用的是LF(Lax-Friedrichs)格式結(jié)合9 階WENO(weighted essentially non-oscillatory)模板的方法求解[36];輸運(yùn)通量的空間導(dǎo)數(shù)項(xiàng)?FD/?x和?GD/?y采用了10 階中心差分的格式進(jìn)行求解; 而化學(xué)反應(yīng)源項(xiàng)S和解U的時(shí)間導(dǎo)數(shù)項(xiàng)?U/?t則采用了3 階Runge-Kutta 方法進(jìn)行求解.上述控制方程及其計(jì)算方法在以前的研究中已經(jīng)得到了廣泛應(yīng)用[31-33].

    1.2 計(jì)算條件與驗(yàn)證

    本文模擬的激波與預(yù)混火焰界面相互作用的計(jì)算構(gòu)型如圖1(a) 所示,計(jì)算域?yàn)槎S矩形形狀.初始時(shí)刻,平面入射激波(incident shock wave,ISW) 位于計(jì)算域距下邊界為l1的位置并開始向上運(yùn)動,在運(yùn)動過程中與波長為λ,初始振幅為a0的正弦形預(yù)混火焰界面(flame interface,FI)發(fā)生相互作用.令火焰界面的下方為密度較大的未燃預(yù)混氣(綠色部分),火焰上方為密度較小的等壓燃燒條件下的已燃?xì)怏w(藍(lán)色部分),則激波從大密度側(cè)向小密度側(cè)越過火焰界面后會形成繼續(xù)上行的激波和向下反射的稀疏波.上行激波行至計(jì)算域頂端的剛性壁面后會反射下行,下行激波則再次和火焰界面發(fā)生相互作用,作用后的激波會繼續(xù)反射向頂端壁面方向運(yùn)動,經(jīng)過再次反射后又會下行與界面發(fā)生作用.這個(gè)過程會一直進(jìn)行下去,直至反射激波強(qiáng)度逐漸減小至激波不再明顯為止.

    圖1 正弦預(yù)混火焰界面的(a)初始形態(tài)和(b)典型的發(fā)展過程.ISW 為入射激波,FI 為火焰界面Fig.1 (a)Initially sinusoidal pattern and(b)developed pattern of flame interface induced by shock wave.ISW is incident shock wave,FI is flame interface

    令初始的正弦形火焰界面數(shù)學(xué)形式如下

    這里a(x)代表正弦火焰界面的流向(y方向)振幅沿x方向的變化.從界面的角度看,在入射激波的作用下,火焰界面經(jīng)RM 不穩(wěn)定過程發(fā)生變化,典型的形態(tài)如圖1(b)所示.可以看到,二維正弦形界面在激波的單次作用下演變?yōu)椤搬?spike)?泡(bubble)”結(jié)構(gòu),其中釘結(jié)構(gòu)比較復(fù)雜,它還包含了其上的“帽(cap)”結(jié)構(gòu)和兩側(cè)的“渦(vortex)” 結(jié)構(gòu).定義“帽” 的頂部到“泡” 的底部的流向長度為界面混合區(qū)寬度WMZ,則該量是評價(jià)界面RM 不穩(wěn)定發(fā)展過程的重要物理量[12,24].

    本文計(jì)算中,選取計(jì)算域流向長度為l1+l2=0.6 m,其中l(wèi)1=0.1 m 以保證整個(gè)計(jì)算過程中入射激波與界面作用后形成的下行稀疏波不會走出下邊界而影響邊界條件; 計(jì)算域的橫向長度為正弦界面波長λ=0.012 8 m.計(jì)算域的左右邊界為周期性邊界條件,上邊界為絕熱剛性無滑移的壁面條件,下邊界為零梯度邊界條件.假設(shè)本文氣體符合量熱完全氣體的假設(shè),并令初始時(shí)刻波前未燃?xì)怏w(圖1(a)中綠色區(qū)域)的初始溫度和初始壓力分別為T0和p0,則波后未燃?xì)?圖1(a)中紅色區(qū)域)的狀態(tài)可由激波關(guān)系式求出.本文界面初始形態(tài)參數(shù)、預(yù)混氣的熱力學(xué)和動力學(xué)參數(shù)的取值[33]均列入表1.

    表1 計(jì)算參數(shù)[33]Table 1 Computational parameters[33]

    此外,在本文計(jì)算中設(shè)初始平面入射激波Mach數(shù)為M0,則根據(jù)一維激波與界面相互作用的數(shù)值模擬,可以確定入射激波作用后的界面Atwood 數(shù)A=(ρ1?ρ2)/(ρ1+ρ2) (這里ρ1和ρ2分別代表界面兩側(cè)重和輕介質(zhì)的密度)和界面兩側(cè)流向速度差?V,以及第一次反射激波作用后的界面Atwood 數(shù)A1和界面兩側(cè)的流向速度差?V1,這些與一維流動有關(guān)的參數(shù)值也列入表1.

    為研究不同預(yù)混氣的化學(xué)活性對RM 不穩(wěn)定發(fā)展的影響,用不同活化能的數(shù)值來表征不同的化學(xué)反應(yīng)活性,活化能的變化的范圍列入表1 中,此外,還計(jì)算了化學(xué)凍結(jié)流條件下的預(yù)混火焰界面的發(fā)展過程.化學(xué)凍結(jié)流是一種保留了預(yù)混火焰熱力學(xué)狀態(tài)但將化學(xué)動力學(xué)過程進(jìn)行凍結(jié)(設(shè)Ea無窮大,即反應(yīng)速率為零)的理想狀態(tài),考慮的目的是為了通過與反應(yīng)流計(jì)算結(jié)果的對比來分析活化能的影響.在本文計(jì)算中,計(jì)算了7 種不同活化能(Ea/RT0=30.0,31.2,31.6,32.0,32.5,38.2,50.0)的反應(yīng)流以及凍結(jié)流(Ea/RT0=∞)的界面混合區(qū)發(fā)展過程.但在后面的分析中,主要是以Ea/RT0=30.0,38.2,50.0 和凍結(jié)流等4 種典型條件的結(jié)果為主,其中活化能Ea/RT0=38.2的算例為基準(zhǔn)算例,其活化能取值來自文獻(xiàn)[37].

    圖2 不同網(wǎng)格尺寸下計(jì)算的混合區(qū)寬度隨時(shí)間的變化Fig.2 Time histories of width of mixing zone with different grid resolutions

    為考察本文計(jì)算結(jié)果的有效性,首先針對Ea=38.2RT0條件下的激波與預(yù)混火焰相互作用過程的數(shù)值模擬進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證.本文采用了均勻正交的正方形網(wǎng)格來劃分計(jì)算區(qū)域,選取了三種不同的網(wǎng)格尺寸(?x=?y=0.1 mm,0.05 mm 和0.025 mm)進(jìn)行激波與正弦形預(yù)混火焰界面相互作用過程的數(shù)值模擬.圖2 給出了3 種網(wǎng)格尺寸下計(jì)算得到的火焰界面發(fā)展過程中混合區(qū)寬度(WMZ,定義見圖1(b))隨時(shí)間的變化過程.圖中時(shí)間范圍I,II 分別為入射激波和第一次反射激波作用后WMZ的變化階段,可以發(fā)現(xiàn)階段II 的初期存在一個(gè)短暫的WMZ下降的過程,這是由于反射激波掠過火焰界面時(shí)對界面的壓縮所導(dǎo)致的.圖中結(jié)果表明,在階段I 和II 中,3 種網(wǎng)格尺寸計(jì)算的結(jié)果均相互吻合較好,但在階段II 之后的后繼反射激波與火焰界面的作用過程中,?x=?y=0.1 mm 尺寸下計(jì)算的結(jié)果要明顯低于?x=?y=0.05 mm 和0.025 mm 兩種尺寸下計(jì)算的結(jié)果,盡管后兩者的計(jì)算結(jié)果之間也略有差別.階段II 后計(jì)算結(jié)果的差異不僅體現(xiàn)了計(jì)算網(wǎng)格尺寸的影響,實(shí)際上也表明火焰界面的失穩(wěn)向著小尺度結(jié)構(gòu)甚至是湍流混合的方向發(fā)展,這超出了本文采用的控制方程和網(wǎng)格尺寸所研究的范圍,因此,本文研究僅考慮階段I 和II 的相互作用過程.根據(jù)網(wǎng)格無關(guān)性測試的結(jié)果,本文采用?x=?y=0.05 mm 的網(wǎng)格尺寸對所有算例開展數(shù)值模擬,每一個(gè)算例總網(wǎng)格數(shù)均為307.2 萬個(gè),為提高計(jì)算效率采用了基于MPI的并行策略編制程序開展計(jì)算.

    2 結(jié)果與討論

    2.1 階段I 的界面演化

    圖3 不同化學(xué)活性下混合區(qū)寬度隨時(shí)間的變化. SI 為入射激波作用階段混合區(qū)增長速率,SII 為反射激波作用階段混合區(qū)增長速率Fig.3 Time histories of width of mixing zone with different chemical activities. SI is growth rate of mixing zone at incident shock wave stage,SII is growth rate of mixing zone at reflected shock wave stage

    圖3 首先給出了3 個(gè)典型活化能(Ea=30.0RT0,38.2RT0和50.0RT0) 的反應(yīng)流以及凍結(jié)流下反應(yīng)區(qū)寬度(WMZ) 隨時(shí)間的變化.可以看到,階段I 中四個(gè)算例在經(jīng)歷了短暫的相同的發(fā)展過程(t< 0.1 ms)之后,其WMZ的增長快慢開始不同,活化能最小(Ea=30.0RT0)的算例增長最快,而活化能最大(Ea=50.0RT0) 以及凍結(jié)流的算例增長最慢.進(jìn)一步觀察可以發(fā)現(xiàn),階段I 的發(fā)展過程中盡管不同化學(xué)反應(yīng)活性下的界面混合區(qū)增長快慢不同,但均呈現(xiàn)線性增長的形式,通過對計(jì)算結(jié)果的線性擬合可以確定其WMZ的增長速率SI.

    在進(jìn)一步定量研究SI之前,首先考察了火焰界面在階段I 的演化規(guī)律和混合區(qū)寬度的增長行為.圖4 顯示了階段I 中4 個(gè)典型時(shí)刻(t=0.2,0.3,0.4 和0.5 ms)下,用反應(yīng)物質(zhì)量分?jǐn)?shù)(Y)標(biāo)記的界面形態(tài)的變化,每個(gè)時(shí)刻下的4 個(gè)子圖從左到右分別代表了活化能為Ea=30.0RT0,38.2RT0和50.0RT0的反應(yīng)流以及凍結(jié)流等4 種不同反應(yīng)活性的情況.在t=0.2 ms時(shí)(圖4(a)),4 種條件下的界面結(jié)構(gòu)均呈現(xiàn)出“釘?泡” 結(jié)構(gòu),釘結(jié)構(gòu)的上方和兩側(cè)也分別呈現(xiàn)出“帽”結(jié)構(gòu)和渦結(jié)構(gòu)(見圖1(b)的說明).但是,當(dāng)反應(yīng)活性強(qiáng)(Ea=30.0RT0)時(shí),渦結(jié)構(gòu)內(nèi)反應(yīng)物幾乎完全消耗,這說明化學(xué)反應(yīng)的作用首先是發(fā)生在渦結(jié)構(gòu)中的,其原因主要是因?yàn)镽M 不穩(wěn)定過程中形成的渦結(jié)構(gòu)使得界面兩側(cè)的未燃?xì)馀c已燃?xì)膺M(jìn)行了充分的摻混,這種強(qiáng)化的傳熱和傳質(zhì)過程進(jìn)一步促進(jìn)了化學(xué)反應(yīng)的進(jìn)行.該時(shí)刻下的結(jié)果還表明,不同反應(yīng)活性下混合區(qū)寬度的變化差異不大.當(dāng)t=0.3 ms 時(shí)(圖4(b)),可以發(fā)現(xiàn)隨著反應(yīng)活性的增強(qiáng),渦結(jié)構(gòu)內(nèi)的混合過程逐漸被化學(xué)反應(yīng)過程所替代,其中,Ea=30.0RT0條件下的反應(yīng)流其渦結(jié)構(gòu)不再明顯,界面上方形態(tài)主要以“釘?帽” 結(jié)構(gòu)為主并明顯向上移動,從而導(dǎo)致其混合區(qū)寬度也較其他幾種情況更長.當(dāng)t=0.4 ms時(shí)(圖4(c)),除了界面上方的“釘?帽?渦” 結(jié)構(gòu)繼續(xù)發(fā)展外,還可以看到界面下方的泡結(jié)構(gòu)隨著反應(yīng)活性的增強(qiáng)也逐漸向下移動,活性越強(qiáng),“泡”結(jié)構(gòu)向下發(fā)展越快,因此導(dǎo)致整個(gè)混合區(qū)寬度隨著活性的增強(qiáng)而逐漸變大,進(jìn)而導(dǎo)致圖3 中階段I 的WMZ變化規(guī)律.到t=0.5 ms 時(shí)(圖4(d)),可以看到,反應(yīng)活性最強(qiáng)(Ea=30.0RT0) 的算例,其“帽” 結(jié)構(gòu)已經(jīng)由化學(xué)反應(yīng)所消耗,因而只留下釘結(jié)構(gòu),這也導(dǎo)致其反應(yīng)區(qū)寬度出現(xiàn)了突然的下降(見圖3);而活性次強(qiáng)(Ea=38.2RT0)的算例,其渦結(jié)構(gòu)不再明顯,但形成的“釘?帽”結(jié)構(gòu)開始使界面混合區(qū)明顯向上發(fā)展,這與Ea=30.0RT0的算例在t=0.3 ms 時(shí)的結(jié)構(gòu)類似.此外,隨著反應(yīng)活性的增強(qiáng),界面下方的泡結(jié)構(gòu)也進(jìn)一步向下發(fā)展.

    圖4 階段I 中不同時(shí)刻的界面形態(tài)Fig.4 Patterns of flame interface at the different instants for stage I

    根據(jù)以上的界面形態(tài)的演化過程,可以發(fā)現(xiàn),階段I 中混合區(qū)寬度隨時(shí)間的變化是由“釘?帽”結(jié)構(gòu)向上發(fā)展以及泡結(jié)構(gòu)向下發(fā)展所共同決定的.化學(xué)反應(yīng)活性對這一過程的影響十分顯著,化學(xué)反應(yīng)越強(qiáng)(活化能越小),“釘?帽”向上發(fā)展以及“泡”向下發(fā)展越顯著,因此使得反應(yīng)區(qū)寬度越大.

    2.2 階段I 混合區(qū)寬度預(yù)測

    2.1 節(jié)首先觀察了不同反應(yīng)活性條件下界面形態(tài)及界面混合區(qū)寬度隨時(shí)間的演化規(guī)律,發(fā)現(xiàn)化學(xué)反應(yīng)對界面上方的“釘?帽”結(jié)構(gòu)以及下方的“泡”結(jié)構(gòu)均有顯著影響.為發(fā)展定量預(yù)測階段I 的混合區(qū)增長速率SI的模型,圖5 進(jìn)一步考察了t=0.4 ms 時(shí)的界面結(jié)構(gòu)及流場特征.圖5(a) 將四種反應(yīng)活性條件下界面輪廓(用Y=0.9 的等值線表示)進(jìn)行了疊加對比,可以發(fā)現(xiàn),界面下方“泡”結(jié)構(gòu)下緣隨著反應(yīng)活性的增加不斷下移,這種下移符合預(yù)混火焰?zhèn)鞑サ奶卣?反應(yīng)活性越大,化學(xué)反應(yīng)速率越大,則火焰?zhèn)鞑ニ俣纫苍酱?因此可以得出結(jié)論:反應(yīng)性預(yù)混火焰界面混合區(qū)寬度的增長其部分原因是由預(yù)混火焰?zhèn)鞑ニ暙I(xiàn)的.根據(jù)經(jīng)典的Frank-Kamenetsky 層流火焰?zhèn)鞑ダ碚揫38],層流預(yù)混火焰的傳播速度可以表達(dá)為

    式中,Sb為“泡” 結(jié)構(gòu)的向下發(fā)展速度,也代表了層流預(yù)混火焰的傳播速度; α 為熱輸運(yùn)系數(shù),α=k0Tn,n=0.7;k0,β,Ea和R的含義及取值見表1;T代表入射波后的已燃區(qū)溫度,本文取T=7.0T0.

    此外,觀察界面上方“釘?帽” 結(jié)構(gòu),可以看出,反應(yīng)活性越強(qiáng),則“帽” 結(jié)構(gòu)的上緣越向上移,但其結(jié)構(gòu)較為復(fù)雜,并不符合預(yù)混火焰的傳播規(guī)律.為此,圖5(b) 進(jìn)一步給出了界面附近的渦量幅值|ω|=|?0.5(?u/?y??v′/?x)|以及速度矢量V=(u,v′)的分布,這里v′=v?vp,vp為入射激波波后流場的平均速度.首先觀察渦量幅值分布可以看出,隨著反應(yīng)活性的增強(qiáng),流場的渦量明顯減弱,這說明化學(xué)反應(yīng)的加入可以抑制渦量的形成,這與之前研究得到的結(jié)論相一致[30].接下來觀察速度矢量的分布可以看到,當(dāng)化學(xué)反應(yīng)較強(qiáng)(如Ea=30.0RT0) 時(shí),釘結(jié)構(gòu)兩側(cè)的渦量變?nèi)?漩渦結(jié)構(gòu)不明顯,因此對入射激波波后流動速度抑制不夠明顯,因此“釘?帽” 結(jié)構(gòu)上方的流體流向速度較大.反之,當(dāng)化學(xué)反應(yīng)較弱(如Ea=50.0RT0) 或化學(xué)反應(yīng)凍結(jié)時(shí),由于釘兩側(cè)存在較強(qiáng)的互為反向的渦旋,渦旋外側(cè)誘導(dǎo)的向下的運(yùn)動速度實(shí)際上抵消了波后流體向上的運(yùn)動.

    圖5 t=0.4 ms 時(shí)不同化學(xué)活性條件下(a)反應(yīng)物質(zhì)量分?jǐn)?shù)Y=0.9 的等值線對比,(b)流場渦量幅度和速度矢量分布對比Fig.5 (a)Comparisons of(a)mass fraction contours with Y=0.9 and of(b)vorticity amplitude and velocity vector among cases with different chemical activities at t=0.4 ms

    為驗(yàn)證上述分析,圖6 給出了計(jì)算域中心(x=λ/2) 處由于渦的旋轉(zhuǎn)所誘導(dǎo)的界面兩側(cè)的速度差?Va(本文稱誘導(dǎo)速度)沿流向的分布,該速度差由界面上、下側(cè)中心處流向速度的差值所確定.圖6 的結(jié)果表明,雖然在較弱活性(Ea=50.0RT0)或凍結(jié)流的條件下,較強(qiáng)的漩渦對在其之間(y=0.031 6 m 左右)的流體所誘導(dǎo)的流向速度差較大,但在“釘?帽” 結(jié)構(gòu)的上方(y>0.032 5 m)其速度差迅速變小至反向流動,說明漩渦外側(cè)的向下流動對速度差的減小起著很大作用.相反,對活性較強(qiáng)的反應(yīng)流(Ea=30.0RT0),由于漩渦強(qiáng)度弱,漩渦外側(cè)誘導(dǎo)的反向流動速度較小,因此渦對之間形成的正向速度差可以一直維持.通過以上分析可知,界面“釘?帽” 結(jié)構(gòu)受化學(xué)反應(yīng)的影響比較復(fù)雜,但是較強(qiáng)的反應(yīng)活性導(dǎo)致了其上方的具有較大的誘導(dǎo)速度差,進(jìn)而使得“釘?帽” 結(jié)構(gòu)向上發(fā)展得更快,由此可將其向上的發(fā)展速度表達(dá)為

    圖6 t=0.4 ms 時(shí)不同化學(xué)活性條件下軸線上流向誘導(dǎo)速度分布Fig.6 Distributions of streamwise induced velocity at centerline for cases with different chemical activities at t=0.4 ms

    式中?V代表一維構(gòu)型下入射激波波后的界面兩側(cè)流向速度差,其值參見表1,δ1為階段I 的經(jīng)驗(yàn)常數(shù).式(14)等號右端exp(?Ea/RT)反映了化學(xué)反應(yīng)活性對界面上緣“釘?帽” 結(jié)構(gòu)發(fā)展的影響,反應(yīng)活性越強(qiáng),該項(xiàng)越大,則Ss越大,因此符合圖5(b)和圖6 的分析.

    本文假設(shè)所考慮的流場具有層流流動特性,因此化學(xué)反應(yīng)和湍流的相互作用可以忽略,故綜合上述分析,可以采用線性疊加的方式給出階段I 中火焰界面的混合區(qū)寬度增長速率SI的表達(dá)式

    式中,SRM代表RM 不穩(wěn)定導(dǎo)致的混合區(qū)寬度增長速率,Sc代表化學(xué)反應(yīng)導(dǎo)致的混合區(qū)寬度增長速率.前者可以采用傳統(tǒng)的各類基于RM 不穩(wěn)定的理論模型進(jìn)行預(yù)測[12-13],后者由于層流流動的假設(shè)仍可采用如下線性疊加的方法進(jìn)行求解

    式中Sb和Ss分別由式(13)和式(14)進(jìn)行求解,其中涉及到的變量僅為活化能Ea,表明了化學(xué)反應(yīng)活性的影響,而其余各常量均可由初始條件和一維激波界面作用得到的參數(shù)來確定.

    注意到式(14)中δ1為經(jīng)驗(yàn)參數(shù),為確定該值,本文對圖3 數(shù)值模擬得到的結(jié)果進(jìn)行線性擬合可以得到SI,令凍結(jié)流的SI=SRM(即凍結(jié)流的WMZ增長速率不考慮化學(xué)反應(yīng),僅由RM 不穩(wěn)定所導(dǎo)致),則由式(15)可以求得數(shù)值模擬條件下的不同Ea的Sc值.圖7 給出了采用數(shù)值模擬得到的Sc隨Ea的變化關(guān)系,可以看出,隨著活化能的減小,Sc呈指數(shù)增大.進(jìn)一步采用式(13)、式(14)和式(16)進(jìn)行計(jì)算,發(fā)現(xiàn)當(dāng)δ1=9.3 時(shí),利用式(16)預(yù)測模型得到的Sc與數(shù)值模擬得到的Sc吻合良好.圖7 中還給出預(yù)測模型中Sb和Ss隨活化能的變化,可以發(fā)現(xiàn),隨著活化能的減小,火焰界面上方“釘?帽” 結(jié)構(gòu)的增長速度比界面下方“泡”的增長速度要更快.

    圖7 數(shù)值模擬與理論預(yù)測的Sc 的比較Fig.7 Comparisons of Sc between numerical and predicted results

    2.3 階段II 界面演化及混合區(qū)寬度分析

    觀察圖3 中階段II(即第一次反射激波作用后的火焰界面發(fā)展階段)可以看出,反射激波作用后,不同化學(xué)反應(yīng)活性的反應(yīng)流和凍結(jié)流條件下,界面混合區(qū)寬度的增長速率(用SII表示)趨于一致,其增長速率也明顯大于階段I 的增長速率.為了解WMZ的增長規(guī)律,圖8 顯示了階段II 的兩個(gè)典型時(shí)刻(t=0.7 ms和0.8 ms)下,不同化學(xué)反應(yīng)活性下界面的發(fā)展形態(tài).圖8(a)的結(jié)果(t=0.7 ms)表明,反射激波與火焰界面作用后,混合區(qū)寬度仍然沿流向持續(xù)增長,當(dāng)反應(yīng)活性較強(qiáng)時(shí),界面呈現(xiàn)出沿流向的細(xì)長結(jié)構(gòu),隨著反應(yīng)活性的下降,界面形狀逐漸變得復(fù)雜,已燃?xì)夂臀慈細(xì)庀嗷交煨纬傻亩S混合區(qū)變得越來越明顯.然而從混合區(qū)寬度的變化來看,不同活性下的混合區(qū)寬度變化并不明顯,即,隨著反應(yīng)活性的增強(qiáng),界面上方的“釘?帽”結(jié)構(gòu)和下方的“泡結(jié)構(gòu)”均同時(shí)下移,表明此時(shí)界面上下存在同樣的預(yù)混火焰?zhèn)鞑ミ^程,兩者傳播方向相同,因此對界面上、下方發(fā)展的影響相互抵消,從而導(dǎo)致混合區(qū)寬度對反應(yīng)活性的依賴不再變得明顯.圖8(b)的結(jié)果(t=0.8 ms)也顯示了類似的規(guī)律.需要注意的是凍結(jié)流條件下,界面形態(tài)開始失去對稱性,這表明到第一次反射激波后期界面的小尺度流動開始增加,界面逐漸向湍流混合的階段發(fā)展,因此,本文的研究僅限于第一次反射激波作用后的階段II.后繼的二次及多次激波與火焰界面的相互作用不在本文研究的范疇.

    圖8 階段II 中不同時(shí)刻的界面形態(tài)Fig.8 Patterns of flame interface at the different instants for stage II

    圖9 進(jìn)一步給出階段II 中t=0.7 ms 時(shí),不同反應(yīng)活性條件下中心線上流向誘導(dǎo)速度?Va沿流向的分布.由結(jié)果也可以看出界面下方反射激波作用之后的區(qū)域(y=0.39 ~0.41 m) 和界面上方的區(qū)域(y> 0.45 m)誘導(dǎo)速度均趨于零,這說明入射激波和反射激波對火焰界面附近的渦旋運(yùn)動的影響相互抵消,因而界面混合區(qū)的發(fā)展不再考慮渦旋運(yùn)動對化學(xué)反應(yīng)的影響,這一點(diǎn)與階段I 中界面混合區(qū)的發(fā)展有所不同.

    圖9 t=0.7 ms 時(shí)不同化學(xué)活性條件下中心線上流向誘導(dǎo)速度分布Fig.9 Distributions of streamwise induced velocity at centerline for cases with different chemical activities at t=0.7 ms

    采用與式(15)相同的形式表達(dá)階段II 的混合區(qū)增長速率SII

    式中變量符號與式(15)相同,下標(biāo)II 代表階段II.根據(jù)前面分析,由于渦旋運(yùn)動不再對火焰界面產(chǎn)生額外影響,界面下方“泡”結(jié)構(gòu)和上方的“釘?帽”結(jié)構(gòu)僅存在相似的火焰向下傳播的過程,因此導(dǎo)致兩者作用相互抵消,即Sc,II=0.因此,混合區(qū)寬度的變化僅由惰性RM 不穩(wěn)定過程所決定,采用反射波后的Mikaelian 模型[8]對其增長速度進(jìn)行預(yù)測,式(17)可表達(dá)為

    式中δ2為階段II 的經(jīng)驗(yàn)常數(shù),?V1和A1的取值見表1.由7 種不同活化能條件下的反應(yīng)流以及凍結(jié)流的數(shù)值模擬結(jié)果可以擬合得到不同的SII值,由此得到的平均值為149.1 m/s,標(biāo)準(zhǔn)偏差為3.2 m/s,將平均增長速度(149.1 m/s)代入式(17)可以得出δ2=0.79.需要說明的是,原始的Mikaelian 模型[8]中經(jīng)驗(yàn)常數(shù)δ2=0.28,但進(jìn)一步的研究也表明[12],δ2的數(shù)值依賴于界面的形態(tài)與結(jié)構(gòu),比如對惰性界面和密度梯度很大的界面,二維和三維計(jì)算采用的δ2分別可取0.55和0.84.本文δ2取0.79 反映了預(yù)混火焰界面與原始Mikaelian 模型中不存在溫差的界面是有所區(qū)別的.

    3 結(jié)論

    采用帶單步反應(yīng)的Navier-Stokes 方程,對入射激波及其反射激波與正弦形狀的預(yù)混火焰界面相互作用過程進(jìn)行了二維數(shù)值模擬,考察了介質(zhì)反應(yīng)活性對界面形態(tài)演化和混合區(qū)增長過程的影響.得到的主要結(jié)論如下:

    (1)正弦形預(yù)混火焰界面在入射激波作用后沿流向不斷拉伸,除常規(guī)RM 不穩(wěn)定誘導(dǎo)的界面混合區(qū)增長外,與化學(xué)反應(yīng)相關(guān)的預(yù)混火焰?zhèn)鞑缑嫦路降摹芭荨苯Y(jié)構(gòu)的發(fā)展起到了促進(jìn)作用,而化學(xué)反應(yīng)與渦結(jié)構(gòu)的相互作用減弱了渦的強(qiáng)度,促進(jìn)了界面上方的流體流向運(yùn)動進(jìn)而使“釘?帽”結(jié)構(gòu)沿流向不斷發(fā)展.化學(xué)活性越強(qiáng),界面的流向增長越快.

    (2) 沿流向發(fā)展的火焰界面在反射激波的作用下,沿流向繼續(xù)發(fā)展,其中界面下方的“泡” 結(jié)構(gòu)和上方“釘?帽” 結(jié)構(gòu)均受到火焰向下傳播的影響,故化學(xué)反應(yīng)的作用抵消,從而導(dǎo)致反射激波作用后界面的發(fā)展基本不依賴于介質(zhì)的化學(xué)反應(yīng)活性.

    (3)根據(jù)上述界面演化機(jī)理,一方面提出了入射激波作用階段,化學(xué)反應(yīng)導(dǎo)致的界面增長速率Sc的表達(dá)式,該表達(dá)式僅依賴于介質(zhì)的反應(yīng)活化能和初始條件,具有較好的通用性;另一方面指出仍可以采用Mikaelian 的模型預(yù)測反射激波作用下火焰界面的混合區(qū)增長速率,但模型常數(shù)要進(jìn)行修正.

    猜你喜歡
    混合區(qū)激波火焰
    支持虛擬車輛輔助假名更新的混合區(qū)位置隱私保護(hù)方案
    《火焰》
    最亮的火焰
    一種基于聚類分析的二維激波模式識別算法
    基于HIFiRE-2超燃發(fā)動機(jī)內(nèi)流道的激波邊界層干擾分析
    漂在水上的火焰
    斜激波入射V形鈍前緣溢流口激波干擾研究
    美國濱海核電廠溫排水混合區(qū)的設(shè)置及啟示
    適于可壓縮多尺度流動的緊致型激波捕捉格式
    考慮邊界反射的河流離岸排放污染混合區(qū)計(jì)算方法
    国产黄色免费在线视频| 亚洲国产精品一区三区| 亚洲美女视频黄频| 在线观看三级黄色| 七月丁香在线播放| 在线观看免费日韩欧美大片 | 午夜激情久久久久久久| 久久久久人妻精品一区果冻| 久久久久人妻精品一区果冻| 亚洲欧美清纯卡通| 亚洲欧美日韩另类电影网站| 99热全是精品| 嘟嘟电影网在线观看| 亚洲三级黄色毛片| 欧美日本中文国产一区发布| 精品国产一区二区三区久久久樱花| 久久久亚洲精品成人影院| 久久99一区二区三区| 九九在线视频观看精品| 久久99一区二区三区| 欧美人与性动交α欧美精品济南到 | 久久久精品免费免费高清| 亚洲,欧美,日韩| 国产精品不卡视频一区二区| 伦理电影大哥的女人| 交换朋友夫妻互换小说| 国产 一区精品| 亚洲美女视频黄频| 亚洲欧美日韩卡通动漫| 成年av动漫网址| 免费看光身美女| 免费人妻精品一区二区三区视频| 亚洲欧洲精品一区二区精品久久久 | 18禁动态无遮挡网站| 欧美老熟妇乱子伦牲交| 亚洲美女黄色视频免费看| av天堂久久9| 亚洲欧美一区二区三区国产| av在线app专区| 亚洲美女搞黄在线观看| 欧美日韩综合久久久久久| 少妇丰满av| 久久综合国产亚洲精品| 欧美成人精品欧美一级黄| 伦精品一区二区三区| 欧美亚洲 丝袜 人妻 在线| 国产亚洲午夜精品一区二区久久| 在线免费观看不下载黄p国产| 欧美97在线视频| 寂寞人妻少妇视频99o| 亚洲一级一片aⅴ在线观看| 亚洲国产最新在线播放| 久久国产精品男人的天堂亚洲 | 色吧在线观看| 人妻人人澡人人爽人人| 大又大粗又爽又黄少妇毛片口| 午夜福利在线观看免费完整高清在| 久久久久久久久大av| 午夜老司机福利剧场| 精品人妻一区二区三区麻豆| 午夜福利网站1000一区二区三区| 视频区图区小说| 国精品久久久久久国模美| 久久久精品区二区三区| 欧美激情国产日韩精品一区| 视频在线观看一区二区三区| 欧美少妇被猛烈插入视频| 最近的中文字幕免费完整| 国产淫语在线视频| 日产精品乱码卡一卡2卡三| av一本久久久久| 熟女人妻精品中文字幕| 国产av一区二区精品久久| 久久热精品热| 丰满乱子伦码专区| 国产一级毛片在线| 九九在线视频观看精品| 97在线视频观看| 最近2019中文字幕mv第一页| 久热这里只有精品99| 美女大奶头黄色视频| 国产精品欧美亚洲77777| 精品人妻熟女av久视频| 99re6热这里在线精品视频| 日本色播在线视频| 性色avwww在线观看| 午夜福利视频在线观看免费| 免费观看在线日韩| 大码成人一级视频| 简卡轻食公司| 亚洲在久久综合| 国产一区二区三区av在线| xxxhd国产人妻xxx| 精品亚洲成国产av| 精品人妻偷拍中文字幕| 亚洲精品国产色婷婷电影| 少妇熟女欧美另类| 日韩熟女老妇一区二区性免费视频| 草草在线视频免费看| 18禁在线无遮挡免费观看视频| 久久精品国产鲁丝片午夜精品| 久久国产精品大桥未久av| 亚洲av.av天堂| 亚洲精品av麻豆狂野| 精品99又大又爽又粗少妇毛片| 久久女婷五月综合色啪小说| 免费播放大片免费观看视频在线观看| 日韩亚洲欧美综合| 亚洲,欧美,日韩| 26uuu在线亚洲综合色| 大片免费播放器 马上看| 久久精品国产a三级三级三级| 久久精品人人爽人人爽视色| 99热国产这里只有精品6| videossex国产| 91成人精品电影| 国产女主播在线喷水免费视频网站| 国产精品免费大片| 亚洲国产成人一精品久久久| .国产精品久久| 国产精品人妻久久久影院| 五月天丁香电影| 精品人妻在线不人妻| 大香蕉久久成人网| videossex国产| 男人添女人高潮全过程视频| 男女啪啪激烈高潮av片| 国产免费视频播放在线视频| videos熟女内射| 日本黄大片高清| 少妇的逼好多水| 在线免费观看不下载黄p国产| 亚洲精品一区蜜桃| 日本黄大片高清| 久久ye,这里只有精品| 亚洲成人一二三区av| 性高湖久久久久久久久免费观看| 色吧在线观看| 婷婷色av中文字幕| 国产成人精品久久久久久| 男女国产视频网站| 乱人伦中国视频| 两个人的视频大全免费| 视频中文字幕在线观看| 久久久久久久久久久丰满| 三级国产精品欧美在线观看| a级毛片免费高清观看在线播放| 亚洲成人一二三区av| 亚洲人成77777在线视频| 国产精品人妻久久久影院| 一级黄片播放器| 九草在线视频观看| √禁漫天堂资源中文www| 国产69精品久久久久777片| 国产精品.久久久| 久久久久视频综合| 999精品在线视频| 久久女婷五月综合色啪小说| 少妇的逼好多水| 国产精品一二三区在线看| 亚洲中文av在线| 亚洲精品色激情综合| 最近的中文字幕免费完整| 欧美日本中文国产一区发布| 一本—道久久a久久精品蜜桃钙片| 日韩 亚洲 欧美在线| 精品少妇黑人巨大在线播放| 国产伦精品一区二区三区视频9| 高清欧美精品videossex| 午夜视频国产福利| 日日撸夜夜添| 麻豆精品久久久久久蜜桃| videosex国产| 国产亚洲欧美精品永久| 亚洲av综合色区一区| 色网站视频免费| 一级毛片电影观看| 乱码一卡2卡4卡精品| 久久久亚洲精品成人影院| 大香蕉97超碰在线| 午夜视频国产福利| 少妇人妻久久综合中文| 九色亚洲精品在线播放| 国产有黄有色有爽视频| 免费观看性生交大片5| 久久久国产欧美日韩av| 蜜臀久久99精品久久宅男| 精品少妇内射三级| 久热这里只有精品99| 日韩伦理黄色片| 色婷婷久久久亚洲欧美| 色5月婷婷丁香| av播播在线观看一区| 18禁在线播放成人免费| 人妻少妇偷人精品九色| 在线观看免费视频网站a站| 亚洲第一区二区三区不卡| 桃花免费在线播放| 久久韩国三级中文字幕| 亚洲国产成人一精品久久久| 国产在线视频一区二区| 男女边摸边吃奶| 国产男女内射视频| 春色校园在线视频观看| 中文乱码字字幕精品一区二区三区| 精品亚洲乱码少妇综合久久| 国产欧美亚洲国产| av在线观看视频网站免费| 国产一区二区三区av在线| 亚洲精品国产色婷婷电影| 久久99热这里只频精品6学生| 在线免费观看不下载黄p国产| 久久久欧美国产精品| 一级,二级,三级黄色视频| 一级毛片电影观看| 人妻系列 视频| 日韩在线高清观看一区二区三区| 这个男人来自地球电影免费观看 | 成人综合一区亚洲| 国产精品秋霞免费鲁丝片| 亚洲国产精品成人久久小说| 国产探花极品一区二区| 色哟哟·www| 三上悠亚av全集在线观看| 美女大奶头黄色视频| 毛片一级片免费看久久久久| 精品人妻一区二区三区麻豆| 大话2 男鬼变身卡| 亚洲精品视频女| 熟女av电影| 成人毛片60女人毛片免费| 啦啦啦在线观看免费高清www| 最近手机中文字幕大全| 日本91视频免费播放| 国产精品一区二区三区四区免费观看| 哪个播放器可以免费观看大片| 亚洲成人手机| 又黄又爽又刺激的免费视频.| 欧美xxⅹ黑人| 国产精品免费大片| 高清欧美精品videossex| 国产欧美日韩一区二区三区在线 | 亚洲,一卡二卡三卡| 国产一区有黄有色的免费视频| 中文天堂在线官网| 免费av中文字幕在线| 国产成人aa在线观看| 国产日韩欧美在线精品| 精品国产露脸久久av麻豆| 日韩av免费高清视频| 蜜桃国产av成人99| 成人免费观看视频高清| 久久99精品国语久久久| 精品久久久精品久久久| 亚洲人成网站在线观看播放| 汤姆久久久久久久影院中文字幕| 欧美日韩视频高清一区二区三区二| 一级黄片播放器| 国产在线视频一区二区| 免费播放大片免费观看视频在线观看| 日日爽夜夜爽网站| 亚洲国产精品国产精品| 插逼视频在线观看| 2018国产大陆天天弄谢| 97精品久久久久久久久久精品| kizo精华| 亚洲天堂av无毛| 男男h啪啪无遮挡| 97在线人人人人妻| 永久免费av网站大全| 欧美另类一区| 亚洲精品中文字幕在线视频| 岛国毛片在线播放| 日韩不卡一区二区三区视频在线| 人妻人人澡人人爽人人| 欧美97在线视频| 欧美日韩成人在线一区二区| 日日爽夜夜爽网站| 国产乱人偷精品视频| 国产免费一级a男人的天堂| 男女高潮啪啪啪动态图| 熟妇人妻不卡中文字幕| av.在线天堂| 久久ye,这里只有精品| 纯流量卡能插随身wifi吗| 日本黄色日本黄色录像| 3wmmmm亚洲av在线观看| 国产av国产精品国产| 飞空精品影院首页| 欧美xxⅹ黑人| 一级二级三级毛片免费看| av线在线观看网站| 亚洲一区二区三区欧美精品| 91在线精品国自产拍蜜月| 欧美成人午夜免费资源| videosex国产| 大码成人一级视频| 国产视频内射| 人人妻人人澡人人爽人人夜夜| 国产一区亚洲一区在线观看| 97精品久久久久久久久久精品| 久久鲁丝午夜福利片| 日日撸夜夜添| 下体分泌物呈黄色| 亚洲久久久国产精品| 美女国产视频在线观看| 中文乱码字字幕精品一区二区三区| 97超碰精品成人国产| 亚洲欧美成人综合另类久久久| 在现免费观看毛片| 赤兔流量卡办理| 国产熟女午夜一区二区三区 | 欧美日本中文国产一区发布| 日韩中文字幕视频在线看片| 久久人人爽人人片av| 五月开心婷婷网| 国产精品免费大片| 午夜精品国产一区二区电影| 国产欧美日韩一区二区三区在线 | 内地一区二区视频在线| 欧美人与性动交α欧美精品济南到 | 婷婷色麻豆天堂久久| 中文乱码字字幕精品一区二区三区| 99视频精品全部免费 在线| 午夜福利视频精品| .国产精品久久| 免费观看无遮挡的男女| 日韩熟女老妇一区二区性免费视频| 日本av手机在线免费观看| 亚洲精品一二三| 九色成人免费人妻av| 国产成人精品一,二区| 国产欧美日韩综合在线一区二区| 国产69精品久久久久777片| freevideosex欧美| 国产亚洲一区二区精品| 只有这里有精品99| 亚洲精品成人av观看孕妇| 亚洲av欧美aⅴ国产| 免费日韩欧美在线观看| 亚洲欧美精品自产自拍| 精品一区二区三区视频在线| 在线看a的网站| 插阴视频在线观看视频| 久久久久久久久久人人人人人人| 午夜免费鲁丝| 久久鲁丝午夜福利片| 日韩成人av中文字幕在线观看| 日日摸夜夜添夜夜添av毛片| 中文字幕最新亚洲高清| 国产毛片在线视频| 国产黄片视频在线免费观看| 欧美+日韩+精品| 丰满迷人的少妇在线观看| 亚洲三级黄色毛片| 夫妻午夜视频| 80岁老熟妇乱子伦牲交| 精品久久国产蜜桃| 久久久精品免费免费高清| 一级a做视频免费观看| 一级毛片aaaaaa免费看小| 精品午夜福利在线看| 老女人水多毛片| 十分钟在线观看高清视频www| 新久久久久国产一级毛片| 久久久久久伊人网av| av黄色大香蕉| 18禁裸乳无遮挡动漫免费视频| 亚洲少妇的诱惑av| 五月玫瑰六月丁香| 国产色爽女视频免费观看| 成人国语在线视频| 乱人伦中国视频| 如日韩欧美国产精品一区二区三区 | 免费人成在线观看视频色| 大码成人一级视频| 久久久精品免费免费高清| 免费av中文字幕在线| 国产午夜精品一二区理论片| 国产精品欧美亚洲77777| 人妻 亚洲 视频| 爱豆传媒免费全集在线观看| 国产成人精品一,二区| 国产男人的电影天堂91| 午夜免费观看性视频| 中文精品一卡2卡3卡4更新| 国产精品麻豆人妻色哟哟久久| 日本vs欧美在线观看视频| 精品99又大又爽又粗少妇毛片| 男人操女人黄网站| 国产视频内射| 欧美日韩精品成人综合77777| 亚洲欧美一区二区三区国产| 亚洲av国产av综合av卡| 蜜桃国产av成人99| 精品国产一区二区久久| 精品亚洲成a人片在线观看| 亚洲精品乱久久久久久| 亚洲精品久久成人aⅴ小说 | 久久国产精品男人的天堂亚洲 | 国产成人aa在线观看| 日韩av在线免费看完整版不卡| 一边摸一边做爽爽视频免费| 欧美日韩成人在线一区二区| 国产一级毛片在线| 国产一区二区在线观看av| 精品亚洲成a人片在线观看| 国产欧美另类精品又又久久亚洲欧美| 午夜免费观看性视频| 国产亚洲最大av| 在线播放无遮挡| 18禁在线播放成人免费| 少妇高潮的动态图| 美女cb高潮喷水在线观看| 男女免费视频国产| 国产 精品1| 日产精品乱码卡一卡2卡三| 欧美3d第一页| 国产伦精品一区二区三区视频9| videos熟女内射| 在线观看三级黄色| 亚洲精品一二三| xxxhd国产人妻xxx| 人体艺术视频欧美日本| 曰老女人黄片| 国产午夜精品久久久久久一区二区三区| 91久久精品国产一区二区成人| 日韩,欧美,国产一区二区三区| 亚洲,一卡二卡三卡| .国产精品久久| 欧美 亚洲 国产 日韩一| 久久人妻熟女aⅴ| av女优亚洲男人天堂| 亚洲av国产av综合av卡| 亚洲国产成人一精品久久久| 亚洲综合色惰| 日韩一区二区视频免费看| 亚洲av电影在线观看一区二区三区| 伊人久久精品亚洲午夜| 一级毛片我不卡| 在线看a的网站| 少妇 在线观看| 欧美日韩视频精品一区| 精品人妻一区二区三区麻豆| 久久久国产一区二区| 97超碰精品成人国产| 久久精品久久精品一区二区三区| 黑丝袜美女国产一区| videossex国产| 国产免费一级a男人的天堂| 亚洲国产色片| 亚洲欧美精品自产自拍| 国内精品宾馆在线| 亚洲精品,欧美精品| 狠狠精品人妻久久久久久综合| 女人精品久久久久毛片| 亚洲精品久久成人aⅴ小说 | 久久毛片免费看一区二区三区| 观看美女的网站| 蜜桃在线观看..| 国产淫语在线视频| 日韩制服骚丝袜av| 午夜精品国产一区二区电影| 日韩伦理黄色片| 亚洲精品国产av蜜桃| 亚洲av国产av综合av卡| 午夜激情av网站| 满18在线观看网站| 一区二区三区免费毛片| 日韩在线高清观看一区二区三区| 大片电影免费在线观看免费| 国产亚洲av片在线观看秒播厂| 汤姆久久久久久久影院中文字幕| 曰老女人黄片| 免费高清在线观看日韩| 高清不卡的av网站| 欧美老熟妇乱子伦牲交| 亚洲国产欧美在线一区| 中文字幕人妻熟人妻熟丝袜美| 日本午夜av视频| 少妇熟女欧美另类| 亚洲av国产av综合av卡| 国国产精品蜜臀av免费| 啦啦啦啦在线视频资源| 男女啪啪激烈高潮av片| 亚洲欧美清纯卡通| 亚洲人成77777在线视频| 搡老乐熟女国产| 两个人免费观看高清视频| 亚州av有码| 乱人伦中国视频| 国产 一区精品| 色婷婷久久久亚洲欧美| 欧美日韩一区二区视频在线观看视频在线| 日韩中文字幕视频在线看片| 丁香六月天网| 国产国语露脸激情在线看| 在现免费观看毛片| 精品国产露脸久久av麻豆| 丝袜美足系列| 亚洲av二区三区四区| 欧美日韩视频高清一区二区三区二| 亚洲国产毛片av蜜桃av| 国产精品人妻久久久影院| 精品亚洲乱码少妇综合久久| 国产欧美另类精品又又久久亚洲欧美| 日韩精品免费视频一区二区三区 | 22中文网久久字幕| 在线免费观看不下载黄p国产| 国产永久视频网站| 亚洲人与动物交配视频| 又大又黄又爽视频免费| av专区在线播放| 日日爽夜夜爽网站| 欧美xxⅹ黑人| 美女视频免费永久观看网站| 成人国产av品久久久| 男女无遮挡免费网站观看| 天美传媒精品一区二区| 国产成人精品久久久久久| 免费看av在线观看网站| 亚洲精品色激情综合| 亚洲,欧美,日韩| 我要看黄色一级片免费的| 汤姆久久久久久久影院中文字幕| 亚洲少妇的诱惑av| 成人漫画全彩无遮挡| 人体艺术视频欧美日本| 久久99热6这里只有精品| 亚洲av.av天堂| 国产精品不卡视频一区二区| 一级毛片我不卡| 只有这里有精品99| 少妇精品久久久久久久| 亚洲情色 制服丝袜| 中文字幕人妻丝袜制服| 亚洲精品中文字幕在线视频| 哪个播放器可以免费观看大片| 内地一区二区视频在线| 最近最新中文字幕免费大全7| 在线观看免费视频网站a站| 亚洲精品av麻豆狂野| 欧美精品一区二区大全| 蜜桃久久精品国产亚洲av| 国产精品欧美亚洲77777| 成人亚洲精品一区在线观看| 亚洲精品第二区| 日本vs欧美在线观看视频| 91精品一卡2卡3卡4卡| 国产在视频线精品| 国产黄频视频在线观看| 在线观看一区二区三区激情| 久久精品久久精品一区二区三区| 亚洲精品aⅴ在线观看| 欧美日韩亚洲高清精品| √禁漫天堂资源中文www| 人人妻人人澡人人爽人人夜夜| 免费黄色在线免费观看| 色吧在线观看| 免费看光身美女| 曰老女人黄片| 亚洲国产av影院在线观看| 秋霞伦理黄片| 热re99久久国产66热| 18在线观看网站| 中文字幕人妻熟人妻熟丝袜美| 黑人欧美特级aaaaaa片| 在线看a的网站| 18禁观看日本| 狂野欧美激情性xxxx在线观看| 九九在线视频观看精品| 蜜桃国产av成人99| 午夜久久久在线观看| 看十八女毛片水多多多| 国产乱来视频区| 国产一区亚洲一区在线观看| 国产成人精品福利久久| 国产成人aa在线观看| 亚洲av.av天堂| 一个人看视频在线观看www免费| 国产精品久久久久久精品电影小说| 少妇人妻 视频| 热99久久久久精品小说推荐| 日韩免费高清中文字幕av| 一区在线观看完整版| 精品人妻在线不人妻| 最后的刺客免费高清国语| 日本av免费视频播放| h视频一区二区三区| 亚洲成人手机| 成人手机av| 精品一品国产午夜福利视频| 国模一区二区三区四区视频| 制服诱惑二区| 成人午夜精彩视频在线观看| 国产免费一区二区三区四区乱码| 亚洲美女黄色视频免费看| 一级毛片aaaaaa免费看小| 99九九在线精品视频| 亚洲熟女精品中文字幕| 日韩欧美一区视频在线观看| 国产免费一区二区三区四区乱码| 国产精品成人在线| 在线观看免费视频网站a站| 一级毛片我不卡| 精品久久久久久久久亚洲| 日韩欧美一区视频在线观看| 一本久久精品| 少妇丰满av| 国产伦理片在线播放av一区| 美女国产高潮福利片在线看| 少妇人妻久久综合中文| 伊人久久精品亚洲午夜| 在线免费观看不下载黄p国产| 精品久久蜜臀av无|