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

    磁控條件下激波沖擊三角形氣柱過程的數(shù)值研究?

    2018-11-28 10:40:14董國丹張煥好林震亞秦建華陳志華郭則慶沙莎
    物理學(xué)報 2018年20期
    關(guān)鍵詞:氣柱渦量重質(zhì)

    董國丹 張煥好 林震亞 秦建華 陳志華 郭則慶 沙莎

    1)(南京理工大學(xué),瞬態(tài)物理國家重點實驗室,南京 210094)

    2)(北京電子工程總體研究所,北京 100854)

    本文基于磁流體動力學(xué)方程組,在保證磁場散度為零的條件下,采用CTU+CT(corner transport upwind+constrained transport)算法,對有無磁場控制下激波與重質(zhì)或輕質(zhì)三角形氣柱相互作用過程進(jìn)行數(shù)值研究.結(jié)果表明:無論有無磁場,兩氣柱在激波沖擊下均具有完全不同的波系結(jié)構(gòu)和射流現(xiàn)象.其中,入射激波與重氣柱發(fā)生常規(guī)折射,形成介質(zhì)射流,而與輕氣柱作用則發(fā)生非常規(guī)折射,形成反相空氣射流.無磁場時,氣柱在激波沖擊下,產(chǎn)生Richtmyer-Meshkov和Kelvin-Helmholtz不穩(wěn)定性,界面出現(xiàn)次級渦序列,重氣柱上下角卷起形成主渦對,輕氣柱空氣射流穿過下游界面后形成偶極子渦.施加橫向磁場后,次級渦序列、主渦對以及偶極子渦均消失.進(jìn)一步研究表明,在磁場作用下,洛倫茲力將不穩(wěn)定性誘導(dǎo)產(chǎn)生的渦量向界面兩側(cè)的Alfvén波上輸運,減少界面渦量沉積,抑制界面卷起失穩(wěn).最終,渦量沿界面兩側(cè)形成相互遠(yuǎn)離的渦層,界面不穩(wěn)定性得到控制.此外,定量分析表明磁場能加快兩氣柱上游界面的運動,抑制下游界面的運動,且對輕氣柱的控制效果更好.

    1 引 言

    激波與不同流體分界面相互作用的過程中,誘導(dǎo)界面微弱擾動從線性發(fā)展為非線性,此現(xiàn)象廣泛存在于兵器發(fā)射、航空、航天、天體力學(xué)、地質(zhì)物理、核能以及化工等領(lǐng)域,其蘊含了湍流轉(zhuǎn)捩以及激波與渦相互作用等復(fù)雜的流體物理現(xiàn)象[1?3].Richtmyer基于Taylor線性理論,對激波沖擊界面、誘導(dǎo)界面失穩(wěn)的現(xiàn)象進(jìn)行研究,提出了預(yù)測界面增長的脈沖加速模型,該預(yù)測隨后被Meshkov的激波管實驗所證實.因此激波沖擊氣體分界面,進(jìn)而誘導(dǎo)界面失穩(wěn)的現(xiàn)象被稱為Richtmyer-Meshkov(R-M)不穩(wěn)定性[4,5].

    自R-M不穩(wěn)定性提出以來,人們對其進(jìn)行了大量的研究.Rudinger和Somers[6]實驗研究了激波分別沖擊小尺度H2,He與SF6氣柱的宏觀運動過程,提出了一個預(yù)測界面渦量和速度的簡單理論模型.隨后,Haas和Sturtevant[7]對激波與球形和柱形He,R22氣柱進(jìn)行研究,他們結(jié)合光學(xué)知識進(jìn)行分析,發(fā)現(xiàn)激波經(jīng)過重質(zhì)R22氣柱后,在氣柱內(nèi)部匯聚,而后經(jīng)過輕質(zhì)He氣柱后發(fā)散,但對于激波與兩種介質(zhì)氣柱相互作用過程中具體波系發(fā)展沒有給出詳細(xì)的說明.Layes等[8,9]采用高速紋影圖像技術(shù)對激波沖擊輕質(zhì)He、重質(zhì)Kr以及等質(zhì)N2氣柱進(jìn)行了實驗研究,發(fā)現(xiàn)在輕氣柱上游頂點會形成反向空氣射流,重氣柱下游界面則會形成介質(zhì)射流,而等質(zhì)N2氣柱界面形狀無明顯變化,但是對于兩種完全不同射流形成機理沒有給出明確的解釋.Ranjan等[10,11]分析激波與不同馬赫數(shù)He球形氣泡作用過程中的復(fù)雜波系、界面發(fā)展和渦量產(chǎn)生等復(fù)雜物理現(xiàn)象,發(fā)現(xiàn)作用過程中會產(chǎn)生不同的折射現(xiàn)象,這與激波入射角和介質(zhì)有關(guān).

    早期的研究多局限于球形和柱形界面,但是激波與界面作用過程中的折射現(xiàn)象與入射角密切相關(guān).近年來,人們開始對不同幾何外形氣柱的波系結(jié)構(gòu)和界面發(fā)展進(jìn)行深入研究.Zhai等[12]實驗研究了激波與正方形、長方形,三角形、菱形輕質(zhì)(He)氣柱作用過程中復(fù)雜波系結(jié)構(gòu)以及界面擾動發(fā)展過程,發(fā)現(xiàn)激波與輕氣柱發(fā)生非常規(guī)折射,并出現(xiàn)非常規(guī)折射之間的轉(zhuǎn)換.Luo等[13]對激波與正方形、長方形、三角形、菱形重質(zhì)(SF6)氣柱進(jìn)行研究,發(fā)現(xiàn)不同幾何外形氣柱具有不同的波系結(jié)構(gòu)和射流現(xiàn)象.但關(guān)于重和輕三角形界面不穩(wěn)定性仍需更加詳細(xì)的對比研究,特別是輕氣柱后期出現(xiàn)的偶極子渦.Dong等[14]對激波與120?及60?V形界面作用過程進(jìn)行實驗研究,結(jié)果表明不同激波角能誘導(dǎo)不同強度的斜壓機制,最終誘導(dǎo)界面產(chǎn)生不同的渦量運動速度,且在60?界面時所形成的渦更多更清晰.沙莎等[15?17]詳細(xì)研究了激波與球形R22氣柱及梯形SF6氣柱作用中射流形成機理與及界面演化過程,發(fā)現(xiàn)折射激波聚焦于氣柱左側(cè)內(nèi)部而引導(dǎo)射流,同時氣柱上下兩側(cè)大量環(huán)境氣體被卷吸進(jìn)渦核中.

    此外,因Rayleigh-Taylor(R-T)不穩(wěn)定性常伴隨著R-M不穩(wěn)定性出現(xiàn)[18],它們在天體物理、大氣物理、慣性約束核聚變和高速飛行器等領(lǐng)域同樣有著非常重要而廣泛的應(yīng)用和研究意義[19].然而,這些領(lǐng)域中物質(zhì)多呈現(xiàn)第四態(tài)(即等離子狀態(tài)),所以利用磁場來控制界面不穩(wěn)定性的研究成為近年來的熱點和難點.李源和羅喜勝[20]研究了磁場對R-T不穩(wěn)定性的作用,發(fā)現(xiàn)磁場的非線性項對R-T不穩(wěn)定性有較大影響.磁流體動力學(xué)(magnetohydrodynamics,MHD)將流體力學(xué)和電磁學(xué)結(jié)合以描述導(dǎo)電流體在電磁場中的運動,因此本文基于MHD方程組研究磁場對R-M不穩(wěn)定性的影響.R-M不穩(wěn)定性中MHD效應(yīng)的研究起源于磁頂層動力學(xué)[21].實驗中要獲得穩(wěn)定的等離子體比較困難,因此數(shù)值模擬在研究MHD效應(yīng)對R-M不穩(wěn)定性的影響中具有十分重要的價值.Samtaney[22]通過對非線性理想磁流體進(jìn)行二維數(shù)值模擬,研究了MHD效應(yīng)對斜平面R-M不穩(wěn)定性的影響,并認(rèn)為磁場可有效控制界面R-M不穩(wěn)定性,且當(dāng)磁場方向與激波運動方向一致時,磁場作用于界面斜壓渦量,使界面渦量向周圍擴(kuò)散,減少界面上渦量的堆積,從而抑制界面的不穩(wěn)定性.Wheatley等[23?25]基于不可壓線性MHD方程組研究發(fā)現(xiàn)當(dāng)磁場方向與激波的運動方向垂直時,磁場通過干擾間斷處激波折射過程,減少界面渦量沉積從而抑制不穩(wěn)定性的發(fā)展.隨后,Sano等[26]對平行磁場控制R-M不穩(wěn)定性的臨界磁場強度進(jìn)行研究,發(fā)現(xiàn)臨界磁場強度與入射激波馬赫數(shù)密切相關(guān),并給出了MHD中控制R-M不穩(wěn)定性臨界條件的公式.Cao等[27]通過理論推導(dǎo),發(fā)現(xiàn)在橫向磁場下,洛倫茲力作用于界面,抑制界面不穩(wěn)定性的增長.Mostert等[28]研究了磁場對球形和柱形爆炸中R-M不穩(wěn)定性的影響.Lin等[29,30]研究了非理想環(huán)境下磁場對重質(zhì)氣團(tuán)爆炸的影響,發(fā)現(xiàn)磁場可有效控制R-M不穩(wěn)定性,且在磁場作用下整個高密度氣團(tuán)會被壓縮,同時電阻和雙極擴(kuò)散效應(yīng)會衰弱磁場的作用,但雙極擴(kuò)散效應(yīng)還會增大磁壓力作用范圍.

    上述實驗研究、數(shù)值模擬和理論預(yù)測都表明斜壓渦量的積累是R-M不穩(wěn)定性產(chǎn)生的主要原因,斜壓渦量是因密度和壓力梯度的不一致(斜壓效應(yīng))而產(chǎn)生.因此對重和輕氣柱R-M不穩(wěn)定性的對比研究十分必要.Dong等[14]的研究證明60?入射角時,界面產(chǎn)生的渦對更多更清晰,因此本文選取封閉的正三角形氣柱進(jìn)行研究.目前MHD效應(yīng)對R-M不穩(wěn)定性的研究多局限于單模重質(zhì)界面,對于封閉重質(zhì)和輕氣柱的研究仍欠缺.此外,雖然Samtaney[22]對法向磁場作用下MHD效應(yīng)、對45?斜平面重質(zhì)界面R-M不穩(wěn)定性進(jìn)行了研究,但是在橫向磁場作用下MHD效應(yīng)對重質(zhì)和輕質(zhì)兩種界面的作用仍需完善.因此,為了更好地闡明不同介質(zhì)氣柱波系的發(fā)展和演化以及磁控下氣柱形態(tài)演化,本文選取正三角形重質(zhì)和輕氣柱進(jìn)行研究.采用CTU+CT(corner transport upwind+constrained transport)算法求解MHD方程組[29?31],其中CTU算法用于計算多維積分,CT算法用于保證磁場散度為零.分別對有無磁場情況下激波沖擊重質(zhì)與輕質(zhì)三角形氣柱的過程進(jìn)行數(shù)值研究,詳細(xì)分析了兩種氣柱具體波系發(fā)展以及界面演化,并對磁場控制界面不穩(wěn)定性發(fā)展的機理進(jìn)行探討.最后,定量分析了磁場控制作用下重質(zhì)和輕質(zhì)兩種界面不穩(wěn)定性的發(fā)展規(guī)律.

    2 計算方法和模型

    本文基于非理想MHD方程組[29,30],采用非分裂的 CTU+CT[31,32]算法進(jìn)行求解. 其中,CTU(corner transport upwind)算法基于PPM(piecewise parabolic method)對守恒量進(jìn)行重構(gòu),是一種迎風(fēng)格式的有限體積積分算子.另外,為保證磁場散度為零,在計算Godunov通量時,結(jié)合用于計算電場的CT(constrained transport)算法進(jìn)行通量重構(gòu).

    圖1為平面激波沖擊三角形氣柱的計算模型.計算域為[?0.01,0.29]× [?0.06,0.06],其長和寬分別為D=0.3 m,H=0.12 m,等邊三角形氣柱的邊長為L=0.07 m.初始時,三角形左頂點與左邊界的水平距離為0.01 m,平面入射激波位于左邊界上,并隨后自左向右傳播.重與輕質(zhì)氣柱內(nèi)分別充滿R22和75%(質(zhì)量分?jǐn)?shù),下同)He+25%Air的氣體,氣柱外為空氣,內(nèi)外壓力均為1 atm(1 atm=1.013×105Pa).相應(yīng)氣體參數(shù)參照表1.

    為了更好地反映磁場對氣柱變形的控制,計算域的上、下及右邊界設(shè)為固壁.計算域采用均勻分布的笛卡爾網(wǎng)格,經(jīng)網(wǎng)格無關(guān)性測試后,總網(wǎng)格數(shù)為1500×600.本文中,設(shè)初始橫向磁場強度分別為B=0 T和B=0.01 T,分別對應(yīng)有控和無控條件,其方向平行于y軸.另外,為了使氣體受磁場影響,假設(shè)氣體已經(jīng)電離.根據(jù)文獻(xiàn)[30,31],取等離子體初始電導(dǎo)率為107S/m,熱導(dǎo)率為1.4 W/(m·K?1),黏性系數(shù)為3.72×10?5Pa/s,霍爾系數(shù)為10?7m3/C,雙極擴(kuò)散系數(shù)隨粒子數(shù)密度ρ的增加而減小,因此雙極擴(kuò)散系數(shù)為(10?5/ρ)× 10?5m2·Pa·s?1.

    圖1 模型示意圖Fig.1.Schematic of computational model.

    圖2 激波沖擊重質(zhì)和三角形氣柱過程中計算與實驗[12,13]紋影圖的對比(ds,繞射激波;s1,激波1;tp,三波點;ts,透射激波;ms,馬赫桿;sts,二次透射激波) (a)重氣柱;(b)輕氣柱Fig.2.Comparison of numerical and experimental[12,13]schlieren images of the interactions between shock waves and heavy and light triangular cylinders(ds,dif f racted shock;s1,shock wave1;tp,triple point;ts,transmitted shock;ms,Mach stem;sts,secondary transmitted shock):(a)Heavy cylinders;(b)light cylinders.

    表1 氣體參數(shù)Table 1.Gas parameters

    3 結(jié)果與討論

    3.1 算例驗證

    圖2為激波沖擊重質(zhì)(a)與輕質(zhì)(b)三角形氣柱的計算紋影圖(右)與實驗紋影圖(左)[12,13]的對比,驗證算例的條件均與相應(yīng)的實驗條件一致.可見,本文數(shù)值方法捕捉到的入射、反射、折射和透射等復(fù)雜波系以及氣柱形態(tài)的演化均與相應(yīng)的實驗結(jié)果符合.另外,本文結(jié)果能更清晰地反映出重氣柱內(nèi)三波點處的滑移線以及輕氣柱外三波點處的滑移線.

    3.2 無磁場時重氣柱

    圖3為無磁場控制(B=0 T)下,激波與重氣柱作用過程的計算紋影圖.圖4為不同時刻流場壓力與速度分布圖,其中A為氣柱左頂點,B為繞射激波聚焦點,C1,C2為氣柱上下角.由圖3(a)可見,入射激波is沖擊氣柱時,其中間段因入射角為零,在氣柱左頂點發(fā)生透射,在氣柱內(nèi)產(chǎn)生透射激波ts.且因重質(zhì)氣體聲阻抗大,ts的傳播速度小于is.同時,入射激波is的上下兩段分別在氣柱上下界面發(fā)生正規(guī)折射(regular refraction,RRR),形成弧形反射激波crs與折射激波ras.隨后,氣柱內(nèi)ras和ts相互作用產(chǎn)生激波s1和兩個三波點tp以及三波點處的滑移線.

    入射激波is與弧形反射激波crs向下游傳播的過程中,持續(xù)與氣柱上下界面發(fā)生作用,誘導(dǎo)界面出現(xiàn)R-M不穩(wěn)定性(圖3(b)—(e)).另外,因重質(zhì)氣體的阻擋作用,氣柱外側(cè)氣體波后速度大于氣柱內(nèi)波后速度(圖4(a)),K-H不穩(wěn)定性開始出現(xiàn).因此在斜壓效應(yīng)與K-H不穩(wěn)定性的共同作用下,氣柱左頂點卷起形成向外翻轉(zhuǎn)的渦對,如圖3(b)所示.因此氣柱上下界面在R-M與K-H不穩(wěn)定性的共同作用下,開始失穩(wěn)并逐漸卷起形成珠狀小渦序列(圖 3(b)—(f)).

    當(dāng)is繞過C1與C2后,形成繞射激波ds,其一端在下游界面與ras相連,另一端則與is相連(圖3(b)),且繞射過程中還會形成向上游傳播的膨脹扇ps.同時,折射激波ras作用于下游界面,在氣柱內(nèi)產(chǎn)生反射稀疏波rrs(圖3(b)).由于C1與C2處存在較大的剪切速度(圖4(b)),因此氣柱上下角在斜壓效應(yīng)和K-H不穩(wěn)定的共同作用下卷起形成渦對,并最終發(fā)展成為兩個主渦(圖3(b)—(g)).

    隨著ts向下游傳播(圖3(a)—(c)),其長度不斷增大,同時折射激波ras不斷縮短并在t=230μs時消失.當(dāng)ts穿過下游界面時,分別形成向下游與上游傳播的弧形透射激波cts與弱反射激波ws(圖3(d)).此時rrs穿過s1,并繼續(xù)向氣柱中心運動(圖3(d)),同時其波后產(chǎn)生局部低壓區(qū)(圖4(b)).

    t=240μs時,上下兩條繞射激波ds在中心軸上聚焦并反射,形成反射激波drs(圖3(e)),同時cts被分成上下兩段(圖3(e)).繞射激波ds聚焦(圖3(e))時,聚焦點B后產(chǎn)生局部高壓區(qū),對下游界面造成瞬間的強沖擊,導(dǎo)致碰撞區(qū)界面的流場速度瞬間上升,誘導(dǎo)氣柱界面的局部變形(圖3(f)),并逐漸形成一道向下游傳播的蘑菇狀射流(圖3(g)).

    圖3 無磁場時,激波與重氣柱作用過程的計算紋影圖 (brs,尾壁反射激波;crs,弧形反射激波;cts,弧形透射激波;is,入射激波;ps,膨脹扇;ras,折射激波;rrs,反射稀疏波;ws,弱反射激波) (a)t=110μs;(b)t=190μs;(c)t=220μs;(d)t=230μs;(e)t=240μs;(f)t=330μs;(g)t=1080μs;(h)t=1900μsFig.3.Schlieren image sequences of interaction between shock wave and heavy cylinder without a magnetic f i eld(brs,back wall ref l ected shock;crs,curved ref l ected shock;cts,curved transmitted wave;is,incident shock;ps,expansion fan;ras,refracted shock;rrs,rarefaction wave;ws,weak shock):(a)t=110μs;(b)t=190μs;(c)t=220μs;(d)t=230μs;(e)t=240μs;(f)t=330μs;(g)t=1080μs;(h)t=1900μs.

    圖4 不同時刻壓力(左)和速度(右)云圖 (a)t=90μs;(b)t=240μsFig.4.Cloud chat of pressure(left)and velocity(right)at dif f erent time:(a)t=90μs;(b)t=240μs.

    由圖3(f)可知,氣柱內(nèi)的s1,rrs與ws在氣柱上下界面間來回碰撞,在氣柱內(nèi)形成復(fù)雜的波系結(jié)構(gòu).期間,復(fù)雜波系與界面上小渦序列多次相互作用,加速了界面的失穩(wěn)而形成復(fù)雜的渦結(jié)構(gòu)(圖3(f),(g)),從而加劇了重質(zhì)氣體與環(huán)境氣體的混合.此外,在t=1080μs后 (圖3(g)),入射激波在右邊界反射而回傳的尾壁反射激波brs到達(dá)氣柱尾部,隨后再次與氣柱發(fā)生作用,加劇了氣柱界面的失穩(wěn),使氣柱形狀更為復(fù)雜,如圖3(h)所示.

    3.3 無磁場時輕氣柱

    圖5為無磁場時,激波與輕氣柱作用過程的計算紋影圖.圖6為不同時刻,流場速度分量云圖.由圖5(a)可知,入射激波is沖擊氣柱后,在氣柱內(nèi)產(chǎn)生透射激波ts1.由于輕氣柱聲阻抗小,ts1的傳播速度大于is,在氣柱外產(chǎn)生自由前導(dǎo)激波fps.該自由前導(dǎo)激波fps在界面外與入射激波is發(fā)生自由前導(dǎo)折射(free precursor refraction,FPR)[13,34],伴隨著一系列反射稀疏波rrw,形成向上游傳播的弧形反射激波crs、與界面相連的激波sk以及滑移線sl.

    以復(fù)合改性植物膠為稠化劑,有機絡(luò)合物為交聯(lián)劑,引入抑制劑進(jìn)一步控制緩交時間,并添加少量表面活性劑提高解堵后破膠液的排液能力,形成了暫堵壓井膠塞配方,見表1。

    因輕氣柱聲阻抗小,is沖擊氣柱左頂點后,壓縮頂點向下游運動,并在頂點處形成一個向下游運動的高速區(qū),如圖6(a)所示.該高速區(qū)持續(xù)作用在氣柱上游頂點上,最終形成一個反相空氣射流結(jié)構(gòu)(圖 5(a)—(h)).

    ts1穿過下游界面時,形成向下游傳播的二次透射激波ts2與向上游傳播的弱激波ws,此時fps繞過氣柱上下角與ts2相連(圖5(b)).在t=170 μs時,弱激波ws到達(dá)氣柱左頂點形成的空氣射流渦環(huán)處,并與射流渦環(huán)相互作用.但因射流渦環(huán)上下渦核的旋向相反,弱激波被剪斷成3段(圖5(c)).隨著fps離開界面,FPR會轉(zhuǎn)變?yōu)橛蒮ps,is,激波sn,馬赫桿ms以及三波點tp組成的雙馮·諾依曼折射(twin von Neumann refraction,TNR)[13,34](圖5(d)).

    t=190μs時,反射激波crs到達(dá)計算域的上下邊界,并發(fā)生碰撞反射,形成的反射激波bs向氣柱中心軸運動.在t=260μs時,反射激波再次跟氣柱作用,但因氣柱內(nèi)外聲阻抗不同,壁面反射激波bs被剪成3段(圖5(e)),其中氣柱內(nèi)段傳播較快.氣柱內(nèi)復(fù)雜波系在t=800μs時基本耗散.t=1050μs時,入射激波在右邊界反射而回傳的尾壁反射激波brs沖擊氣柱.由此可見,相對于重氣柱,激波與輕氣柱作用過程中激波結(jié)構(gòu)及界面演變完全不一樣,重氣柱內(nèi)的激波呈匯聚趨勢,而輕氣柱內(nèi)激波發(fā)散.

    由圖5(d)—(f)可知,氣柱左頂點處的射流結(jié)構(gòu)形成后,向下游運動的同時不斷將氣柱上下界面卷入到射流的反向渦環(huán)內(nèi),使渦環(huán)尺度持續(xù)增大同時,在反相空氣射流的拖動下,界面在y方向產(chǎn)生較大的剪切速度(圖6(b)),因此在K-H不穩(wěn)定作用下,射流尾部氣柱界面開始失穩(wěn)并卷起渦串(圖5(e),(f)).隨著界面的發(fā)展,射流渦前端面與下游界面的距離越來越近(圖5(d)—(f)),并在t=800μs時穿透下游界面,形成向下游運動的偶極子渦[33],同時下游界面幾乎被分割成對稱的兩半,如圖5(f)所示.隨后,尾壁反射激波brs沖擊氣柱,進(jìn)一步加強了氣柱界面的失穩(wěn),使界面上的大渦結(jié)構(gòu)破碎形成復(fù)雜的湍流渦結(jié)構(gòu)(圖5(h)).

    圖5 無磁場情況下,激波與輕氣柱作用過程的密度紋影圖(bs,壁面反射激波;fps,自由前導(dǎo)激波;rrw,反射稀疏波;sk,激波k;sl,滑移線;sn,激波n;ts1,一次透射激波;ts2,二次透射激波) (a)t=110μs;(b)t=150μs;(c)t=170μs;(d)t=190μs;(e)t=260μs;(f)t=800μs;(g)t=1050μs;(h)t=1900μsFig.5.Density schlieren image sequences interaction between shock wave and light cylinder without a magnetic fi eld(bs,wall re fl ected shock wave;fps,free precursor shock;rrw,re fl ected rare shock;sk,shock wave k;sl,slip line;sn,shock wave n;ts1, fi rstly transmitted shock;ts2,secondary transmitted shock):(a)t=110μs;(b)t=150μs;(c)t=170μs;(d)t=190μs;(e)t=260μs;(f)t=800μs;(g)t=1050μs;(h)t=1900μs.

    圖6 t=90μs和t=600μs速度分量云圖 (a)x方向的速度分量;(b)y方向的速度分量Fig.6.Velocity components cloud chart at t=90 μs and t=600 μs:(a)Velocity component of x direction;(b)velocity component of y direction.

    圖7 為兩氣柱中軸線上壓力和速度曲線圖,其中實線為重氣柱,虛線為輕氣柱.結(jié)合圖3與圖7可知,重氣柱相當(dāng)于一個凸透鏡,使激波在內(nèi)部匯聚,形成高壓區(qū);而輕氣柱相當(dāng)于一個凹透鏡,使激波發(fā)散,因而輕氣柱內(nèi)部壓力不高.對于重氣柱,繞射激波ds聚焦時(t=240μs),在中軸線上產(chǎn)生320 kPa的高壓,膨脹波rrs氣柱內(nèi)產(chǎn)生175 kPa的相對低壓(圖7(a)),氣柱界面兩邊極大的壓差導(dǎo)致下游界面形成介質(zhì)射流.而輕氣柱,其左頂點在激波撞擊后產(chǎn)生高速向下游運動,速度為280 m/s,從而導(dǎo)致反向空氣射流的形成.此外,重氣柱內(nèi)波系匯聚且持續(xù)與氣柱界面作用,所以R-M不穩(wěn)定是重氣柱失穩(wěn)的主要原因,這也與Luo等[13]的結(jié)果一致.對于輕氣柱,雖然波系在氣柱內(nèi)發(fā)射,但入射激波作用下產(chǎn)生的高速射流頭部會在界面兩側(cè)產(chǎn)生較大的剪切效應(yīng),因此K-H不穩(wěn)定性在輕氣柱失穩(wěn)過程中占主導(dǎo)作用.

    圖7 不同時刻,氣柱對稱軸上壓力和速度沿軸線的分布(H重氣柱;L輕氣柱) (a)壓力;(b)速度Fig.7.Distributions of pressure and velocity along symmetric axis at dif f erent time(H,heavy cylinders;L,light cylinders):(a)Pressure;(b)velocity.

    3.4 有磁場時界面演化

    圖8 為加入B=0.01 T橫向磁場后,激波與兩種介質(zhì)三角形氣柱作用過程的計算紋影圖.其中上圖為重氣柱,下圖為輕氣柱.在磁場作用下,流場波系結(jié)構(gòu)及演化與無磁場時一致(圖3和圖5),因而不再綴述且無磁場時的數(shù)值結(jié)果與相應(yīng)的實驗結(jié)果符合,可以保證本文算法對磁場作用下激波與界面作用過程模擬的準(zhǔn)確性.由圖8可知,在磁場作用下,重氣柱除了上下角因斜壓效應(yīng)而卷起形成的渦以及下游邊界中心形成的射流結(jié)構(gòu)外,其余界面保持光滑,射流渦的表面也變得光滑(圖8(a1)—(d1));輕氣柱空氣射流頭部變得光滑,偶極子渦對也消失(圖8(a2)—(d2)).

    圖8 磁場控制下(B=0.01 T),激波與重質(zhì)(上)和輕質(zhì)(下)氣柱作用過程的密度紋影圖 (a1)t=240μs;(b1)t=400μs;(c1)t=1200 μs;(d1)t=1900 μs;(a2)t=150 μs;(b2)t=600 μs;(c2)t=1050 μs;(d2)t=1900 μsFig.8.Density schlieren image sequences of the shock wave interaction with heavy(upper)and light(lower)cylinder under magnetic f i eld control(B=0.01 T):(a1)t=240μs;(b1)t=400μs;(c1)t=1200μs;(d1)t=1900μs;(a2)t=150μs;(b2)t=600μs;(c2)t=1050μs;(d2)t=1900μs.

    圖9(a)和圖9(b)為有無磁場情況下,氣柱界面渦量的分布云圖,其中上圖為重氣柱 (400μs),下圖為輕氣柱(200μs).重氣柱密度梯度垂直界面指向下游,壓力梯度垂直波陣面指向上游,在斜壓效應(yīng)的作用下,重氣柱上界面附近流體順時針旋轉(zhuǎn)生成負(fù)渦量,相應(yīng)的下界面處流體逆時針旋轉(zhuǎn)生成正渦量,而輕氣柱則相反,具體渦旋方向如圖1所示.由圖9可知,無磁場時,斜壓渦量沉積在界面上(圖9(a)),并不斷誘導(dǎo)界面卷起小渦序列進(jìn)而失穩(wěn).加入橫向磁場后,渦量不再沉積于界面,而是在界面兩側(cè)形成相互遠(yuǎn)離的渦層(圖9(b)).

    圖9 渦量和洛倫茲力云圖(上圖:重氣柱,t=400μs;下圖:輕氣柱,t=200μs)(a)Vorticity(B=0 T);(b)vorticity(B=0.01 T);(c)Lorentz force(y);(d)Lorentz force(x)Fig.9.Vorticities and Lorentz force distribution(upper,heavy cylinder,t=400μs;lower,light cylinder t=200μs):(a)Vorticity(B=0 T);(b)vorticity(B=0.01 T);(c)Lorentz force(y);(d)Lorentz force(x).

    圖10 兩氣柱界面渦旋方向和洛倫茲力f的方向 (a)重氣柱;(b)輕氣柱Fig.10.Directions of rotating vorticities and Lorentz forces f on the cylinder interfaces:(a)Heavy cylinder;(b)light cylinder.

    為了進(jìn)一步從理論上分析磁場作用下渦量不能沉積在界面的機理,以間斷面為參考坐標(biāo)系,考慮動量和切向電場的守恒關(guān)系式[35]:

    其中方括號代表界面兩側(cè)的差值;下標(biāo)n和t分別代表法向和切向分量;ρ代表密度;Bn,Bt分別代表法向和切向磁場強度;vn,vt分別代表法向和切向速度.連續(xù)間斷面上vn=0,且界面上Bn=0,由(1),(2)式可得[vt]=[Bt]=0.因此連續(xù)間斷面上不能存在渦層和電流層[26].盡管激波沖擊氣柱的瞬間會在氣柱界面上產(chǎn)生渦量,隨后該渦量也會在磁場作用下遠(yuǎn)離界面.

    Samtaney等[22?25]對界面附近激波折射過程進(jìn)行詳細(xì)的研究后發(fā)現(xiàn),在磁場作用下,激波折射過程會產(chǎn)生MHD波.Samtaney[22]對平面激波沖擊45?斜平面的研究表明,激波沖擊界面產(chǎn)生的反射激波在磁場作用下會轉(zhuǎn)變?yōu)橐粚β瓷浜涂旆瓷浯怕暡?一種MHD波),透射激波轉(zhuǎn)變?yōu)橐粚β干浜涂焱干浯怕暡?另外,渦量在磁場作用下不再沉積于界面,而是分布在慢反射和慢透射MHD波上.本文兩種介質(zhì)封閉三角形界面較Wheatley和Samtaney的單模重質(zhì)界面更為復(fù)雜,產(chǎn)生的波系結(jié)構(gòu)和射流現(xiàn)象也復(fù)雜多變,因而更具有研究前景和價值.從圖11的磁力線和密度分布圖(100μs)可知,在激波干擾下,界面附近磁力線扭曲.一方面扭曲的磁力線會產(chǎn)生阿爾文波(一種慢MHD橫波).另一方面扭曲的磁力線會增大平行于界面的磁分量,從而在界面兩側(cè)產(chǎn)生較大的洛倫茲力[27].最終界面渦量在洛倫茲力的作用下形成相互遠(yuǎn)離的渦層(圖9(b)),該渦層將附于Alfvén波上,并隨著Alfvén波逐漸遠(yuǎn)離界面,最終界面因缺少渦量沉積而穩(wěn)定.

    圖11 兩氣柱磁力線和密度分布圖 (a)重氣柱;(b)輕氣柱Fig.11.Spatial distribution of magnetic f i eld lines and density of two gas cylinders:(a)Heavy cylinder;(b)light cylinder.

    圖12 激波與氣柱作用過程中氣柱上下游界面運動情況 (a)重氣柱;(b)輕氣柱Fig.12.Movement of upper and lower interfaces during the interaction of a planar shock wave and cylinder:(a)Heavy cylinder;(b)light cylinder.

    圖12 為氣柱界面與左邊界距離隨時間的變化曲線.無磁場控制時,兩氣柱上游界面(Lx)均以恒定的速度向下游運動[13,14],但因輕質(zhì)氣體的聲阻抗小,其速度更快.下游界面(Rx)則在入射激波到達(dá)前靜止,隨后同樣以恒定的速度向下游運動[12,13].加入磁場后(BLx與BRx),BLx與Lx基本重合.然而,磁場會阻止氣柱上下角處渦結(jié)構(gòu)的進(jìn)一步卷起,因此在上下角卷起渦結(jié)構(gòu)后(t>400μs),BRx與Rx的差距逐漸加大.線圖中拐點(t≈1100μs)為尾壁反射激波沖擊氣柱時刻,可見磁場使拐點后下游界面基本不動.由此可知,磁場可有效控制氣柱的運動,且對輕氣柱界面運動具有更好的控制效果.

    圖13和圖14則分別為重與輕氣柱縱向和橫向尺寸隨時間的變化,圖中拐點(t≈1100μs)為尾壁反射激波沖擊氣柱所致.對于重氣柱(圖13),磁場可有效抑制界面兩個特征尺度的變形.但對于輕氣柱(圖14),磁場能有效控制縱向尺寸L的變化,卻并不能控制橫向尺寸H的變化.由此可知,磁場對重氣柱變形的控制效果更好.

    圖13 激波與重氣柱作用過程中長度L和高度H的變化情況 (a)縱向長度L;(b)橫向高度HFig.13.Development of length L and height H of the triangular interface during the interaction between shock waves and heavy cylinder:(a)Length L;(b)height H.

    4 結(jié) 論

    本文基于非理想MHD方程組,為保證磁場散度為零,采用非分裂的CTU+CT算法分別對有無磁場情況下激波沖擊重質(zhì)和輕質(zhì)三角形氣柱的過程進(jìn)行數(shù)值研究.得到以下結(jié)論.

    無論有無磁場,激波與重質(zhì)或輕氣柱相互作用的過程中均具有不同的波系結(jié)構(gòu)和射流現(xiàn)象.具體而言,入射激波與重氣柱作用時發(fā)生常規(guī)折射,氣柱相當(dāng)于一個凸透鏡,使復(fù)雜波系在氣柱內(nèi)匯聚;入射激波與輕氣柱作用時則發(fā)生非常規(guī)折射,氣柱相當(dāng)于一個凹透鏡,使復(fù)雜波系經(jīng)過氣柱后發(fā)散.另外,繞射激波在重氣柱中心軸上聚焦產(chǎn)生局部高壓,誘導(dǎo)重氣柱下游界面中心處出現(xiàn)介質(zhì)射流;而入射激波與輕氣柱作用時,輕氣柱左頂點卷起形成空氣射流.

    無磁場時,重氣柱上下角卷起形成主渦對,輕氣柱空氣射流會穿過下游界面形成相互分離的偶極子渦.同時,兩氣柱界面均出現(xiàn)大量次級渦序列.此外,對于重氣柱,因氣柱內(nèi)波系匯聚并持續(xù)與界面作用,R-M不穩(wěn)定性在其失穩(wěn)過程中占主導(dǎo)地位;而對于輕氣柱,后期氣柱內(nèi)波系耗散,射流頭部運動較快,K-H不穩(wěn)定性在失穩(wěn)過程中占主導(dǎo)地位.

    施加磁場后,洛倫茲力作用于界面渦量,將渦量向界面兩側(cè)輸運,減少界面渦量沉積,抑制界面卷起失穩(wěn).最終界面兩側(cè)會形成兩個相互遠(yuǎn)離的渦層,該渦層附于Alfvén波上,并隨著Alfvén波遠(yuǎn)離界面.最終,在磁場作用下,重氣柱的主渦對、輕氣柱的偶極子渦以及兩氣柱界面次級渦序列也消失,界面不穩(wěn)定性得到控制.最后定量分析表明,磁場能有效控制氣柱界面的運動,且對輕氣柱界面運動具有更好的控制效果.此外,磁場可控制重氣柱縱向與橫向兩個特征尺寸的增長,卻不能控制輕氣柱的橫向尺寸的增長.

    猜你喜歡
    氣柱渦量重質(zhì)
    理想氣體與液柱類問題的處理方法探究
    預(yù)壓縮氣墊包裝系統(tǒng)靜力及動力學(xué)特性研究
    包裝工程(2024年1期)2024-01-20 06:17:38
    激波誘導(dǎo)雙層氣柱演化的偏心效應(yīng)研究
    氣體物理(2022年2期)2022-03-31 12:49:16
    含沙空化對軸流泵內(nèi)渦量分布的影響
    國家公園建設(shè)重質(zhì)不重量
    綠色中國(2019年17期)2019-11-26 07:04:42
    自由表面渦流動現(xiàn)象的數(shù)值模擬
    重質(zhì)高酸原油高效破乳劑研究
    航態(tài)對大型船舶甲板氣流場的影響
    The application of numerical simulation of delta wing with blunt leading edge using RANS/LES hybrid method
    重質(zhì)純堿不同生產(chǎn)工藝的對比
    ponron亚洲| 97人妻精品一区二区三区麻豆| 亚洲精品日韩在线中文字幕| 最近中文字幕2019免费版| 国产永久视频网站| 久久久久久久久久人人人人人人| 国产午夜精品一二区理论片| 国产在线一区二区三区精| 又爽又黄a免费视频| 伊人久久国产一区二区| 中文字幕人妻熟人妻熟丝袜美| 国产色爽女视频免费观看| 十八禁国产超污无遮挡网站| 精品久久久久久久久久久久久| 国产淫语在线视频| 精品国内亚洲2022精品成人| 亚洲图色成人| 精品人妻一区二区三区麻豆| 亚洲国产成人一精品久久久| 国产精品伦人一区二区| 九九爱精品视频在线观看| 亚洲一区高清亚洲精品| 韩国高清视频一区二区三区| 亚洲欧美中文字幕日韩二区| 在线 av 中文字幕| 日韩av在线大香蕉| 国产一区二区亚洲精品在线观看| 一级毛片 在线播放| 亚洲av成人av| 色吧在线观看| 丰满乱子伦码专区| 免费观看无遮挡的男女| 精品久久久精品久久久| 人妻少妇偷人精品九色| 干丝袜人妻中文字幕| 伊人久久国产一区二区| 日韩一区二区视频免费看| 十八禁网站网址无遮挡 | 亚洲精品国产av蜜桃| 日韩 亚洲 欧美在线| 亚洲va在线va天堂va国产| 亚洲欧美清纯卡通| 99热6这里只有精品| 久久精品久久久久久久性| 亚洲成人一二三区av| 少妇人妻精品综合一区二区| 亚洲欧美清纯卡通| 中文乱码字字幕精品一区二区三区 | 3wmmmm亚洲av在线观看| 久久久成人免费电影| 久久人人爽人人爽人人片va| 国产亚洲91精品色在线| 哪个播放器可以免费观看大片| 日韩一本色道免费dvd| 男女啪啪激烈高潮av片| av在线蜜桃| 久久久午夜欧美精品| 亚洲欧美精品专区久久| 亚洲国产av新网站| 亚洲欧洲日产国产| 国精品久久久久久国模美| 亚洲最大成人手机在线| 日本欧美国产在线视频| 啦啦啦啦在线视频资源| av线在线观看网站| 99久久精品一区二区三区| 白带黄色成豆腐渣| 日韩欧美 国产精品| 最近中文字幕2019免费版| av福利片在线观看| 22中文网久久字幕| 午夜福利高清视频| 纵有疾风起免费观看全集完整版 | 男的添女的下面高潮视频| 久久久久精品性色| 一夜夜www| 精品国产露脸久久av麻豆 | 在线免费观看不下载黄p国产| 国产爱豆传媒在线观看| 久久久久久久国产电影| 精品久久久久久电影网| 少妇高潮的动态图| 成人鲁丝片一二三区免费| 综合色av麻豆| 久久久午夜欧美精品| 成人午夜精彩视频在线观看| 午夜激情欧美在线| 午夜福利在线在线| 啦啦啦啦在线视频资源| www.色视频.com| 丝瓜视频免费看黄片| 亚洲av一区综合| 天堂中文最新版在线下载 | 亚洲激情五月婷婷啪啪| 嘟嘟电影网在线观看| 中国国产av一级| av在线蜜桃| 可以在线观看毛片的网站| 日韩欧美精品免费久久| 亚洲精品成人久久久久久| 国产激情偷乱视频一区二区| 天天躁日日操中文字幕| 亚洲精品自拍成人| 九草在线视频观看| 九九在线视频观看精品| 日日啪夜夜撸| 国产亚洲5aaaaa淫片| 精品国产三级普通话版| 亚洲内射少妇av| 91在线精品国自产拍蜜月| 中国国产av一级| 精品一区在线观看国产| 免费观看性生交大片5| 一个人观看的视频www高清免费观看| 国产老妇伦熟女老妇高清| 观看免费一级毛片| 免费看光身美女| 边亲边吃奶的免费视频| 午夜福利网站1000一区二区三区| av在线亚洲专区| 国产午夜精品一二区理论片| 亚洲,欧美,日韩| 国产乱人偷精品视频| 午夜激情福利司机影院| 超碰97精品在线观看| 成人亚洲精品av一区二区| 国产亚洲av片在线观看秒播厂 | 蜜桃久久精品国产亚洲av| 老司机影院成人| 国产老妇伦熟女老妇高清| 日本黄色片子视频| 久久久久久久亚洲中文字幕| 国产亚洲精品av在线| 人人妻人人澡人人爽人人夜夜 | 高清午夜精品一区二区三区| 伊人久久精品亚洲午夜| 国产伦在线观看视频一区| 99视频精品全部免费 在线| 91av网一区二区| 男女边吃奶边做爰视频| 欧美日韩国产mv在线观看视频 | 婷婷六月久久综合丁香| 听说在线观看完整版免费高清| 午夜免费激情av| 欧美高清性xxxxhd video| 日韩成人伦理影院| 麻豆久久精品国产亚洲av| 深爱激情五月婷婷| 中国国产av一级| 日韩欧美精品免费久久| 久久久a久久爽久久v久久| 欧美激情久久久久久爽电影| 寂寞人妻少妇视频99o| 26uuu在线亚洲综合色| 色吧在线观看| 成人漫画全彩无遮挡| 欧美成人a在线观看| 丝袜美腿在线中文| 国产一级毛片七仙女欲春2| 国产精品一区二区性色av| 熟女电影av网| 全区人妻精品视频| 亚洲精品自拍成人| 亚洲伊人久久精品综合| 色吧在线观看| 亚洲熟妇中文字幕五十中出| 爱豆传媒免费全集在线观看| 菩萨蛮人人尽说江南好唐韦庄| 亚洲自偷自拍三级| 国产成人精品婷婷| 亚洲内射少妇av| 一本久久精品| av线在线观看网站| 亚洲性久久影院| 国产高潮美女av| 日本与韩国留学比较| 亚洲成人中文字幕在线播放| 97超视频在线观看视频| 国产av国产精品国产| 国产成人aa在线观看| 夜夜爽夜夜爽视频| 人体艺术视频欧美日本| 久久这里只有精品中国| 欧美一区二区亚洲| 成人性生交大片免费视频hd| 国产在线一区二区三区精| 久久久久国产网址| 神马国产精品三级电影在线观看| 午夜激情久久久久久久| 国产在视频线精品| 别揉我奶头 嗯啊视频| 国产成人freesex在线| 亚洲精品视频女| 日韩中字成人| 九草在线视频观看| 久久久午夜欧美精品| 国产成人精品婷婷| 免费观看的影片在线观看| 插阴视频在线观看视频| av黄色大香蕉| 搞女人的毛片| 永久免费av网站大全| 国内揄拍国产精品人妻在线| 国产高清有码在线观看视频| 亚洲一区高清亚洲精品| 18+在线观看网站| freevideosex欧美| 男女边吃奶边做爰视频| videossex国产| 精品人妻一区二区三区麻豆| 九草在线视频观看| 日韩伦理黄色片| 视频中文字幕在线观看| 日韩 亚洲 欧美在线| 日韩电影二区| 22中文网久久字幕| 国产黄a三级三级三级人| 国产精品久久久久久精品电影小说 | 午夜福利网站1000一区二区三区| 久久午夜福利片| 国产在线男女| 国产一区亚洲一区在线观看| 亚洲图色成人| 91av网一区二区| 好男人在线观看高清免费视频| 伊人久久精品亚洲午夜| 搡老乐熟女国产| 精品人妻一区二区三区麻豆| 精品久久久久久久久亚洲| 好男人在线观看高清免费视频| 免费看光身美女| 国产黄a三级三级三级人| 亚洲美女搞黄在线观看| 国产乱人视频| 亚洲四区av| 不卡视频在线观看欧美| 亚洲18禁久久av| 极品教师在线视频| 国产高潮美女av| 欧美成人a在线观看| 国产精品不卡视频一区二区| av免费在线看不卡| 免费看av在线观看网站| av天堂中文字幕网| 一级毛片 在线播放| 国产片特级美女逼逼视频| 国产亚洲午夜精品一区二区久久 | 亚洲精品久久久久久婷婷小说| 亚洲自偷自拍三级| 亚洲av免费高清在线观看| 国产免费又黄又爽又色| 人人妻人人澡人人爽人人夜夜 | 亚洲国产最新在线播放| 精品人妻一区二区三区麻豆| 国产高清不卡午夜福利| 亚洲国产av新网站| 99久久精品国产国产毛片| 精品久久久噜噜| 久久精品人妻少妇| 亚洲精品国产成人久久av| 激情五月婷婷亚洲| 日本黄色片子视频| 国产精品99久久久久久久久| 国内少妇人妻偷人精品xxx网站| 97人妻精品一区二区三区麻豆| 日韩欧美三级三区| 麻豆精品久久久久久蜜桃| 欧美不卡视频在线免费观看| 免费av毛片视频| 啦啦啦中文免费视频观看日本| 免费黄网站久久成人精品| 99久久九九国产精品国产免费| 26uuu在线亚洲综合色| 搡老妇女老女人老熟妇| 婷婷色av中文字幕| 亚洲色图av天堂| 日本色播在线视频| 午夜亚洲福利在线播放| 国产午夜精品论理片| 亚洲精品成人久久久久久| 毛片一级片免费看久久久久| 中文字幕人妻熟人妻熟丝袜美| 肉色欧美久久久久久久蜜桃 | 国产乱人视频| 亚洲欧洲日产国产| 禁无遮挡网站| 秋霞伦理黄片| 日韩欧美精品v在线| 女人十人毛片免费观看3o分钟| 欧美3d第一页| 国产伦精品一区二区三区四那| 久久99精品国语久久久| 大片免费播放器 马上看| 国产乱人偷精品视频| 免费观看a级毛片全部| 熟女人妻精品中文字幕| 国产熟女欧美一区二区| 久久综合国产亚洲精品| 日韩一区二区视频免费看| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 欧美3d第一页| 久久热精品热| 国内精品宾馆在线| av福利片在线观看| 亚洲欧美日韩卡通动漫| 日韩成人伦理影院| 亚洲国产欧美人成| 一级爰片在线观看| 日韩制服骚丝袜av| 亚洲av不卡在线观看| 汤姆久久久久久久影院中文字幕 | 国产精品一区二区在线观看99 | 亚洲人成网站在线观看播放| 国产一区亚洲一区在线观看| 伊人久久国产一区二区| eeuss影院久久| 国产69精品久久久久777片| 亚洲精品aⅴ在线观看| 久久久久九九精品影院| 三级经典国产精品| 久久99热这里只频精品6学生| 久久这里有精品视频免费| 爱豆传媒免费全集在线观看| 精品亚洲乱码少妇综合久久| 久久精品国产自在天天线| 成人鲁丝片一二三区免费| 午夜日本视频在线| av女优亚洲男人天堂| 亚洲电影在线观看av| 欧美激情在线99| 伦精品一区二区三区| 噜噜噜噜噜久久久久久91| 午夜免费观看性视频| 国产欧美另类精品又又久久亚洲欧美| 国产乱人偷精品视频| 最新中文字幕久久久久| 欧美 日韩 精品 国产| 街头女战士在线观看网站| 日韩 亚洲 欧美在线| 国产精品一区二区性色av| 美女cb高潮喷水在线观看| 亚洲成人av在线免费| 一个人免费在线观看电影| 国产伦理片在线播放av一区| 综合色av麻豆| 成人午夜精彩视频在线观看| 精品久久久久久久久av| 天堂网av新在线| av免费在线看不卡| 国产在视频线在精品| 久久这里只有精品中国| 久久久久久久国产电影| 美女主播在线视频| 91精品国产九色| 国产午夜精品久久久久久一区二区三区| 亚洲精品成人av观看孕妇| 午夜福利高清视频| 久久久久九九精品影院| 国产精品不卡视频一区二区| av播播在线观看一区| 国产综合精华液| 啦啦啦啦在线视频资源| 色网站视频免费| 国产精品人妻久久久影院| 成人二区视频| 国产免费福利视频在线观看| 天美传媒精品一区二区| 成年女人在线观看亚洲视频 | 亚洲国产欧美在线一区| 亚洲成人一二三区av| 一区二区三区高清视频在线| 日日摸夜夜添夜夜爱| 免费人成在线观看视频色| 亚洲成人精品中文字幕电影| 毛片女人毛片| 日韩在线高清观看一区二区三区| 国产 亚洲一区二区三区 | 亚洲国产色片| 看免费成人av毛片| 国产黄a三级三级三级人| 久久久久久久大尺度免费视频| 国产成人精品久久久久久| 日本av手机在线免费观看| 国产亚洲av片在线观看秒播厂 | 欧美性感艳星| videos熟女内射| 欧美日韩在线观看h| 色播亚洲综合网| 国产淫语在线视频| 少妇被粗大猛烈的视频| 老司机影院成人| 免费无遮挡裸体视频| 国产精品日韩av在线免费观看| 蜜桃亚洲精品一区二区三区| av.在线天堂| 国国产精品蜜臀av免费| 成人二区视频| 精品熟女少妇av免费看| 听说在线观看完整版免费高清| 日本三级黄在线观看| 日韩欧美 国产精品| 久久97久久精品| 精品久久久久久久人妻蜜臀av| 汤姆久久久久久久影院中文字幕 | 18禁裸乳无遮挡免费网站照片| 亚洲丝袜综合中文字幕| 亚洲欧美精品专区久久| 男女视频在线观看网站免费| 最近最新中文字幕大全电影3| 午夜福利网站1000一区二区三区| 97人妻精品一区二区三区麻豆| 午夜亚洲福利在线播放| 欧美一级a爱片免费观看看| 国内精品美女久久久久久| 有码 亚洲区| 18禁动态无遮挡网站| 2022亚洲国产成人精品| 性色avwww在线观看| 麻豆成人av视频| 一级片'在线观看视频| 成年av动漫网址| 蜜臀久久99精品久久宅男| 一区二区三区乱码不卡18| 欧美潮喷喷水| 色综合亚洲欧美另类图片| 麻豆精品久久久久久蜜桃| 亚洲成人av在线免费| 国国产精品蜜臀av免费| 国产91av在线免费观看| 熟女电影av网| 亚洲aⅴ乱码一区二区在线播放| 国产爱豆传媒在线观看| 久久精品综合一区二区三区| 免费看av在线观看网站| 亚洲成人av在线免费| 最近中文字幕高清免费大全6| 国产欧美日韩精品一区二区| 我的女老师完整版在线观看| 极品教师在线视频| 国产黄频视频在线观看| 亚洲成人久久爱视频| 大陆偷拍与自拍| 亚洲av一区综合| 国产v大片淫在线免费观看| 国产一级毛片七仙女欲春2| 男人舔奶头视频| 欧美日韩视频高清一区二区三区二| 欧美 日韩 精品 国产| 精品国产露脸久久av麻豆 | 97热精品久久久久久| 日韩视频在线欧美| av在线观看视频网站免费| 欧美xxxx黑人xx丫x性爽| 午夜免费男女啪啪视频观看| or卡值多少钱| 看免费成人av毛片| 插阴视频在线观看视频| 成人毛片a级毛片在线播放| 亚洲精品第二区| 免费观看a级毛片全部| 亚洲国产精品国产精品| 高清av免费在线| 中文资源天堂在线| 卡戴珊不雅视频在线播放| 亚洲av二区三区四区| 一区二区三区乱码不卡18| 精品一区二区三区人妻视频| 少妇熟女欧美另类| 亚洲成人中文字幕在线播放| 尤物成人国产欧美一区二区三区| 午夜精品在线福利| ponron亚洲| 欧美3d第一页| 久久久久久伊人网av| 精品国内亚洲2022精品成人| 午夜亚洲福利在线播放| 在线天堂最新版资源| 国产黄色视频一区二区在线观看| 亚洲精品日韩av片在线观看| 乱码一卡2卡4卡精品| 国国产精品蜜臀av免费| www.av在线官网国产| 插阴视频在线观看视频| 亚洲av国产av综合av卡| 欧美高清性xxxxhd video| 国产男人的电影天堂91| 欧美激情国产日韩精品一区| 插阴视频在线观看视频| 亚洲av中文av极速乱| 国产精品久久久久久久久免| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 99久久精品一区二区三区| 午夜亚洲福利在线播放| 亚洲电影在线观看av| 亚洲精品成人av观看孕妇| 久久久久精品久久久久真实原创| 大又大粗又爽又黄少妇毛片口| 免费看美女性在线毛片视频| 熟女电影av网| 国产乱来视频区| 最新中文字幕久久久久| 大又大粗又爽又黄少妇毛片口| 一级毛片久久久久久久久女| 国产精品久久久久久久电影| 国产精品国产三级国产专区5o| 搡女人真爽免费视频火全软件| 国产探花在线观看一区二区| 免费观看性生交大片5| 欧美三级亚洲精品| 成人鲁丝片一二三区免费| 亚洲av.av天堂| 亚洲av中文字字幕乱码综合| 在线观看免费高清a一片| 亚洲欧美精品自产自拍| 成人二区视频| 五月玫瑰六月丁香| 久久精品国产亚洲av涩爱| 能在线免费看毛片的网站| 午夜日本视频在线| 直男gayav资源| 夫妻午夜视频| 精品酒店卫生间| 亚洲国产色片| 少妇丰满av| 成年女人在线观看亚洲视频 | 特大巨黑吊av在线直播| 欧美精品国产亚洲| 免费看不卡的av| 性插视频无遮挡在线免费观看| 肉色欧美久久久久久久蜜桃 | 久久久久免费精品人妻一区二区| videossex国产| freevideosex欧美| 直男gayav资源| 五月玫瑰六月丁香| 天天一区二区日本电影三级| 亚洲精品日本国产第一区| 亚洲成人中文字幕在线播放| 亚洲三级黄色毛片| 久久97久久精品| 欧美日韩一区二区视频在线观看视频在线 | 精品久久久久久久久亚洲| 日韩视频在线欧美| 亚洲欧美日韩卡通动漫| 男的添女的下面高潮视频| 欧美成人精品欧美一级黄| 亚洲av成人av| 国产一区二区三区综合在线观看 | 色吧在线观看| 伦精品一区二区三区| 国产亚洲av嫩草精品影院| 亚洲欧美日韩无卡精品| 日本一本二区三区精品| 国产黄a三级三级三级人| 欧美高清成人免费视频www| 淫秽高清视频在线观看| 国产又色又爽无遮挡免| 亚洲最大成人av| 久久亚洲国产成人精品v| 日本免费a在线| 亚洲最大成人中文| 日韩av在线大香蕉| 精品不卡国产一区二区三区| 在现免费观看毛片| 亚洲经典国产精华液单| 午夜老司机福利剧场| 永久网站在线| 51国产日韩欧美| 国产精品福利在线免费观看| 中文字幕亚洲精品专区| 精品一区二区免费观看| 又黄又爽又刺激的免费视频.| 黑人高潮一二区| 午夜视频国产福利| 成人漫画全彩无遮挡| 成人高潮视频无遮挡免费网站| 嫩草影院新地址| 欧美日韩综合久久久久久| 老司机影院成人| 午夜福利在线在线| 欧美 日韩 精品 国产| 国产成人免费观看mmmm| 久久精品熟女亚洲av麻豆精品 | 人妻系列 视频| 2022亚洲国产成人精品| 成人亚洲精品av一区二区| 久久久久久久午夜电影| 精品国产一区二区三区久久久樱花 | 国产午夜精品久久久久久一区二区三区| 国产成人a区在线观看| 免费观看无遮挡的男女| 亚洲精品日韩在线中文字幕| 真实男女啪啪啪动态图| 日本爱情动作片www.在线观看| 国产爱豆传媒在线观看| 精品久久久久久久久久久久久| 在线观看美女被高潮喷水网站| 纵有疾风起免费观看全集完整版 | 亚洲欧美日韩无卡精品| 联通29元200g的流量卡| 欧美xxxx黑人xx丫x性爽| 中文欧美无线码| 可以在线观看毛片的网站| 97热精品久久久久久| 国产 一区精品| 成人一区二区视频在线观看| 好男人视频免费观看在线| 亚洲成人精品中文字幕电影| 一级毛片aaaaaa免费看小| 欧美xxxx黑人xx丫x性爽| 精品人妻视频免费看| 亚洲av中文字字幕乱码综合| 亚洲欧美清纯卡通| 日韩伦理黄色片| 日韩成人伦理影院|