不同初始磁場對激波沖擊R22重氣柱過程影響的數值模擬*

2018-03-07 04:56:15林震亞郭則慶張煥好陳志華

爆炸與沖擊 2018年2期

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

(南京理工大學瞬態物理國家重點實驗室，江蘇南京 210094)

在激波穿過物質界面過程中，由于激波區域存在較強的壓力梯度，而物質界面區域存在較強的密度梯度，兩者的不平衡使激波作用下的物質界面失穩形成不同尺度的渦，并最終向湍流轉捩，從而增強了界面兩側流體的混合[1]，同時在界面處還會出現Richtmyer-Meshkov(RM)不穩定與Kelvin-Helmholtz(KH)不穩定現象，因而激波與界面相互作用是一個復雜的多尺度、強非線性的物理問題。此類現象廣泛存在于超燃沖壓發動機的燃料混合、慣性約束核聚變、水中炸藥爆炸以及天體物理中的超新星爆炸等領域，因此有著重要而廣泛的研究背景。

由于激波-界面相互作用問題在學術和工程領域的研究價值，對激波與不同密度的氣體界面相互作用問題有了大量研究。Haas等[2]對弱激波與R22重氣柱、氣泡以及He輕氣柱、氣泡的作用過程進行了實驗研究，對氣柱、氣泡變形過程進行了分析討論。Tomkins等[3]對激波與SF6重氣柱的作用過程進行了實驗研究，研究分析了兩種介質的混合機理。范美如等[4-5]對激波與矩形、橢圓、菱形以及兩種三角形五種不同形狀的SF6氣柱作用過程進行了數值研究，對比分析了這幾種形狀界面的波系、渦量以及氣體界面的演變。王顯圣等[6]數值研究了入射激波以及反射激波與SF6重氣柱作用過程。Si等[7]利用高速紋影技術對平面入射激波以及反射激波與He輕氣泡以及SF6重氣泡作用過程進行了實驗研究，結果清晰顯示了激波誘導氣泡變形的過程，并對整個過程中氣泡尺寸變化以及流場環量變化進行了定量分析。

非理想情況下關于RM不穩定性控制方法方面的研究尚處于起步階段，研究表明，合適的磁場對氣體界面RM不穩定性的產生和發展具有一定的抑制作用，研究仍主要集中在理論研究方面。Chandrasekhar[8]和Wheatley等[9-10]采用線性理論方法，推導了垂直或平行于界面的磁場對不可壓縮流體界面不穩定性的影響，得出平行于物質界面的磁場能夠抑制Rayleigh-Taylor(RT)和RM不穩定性的結論。Cao等[11]同時考慮了初始平行于界面的磁場和剪切流對界面不穩定性的影響，發現磁場能夠抑制界面不穩定性，而剪切流將會促進界面不穩定性的發展，兩者之間存在競爭機制。由于線性理論所得到的平面流場中的RM不穩定性與垂直于平面的磁場無關，這與實際不符，為了研究此種情況下磁場的效應，Khan等[12]在勢流理論的基礎上，建立了新的模型，分析發現在非線性效應下，磁場能夠影響界面不穩定性的發展。Shin等[13]對不同方向磁場中平面激波與球形氣泡相互作用問題進行了數值模擬，發現在界面發展初期，磁場效應很小，而到了發展后期，即使是弱磁場，對界面形狀的影響也很明顯，磁場越強，界面內外物質混合率越低，后期湍流發展被抑制程度越強。李源等[14]理論分析了磁場中非理想流體的RT不穩定性氣泡的演化過程，在與磁場垂直的平面中，綜合考慮黏性和表面張力的影響，推導了描述二維非理想磁流體RT不穩定性氣泡運動的控制方程組，分析了流體黏性、表面張力和磁場對氣泡界面不穩定性發展的影響。

目前，關于RM不穩定性的研究主要用于理論和機理分析，由于高維可壓方程沒有解析解，因此理論研究主要針對不可壓縮MHD方程；對于機理分析，不論是實驗還是數值模擬，目前的分析重心在激波、阿爾文波等物理現象上。由于非理想情況下無法計算阿爾文波，同時激波與磁場并無直接關聯，因此本文中側重于分析渦量、磁壓力與磁能量。本文中，基于MHD方程，采用CTU+CT算法，以渦量為切入點討論不同磁感應強度對激波沖擊R22重質氣柱過程中界面不穩定性的影響，揭示磁場對不穩定性控制的機理。

1 數值方法與計算模型

