• <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)工藝的對比
    不卡一级毛片| 久9热在线精品视频| 我的老师免费观看完整版| 免费黄网站久久成人精品| 亚洲人与动物交配视频| 国模一区二区三区四区视频| 亚洲无线在线观看| 免费电影在线观看免费观看| 国产免费av片在线观看野外av| 999久久久精品免费观看国产| 99久久九九国产精品国产免费| 国产亚洲av嫩草精品影院| 国内揄拍国产精品人妻在线| 搡女人真爽免费视频火全软件 | 免费搜索国产男女视频| 69人妻影院| 毛片一级片免费看久久久久 | 午夜福利视频1000在线观看| 国内揄拍国产精品人妻在线| 亚洲乱码一区二区免费版| 精品福利观看| 久久久久久久久久成人| 日韩大尺度精品在线看网址| 99热这里只有是精品在线观看| 俺也久久电影网| 长腿黑丝高跟| 亚洲综合色惰| 色噜噜av男人的天堂激情| 国产精品久久久久久亚洲av鲁大| 午夜福利在线观看免费完整高清在 | 久久欧美精品欧美久久欧美| 小说图片视频综合网站| 俄罗斯特黄特色一大片| 成年版毛片免费区| 国产精品久久电影中文字幕| 欧美成人一区二区免费高清观看| 亚洲精品粉嫩美女一区| 99热网站在线观看| 日本五十路高清| 日韩精品青青久久久久久| 国产熟女欧美一区二区| 国产aⅴ精品一区二区三区波| 精品无人区乱码1区二区| 99riav亚洲国产免费| 狂野欧美激情性xxxx在线观看| 18禁黄网站禁片免费观看直播| 日韩精品中文字幕看吧| 国产在线精品亚洲第一网站| 亚洲人与动物交配视频| 亚洲va在线va天堂va国产| 91久久精品国产一区二区成人| 欧美一区二区国产精品久久精品| 日韩精品青青久久久久久| 日本熟妇午夜| 国产成人一区二区在线| 午夜福利视频1000在线观看| 99久久精品热视频| 午夜免费男女啪啪视频观看 | 成人三级黄色视频| 99久久久亚洲精品蜜臀av| 国产精品98久久久久久宅男小说| 看黄色毛片网站| 精品久久久久久久久久免费视频| 真实男女啪啪啪动态图| 麻豆av噜噜一区二区三区| 国产精品美女特级片免费视频播放器| 国产亚洲av嫩草精品影院| 美女cb高潮喷水在线观看| 国模一区二区三区四区视频| 夜夜爽天天搞| 偷拍熟女少妇极品色| 99在线人妻在线中文字幕| 欧美日韩精品成人综合77777| 变态另类成人亚洲欧美熟女| 国产精品乱码一区二三区的特点| 极品教师在线视频| 国产成人aa在线观看| 亚洲av免费高清在线观看| 性插视频无遮挡在线免费观看| 免费观看人在逋| 亚洲国产欧美人成| eeuss影院久久| 欧美一区二区精品小视频在线| 少妇高潮的动态图| 免费观看精品视频网站| 亚洲av不卡在线观看| 成年免费大片在线观看| 亚洲欧美激情综合另类| 久久午夜福利片| 国产伦人伦偷精品视频| 色综合站精品国产| 国产男人的电影天堂91| 九九久久精品国产亚洲av麻豆| 国产亚洲精品久久久久久毛片| 亚洲成av人片在线播放无| 国产aⅴ精品一区二区三区波| 人妻少妇偷人精品九色| 狠狠狠狠99中文字幕| 日韩欧美一区二区三区在线观看| 日本黄色视频三级网站网址| 一个人看视频在线观看www免费| 国产综合懂色| 国产淫片久久久久久久久| 国产精品三级大全| 床上黄色一级片| 蜜桃久久精品国产亚洲av| 欧美色视频一区免费| 亚洲va在线va天堂va国产| 成人无遮挡网站| 毛片一级片免费看久久久久 | 亚洲av美国av| 自拍偷自拍亚洲精品老妇| 亚洲中文字幕一区二区三区有码在线看| 亚洲无线观看免费| 真实男女啪啪啪动态图| 国产精品久久电影中文字幕| 亚洲三级黄色毛片| 亚洲人成网站高清观看| 久久久久久久久久久丰满 | 成年女人毛片免费观看观看9| 午夜爱爱视频在线播放| 欧美色视频一区免费| 99久国产av精品| 国产淫片久久久久久久久| 久久九九热精品免费| 在线免费观看的www视频| 成人毛片a级毛片在线播放| 日本成人三级电影网站| 国产欧美日韩精品一区二区| 免费一级毛片在线播放高清视频| 国产成人a区在线观看| 国产在线精品亚洲第一网站| 窝窝影院91人妻| 日日夜夜操网爽| 十八禁国产超污无遮挡网站| av在线老鸭窝| 欧美极品一区二区三区四区| 欧美高清性xxxxhd video| 国内久久婷婷六月综合欲色啪| 高清在线国产一区| 老师上课跳d突然被开到最大视频| 国产在线精品亚洲第一网站| 淫妇啪啪啪对白视频| 午夜精品一区二区三区免费看| 国产综合懂色| 亚洲三级黄色毛片| 国产亚洲精品综合一区在线观看| 中文资源天堂在线| 成人永久免费在线观看视频| 内射极品少妇av片p| 国产毛片a区久久久久| 免费无遮挡裸体视频| 亚洲精品一卡2卡三卡4卡5卡| 在线天堂最新版资源| 亚洲 国产 在线| 国产真实伦视频高清在线观看 | 听说在线观看完整版免费高清| 在线免费观看的www视频| 天堂√8在线中文| 久久久精品欧美日韩精品| 一区二区三区四区激情视频 | 啦啦啦啦在线视频资源| 夜夜夜夜夜久久久久| 啦啦啦啦在线视频资源| 久99久视频精品免费| 亚洲欧美清纯卡通| 亚洲五月天丁香| 日本 欧美在线| 日本欧美国产在线视频| a在线观看视频网站| 国产亚洲精品久久久com| 精品欧美国产一区二区三| 乱系列少妇在线播放| 一级黄片播放器| 在线看三级毛片| 成人二区视频| 老师上课跳d突然被开到最大视频| 欧美绝顶高潮抽搐喷水| 22中文网久久字幕| 无遮挡黄片免费观看| 色噜噜av男人的天堂激情| 国国产精品蜜臀av免费| 超碰av人人做人人爽久久| 亚洲狠狠婷婷综合久久图片| 俄罗斯特黄特色一大片| 国产高清有码在线观看视频| 少妇猛男粗大的猛烈进出视频 | 综合色av麻豆| 欧美性感艳星| 俄罗斯特黄特色一大片| 在线观看美女被高潮喷水网站| 我要看日韩黄色一级片| 欧美黑人巨大hd| 一本久久中文字幕| 最近在线观看免费完整版| 在线观看一区二区三区| 午夜免费成人在线视频| 日韩中字成人| 香蕉av资源在线| 久久精品综合一区二区三区| 观看美女的网站| 精品久久久久久久末码| 在现免费观看毛片| 麻豆av噜噜一区二区三区| 亚洲一区高清亚洲精品| 精品人妻一区二区三区麻豆 | 超碰av人人做人人爽久久| 高清毛片免费观看视频网站| 中文字幕精品亚洲无线码一区| 久久午夜亚洲精品久久| 热99re8久久精品国产| 很黄的视频免费| 校园春色视频在线观看| 国产在线男女| 在线a可以看的网站| 国内少妇人妻偷人精品xxx网站| 久久精品国产清高在天天线| 91在线精品国自产拍蜜月| 午夜爱爱视频在线播放| 夜夜夜夜夜久久久久| 极品教师在线视频| 老女人水多毛片| 波多野结衣巨乳人妻| 国产精品亚洲美女久久久| 国产乱人伦免费视频| 午夜福利高清视频| 亚洲精品日韩av片在线观看| eeuss影院久久| 亚洲精品色激情综合| 成年女人永久免费观看视频| 久久久精品欧美日韩精品| 亚洲国产精品成人综合色| 一本一本综合久久| 国产亚洲欧美98| 精品人妻一区二区三区麻豆 | 久久6这里有精品| 一边摸一边抽搐一进一小说| 欧美xxxx性猛交bbbb| 99久国产av精品| 成人综合一区亚洲| 日日摸夜夜添夜夜添av毛片 | 天堂动漫精品| a级毛片a级免费在线| 岛国在线免费视频观看| 男人和女人高潮做爰伦理| 成人鲁丝片一二三区免费| 国产大屁股一区二区在线视频| 在线播放无遮挡| 在线观看66精品国产| 日韩欧美国产一区二区入口| 18禁黄网站禁片免费观看直播| 欧美黑人欧美精品刺激| 午夜亚洲福利在线播放| 国产成人影院久久av| 国产日本99.免费观看| 日本在线视频免费播放| 搞女人的毛片| 老熟妇仑乱视频hdxx| 男女做爰动态图高潮gif福利片| 国产男人的电影天堂91| 亚洲一级一片aⅴ在线观看| 我要搜黄色片| 99久久无色码亚洲精品果冻| 国内精品久久久久久久电影| 又爽又黄a免费视频| 成人精品一区二区免费| 91狼人影院| 欧美日韩瑟瑟在线播放| 97碰自拍视频| 美女黄网站色视频| 99九九线精品视频在线观看视频| 最后的刺客免费高清国语| 精品午夜福利视频在线观看一区| 在线国产一区二区在线| 99视频精品全部免费 在线| 国产爱豆传媒在线观看| 日韩精品有码人妻一区| 18禁在线播放成人免费| 国模一区二区三区四区视频| 欧美丝袜亚洲另类 | 成人毛片a级毛片在线播放| 精品一区二区三区视频在线观看免费| 国产一区二区在线观看日韩| 偷拍熟女少妇极品色| 欧美性感艳星| 久久久色成人| 国产v大片淫在线免费观看| 国产探花在线观看一区二区| 人人妻人人看人人澡| 日韩,欧美,国产一区二区三区 | 尤物成人国产欧美一区二区三区| 日日夜夜操网爽| 在线观看免费视频日本深夜| 日韩中字成人| 中文字幕熟女人妻在线| 久久国内精品自在自线图片| 高清在线国产一区| 成人一区二区视频在线观看| 两性午夜刺激爽爽歪歪视频在线观看| 动漫黄色视频在线观看| 国语自产精品视频在线第100页| 老熟妇仑乱视频hdxx| 免费不卡的大黄色大毛片视频在线观看 | 搡女人真爽免费视频火全软件 | 99热这里只有是精品50| 中文亚洲av片在线观看爽| 热99在线观看视频| 欧美成人a在线观看| 十八禁国产超污无遮挡网站| 日韩中文字幕欧美一区二区| 国产 一区精品| 乱人视频在线观看| 男女视频在线观看网站免费| 给我免费播放毛片高清在线观看| 国产中年淑女户外野战色| 99热这里只有精品一区| 国产精品99久久久久久久久| 国产高潮美女av| www.色视频.com| 一级毛片久久久久久久久女| 亚洲成人久久爱视频| 99热这里只有是精品50| 久久香蕉精品热| 国产高清激情床上av| 露出奶头的视频| 麻豆国产av国片精品| 国产aⅴ精品一区二区三区波| 日韩在线高清观看一区二区三区 | 在线观看av片永久免费下载| 哪里可以看免费的av片| 18禁裸乳无遮挡免费网站照片| 免费黄网站久久成人精品| 伊人久久精品亚洲午夜| 亚洲真实伦在线观看| 亚洲自偷自拍三级| 国产亚洲精品久久久com| 精品久久久噜噜| 一个人看视频在线观看www免费| 日韩欧美 国产精品| 他把我摸到了高潮在线观看| 久久久精品欧美日韩精品| 午夜免费激情av| 国产精品久久久久久久久免| 亚洲av成人精品一区久久| 三级男女做爰猛烈吃奶摸视频| 自拍偷自拍亚洲精品老妇| 麻豆国产97在线/欧美| 俺也久久电影网| 亚洲久久久久久中文字幕| 日韩在线高清观看一区二区三区 | 午夜福利18| 国产亚洲91精品色在线| 午夜亚洲福利在线播放| 精品无人区乱码1区二区| 国内精品久久久久精免费| 特级一级黄色大片| 神马国产精品三级电影在线观看| 五月玫瑰六月丁香| 欧美黑人巨大hd| 久久久久久久久久久丰满 | 亚洲第一区二区三区不卡| www日本黄色视频网| 亚洲av日韩精品久久久久久密| 99riav亚洲国产免费| 日本三级黄在线观看| 欧美高清成人免费视频www| 99视频精品全部免费 在线| 在线观看一区二区三区| АⅤ资源中文在线天堂| 亚洲国产色片| 午夜精品一区二区三区免费看| 久久久久久九九精品二区国产| 九色成人免费人妻av| 国模一区二区三区四区视频| 老司机午夜福利在线观看视频| 国产在线精品亚洲第一网站| 我的老师免费观看完整版| 99久久久亚洲精品蜜臀av| 久久人人精品亚洲av| 欧美一区二区国产精品久久精品| 91久久精品国产一区二区成人| 天天躁日日操中文字幕| 免费观看的影片在线观看| 成人二区视频| 欧美激情在线99| 国内精品宾馆在线| 久久草成人影院| 观看免费一级毛片| 精品久久久久久久久久免费视频| 国产精品野战在线观看| 婷婷亚洲欧美| 国产男靠女视频免费网站| 国产精品国产三级国产av玫瑰| 狂野欧美激情性xxxx在线观看| 69av精品久久久久久| 亚洲五月天丁香| 熟女人妻精品中文字幕| 久久精品国产鲁丝片午夜精品 | av专区在线播放| 国产精品亚洲一级av第二区| 欧美bdsm另类| 一级毛片久久久久久久久女| 桃红色精品国产亚洲av| 精品午夜福利视频在线观看一区| 欧美高清成人免费视频www| 国产免费男女视频| 国内精品美女久久久久久| 两个人的视频大全免费| 亚洲综合色惰| 精品久久久久久久末码| 一个人免费在线观看电影| 国内精品宾馆在线| 中文字幕免费在线视频6| 婷婷精品国产亚洲av| 五月玫瑰六月丁香| 99久久无色码亚洲精品果冻| 男人狂女人下面高潮的视频| 日韩在线高清观看一区二区三区 | 成人二区视频| 国产精品人妻久久久久久| 1000部很黄的大片| 99热网站在线观看| 免费搜索国产男女视频| 亚洲欧美精品综合久久99| 国产精品98久久久久久宅男小说| 精品人妻视频免费看| 97人妻精品一区二区三区麻豆| 亚洲国产欧美人成| 嫁个100分男人电影在线观看| 国产v大片淫在线免费观看| av黄色大香蕉| 少妇高潮的动态图| 色噜噜av男人的天堂激情| 嫩草影院新地址| 欧美色视频一区免费| 又黄又爽又刺激的免费视频.| 在线a可以看的网站| 久久国产乱子免费精品| 最近中文字幕高清免费大全6 | or卡值多少钱| 色精品久久人妻99蜜桃| 一级av片app| 欧美一级a爱片免费观看看| 国产高潮美女av| 九色成人免费人妻av| videossex国产| 观看免费一级毛片| 国产精品98久久久久久宅男小说| 国内揄拍国产精品人妻在线| 少妇猛男粗大的猛烈进出视频 | 久久人妻av系列| 国产高潮美女av| 亚洲精品国产成人久久av| 欧美xxxx性猛交bbbb| 国产激情偷乱视频一区二区| 在线观看一区二区三区| 黄色视频,在线免费观看| www.www免费av| 蜜桃久久精品国产亚洲av| 亚洲性久久影院| 午夜福利18| 在线免费观看的www视频| 成人亚洲精品av一区二区| 三级男女做爰猛烈吃奶摸视频| 久久人人爽人人爽人人片va| 国产日本99.免费观看| 伦理电影大哥的女人| 亚洲狠狠婷婷综合久久图片| 免费黄网站久久成人精品| 国产午夜福利久久久久久| 国产一区二区亚洲精品在线观看| 美女高潮的动态| 亚洲欧美日韩东京热| 久久精品国产清高在天天线| 亚洲人成伊人成综合网2020| 亚洲国产精品合色在线| 亚洲乱码一区二区免费版| 久久久久久久久久黄片| 91麻豆av在线| 午夜激情福利司机影院| 一本精品99久久精品77| 日本黄大片高清| 天堂动漫精品| 亚洲av免费高清在线观看| 成人美女网站在线观看视频| 日本黄色视频三级网站网址| 91麻豆精品激情在线观看国产| 久久天躁狠狠躁夜夜2o2o| 亚洲精品乱码久久久v下载方式| 白带黄色成豆腐渣| 2021天堂中文幕一二区在线观| 国产高清视频在线观看网站| 给我免费播放毛片高清在线观看| 直男gayav资源| 国产乱人伦免费视频| 他把我摸到了高潮在线观看| 99国产精品一区二区蜜桃av| 欧美另类亚洲清纯唯美| 国产精品一区二区三区四区免费观看 | 99热精品在线国产| 久久精品91蜜桃| 97超视频在线观看视频| 亚洲av免费高清在线观看| 啦啦啦韩国在线观看视频| 国产高清视频在线播放一区| a级毛片免费高清观看在线播放| 午夜亚洲福利在线播放| 欧美高清性xxxxhd video| 伦理电影大哥的女人| 久久人人精品亚洲av| 一区二区三区四区激情视频 | 久久精品国产亚洲av天美| 在线播放无遮挡| av天堂在线播放| 久久久久久久精品吃奶| 尤物成人国产欧美一区二区三区| 不卡一级毛片| 国国产精品蜜臀av免费| 精品午夜福利视频在线观看一区| 色哟哟哟哟哟哟| 国产精品三级大全| 99热6这里只有精品| 国产精品野战在线观看| 国产极品精品免费视频能看的| 亚洲国产欧洲综合997久久,| 在线免费十八禁| 国产精品电影一区二区三区| 亚洲精品日韩av片在线观看| 日本 欧美在线| 日韩欧美国产在线观看| 国产免费男女视频| 精品国产三级普通话版| 国产亚洲av嫩草精品影院| 真实男女啪啪啪动态图| 国产精品综合久久久久久久免费| 国产精品三级大全| aaaaa片日本免费| 亚洲av中文字字幕乱码综合| 亚洲国产精品sss在线观看| 亚洲五月天丁香| 国产伦精品一区二区三区四那| 国产一区二区三区av在线 | 亚洲欧美清纯卡通| 99热精品在线国产| 亚洲自拍偷在线| 成人无遮挡网站| 精品久久久久久久人妻蜜臀av| 色在线成人网| 一夜夜www| 99久国产av精品| 日本黄色视频三级网站网址| 国产精品爽爽va在线观看网站| 精品午夜福利视频在线观看一区| 两个人的视频大全免费| 夜夜爽天天搞| h日本视频在线播放| 精品一区二区免费观看| 美女cb高潮喷水在线观看| 亚洲av成人精品一区久久| 国产女主播在线喷水免费视频网站 | 无人区码免费观看不卡| 午夜日韩欧美国产| 久久精品国产亚洲av香蕉五月| 好男人在线观看高清免费视频| 国产亚洲精品av在线| 亚洲精品日韩av片在线观看| 色播亚洲综合网| 国产在线精品亚洲第一网站| 日本熟妇午夜| 成年女人永久免费观看视频| 欧美另类亚洲清纯唯美| 精品一区二区三区人妻视频| a级毛片a级免费在线| 蜜桃亚洲精品一区二区三区| 亚洲第一区二区三区不卡| www.色视频.com| 一本精品99久久精品77| 久久国产乱子免费精品| 欧美成人一区二区免费高清观看| 又黄又爽又刺激的免费视频.| 国产精品无大码| av女优亚洲男人天堂| 午夜福利在线观看免费完整高清在 | 国内久久婷婷六月综合欲色啪| 久久热精品热| 中国美女看黄片| 亚洲午夜理论影院| 日韩精品有码人妻一区| 麻豆国产av国片精品| 99热精品在线国产| 天堂动漫精品| 老司机福利观看| 人妻夜夜爽99麻豆av| 欧美日韩乱码在线| 色综合婷婷激情| 久久久精品欧美日韩精品| 又黄又爽又刺激的免费视频.| 少妇人妻一区二区三区视频| 三级毛片av免费| 欧美日韩国产亚洲二区| 制服丝袜大香蕉在线| 国产精品电影一区二区三区| eeuss影院久久| 麻豆成人午夜福利视频| av.在线天堂| 国产淫片久久久久久久久| 少妇高潮的动态图| 色哟哟·www| 99热这里只有是精品在线观看| 波多野结衣高清无吗| 日本一二三区视频观看| 欧美激情久久久久久爽电影| 中文字幕熟女人妻在线| 亚洲精品乱码久久久v下载方式| 婷婷色综合大香蕉|