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

    磁場對激波沖擊R22重氣柱作用過程影響的數(shù)值模擬*

    2017-07-31 21:03:15林震亞張煥好陳志華
    爆炸與沖擊 2017年4期
    關(guān)鍵詞:氣柱不穩(wěn)定性激波

    林震亞,張煥好,陳志華,劉 迎

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

    磁場對激波沖擊R22重氣柱作用過程影響的數(shù)值模擬*

    林震亞,張煥好,陳志華,劉 迎

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

    為研究平面入射激波與磁化R22重質(zhì)圓形氣柱的作用過程,首先通過數(shù)值方法得到了不同初始條件下激波誘導(dǎo)R22氣柱的Kelvin-Helmholtz(KH)及Richtmyer-Meshkov(RM)不穩(wěn)定性導(dǎo)致的重氣柱變形過程,并詳細(xì)討論了不同情況下透射激波在氣柱內(nèi)聚焦誘導(dǎo)射流的過程;然后在加入磁場的情況下,采用CTU+CT算法進(jìn)行數(shù)值模擬,以保證數(shù)值結(jié)果滿足任意時(shí)刻磁場的散度為零。計(jì)算結(jié)果表明:磁場對激波誘導(dǎo)R22氣柱不穩(wěn)定性具有抑制作用;法向磁場和流向磁場都可以很好地抑制RM不穩(wěn)定性;對于KH不穩(wěn)定性,法向磁場的控制效果更好,不僅可以抑制界面上渦串的卷起,還可以阻止主渦的發(fā)展,而流向磁場做不到后者;磁場對射流影響不大,射流處的磁能量可以一定程度上抑制射流的衰減,同時(shí)法向磁場可以減小聚焦時(shí)壓力及速度峰值。

    磁化氣體;磁流體動力學(xué);磁場;不穩(wěn)定性;激波

    激波與兩種不同密度氣體界面相互作用時(shí),會在界面附近發(fā)生反射、折射及繞射等現(xiàn)象,同時(shí)在界面處還會產(chǎn)生Richtmyer-Meshkov(RM)和Kelvin-Helmholtz(KH)不穩(wěn)定性,導(dǎo)致界面扭曲變形。該現(xiàn)象廣泛存在于超燃沖壓發(fā)動機(jī)中的燃料混合、慣性約束核聚變、水中炸藥爆炸及天體物理中的超新星爆炸等領(lǐng)域,因而具有重要的研究價(jià)值。

    目前,對激波與氣體界面相互作用的研究主要集中以下兩方面:一是氣體界面在RM不穩(wěn)定性與KH不穩(wěn)定性作用下的演變過程及其相關(guān)復(fù)雜流動機(jī)理;二是界面演化后期的失穩(wěn)、破碎及湍流混合問題。K.A.Meyer等[1]對激波與氣體界面作用時(shí)的RM不穩(wěn)定性進(jìn)行了數(shù)值模擬,得到了與Meshkov實(shí)驗(yàn)定性一致的結(jié)果;Q.Zhang等[2]研究了二維和三維初始擾動的發(fā)展過程,得到了氣泡、尖釘速度及界面發(fā)展的總體增長率的計(jì)算公式;N.N.Anuchina等[3]對初始擾動幅度下相同入射強(qiáng)度激波與正弦擾動氣體界面的作用過程進(jìn)行了數(shù)值模擬,結(jié)果表明三維數(shù)值計(jì)算得到的擾動增長率高于二維結(jié)果;G.A.Ruev等[4]對馬赫數(shù)為3.5的平面入射激波與正弦擾動氣體界面的作用過程進(jìn)行了數(shù)值模擬,給出了氣體界面RM不穩(wěn)定性的發(fā)展過程及兩種氣體混合區(qū)的寬度變化;J.H.J.Niederhaus等[5]對不同強(qiáng)度平面入射激波與多種不同介質(zhì)的氣柱、氣泡相互作用過程進(jìn)行了大量的數(shù)值研究,重點(diǎn)分析了作用過程中流場環(huán)量的變化和兩種氣體的混合機(jī)理;同時(shí),B.Thornber等[6]對低馬赫數(shù)的平面入射激波與R22方柱作用過程的數(shù)值模擬,結(jié)果與Holder等的實(shí)驗(yàn)結(jié)果完全相吻;B.Hejazialhosseini等[7]對不同馬赫數(shù)下平面入射激波與氦氣氣泡作用過程進(jìn)行了數(shù)值研究,得到了氣泡的變形過程及氣泡表面在KH不穩(wěn)定作用下次級渦的產(chǎn)生過程;O.Schilling等[8]采用數(shù)值方法研究了平面入射激波與Air/R22氣體界面的相互作用過程,并分析了反射激波與分界面作用前后流場結(jié)構(gòu)和混合層演變;B.Tian等[9]利用ALE(arbitrary Lagrangian-Eulerian)和GALE(global arbitrary Lagrangian-Eulerian)方法對多種激波與氣體界面作用問題進(jìn)行了數(shù)值模擬;S.K.Shankar等[10]采用數(shù)值方法模擬了平面入射激波與重氣柱作用過程,結(jié)果表明流場初始壓力和密度梯度是影響兩個(gè)主渦產(chǎn)生的關(guān)鍵因素;C.Bailie等[11]對不同強(qiáng)度平面入射激波與Air/R22正弦小擾動界面的作用過程進(jìn)行了數(shù)值模擬,并分析入射激波強(qiáng)度對作用過程的影響。

    目前,對界面RM不穩(wěn)定性控制方法的研究表明,合適的磁場對氣體界面RM不穩(wěn)定性的產(chǎn)生與發(fā)展具有一定的抑制作用,但仍處在理論研究階段。S.Chandrasekhar等[12]和V.Wheatley等[13-14]采用線性理論方法推導(dǎo)了垂直或平行于界面的磁場對不可壓縮流體界面不穩(wěn)定性的影響,結(jié)果表明:平行于界面的磁場能夠抑制RT和RM不穩(wěn)定性;但是在線性理論下,平面流場中的RM不穩(wěn)定性與垂直于平面的磁場無關(guān),與實(shí)際情況不符。為進(jìn)一步研究磁場效應(yīng),M.Khan等[15]在勢流理論的基礎(chǔ)上建立了新模型,發(fā)現(xiàn)考慮非線性效應(yīng)后,磁場能夠影響界面不穩(wěn)定性的發(fā)展。J.Cao等[16]同時(shí)考慮了初始平行于界面的磁場和剪切流對界面不穩(wěn)定性的影響,發(fā)現(xiàn)磁場能夠抑制界面不穩(wěn)定性,而剪切流將會促進(jìn)界面不穩(wěn)定性的發(fā)展,兩者之間存在競爭機(jī)制。M.S.Shin等[17]對不同方向磁場作用下,平面激波與球形氣泡相互作用過程進(jìn)行了數(shù)值模擬,發(fā)現(xiàn)在發(fā)展初期,磁場效應(yīng)很小,而在發(fā)展后期,即使是弱磁場,其對界面形狀的影響也很明顯。此外,磁場越強(qiáng),界面內(nèi)外物質(zhì)混合率越低,并越能抑制后期湍流的發(fā)展。李源等[18]基于理論分析的方法,綜合考慮黏性和表面張力的影響,考察了磁場中非理想流體中RT不穩(wěn)定性氣泡的演化過程,并推導(dǎo)出描述二維非理想磁流體RT不穩(wěn)定性氣泡運(yùn)動的控制方程組,得到了流體黏性、表面張力和磁場對氣泡界面不穩(wěn)定性發(fā)展的影響。

    綜上所述,目前對磁場控制不穩(wěn)定性的研究主要集中在基于不可壓理想磁流場的解析推導(dǎo)及磁場影響下所產(chǎn)生現(xiàn)象的數(shù)值模擬兩個(gè)方面。本文中,基于磁流體動力學(xué)(magneto-h(huán)ydro-dynamics,MHD)方程,采用CTU+CT算法對激波沖擊磁化R22氣柱進(jìn)行數(shù)值模擬,并以磁壓力為切入點(diǎn)研究磁場對不穩(wěn)定性及射流的影響。

    1 數(shù)值方法及計(jì)算模型

    1.1 數(shù)值方法

    在滿足流體力學(xué)方程的條件下,帶電荷的流體團(tuán)更有利于形成團(tuán)塊,主要原因如下:(1)帶電粒子間的碰撞頻率遠(yuǎn)大于中性粒子之間的碰撞頻率,因此其碰撞平均自由程較短;(2)有磁場存在時(shí),磁場對帶電粒子有側(cè)向約束。洛倫茲力使帶電粒子在垂直于磁場方向只能保持在回轉(zhuǎn)半徑范圍內(nèi)運(yùn)動。等離子體通常比中性粒子有更低的密度和更高的溫度,由于整個(gè)過程在短時(shí)間內(nèi)完成,本文中考慮高磁雷諾數(shù)的情況,使用MHD方程來研究激波沖擊R22氣柱中出現(xiàn)的不穩(wěn)定性:

    式中:ρ為質(zhì)量密度,v為速度矢量,p為壓力,E為能量,μ為磁導(dǎo)率,B為磁感應(yīng)強(qiáng)度,^b為B的單位矢量,Π為包含各向同性與各向異性部分的黏性應(yīng)力張量,分別由系數(shù)ν0和ν‖控制,γ為比熱比,T為溫度,I為單位張量,符號“∶”表示雙點(diǎn)積運(yùn)算。在各向異性的情況下,黏性通量限制在平行于磁場線方向。Q為熱通量,同樣包含各向同性與各向異性的情況,分別由系數(shù)κ0和κ‖控制。在各向異性的情況下,熱通量同樣限制在平行于磁場線方向。同時(shí),感應(yīng)方程包含多種非理想MHD效應(yīng),包括歐姆損耗(由電阻率η控制),霍爾效應(yīng)(由系數(shù)ηH控制)和雙極性擴(kuò)散(由系數(shù)ηAD控制)。

    本文中采用CTU+CT算法對MHD方程組進(jìn)行求解,其中CTU算法是一種用于求解守恒雙曲系統(tǒng)的不分離的二維有限體積算法,其三維變形則由J.Saltzman[19]提出。另外,由于傳統(tǒng)的12-solve CTU算法忽略了MHD方程中始終存在的約束·B=0,不能作為求解MHD方程的有效算法,因此需要構(gòu)建一個(gè)Godunov通量中網(wǎng)格單元邊緣的平均電場算法,通常稱為CT算法,該算法可以描述為以Godunov電場為預(yù)測值,以CT電場為校正值的預(yù)測過程。鑒于在文獻(xiàn)[20]中εcCT已顯示的最佳性能,本文中采用εcCT算法。該算法包含逆風(fēng)偏差(根據(jù)接觸模式),可以將網(wǎng)格對齊的平行平面流動降低到正確的Godunov電動勢。

    1.2 計(jì)算模型

    圖1為激波與R22氣柱相互作用的計(jì)算模型,其中R22氣泡直徑為0.05m,氣柱中心與左邊界的距離為0.05m,與右邊界的距離為0.25m,計(jì)算域高度為0.089m。初始時(shí),氣柱內(nèi)采用R22作為重質(zhì)氣體,外部介質(zhì)為空氣,氣體參數(shù)如表1所示。氣柱內(nèi)外壓力均為101 325Pa,入射激波從左向右傳播,來流馬赫數(shù)Ma=1.22,上、下邊界設(shè)為固體反射邊界,而右邊界設(shè)為出口條件。計(jì)算域網(wǎng)格采用均勻分布的笛卡爾網(wǎng)格,經(jīng)網(wǎng)格收斂性測試后,區(qū)域網(wǎng)格總數(shù)設(shè)置為3 000×890。本文中,設(shè)初始磁場強(qiáng)度為0.01T(β≈2 500),令平行于來流的方向?yàn)閤方向,垂直于來流的方向?yàn)閥方向。

    圖1 計(jì)算模型Fig.1 Computational model

    表1 氣體參數(shù)Table 1 Gas parameters used in this paper

    為了使氣體受磁場影響,首先通過高溫氣體電離,溫度為12 000K。由于含等離子體,電導(dǎo)率、熱導(dǎo)率、黏性系數(shù)、霍爾系數(shù)以及雙極擴(kuò)散系數(shù)由以下方法確定。電導(dǎo)率采用基于玻爾茲曼方程得到的Spitzer公式進(jìn)行計(jì)算,即由玻爾茲曼方程的碰撞積分出發(fā),忽略短程碰撞的影響,考慮剩余的碰撞項(xiàng),得到弱電場情況下由離子與電子所組成的等離子體電導(dǎo)率(σ)為:

    式中:j=m/(2kBT),其中m為質(zhì)量,kB為Boltzmann常數(shù);γ0=0.1Z+0.49;e為電子電荷量;Z為離子電荷數(shù)。

    如果碰撞積分到德拜長度為止,則庫侖對數(shù)項(xiàng)為:

    計(jì)算可得,初始電導(dǎo)率約為107S/m。

    當(dāng)氣體中相當(dāng)部分的原子被電離后,靜電力將對粒子間碰撞起主要作用,若考慮遠(yuǎn)程作用,則離子碰撞將不能視為兩體方式。通過求解玻爾茲曼方程得到一階近似電離氣體熱導(dǎo)率為:

    根據(jù)S.Chapman等[22]的理論可以推導(dǎo)出kc、ν的二階近似值。盡管二階近似值與一階近似值有顯著的差異,但其后各階近似所增加的修正量都是很微小的。假定黏度均由重離子引起,電子黏度很小,由于電子的熱運(yùn)動速度遠(yuǎn)大于離子的熱運(yùn)動速度,因此熱擴(kuò)散效應(yīng)主要是由電子引起。經(jīng)過估算:本文中黏性系數(shù)約為3.72×10-5Pa·s,熱導(dǎo)率約為1.4W/(m·K)。

    由于準(zhǔn)電中性條件,電子和離子在等離子體內(nèi)不可能獨(dú)立運(yùn)動,當(dāng)電子快速離開某等離子體單元時(shí),單元內(nèi)將形成由正電荷產(chǎn)生的電場,該電場阻礙電子繼續(xù)離開所在體積單元,并促使離子更快地離去,這種擴(kuò)散狀態(tài)稱為雙極型狀態(tài)。由于電子的擴(kuò)散系數(shù)遠(yuǎn)大于離子,因此:

    式中:ηi為正離子的熱擴(kuò)散系數(shù),Te、Ti分別為電子、離子溫度。由此可見,雙極擴(kuò)散系數(shù)恒小于電子自由擴(kuò)散系數(shù),而大于離子擴(kuò)散系數(shù),因此雙極性電場極大影響了電子的定向速度。本文中取雙極擴(kuò)散系數(shù)為(10-5/ρ)m2·Pa·s-1。

    Hall效應(yīng)可在磁流體通道中產(chǎn)生橫向洛倫茲力,而且體力的不均勻分布可導(dǎo)致二次流的產(chǎn)生,對流體的運(yùn)動產(chǎn)生重要影響,因此需要考慮霍爾效應(yīng)。此處,霍爾系數(shù)ηH取為10-7m3/C。

    2 結(jié)果與討論

    為了驗(yàn)證數(shù)值方法與計(jì)算模型的可靠性,圖2給出了不加磁場時(shí)入射激波與R22氣柱相互作用過程計(jì)算結(jié)果與文獻(xiàn)[21]中實(shí)驗(yàn)結(jié)果的對比,兩者具有相同的初始條件。由圖2可知:計(jì)算結(jié)果中激波的反射、繞射及發(fā)展過程均與實(shí)驗(yàn)結(jié)果完全相符,兩者氣柱的變形過程也相似;但是,由于初始時(shí)兩種氣體接觸界面采用間斷面形式,因此計(jì)算結(jié)果中接觸界面上的KH不穩(wěn)定性更明顯,從而導(dǎo)致界面上、下兩側(cè)的湍流轉(zhuǎn)捩過程(見圖2(e))與實(shí)驗(yàn)結(jié)果稍有不同。

    圖2無磁場情況下,激波與R22氣柱相互作用過程的計(jì)算紋影(上)與文獻(xiàn)[21]實(shí)驗(yàn)結(jié)果(下)的對比Fig.2 Numerical results of the interaction of shock wave with the R22air column in this paper(upper row),and the experimental result in reference[21](lower row).

    圖3 為后期激波與R22氣柱作用過程的計(jì)算紋影。從圖3可以看出:由于激波對不同密度氣體界面的作用,界面產(chǎn)生了RM不穩(wěn)定性;在作用初期,隨機(jī)擾動在氣體界面發(fā)生,擾動振幅在界面線性增長;隨著變形的加劇,振幅經(jīng)歷非線性增長,“氣泡”與“尖釘”結(jié)構(gòu)出現(xiàn);隨著運(yùn)動的繼續(xù)進(jìn)行,在KH不穩(wěn)定性的影響下,“尖釘”破碎,“氣泡”被拉長,最終形成湍流混合狀態(tài)。

    圖3后期激波與R22氣柱作用過程計(jì)算紋影Fig.3 Numerical results of the interaction of shock wave with the R22air column in later period.

    圖4 為加入均勻磁場后,激波與R22氣柱相互作用過程的密度紋影,其中上方子圖的磁場方向垂直于來流方向(即磁場方向沿y軸正向),而下方子圖中磁場方向平行于來流方向。由圖4可知,與無磁場情況(見圖2)相比,磁作用力能明顯抑制界面處的KH不穩(wěn)定性,且垂直于來流方向的磁作用力對不穩(wěn)定性的控制較平行于來流方向的磁作用力更明顯。為了更清晰地描述此現(xiàn)象,將動量方程中與磁場相關(guān)的項(xiàng)拆分為磁壓力項(xiàng)(B·B)/(2μ)與磁張力項(xiàng)·(BB)/μ兩部分,這兩項(xiàng)均為二階張量。然而,由于計(jì)算結(jié)果中磁壓力值遠(yuǎn)大于磁張力值,下文只對磁壓力的變化情況進(jìn)行討論。

    由式(4)可知,垂直于速度方向的磁場分量梯度較大,磁場變化較快,該效應(yīng)會引起磁場在某些位置(如垂直于磁場方向的界面處)的集中。因此,隨著時(shí)間的推移,氣泡界面處的磁場強(qiáng)度不斷增大,妨礙輕質(zhì)氣體和重質(zhì)氣體的相互混合。由圖4可知:當(dāng)初始磁場垂直于來流方向時(shí),磁壓力作用在x軸方向,可以有效阻止氣柱兩側(cè)的卷起;在t=450μs之后,由于斜壓作用加劇,氣柱尾部界面開處始形成渦結(jié)構(gòu);然而,沿x軸方向的磁壓力抑制了不穩(wěn)定性的進(jìn)一步發(fā)展,最終在氣柱尾部界面附近形成了上、下兩個(gè)突起,如圖4(f)所示。當(dāng)磁場平行于來流方向時(shí),磁壓力作用在y軸方向,因上側(cè)界面處x軸向速度的差異,在界面處卷起渦串,如圖4(c)和圖4(d)所示;同時(shí),由于斜壓效應(yīng),氣柱上、下兩側(cè)卷起兩個(gè)渦環(huán),如圖4(e)所示;隨后,沿y軸方向的磁壓力阻止了主渦的進(jìn)一步發(fā)展,最終在主渦位置形成“尖釘”狀結(jié)構(gòu),見圖4(f)。由于磁化氣柱的熱導(dǎo)率遠(yuǎn)大于空氣的,界面處氣體混合速度較無磁場時(shí)更快。

    圖4 B=0.01T時(shí)不同磁場方向?qū)げㄅcR22氣柱作用過程的影響Fig.4 Numerical results of the interaction process of shock wave and R22air column when B=0.01T

    圖5和圖6分別給出了磁場方向垂直或平行于來流方向時(shí)氣柱對稱軸上壓力(p)與x方向速度(u)的分布。由圖5和圖6可知,與無磁場的情況[23]相比,初始壓力與速度變化不大。當(dāng)t=200μs時(shí),透射激波在x軸上發(fā)生聚焦,并且在匯聚處產(chǎn)生局部高壓區(qū):初始磁場沿y軸正方向時(shí),壓力曲線峰值pm=700kPa,初始磁場沿x軸正方向時(shí),pm=760kPa,均小于無磁場情況。此時(shí),x方向速度同樣處于峰值,但施加沿x軸方向磁場時(shí)的速度(u=153m/s)遠(yuǎn)大于無磁場(u=100m/s)及施加沿y軸方向磁場時(shí)的速度(u=125m/s)。隨著相交點(diǎn)向下游移動,聚焦點(diǎn)的壓力峰值不斷下降,速度逐漸增大,并在t=250μs時(shí)達(dá)到最大值(u≈170m/s)。隨后,在流體阻力作用下,速度不斷衰減。綜上所述,磁場的加入對壓力影響不大,但卻能提升透射激波聚焦時(shí)速度峰值,而磁場沿x軸方向時(shí)影響更劇烈。此外,在激波衰減過程中,沿y方向磁場較沿x方向磁場更能減緩射流的衰減。

    圖5 初始磁場沿y軸正方向時(shí)入射激波誘導(dǎo)射流產(chǎn)生過程Fig.5 Generation process of jet induced by incident shock wave when initial magnetic field is along yaxis

    圖6 初始磁場沿x軸正方向時(shí)入射激波誘導(dǎo)射流產(chǎn)生過程Fig.6 Generation process of jet induced by incident shock wave when initial magnetic field is along xaxis

    圖7 給出了磁場沿y軸正方向或x軸正方向時(shí)激波與R22氣柱作用過程的計(jì)算紋影。由圖7可知,磁場能夠較好抑制RM不穩(wěn)定性。當(dāng)初始磁場垂直于來流方向時(shí),界面處的橫向磁壓力對氣體混合的抑制效果較好,同時(shí)氣柱截面被壓扁拉長,兩個(gè)主渦在磁場作用下受到抑制,最終被拉成長條狀;當(dāng)初始磁場平行于來流方向時(shí),后期氣柱寬度小于磁場垂直于來流方向的情形,雖然形成兩個(gè)主渦,但磁場抑制了氣柱的破碎及次級渦的生成。由此可見,加入磁場后,氣柱界面處沒有在KH不穩(wěn)定性的作用下持續(xù)失穩(wěn)并卷起“珠狀”次級渦,氣柱頂端也未出現(xiàn)RM不穩(wěn)定性,因此,兩種方向的初始磁場對不穩(wěn)定性都有較好的抑制作用。

    為更清楚地看到磁場對氣柱界面的影響,圖8給出了初始磁場沿y軸正方向時(shí)不同時(shí)刻磁壓力和磁能量的分布。由圖8可知,初始時(shí)刻,透射激波在氣柱內(nèi)部向下游傳播,圓弧反射激波向上游傳播,磁壓力較大處集中在氣柱界面處,磁能量較大處集中在界面處及透射激波尾端。當(dāng)t=165μs時(shí),x方向磁壓力沿x軸對稱分布,而y方向磁壓力沿x軸反對稱分布;同時(shí)x方向磁壓力對氣體界面內(nèi)側(cè)產(chǎn)生向右拉力,對氣體界面外側(cè)則產(chǎn)生向左拉力,而y方向磁壓力對氣體界面內(nèi)側(cè)產(chǎn)生向氣柱中心的拉力,而對氣體界面外側(cè)產(chǎn)生向外的拉力。另外,當(dāng)入射激波與氣柱作用時(shí),因氣柱內(nèi)外氣體聲阻抗的不同,氣柱內(nèi)外會產(chǎn)生速度差。氣柱外部界面處,空氣x方向速度大于界面速度,而氣柱內(nèi)部界面處,空氣x方向速度小于當(dāng)前界面速度。由于磁場變化率與速度差密切相關(guān),因此這兩處的速度梯度會帶來磁場的變化,由此產(chǎn)生的磁壓力進(jìn)一步抑制氣柱界面受內(nèi)外剪切力的影響(同理,對空氣y方向速度,也有類似結(jié)論)。由圖7(c)和圖7(d)可知,發(fā)展后期,磁能量與磁壓力分布類似,包裹在界面處。綜上所述,當(dāng)初始磁場垂直于來流方向時(shí),磁壓力始終包裹著界面,阻止界面受周圍環(huán)境的進(jìn)一步影響,從而防止界面上渦量的生成和擴(kuò)散,抑制界面的不穩(wěn)定性。

    圖7 后期激波與R22氣柱作用過程密度紋影Fig.7 Numerical results of interaction process of shock wave with R22air column in later period

    圖8 初始磁場沿y軸正方向時(shí)磁壓力及磁能量的分布(每幅子圖中,左上為磁壓力,左下為磁能量,中間和右側(cè)分別為磁壓力在x、y方向的分量)Fig.8 Magnetic pressure and magnetic energy when initial magnetic field is along yaxis(upper left:magnetic pressure;lower left:magnetic energy;middle:magnetic pressure along xaxis;right:magnetic pressure along yaxis)

    圖9為初始磁場沿y軸正方向時(shí)磁壓力(pm)和磁能量(Em)沿橫向(x方向)及縱向(y方向)的分布情況。由圖9可知:當(dāng)透射激波在氣柱內(nèi)傳播時(shí),磁能量在兩個(gè)位置取值較大,一處是氣柱左側(cè)界面處,另一處位于透射激波尾部,與之前分析相符;界面處的磁壓力遠(yuǎn)大于其他位置。由圖9可知,氣柱上、下界面處的磁壓力及磁能量遠(yuǎn)大于左、右界面處。圖9(c)顯示:當(dāng)透射及入射激波先后聚焦后,會在氣柱尾端射流位置形成局部高磁壓力區(qū)域,同時(shí),該區(qū)域磁能量很低,接近于零。由于圖9(c)和圖9(d)中射流位置的磁壓力及磁能量很低,因此磁場對射流并無明顯影響;但是射流前、后的磁能量可以一定程度上抑制射流的衰減。此外,圖9中不同時(shí)刻磁壓力和磁能量分布曲線的形狀無明顯變化,但上、下界面處磁壓力和磁能量的峰值隨時(shí)間不斷增加,這意味著界面處磁場強(qiáng)度隨時(shí)間不斷增大。

    圖10為初始磁場沿x軸正方向時(shí)磁壓力及磁能量的分布。由圖10可知,磁壓力與磁能量的分布類似,其峰值均集中在界面附近。與初始磁場沿y軸正方向的情況相似,x方向磁壓力沿x軸對稱分布,y方向磁壓力沿x軸反對稱分布,同時(shí)磁壓力與磁能量的作用區(qū)域更緊貼界面。由于磁壓力較弱且分布不均勻,圖11給出了邊界處磁壓力矢量分布。由圖11(a)可知,當(dāng)磁場平行于來流方向時(shí),初始磁壓力分布于界面處,且與剪切力及外界對其壓力的合力方向相反。隨后,在剪切力的作用下,氣柱表面開始卷起渦,磁壓力的方向隨著界面變形而改變,如圖11(b)所示。在主渦生成時(shí),磁壓力較大的位置依在界面附近,此對渦中心部分影響不大,如圖11(c)和圖11(d)所示。由此可見:當(dāng)渦量大于零時(shí)(渦逆時(shí)針方向),磁壓力指向渦動中心;而當(dāng)渦量小于零時(shí),磁壓力指向外部。

    圖9 初始磁場沿y軸正方向時(shí)的磁壓力及磁能量分布Fig.9 Distribution of magnetic pressure and magnetic energy when initial magnetic field is along yaxis

    圖10 初始磁場沿x軸正方向時(shí)磁壓力及磁能量分布(每幅子圖中,左上為磁壓力,左下為磁能量,中間和右側(cè)分別為磁壓力在x、y方向的分量)Fig.10 Magnetic pressure and magnetic energy when initial magnetic field is along xaxis(upper left:magnetic pressure;lower left:magnetic energy;middle:magnetic pressure along xaxis;right:magnetic pressure along yaxis)

    圖12給出了初始磁場沿x軸正方向時(shí)磁壓力及磁能量沿橫向及縱向的分布。由圖12(a)和圖12(b)可知,初始時(shí)刻,橫軸和縱軸上磁壓力及磁能量均小于初始磁場沿y軸正方向的情況。當(dāng)t=165μs時(shí),橫軸上的磁壓力集中在透射激波尾部,磁能量在頂部界面后迅速減小,接近透射激波時(shí)快速增大并達(dá)到峰值;而縱軸上磁能量和磁壓力的峰值均在上、下界面處。當(dāng)t=200μs時(shí),橫軸上的磁壓力集中在氣柱尾端,磁能量仍然在透射激波尾部,縱軸上的磁壓力不斷增大;此時(shí)由于渦在界面處開始生成,縱軸上磁壓力的分布不對稱。當(dāng)t=370μs時(shí),透射激波聚集,在尾部形成射流,圓弧形二次激波向外傳播,并最終碰到上、下界面發(fā)生反射;此時(shí),橫軸上磁壓力與磁能量最大值位置相同,均在氣柱尾部界面處達(dá)到峰值,而縱軸上的磁壓力峰值則繼續(xù)增大,如圖12(c)所示。當(dāng)兩個(gè)主渦形成后,縱軸處主渦位置的磁能量較為明顯,但此時(shí)磁能量及磁壓力最大位置仍在界面附近;橫軸上磁壓力峰值減小,磁能量繼續(xù)增加,并集中在氣柱尾部,如圖12(d)所示。

    圖11 初始磁場沿x軸正方向時(shí)邊界處磁壓力矢量圖Fig.11 Vector diagram of magnetic pressure on the interface when initial magnetic field is along xaxis

    圖12 初始磁場為x軸正方向時(shí)磁壓力及磁能量分布Fig.12 Distribution of magnetic pressure and magnetic energy when initial magnetic field is along xaxis

    3 結(jié) 論

    基于CTU+CT算法,對磁場環(huán)境下激波沖擊重質(zhì)氣體進(jìn)行數(shù)值模擬,并與不加磁場的情況進(jìn)行對比,所得結(jié)論如下。

    (1)R22氣柱與入射激波作用后,由于KH不穩(wěn)定性,氣柱上、下兩側(cè)不斷卷起最終形成兩個(gè)主渦。同時(shí),由于RM不穩(wěn)定性,后期氣柱左端界面處產(chǎn)生了“尖釘”與“氣泡”狀結(jié)構(gòu)。磁場能夠抑制界面不穩(wěn)定性,平行磁場通過減少界面上渦量的生成減緩不穩(wěn)定性的增長;垂直磁場則通過包裹界面的磁力線壓縮界面,抑制剪切力對界面的作用,使界面處不再卷起,從而抑制界面不穩(wěn)定性。

    (2)垂直磁場和平行磁場都可以很好地抑制RM不穩(wěn)定性,對于KH不穩(wěn)定性,垂直磁場的控制效果更好,不僅可以抑制界面上渦串的卷起,還可以阻止主渦的發(fā)展,而平行磁場做不到后者。磁場對射流影響不大,射流處的磁能量可以在一定程度上抑制射流的衰減,同時(shí),垂直磁場可以減小射流聚焦時(shí)壓力及速度峰值。

    (3)當(dāng)磁場強(qiáng)度較?。é隆? 500)時(shí),初始磁場對激波與R22氣柱作用過程沒有顯著影響,隨后磁壓力及磁能量不斷增加,從而抑制不穩(wěn)定性的發(fā)展。磁壓力及磁能量對界面影響很大,在上、下界面處達(dá)到峰值。同時(shí),磁壓力對界面沒有壓縮作用,而是產(chǎn)生一個(gè)反向作用力減小外界對界面的影響。當(dāng)渦形成后,磁場對渦量影響不大,但可以通過抑制外界的剪切力阻止其進(jìn)一步的發(fā)展。

    [1] Meyer K A,Blewett P J.Numerical investigation of the stability of a shock-accelerated interface between two fluids[J].Physics of Fluids,1972,15(5):753-759.

    [2] Zhang Q,Sohn S I.An analytical nonlinear theory of Richtmyer-Meshkov instability[J].Physics Letters A,1996,212(3):149-155.

    [3] Anuchina N N,Volkov V I,Gordeychuk V A,et al.Numerical simulations of Rayleigh-Taylor and Richtmyer-Meshkov instability using MAH-3code[J].Journal of Computational and Applied Mathematics,2004,168(1/2):11-20.

    [4] Ruev G A,F(xiàn)edorov A V,F(xiàn)omin V M.Development of the Richtmyer-Meshkov instability upon interaction of a diffusion mixing layer of two gases with shock waves[J].Journal of Applied Mechanics and Technical Physics,2005,46(3):307-314.

    [5] Niederhaus J H J,Greenough J A,Oakley J G,et al.A computational parameter study for the three-dimensional shock-bubble interaction[J].Journal of Fluid Mechanics,2008,594:85-124.

    [6] Thornber B,Drikakis D,Youngs D.Large-eddy simulation of multi-component compressible turbulent flows using high resolution methods[J].Computers and Fluids,2008,37(7):867-876.

    [7] Hejazialhosseini B,Rossinelli D,Bergdorf M,et al.High order finite volume methods on wavelet-adapted grids with local time-stepping on multicore architectures for the simulation of shock-bubble interactions.[J].Journal of Computational Physics,2010,229(22):8364-8383.

    [8] Schilling O,Latini M.High-order weno simulations of three-dimensional reshocked Richtmyer-Meshkov instability to late times:Dynamics,dependence on initial conditions,and comparisons to experimental data[J].Acta Mathematica Scientia,2010,30(2):595-620.

    [9] Tian B,Shen W,Jiang S,et al.A global arbitrary Lagrangian-Eulerian method for stratified Richtmyer-Meshkov instability[J].Computers and Fluids,2011,46(1):113-121.

    [10] Shankar S K,Kawai S,Lele S K.Two-dimensional viscous flow simulation of a shock accelerated heavy gas cylinder[J].Physics of Fluids,2011,23(2):131.

    [11] Bailie C,Mcfarland J A,Greenough J A,et al.Effect of incident shock wave strength on the decay of Richtmyer-Meshkov instability-introduced perturbations in the refracted shock wave[J].Shock Waves,2012,22(6):511-519.

    [12] Chandrasekhar S.Hydrodynamic and hydromagnetic stability[M].Dover Publications,1961.

    [13] Wheatley V,Pullin D I,Samtaney R.Stability of an impulsively accelerated density interface in magnetohydrodynamics[J].Physical Review Letters,2005,95(12):125002.

    [14] Wheatley V,Samtaney R,Pullin D I.The magnetohydrodynamic Richtmyer-Meshkov instability:The transverse field case[C]∥18th Australasian Fluid Mechanics Conference.Australasian Fluid Mechanics Society,2012.

    [15] Khan M,Mandal L,Banerjee R,et al.Development of Richtmyer-Meshkov and Rayleigh-Taylor instability in presence of magnetic field[J].Nuclear Instruments &Methods in Physics Research A,2011,653(1):2-6.

    [16] Cao J,Wu Z,Ren H,et al.Effects of shear flow and transverse magnetic field on Richtmyer-Meshkov instability[J].Physics of Plasmas,2008,15(4):445-514.

    [17] Shin M S,Stone J M,Snyder G F.The magnetohydrodynamics of shock-cloud interaction in three dimensions[J].Astrophysical Journal,2008,680(1):336-348.

    [18] 李源,羅喜勝.黏性、表面張力和磁場對Rayleigh-Taylor不穩(wěn)定性氣泡演化影響的理論分析[J].物理學(xué)報(bào),2014,2(8):277-285.Li Yuan,Luo Xisheng.Theoretical analysis of effects of viscosity,surface tension,and magnetic field on the bubble evolution of Rayleigh-Taylor instability[J].Acta Physica Sinica,2014,2(8):277-285.

    [19] Saltzman J.An unsplit 3Dupwind method for hyperbolic conservation laws[J].Journal of Computational Physics,1994,115(1):153-168.

    [20] Gardiner T A,Stone J M.An unsplit Godunov method for ideal MHD via constrained transport[J].Journal of Computational Physics,2005,205(2):509-539.

    [21] Haas J,Sturtevant B.Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[J].Journal of Fluid Mechanics,1987,181(1):41-76.

    [22] Chapman S,Cowling T G.The mathematical theory of non-uniform gases[M].London:Cambridge University Press,1970.

    [23] 沙莎,陳志華,薛大文.激波沖擊R22重氣柱所導(dǎo)致的射流與混合研究[J].物理學(xué)報(bào),2013,62(14):144701.Sha Sha,Chen Zhihua,Xue Dawen.The generation of jet and mixing induced by the interaction of shock wave with R22cylinder[J].Acta Physica Sinica,2013,62(14):144701.

    Influence of magnetic field on interaction of shock wave
    with R22heavy gas column

    Lin Zhenya,Zhang Huanhao,Chen Zhihua,Liu Ying

    (National Key Laboratory of Transient Physics,Nanjing University of Science and Technology,Nanjing210094,Jiangsu,China)

    To study the interaction process of the plane incident shock wave with the magnetized R22 heavy gas column,we numerically simulated the deformation process of the shock-wave-induced R22 heavy gas column resulting from Kelvin-Helmholtz(KH)and Richtmyer-Meshkov(RM)instabilities under different initial conditions,and analyzed the jet focusing and inducing process by the transmitted shock wave.When the magnetic field was taken into consideration,the CTU+CT algorithm satisfying the divergence equation of the magnetic field at any time was adopted in the numerical simulation.The results show that the magnetic field is capable of restraining the instability of the shockwave-induced R22gas column.Both the normal magnetic field(vertical to the flow direction)and the tangential magnetic field(parallel to the flow direction)can inhibit the RM instability.However,the restraining of the normal magnetic field is more effective than that of the tangential one with regard to the KH instability,as it can not only inhibit the vortex train rolling up on the interface but also prevent the bound vortex from developing.Besides,it is found that the magnetic field has little influence on the jet,and the magnetic energy at the jet point can suppress the jet attenuation to some extent while the normal magnetic field can reduce the peak pressure and velocity when the transmitted shock wave is focused.

    magnetization gas;magneto-h(huán)ydro-dynamics;magnetic field;instability;shock wave

    O354.5國標(biāo)學(xué)科代碼:1303520

    A

    10.11883/1001-1455(2017)04-0748-11

    (責(zé)任編輯 王玉鋒)

    2016-03-24;

    2016-08-23

    國家自然科學(xué)基金項(xiàng)目(11272156);中國博士后科學(xué)基金項(xiàng)目(2015M571757)

    林震亞(1990- ),男,博士;通信作者:張煥好,122488989@qq.com。

    猜你喜歡
    氣柱不穩(wěn)定性激波
    預(yù)壓縮氣墊包裝系統(tǒng)靜力及動力學(xué)特性研究
    包裝工程(2024年1期)2024-01-20 06:17:38
    激波誘導(dǎo)雙層氣柱演化的偏心效應(yīng)研究
    氣體物理(2022年2期)2022-03-31 12:49:16
    一種基于聚類分析的二維激波模式識別算法
    基于HIFiRE-2超燃發(fā)動機(jī)內(nèi)流道的激波邊界層干擾分析
    磁控條件下激波沖擊三角形氣柱過程的數(shù)值研究?
    斜激波入射V形鈍前緣溢流口激波干擾研究
    可壓縮Navier-Stokes方程平面Couette-Poiseuille流的線性不穩(wěn)定性
    適于可壓縮多尺度流動的緊致型激波捕捉格式
    增強(qiáng)型體外反搏聯(lián)合中醫(yī)辯證治療不穩(wěn)定性心絞痛療效觀察
    前列地爾治療不穩(wěn)定性心絞痛療效觀察
    欧美中文综合在线视频| 午夜福利在线免费观看网站| 人妻一区二区av| 午夜福利,免费看| 一级毛片黄色毛片免费观看视频| 高清黄色对白视频在线免费看| 亚洲av中文av极速乱| 日本av免费视频播放| 成人国产麻豆网| 成年女人毛片免费观看观看9 | 99热网站在线观看| 99香蕉大伊视频| 日韩精品有码人妻一区| 日韩精品有码人妻一区| 在线观看三级黄色| 有码 亚洲区| 亚洲欧洲精品一区二区精品久久久 | 9热在线视频观看99| 国产黄色视频一区二区在线观看| 亚洲国产成人一精品久久久| 丝袜喷水一区| 亚洲人成网站在线观看播放| 男的添女的下面高潮视频| 色吧在线观看| 精品卡一卡二卡四卡免费| 亚洲精品中文字幕在线视频| 久久久久国产网址| 一区在线观看完整版| 男女国产视频网站| 最近2019中文字幕mv第一页| 蜜桃在线观看..| 成人毛片60女人毛片免费| 99热网站在线观看| 最近手机中文字幕大全| 日本欧美视频一区| 免费观看av网站的网址| 精品久久蜜臀av无| 亚洲一区中文字幕在线| 日韩欧美一区视频在线观看| 天堂中文最新版在线下载| 免费在线观看黄色视频的| 久久青草综合色| 黄网站色视频无遮挡免费观看| 精品一区二区三卡| 电影成人av| 亚洲综合色网址| 9191精品国产免费久久| 韩国精品一区二区三区| 成人亚洲精品一区在线观看| 久久久久久伊人网av| 校园人妻丝袜中文字幕| 啦啦啦啦在线视频资源| 亚洲精品美女久久久久99蜜臀 | 大陆偷拍与自拍| 久久精品国产自在天天线| 国产综合精华液| 久久人人爽人人片av| 成年女人在线观看亚洲视频| 男女啪啪激烈高潮av片| 黑人欧美特级aaaaaa片| 亚洲欧美日韩另类电影网站| 黄片小视频在线播放| 日韩av不卡免费在线播放| 色哟哟·www| 国产老妇伦熟女老妇高清| 少妇精品久久久久久久| 男女下面插进去视频免费观看| 日韩中文字幕视频在线看片| 亚洲国产精品成人久久小说| 少妇被粗大的猛进出69影院| 国产黄频视频在线观看| 久久 成人 亚洲| 飞空精品影院首页| 捣出白浆h1v1| 日韩中字成人| av卡一久久| 久热久热在线精品观看| 精品人妻熟女毛片av久久网站| 国产精品偷伦视频观看了| 亚洲av电影在线观看一区二区三区| 最近2019中文字幕mv第一页| 99国产精品免费福利视频| 91国产中文字幕| 中文字幕最新亚洲高清| 韩国av在线不卡| 国产成人午夜福利电影在线观看| 制服诱惑二区| 黄色毛片三级朝国网站| 亚洲精品一区蜜桃| 日本-黄色视频高清免费观看| 国产成人免费观看mmmm| 下体分泌物呈黄色| 熟女av电影| 香蕉国产在线看| xxxhd国产人妻xxx| 免费黄网站久久成人精品| 久久精品人人爽人人爽视色| 一本色道久久久久久精品综合| 99久久精品国产国产毛片| 欧美成人午夜免费资源| 男女免费视频国产| 国产视频首页在线观看| 久久久久久久国产电影| 99精国产麻豆久久婷婷| 亚洲国产欧美在线一区| 狂野欧美激情性bbbbbb| 欧美激情 高清一区二区三区| 9191精品国产免费久久| 欧美成人午夜免费资源| 满18在线观看网站| 欧美少妇被猛烈插入视频| 久久精品夜色国产| 不卡视频在线观看欧美| 精品国产一区二区久久| 99久久综合免费| 熟妇人妻不卡中文字幕| 亚洲欧美精品自产自拍| 国产亚洲最大av| 亚洲欧美日韩另类电影网站| 日韩熟女老妇一区二区性免费视频| 亚洲av国产av综合av卡| 免费不卡的大黄色大毛片视频在线观看| 啦啦啦在线免费观看视频4| 国产精品免费视频内射| 国产色婷婷99| 午夜福利在线观看免费完整高清在| av福利片在线| 午夜免费鲁丝| 亚洲国产色片| 亚洲国产精品国产精品| 国产探花极品一区二区| 青草久久国产| 人人妻人人添人人爽欧美一区卜| 国产精品久久久久久精品电影小说| 超碰97精品在线观看| 国产日韩欧美在线精品| 国产成人一区二区在线| 男女啪啪激烈高潮av片| 国产国语露脸激情在线看| 成人国语在线视频| 两个人看的免费小视频| 99久久精品国产国产毛片| 18禁动态无遮挡网站| 18禁国产床啪视频网站| 免费观看在线日韩| 欧美日韩一级在线毛片| 免费av中文字幕在线| 高清欧美精品videossex| 如何舔出高潮| 久热这里只有精品99| 国产精品av久久久久免费| 日韩伦理黄色片| 多毛熟女@视频| 少妇 在线观看| 国产成人一区二区在线| 视频区图区小说| 久久国内精品自在自线图片| 国产精品二区激情视频| 一级片'在线观看视频| 久久精品夜色国产| 性高湖久久久久久久久免费观看| 多毛熟女@视频| 亚洲精品中文字幕在线视频| 亚洲第一区二区三区不卡| av网站在线播放免费| 免费在线观看视频国产中文字幕亚洲 | xxxhd国产人妻xxx| 久久精品国产鲁丝片午夜精品| 欧美人与善性xxx| 亚洲国产精品一区二区三区在线| 大话2 男鬼变身卡| 国产国语露脸激情在线看| 精品第一国产精品| 王馨瑶露胸无遮挡在线观看| 日韩制服丝袜自拍偷拍| 国产成人精品福利久久| 日产精品乱码卡一卡2卡三| 在线观看三级黄色| 曰老女人黄片| 97在线人人人人妻| 五月开心婷婷网| 又大又黄又爽视频免费| 91精品国产国语对白视频| 性色avwww在线观看| 下体分泌物呈黄色| 精品一品国产午夜福利视频| 亚洲精品视频女| 两个人看的免费小视频| 美女国产高潮福利片在线看| 国产日韩欧美在线精品| 伊人亚洲综合成人网| 精品一区在线观看国产| 国产精品免费视频内射| 天天躁狠狠躁夜夜躁狠狠躁| 国产有黄有色有爽视频| 色婷婷av一区二区三区视频| 日韩欧美精品免费久久| av线在线观看网站| 日韩,欧美,国产一区二区三区| 亚洲精品视频女| 黑人欧美特级aaaaaa片| 99久久人妻综合| 岛国毛片在线播放| 欧美精品一区二区免费开放| 成人漫画全彩无遮挡| 黄色 视频免费看| 90打野战视频偷拍视频| 国产麻豆69| 国产综合精华液| 免费看av在线观看网站| 国产一区二区三区综合在线观看| 亚洲色图 男人天堂 中文字幕| 多毛熟女@视频| 在线观看www视频免费| 午夜福利,免费看| 日韩伦理黄色片| 中文字幕精品免费在线观看视频| 成年av动漫网址| 亚洲成人av在线免费| 日韩制服骚丝袜av| 国产日韩一区二区三区精品不卡| 在现免费观看毛片| 成年动漫av网址| 亚洲av国产av综合av卡| 成人影院久久| 亚洲精品成人av观看孕妇| 啦啦啦啦在线视频资源| 秋霞伦理黄片| 大香蕉久久成人网| xxxhd国产人妻xxx| 国产亚洲精品第一综合不卡| 美女高潮到喷水免费观看| 老鸭窝网址在线观看| 婷婷色综合www| 久热久热在线精品观看| 午夜福利,免费看| 如日韩欧美国产精品一区二区三区| 搡老乐熟女国产| 最近的中文字幕免费完整| 久久婷婷青草| 这个男人来自地球电影免费观看 | 日本vs欧美在线观看视频| 天天躁日日躁夜夜躁夜夜| 男女午夜视频在线观看| 纵有疾风起免费观看全集完整版| 色吧在线观看| 一区二区三区乱码不卡18| 人妻 亚洲 视频| 亚洲av成人精品一二三区| 欧美日韩av久久| 亚洲经典国产精华液单| 美女xxoo啪啪120秒动态图| av在线老鸭窝| 久久精品久久久久久噜噜老黄| 777久久人妻少妇嫩草av网站| 国产男女超爽视频在线观看| 在线观看美女被高潮喷水网站| 午夜福利视频在线观看免费| 久久精品久久久久久噜噜老黄| 777久久人妻少妇嫩草av网站| 久久国产精品大桥未久av| 国产日韩一区二区三区精品不卡| 欧美在线黄色| 少妇人妻久久综合中文| 中文字幕亚洲精品专区| 大香蕉久久成人网| 狂野欧美激情性bbbbbb| 精品福利永久在线观看| 好男人视频免费观看在线| 1024香蕉在线观看| 午夜久久久在线观看| 国产爽快片一区二区三区| 久久免费观看电影| 久久精品国产综合久久久| 国产探花极品一区二区| 国产日韩一区二区三区精品不卡| 一区二区日韩欧美中文字幕| 久久久久国产一级毛片高清牌| 18禁动态无遮挡网站| 麻豆av在线久日| 中文字幕精品免费在线观看视频| 午夜福利在线免费观看网站| 啦啦啦视频在线资源免费观看| 宅男免费午夜| 大片电影免费在线观看免费| 男女午夜视频在线观看| 精品一区二区三区四区五区乱码 | 天堂8中文在线网| 在线观看一区二区三区激情| 亚洲av.av天堂| 国产精品久久久久久久久免| 久久精品亚洲av国产电影网| 侵犯人妻中文字幕一二三四区| 欧美日韩综合久久久久久| 久久久久久久精品精品| 91精品三级在线观看| 欧美国产精品一级二级三级| 国产淫语在线视频| 日韩中文字幕视频在线看片| 精品亚洲成a人片在线观看| 久久精品国产自在天天线| 亚洲精品成人av观看孕妇| 国产女主播在线喷水免费视频网站| 99九九在线精品视频| 国产一区二区三区综合在线观看| 亚洲情色 制服丝袜| 精品一品国产午夜福利视频| 中文精品一卡2卡3卡4更新| 黑人巨大精品欧美一区二区蜜桃| 最黄视频免费看| 久久久精品区二区三区| 欧美亚洲 丝袜 人妻 在线| 五月天丁香电影| 大陆偷拍与自拍| 午夜福利视频精品| 日韩,欧美,国产一区二区三区| 中国国产av一级| 欧美在线黄色| av又黄又爽大尺度在线免费看| 女人高潮潮喷娇喘18禁视频| 免费播放大片免费观看视频在线观看| 欧美成人精品欧美一级黄| 青草久久国产| 天堂俺去俺来也www色官网| 亚洲一码二码三码区别大吗| 曰老女人黄片| 午夜福利网站1000一区二区三区| 一本大道久久a久久精品| 久久国产精品大桥未久av| av国产精品久久久久影院| 久久毛片免费看一区二区三区| xxx大片免费视频| 丝瓜视频免费看黄片| 欧美日韩精品成人综合77777| 妹子高潮喷水视频| 香蕉国产在线看| 秋霞在线观看毛片| 亚洲av国产av综合av卡| 一级黄片播放器| 久久99一区二区三区| 欧美 亚洲 国产 日韩一| 久久精品人人爽人人爽视色| 毛片一级片免费看久久久久| 久久久久久久久久久免费av| 啦啦啦在线免费观看视频4| 伦精品一区二区三区| 成人18禁高潮啪啪吃奶动态图| 久久这里有精品视频免费| 各种免费的搞黄视频| 亚洲图色成人| 女人精品久久久久毛片| 国产亚洲av片在线观看秒播厂| 日本vs欧美在线观看视频| 亚洲美女搞黄在线观看| 久久精品久久久久久噜噜老黄| 晚上一个人看的免费电影| 丝袜人妻中文字幕| 亚洲欧美日韩另类电影网站| 亚洲人成电影观看| 青春草国产在线视频| 亚洲精品日韩在线中文字幕| 性色av一级| 亚洲人成电影观看| 亚洲精品第二区| 精品国产一区二区三区久久久樱花| 毛片一级片免费看久久久久| 免费高清在线观看日韩| 一本大道久久a久久精品| 叶爱在线成人免费视频播放| 丰满少妇做爰视频| 制服丝袜香蕉在线| 久久精品久久精品一区二区三区| 国产探花极品一区二区| 亚洲色图综合在线观看| 久久精品国产亚洲av涩爱| 青青草视频在线视频观看| 亚洲,欧美,日韩| 亚洲情色 制服丝袜| 黄片小视频在线播放| 丰满迷人的少妇在线观看| 欧美av亚洲av综合av国产av | 久久久久久久久免费视频了| 国产毛片在线视频| 天堂8中文在线网| 亚洲男人天堂网一区| 黑人猛操日本美女一级片| 久久婷婷青草| 国产日韩一区二区三区精品不卡| 国产精品二区激情视频| 亚洲,一卡二卡三卡| 女性生殖器流出的白浆| 午夜福利视频精品| 成人亚洲欧美一区二区av| av电影中文网址| 国产熟女欧美一区二区| 国产精品熟女久久久久浪| 大话2 男鬼变身卡| 成年人午夜在线观看视频| 国产成人精品久久久久久| 高清欧美精品videossex| 免费女性裸体啪啪无遮挡网站| 一本大道久久a久久精品| 久久久国产精品麻豆| 免费观看无遮挡的男女| 亚洲人成77777在线视频| 国产成人免费无遮挡视频| 伦理电影大哥的女人| 免费观看无遮挡的男女| 人妻一区二区av| 精品少妇久久久久久888优播| 欧美激情极品国产一区二区三区| 日韩欧美精品免费久久| 肉色欧美久久久久久久蜜桃| 久久久久网色| 黄频高清免费视频| 成年人午夜在线观看视频| 国产成人精品在线电影| 又黄又粗又硬又大视频| 高清黄色对白视频在线免费看| 国产成人aa在线观看| 亚洲av中文av极速乱| 青春草国产在线视频| 街头女战士在线观看网站| 丝袜在线中文字幕| av片东京热男人的天堂| 国产一区有黄有色的免费视频| 亚洲精品乱久久久久久| 亚洲精品久久成人aⅴ小说| www.精华液| 宅男免费午夜| 亚洲,一卡二卡三卡| 色94色欧美一区二区| 又黄又粗又硬又大视频| 老鸭窝网址在线观看| 青春草国产在线视频| 波多野结衣一区麻豆| 久久久久久人妻| 日韩精品有码人妻一区| 新久久久久国产一级毛片| 最近最新中文字幕免费大全7| 青春草视频在线免费观看| 有码 亚洲区| 久久久欧美国产精品| 搡老乐熟女国产| 精品少妇黑人巨大在线播放| 亚洲成人手机| 日本免费在线观看一区| av女优亚洲男人天堂| 在线观看免费日韩欧美大片| 亚洲欧美一区二区三区国产| 纵有疾风起免费观看全集完整版| 欧美 亚洲 国产 日韩一| 欧美人与善性xxx| 成人午夜精彩视频在线观看| 两性夫妻黄色片| 黄色视频在线播放观看不卡| av有码第一页| 夫妻午夜视频| 青春草视频在线免费观看| 久久久久视频综合| 国产亚洲精品第一综合不卡| 中文字幕色久视频| 丝袜喷水一区| 一本—道久久a久久精品蜜桃钙片| 日日啪夜夜爽| 日韩大片免费观看网站| 伊人久久国产一区二区| 一区在线观看完整版| a 毛片基地| 国产成人精品无人区| 欧美日韩精品网址| www.av在线官网国产| 91在线精品国自产拍蜜月| 精品亚洲乱码少妇综合久久| av在线观看视频网站免费| 2021少妇久久久久久久久久久| 日韩,欧美,国产一区二区三区| 卡戴珊不雅视频在线播放| 在线精品无人区一区二区三| 青青草视频在线视频观看| 亚洲,欧美,日韩| 中国国产av一级| 国产xxxxx性猛交| 十分钟在线观看高清视频www| 欧美日韩亚洲高清精品| 波多野结衣一区麻豆| 飞空精品影院首页| 亚洲国产精品一区二区三区在线| 啦啦啦中文免费视频观看日本| kizo精华| 久久久精品94久久精品| 波多野结衣av一区二区av| 欧美人与性动交α欧美精品济南到 | 精品少妇久久久久久888优播| 午夜福利一区二区在线看| av网站在线播放免费| av在线老鸭窝| 亚洲综合色惰| 乱人伦中国视频| 中文字幕人妻熟女乱码| 另类精品久久| 亚洲色图 男人天堂 中文字幕| 国产乱来视频区| freevideosex欧美| 9热在线视频观看99| 高清黄色对白视频在线免费看| 日韩免费高清中文字幕av| 欧美激情高清一区二区三区 | 免费日韩欧美在线观看| 你懂的网址亚洲精品在线观看| 欧美日韩精品网址| 最黄视频免费看| 国产成人精品一,二区| 中国国产av一级| 色94色欧美一区二区| a级毛片黄视频| 最黄视频免费看| 免费少妇av软件| 丰满迷人的少妇在线观看| 国产av精品麻豆| 水蜜桃什么品种好| 日本wwww免费看| 青春草视频在线免费观看| 色视频在线一区二区三区| 亚洲成人av在线免费| 亚洲欧美一区二区三区国产| 欧美bdsm另类| h视频一区二区三区| 一区二区三区四区激情视频| 欧美精品一区二区大全| 蜜桃国产av成人99| 精品久久蜜臀av无| 欧美日本中文国产一区发布| 又黄又粗又硬又大视频| 成人黄色视频免费在线看| 男女午夜视频在线观看| 亚洲情色 制服丝袜| 国产精品av久久久久免费| 建设人人有责人人尽责人人享有的| 日本wwww免费看| 成人午夜精彩视频在线观看| 久久午夜福利片| 精品福利永久在线观看| 欧美成人午夜免费资源| 欧美另类一区| 久久女婷五月综合色啪小说| 国产乱来视频区| 日日撸夜夜添| 欧美老熟妇乱子伦牲交| 国产免费一区二区三区四区乱码| 国产精品国产三级专区第一集| 国产乱人偷精品视频| 欧美精品av麻豆av| 亚洲成人av在线免费| 99热全是精品| 一级毛片 在线播放| 欧美激情 高清一区二区三区| 最近手机中文字幕大全| 黄片无遮挡物在线观看| 国产色婷婷99| 亚洲,一卡二卡三卡| 一级片免费观看大全| 亚洲国产成人一精品久久久| 免费久久久久久久精品成人欧美视频| 黄色配什么色好看| 五月天丁香电影| 日韩精品有码人妻一区| 考比视频在线观看| 美国免费a级毛片| 国产有黄有色有爽视频| 欧美精品亚洲一区二区| 一区二区三区精品91| 飞空精品影院首页| 一区在线观看完整版| 在线亚洲精品国产二区图片欧美| 午夜日本视频在线| 精品国产超薄肉色丝袜足j| 人人妻人人澡人人看| 欧美av亚洲av综合av国产av | 国产亚洲午夜精品一区二区久久| 成年动漫av网址| 欧美少妇被猛烈插入视频| 欧美变态另类bdsm刘玥| 久久影院123| 久久精品国产亚洲av涩爱| 国产成人免费无遮挡视频| 一级毛片黄色毛片免费观看视频| 色婷婷久久久亚洲欧美| 男女边摸边吃奶| 老汉色∧v一级毛片| 视频区图区小说| 免费黄频网站在线观看国产| 国产精品人妻久久久影院| 热99久久久久精品小说推荐| 我的亚洲天堂| 亚洲av电影在线进入| 精品久久久精品久久久| a级毛片在线看网站| 久久精品久久久久久久性| 久久精品久久久久久噜噜老黄| 亚洲国产毛片av蜜桃av| 欧美少妇被猛烈插入视频| 91精品三级在线观看| 日韩欧美精品免费久久| 边亲边吃奶的免费视频| 丰满饥渴人妻一区二区三| 婷婷色综合大香蕉| 亚洲欧美一区二区三区国产| 免费人妻精品一区二区三区视频| 国产欧美日韩一区二区三区在线| 亚洲美女视频黄频| 国产xxxxx性猛交| 有码 亚洲区| 午夜免费男女啪啪视频观看| 啦啦啦中文免费视频观看日本|