1.1 數值方法

采用非理想磁流體動力(MHD)方程組，模擬激波沖擊R22氣柱的作用過程，具體參見文獻[15]，并選用6解CTU+CT算法[16-17]對MHD方程組進行求解。為了清晰地說明程序運行方式，下面對6解CTU+CT算法的主要過程進行描述。

(1)

在y和z方向上的流量也有相似的表達形式。

第3步，在每個界面通過δt/2步長的橫向流量梯度生成PPM界面狀態，其中流體動力學變量(質量、動量及能量密度)為：

(2)

其中，x方向的MHD動量密度源項與能量密度源項分別為：

(3)

(4)

磁場成份在預測的流量區域通過CT電場生成，磁場界面組份通過斯托克斯循環的積分形式形成：

(5)

磁場y、z方向的組分也可由類似的方法得到。這些磁場橫向組份中的MHD源項起點可通過前文描述的感應方程的方向拆分處理。另外，動量與能量密度的MHD源項通過MHD方程的原始變量形式來計算PPM界面狀態。同樣，y、z方向的界面狀態也由類似的方式得到，只要按照上面的描述改變(x,y,z)和(i,j,k)的周期排列。

第4步，對于每個通過δt/2步長更新的界面狀態，計算相關的流量，給出x方向的界面流量：

(6)

y和z方向的界面流量也有相似的表達方式。

(7)

第6步，從n至n+1步更新結果，流體動力學變量(質量、動量和能量密度)通過標準流量積分關系進行迭代：

(8)

同時，磁場的界面平均組份通過斯托克斯循環積分進行迭代，x方向組份為：

(9)

y和z方向的組分也有相似的表達形式。

這相對簡單的3D綜合算法具有二階精度，且在磁場界面通常組份的演化過程中不包含源項。這種算法對于均勻網格的流動，并包含相關對稱性的問題，可降為2D CTU和1D PPM綜合算法。從實驗中觀察到的下降趨勢，對于該算法的CFL<1/2時穩定。

1.2 計算模型

圖1為激波與R22氣柱相互作用的二維計算模型，為了便于分析不同強度磁場對不穩定性的影響，選取與文獻[15]相同的計算域尺寸、氣體參數、初始條件、邊界條件以及網格劃分，并討論了磁感應強度分別為0.01與0.05 T(即β=2 500，100，其中β是磁感應強度的量綱一量，β-1=B2/(2μ0p0))時對流場的影響。另外，為了使氣體受磁場影響，先將氣體進行電離。由文獻[15]，其初始電導率取106S/m，熱導率取1.4 W/m2，黏性系數取3.72×10-5Pa·s，霍爾系數取6×10-6m3/C，同時，雙極擴散系數D隨粒子數密度ρ的增加而減小，因此雙極擴散系數取ρD=10-4m2·Pa/s。

2 結果與討論

圖2為無磁場時入射激波與R22氣柱相互作用過程的渦量(上)與密度紋影圖(下)。由圖可見，計算結果中激波的反射、繞射及發展過程均與數值模擬結果[18]完全相符，且兩者中氣柱的變形過程相似。

當入射激波與氣柱碰撞后，分別形成一道向上游傳播的反射激波與一道向氣柱內部傳播的圓弧透射激波(見圖2(a))。在激波作用下，氣柱在x軸方向上被不斷壓縮。當入射激波繞過氣柱上下兩側圓弧頂點后出現彎曲，形成繞射激波(見圖2(b))。隨后，入射激波在氣柱尾部聚焦，形成向下游傳播的圓弧二次激波，最終在氣柱尾部聚焦形成局部高壓區并產生射流(見圖2(c)～(f))，此過程與實驗結果[19]相符。另外，由渦量分布可見，當入射激波經過氣柱界面后，誘導界面發生RM不穩定而出現渦量(見圖2(a))，且在激波繞經界面上下兩頂點時，引起該處界面的劇烈失穩而卷起形成小渦串(見圖2(b))。隨后，上下界面處反射形成的反射激波又與氣柱的上下界面發生碰撞，加劇了界面的失穩而形成串級旋渦序列。當繞射激波在中心軸上聚焦并反射后，與氣柱持續作用形成兩個主渦(見圖2(c)～(e))。隨著激波在流場中的來回反射，使流場的波系結構復雜多變，并與氣柱發生相互作用，最終使氣柱界面完全失穩(見圖2(f))。后期，隨著界面左側頂端擾動的發展，逐漸形成尖釘狀與氣泡狀結構。

圖3為施加磁感應強度B=0.01 T后，激波與R22氣柱相互作用過程中的渦量分布圖。可見，由于磁場及熱傳導作用，使氣柱邊界在與入射激波作用后會出現內外兩個渦層，此現象在初始磁場垂直于來流方向(垂直磁場)時尤為明顯。對于垂直磁場，氣柱界面保持穩定且無渦量卷起，此時渦量均勻分布于入射激波作用過的界面處，隨后內側渦層逐漸衰減，外側渦層的渦量較大區域出現在尾部，最終在垂直磁場作用下界面的卷起得到抑制，阻止了不穩定性結構的出現。而對于初始磁場平行于來流方向(水平磁場)時，氣柱界面在上下頂端(即剪切力最大處)開始卷起，并在t=250 μs時形成4個渦。隨后，氣柱上下界面處的兩個渦逐漸混合，最終在尾部形成兩個主渦。由此可知，磁感應強度0.01 T的水平磁場雖未能抑制尾部主渦的形成，但卻很好抑制了界面處次級渦的出現。

圖4為B=0.05 T時的渦量分布。可知，當磁感應強度加大時，渦層分離得更快，且不再依附于氣柱界面上。對于垂直磁場情況，初始時的渦量分布(見圖4(a)～(c))與B=0.01 T時相似(見圖3(a)～(c))，但由于兩渦層出現后，外側渦層在入射激波及透射激波先后聚焦后向右運動，同時沿界面法向方向向外擴張，形成斜向后的塊狀結構，最終與上下壁面碰撞并反射，如圖4(d)～(f)所示。同時，內側渦層則沿y軸向內壓縮，最終變為橢圓形。此時，由于界面處渦量較小，無不穩定結構出現，因此不穩定性得到了有效的抑制。而對于水平磁場情況，物質界面上的渦量分布情況與B=0.01 T時相似。此時，氣柱上下兩端頂點處渦量最大，不同的是氣柱前端渦層迅速分離，形成兩個渦層。其中，外側渦層右端逐漸向氣柱尾端移動，逐漸包裹氣柱，而內側渦層尾部則隨著射流向右運動，最終將內側渦層拉長。在水平磁場作用下，后期氣柱尾部雖有卷起，但并未形成主渦(見圖4(f))，同時很好抑制了界面處的不穩定性。由此可見，當磁感應強度增大時，會加速內外兩個渦層分開，并在隨后的運動中脫離氣柱界面，從而更好地抑制激波與界面作用過程中的不穩定性。

圖5是t=945 μs時不同初始磁場情況下的密度紋影圖。可見，當磁感應強度增大時，對RM不穩定性的控制效果更好。同時，垂直磁場對不穩定性的控制效果優于平行磁場，它不僅阻止了主渦的卷起以及次級渦的生成，同時對界面也有更好的壓縮作用，因此對RM及KH不穩定性起到有效的抑制作用。此外，在磁場較小時，渦層附著于界面上，而隨著磁場的增大，渦層逐漸與界面脫離，從而有效地控制界面處的渦量形成。

表1為不同初始條件下不同時刻的平均渦量。當B=0.05,0.01 T時，平均渦量在透射激波聚焦前逐漸增大，在聚焦時達到局部最大值后略微減小(t=200 μs)，隨后開始繼續增大(t=450 μs)，其增速隨時間而變大。當B=0.05 T時，平均渦量始終增大，且增速隨時間也逐漸增大。然而，在聚焦前，垂直磁場的平均渦量略大于水平磁場，但在聚焦后則剛好相反。此外，隨著磁感應強度增大，雖然能有效抑制不穩定性，但平均渦量明顯大于磁場較小的，甚至大于不加磁場的。由此可知，磁場并不能控制渦量的生成，但可以控制渦量的位置，并抑制界面的卷起。

表1 流場平均渦量Table 1 Average vorticity

為了更清楚地說明渦量變化情況，下面將采用渦量的平方進行分析，它能反映流團的旋轉能量，稱為渦度擬能[20]。表2為不同初始條件下不同時刻的平均渦度擬能。在不施加磁場(B=0)時，平均渦度擬能在t=200 μs內逐漸增大。當入射激波在尾部聚焦后，由于流場進入湍流轉捩階段，渦結構逐漸擴散，因此平均渦度擬能逐漸減小，并最終穩定在1.1×107s-2附近。當B=0.01 T時，平均渦度擬能明顯減小，同時垂直磁場比平行磁場對渦度擬能具有更好的抑制作用。由表2可知，雖然加入磁場后平均渦度擬能變化趨勢與無磁場情況類似，但隨著發展后期流場的演化，加劇界面處的擾動，使平均渦度擬能逐漸增加，因而在t=945 μs時，平均渦度擬能分別增至4.533×106和9.524×106s-2。當B=0.05 T時，平均渦度擬能受到更明顯的控制，在激波與氣柱作用過程的前中期，平行與垂直磁場情況下的平均渦度擬能均小于B=0.01 T時的。由此可見，平均渦度擬能清晰地反映磁場對不穩定性的控制。然而，在演化后期，雖然流場進入湍流轉捩階段，但單渦由于耗散作用逐漸減弱，使平均渦度擬能小于發展前中期，因此渦度擬能并不能反映流場的無序情況。

文獻[15]表明，當B=0.01 T時，磁能量與磁壓力較大區域分布在界面附近，且氣柱上下頂點處的磁壓力及磁能量隨著時間不斷增大。圖6～7分別為B=0.05 T時垂直與水行磁場作用下的磁壓力與磁能量分布曲線。對于垂直磁場，磁壓力和磁能量均增大，而當透射激波在水平軸線上發生聚焦后則稍有減小，隨后，磁壓力和磁能量隨著界面的演化繼續增大。同時，因磁場較小時渦層附著于界面，而磁場較大時渦層與界面分離并逐漸遠離，因此磁壓力及磁能量較大的區域位于渦層處而不是界面附近。由此可知，流場中渦量較大的位置磁壓力及磁能量較大，這也與前面的結果相符。

對于水平磁場，由于初始時刻磁場曲率較小，因此來流方向上的磁壓力較小。當氣柱界面處由于剪切力作用開始卷起后，垂直于來流方向的磁壓力對界面KH不穩定性的繼續發展起到了明顯的抑制作用。在這種情況下，磁壓力及磁能量隨時間不斷增大。對于縱向磁場，磁壓力在氣柱界面處快速脈動，上下界面相同橫坐標處磁壓力并不相同，且數值相差較大，與之相對，磁能量大小基本呈上下對稱分布。

3 結論

基于MHD方程及CTU+CT算法，對不同初始磁場下的激波沖擊重質氣柱過程進行了數值模擬，揭示了磁場抑制界面不穩定性的機理，得到以下結論。

(1) 垂直磁場及水平磁場對不穩定性均具有很好的抑制作用，其中水平磁場能很好地抑制界面的不穩定性，卻不能抑制尾部主渦的形成，而垂直磁場則能同時做到以上兩點。此外，隨著磁感應強度的增大，對不穩定性的抑制效果會更好。

(2) 與不加磁場情況相比，加入磁場后平均渦量并未明顯減小。當B=0.01 T時，平均渦量在透射激波聚焦前逐漸增大，在聚焦時達到局部最大值，隨后略微減小，之后繼續增大，且增速隨時間呈增大趨勢，可見B=0.01 T時，平均渦量略小于不加磁場情況。當B=0.05 T時，平均渦量始終增大，且增速隨時間也逐漸增大。此外，當磁感應強度增大時，雖然能有效抑制界面不穩定性，但平均渦量明顯大于磁場較小(B=0.01 T)的情況，甚至大于不加磁場的情形。因而磁場雖不能控制渦量的生成，但可以控制渦量的位置，并抑制界面的卷起。

(3) 平均渦度擬能可較好地反映出磁場對不穩定性的控制效果。當施加磁場后，平均渦度擬能明顯小于無磁場情況，且隨著磁場的增大而減小。同時，垂直磁場比平行磁場更能降低平均渦度擬能，這與不同方向初始磁場對不穩定性的控制效果相符。此外，由于渦的強度比渦的尺度對平均渦度擬能影響更大，因而平均渦量比平均渦度擬能可以更好地反映流場無序情況。

(4) 對于垂直磁場，磁場較小時的磁壓力與磁能量隨時間不斷增大，而磁場較大時的磁壓力與磁能量則在開始時較大，但在透射激波聚焦后出現減小，隨后繼續增大。對于水平磁場，由于初始時刻垂直于磁場方向的速度較小，此時磁壓力作用不明顯。隨著氣柱界面的卷起，磁壓力迅速增大并抑制不穩定性的發展，在此情況下，磁壓力及磁能量隨著時間始終不斷增大。此外，當磁場較小時，渦層附著于界面，渦層與界面隨著磁場增大而出現分離，并隨時間不斷遠離界面。

[1] BROUILLETTE M. The Richtmyer-Meshkov instability[J]. Annual Review of Fluid Mechanics, 2002,34:445-468.

[2] HAAS J F, STURTEVANT B. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[J]. Journal of Fluid Mechanics, 1987,181:41-76.

[3] TOMKINS C, KUMAR S, ORLICZ G, et al. An experimental investigation of mixing mechanisms in shock-accelerated flow[J]. Journal of Fluid Mechanics, 2008,611:131-150.

[4] 范美如,翟志剛,司廷,等.激波與不同形狀界面作用的數值模擬[J].中國科學:物理學、力學、天文學,2011,41(7):862-869.

FAN Meiru, ZHAI Zhigang, SI Ting, et al. Numerical simulation of interaction with different shape accelerated by a planar shock[J]. Scientia Sinica: Physica, Mechanica & Astronomica, 2011,41(7):862-869.

[5] FAN M R, ZHAI Z G, SI T, et al. Numerical study on the evolution of the shock-accelerated SF 6 interface: Influence of the interface shape[J]. Science China: Physics, Mechanics & Astronomy, 2012,55(2):284-296.

[6] 王顯圣,司廷,羅喜勝,等.反射激波沖擊重氣柱的RM不穩定性數值研究[J].力學學報,2012,44(4):664-672.

WANG Xiansheng, SI Ting, LUO Xisheng, et al. Numerical study on the RM instability of a heavy-gas cylinder interacted with reshock[J]. Chinese Journal of Theoretical and Applied Mechanics, 2012,44(4):664-672.

[7] SI T, ZHAI Z, YANG J, et al. Experimental investigation of reshocked spherical gas interfaces[J]. Physics of Fluids, 2012,24(5):054101.

[8] CHANDRASEKHAR S. Hydrodynamic and hydromagnetic stability[M]. Courier Dover Publications, 2013:441-453.

[9] WHEATLEY V, PULLIN D I, SAMTANEY R. Stability of an impulsively accelerated density interface in magnetohydrodynamics[J]. Physical Review Letters, 2005,95(12):125002.

[10] WHEATLEY V, SAMTANEY R, PULLIN D I. The magnetohydrodynamic Richtmyer-Meshkov instability: The transverse field case[C]∥The 18th Australasian Fluid Mechanics Conference. Australasian Fluid Mechanics Society, 2012.

[11] 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):042102.

[12] KHAN M, MANDAL L, BANERJEE R, et al. Development of Richtmyer-Meshkov and Rayleigh-Taylor instability in the presence of magnetic field[J]. Nuclear Instruments and Methods in Physics Research Section: Accelerators, Spectrometers, Detectors and Associated Equipment, 2011,653(1):2-6.

[13] SHIN M S, STONE J M, SNYDER G F. The magnetohydrodynamics of shock-cloud interaction in three dimensions[J]. The Astrophysical Journal, 2008,680(1):336-348.

[14] 李源,羅喜勝.黏性、表面張力和磁場對Rayleigh-Taylor不穩定性氣泡演化影響的理論分析[J].物理學報,2014,63(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,63(8):277-285.

[15] 林震亞,張煥好,陳志華,等.磁場對激波沖擊R22重氣柱作用過程影響的數值模擬[J].爆炸與沖擊,2017,37(4):748-758.

LIN Zhenya, ZHANG Huanhao, CHEN Zhihua, et al. Influence of magnetic field on interaction of shock wave with R22 heavy gas column[J]. Explosion and Shock Waves, 2017,37(4):748-758.

[16] SALTZMAN J. An unsplit 3D upwind method for hyperbolic conservation laws[J]. Journal of Computational Physics, 1994,115(1):153-168.

[17] GARDINER T A, STONE J M. An unsplit Godunov method for ideal MHD via constrained transport in three dimensions[J]. Journal of Computational Physics, 2008,227(8):4123-4141.

[18] 沙莎,陳志華,薛大文.激波沖擊 R22 重氣柱所導致的射流與混合研究[J].物理學報,2013,62(14):291-300.

SHA Sha, CHEN Zhihua, XUE Dawen. The generation of jet and mixing induced by the interaction of shock wave with R22 cylinder[J]. Acta Physica Sinica, 2013,62(14):291-300.

[19] HAAS J F, STURTEVANT B. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[J]. Journal of Fluid Mechanics, 1987,181:41-76.