吳 宇，施紅輝，王 超，葉 斌，張 珂

（浙江理工大學機械與自動控制學院，杭州310018）

液柱在激波沖擊下RM不穩(wěn)定性和破裂過程的數(shù)值計算

吳 宇，施紅輝，王 超，葉 斌，張 珂

（浙江理工大學機械與自動控制學院，杭州310018）

對液柱受到激波沖擊后在氣流中的變形斷裂過程進行了數(shù)值計算，研究了在該過程中Richtmyer-Meshkov（RM）不穩(wěn)定性的具體表現(xiàn)。應用Fluent軟件，數(shù)值模擬了二維（2D）和三維（3D）液柱在激波馬赫數(shù)為1.10、液柱初始直徑為2.76 mm的情況下，氣/液界面上RM不穩(wěn)定性的演化過程以及液體周圍的流場。計算結(jié)果表明，初始擾動數(shù)目對RM不穩(wěn)定性的影響顯著；液柱在橫截面平面（Z方向）發(fā)生變形失穩(wěn)，導致沿液柱軸向出現(xiàn)變形失穩(wěn)。數(shù)值計算的結(jié)果與已有的實驗結(jié)果吻合較好。

RM不穩(wěn)定性；激波；液柱；數(shù)值計算

0 引 言

不同形態(tài)的液體在激波沖擊下失穩(wěn)與破裂的過程和機理，不但是Richtmyer-Meshkov（RM）不穩(wěn)定性研究中的一個關鍵問題，也是一個可壓縮性湍流混合問題。液柱在激波以及波后氣流沖擊下的變形破裂過程中的具體表現(xiàn)與機理，對于進一步研究RM不穩(wěn)定性有相當重要的意義。2001年，Igra等［1］基于歐拉方程及CIP方法對激波與水柱的相互作用進行了數(shù)值計算，并對該方法作了一定程度的改進，使得該算法能夠正確地計算出氣液界面，而不致使界面處的密度階躍模糊化，采用該方法的數(shù)值計算結(jié)果與一些干涉圖片有較好的吻合；同時他們對水柱與固體圓柱進行了對比，結(jié)果發(fā)現(xiàn)，即使是在激波沖擊界面40μs后，仍然可觀察到液柱和固體圓柱一些差別。2004年，柏勁松等［2］采用可壓縮多介質(zhì)高精度PPM方法對實驗模型進行數(shù)值計算，對10模峰對峰、峰對谷初始振幅為1 mm的擾動模型在0、320、640、960μs時刻的果凍環(huán)進行了計算，結(jié)果發(fā)現(xiàn)：果凍內(nèi)外界面變形特征主要表現(xiàn)為Rayleigh-Taylor不穩(wěn)定性；在爆驅(qū)過程中，果凍內(nèi)界面上的擾動發(fā)生反相，并且初始內(nèi)外界面同相時，初始擾動對后續(xù)不穩(wěn)定性現(xiàn)象的影響更為嚴重。2008年，Chen等［3］提出了一種帶有5分成模型的可壓縮多相求解器，對激波與水柱之間的相互作用進行了研究，并計算了平面激波與水柱的相互作用，以及平面激波與雙排水柱的相互作用；研究了水柱分別為水平和豎直的兩種情況，發(fā)現(xiàn)改變放置方向會導致作用過程有差異；使用流動可視化技術展示了激波與水柱相互作用過程中的復雜結(jié)構(gòu)。王濤等［4］提出了基于多相流Volume Fraction（VOF）方法和逐段拋物線方法（piecewise parabolic method，PPM）的多相流逐段拋物線方法（multi-fluid piecewise parabolic method，MFPPM），并對高壓氣體沖擊作用下的氣體/液體交界面上的Richtmyer-Meshkov（RM）不穩(wěn)定性及其引起的流體混合現(xiàn)象進行了數(shù)值計算，計算結(jié)果表明：流體混合區(qū)和尖釘?shù)陌l(fā)展受初始擾動的影響極大，并且該影響在后期更為顯著；氣泡的寬度隨時間線性增長，也略微的受到初始擾動的影響；網(wǎng)格的尺寸不影響流體混合區(qū)以及尖釘和氣泡的變形演化過程。

2011年，施紅輝等［5］在大型水平激波管中，進行了液柱在激波及激波誘導的高速氣流沖擊下變形破碎的實驗研究。該實驗通過控制液柱初始直徑，以及使用不同馬赫數(shù)的激波，對液柱在沖擊后從壓縮變形到破碎霧化的過程進行了觀察，比較了不同實驗之間的差異；對液柱變形過程中的相關參數(shù)進行了測量分析，發(fā)現(xiàn)液柱迎風面出現(xiàn)RM不穩(wěn)定現(xiàn)象。但在實驗照片中，只能做到測量包括液柱初始直徑、尖釘氣泡大概尺寸、破碎霧化使用的大致時間等參數(shù)，不能得到具體的液柱變形過程的包括密度、速度、壓強等物理特性的具體數(shù)據(jù)。本文在文獻［5］液柱實驗的基礎上，根據(jù)該實驗數(shù)據(jù)，以數(shù)值計算為手段，模擬整個液柱實驗過程，利用計算軟件得到實驗過程中任意時刻液柱及其周圍流場的壓強、速度、密度等物理特性的數(shù)值，以進一步研究液柱實驗中的RM不穩(wěn)定性的影響。

1 二維和三維數(shù)學模型與計算方法

根據(jù)守恒關系，結(jié)合k-ε湍流模型和多相流VOF模型，液柱在激波沖擊下RM不穩(wěn)定性和破裂過程可由三維Navier-Stokes（NS）控制方程組表示。二維與三維控制方程組的區(qū)別在于三維方程組多了Z方向的量。三維的控制方程組可由下式表示：

式中φ為通用因變量，Γφ為輸運系數(shù)，Sφ為源項。對各方程而言，φ、Γφ、Sφ的具體含義見表1。

表1 三維控制方程的通用形式中各通用變量的含義

下面給出兩組2D數(shù)值計算的模型。在第一組數(shù)值計算的數(shù)學模型中，選取2 cm×8 cm的計算區(qū)域，計算網(wǎng)格數(shù)設為168 040個。液柱在激波及波后高速氣流沖擊下的變形是非定常過程，本文對該過程進行了模擬，采用壓力速度耦合的PISO算法，同時對密度、動量、能量方程中的對流項采用三階MUSCL離散格式，對壓力項采用Body Force Weighted離散格式，在VOF計算模型中，對體積分數(shù)則采用QUICK格式。具體計算模型如圖1所示。

圖1 2D計算區(qū)域示意

在2D數(shù)值計算中，如果柱體側(cè)面輪廓是直線，那么可認為圓柱體表面是光滑的，高速氣流會把液柱整體向右推動，液柱的外形保持穩(wěn)定不變。因此，在后續(xù)的計算中，采用了帶有初始擾動的邊界模型，設置波長λ=0.05 cm，振幅a=0.01 cm，在液柱側(cè)面一共有40個擾動波段，邊界的擾動情況圖1所示。

第二組模型同樣選取計算區(qū)域為2 cm×8 cm，由于擾動振幅加大，計算網(wǎng)格數(shù)設為172 336個，設置波長λ=0.4 cm，振幅a=0.05 cm，在液柱側(cè)面，共設置有5個擾動波段，計算參數(shù)設置與第一組相同。在第二組數(shù)值計算模型中，減少了液柱側(cè)面上的擾動數(shù)目，相應的延長了擾動波長，并增大了擾動振幅。

使用Gambit軟件進行網(wǎng)格劃分，具體的三維模型如圖2所示。

圖2 計算區(qū)域三維模型

圖2中，左側(cè)區(qū)域表示高壓區(qū)，圓柱體為液柱。因為要詳細的觀察液柱在激波沖擊作用后的變形演化，所以這里需要對實驗區(qū)域的網(wǎng)格進行細化。高壓氣體區(qū)域中的壓力、密度、溫度等物理特性數(shù)值設置高于右側(cè)的低壓氣體，以使得在計算開始時，在這一物理特性間斷面上，自然的形成向右傳播的平面激波。同時設置高壓氣體的初始速度向右，盡可能的再現(xiàn)激波以及波后氣流的運動狀況。圖2中，模型網(wǎng)格數(shù)量設為167 596個，計算時間步長設為5μs。數(shù)值模擬采用壓力速度耦合的PISO算法，同時對密度、動量、能量中的對流項采用二階迎風離散格式，對壓力項采用Body Force Weighted離散格式，在VOF計算模型中，對體積分數(shù)則采用幾何重構(gòu)（Geo-Reconstruct）格式。

2 結(jié)果分析與討論

當激波與波后氣流對液柱產(chǎn)生作用時，液柱表面會的出現(xiàn)比較復雜的變化。根據(jù)文獻［5］中的實驗結(jié)果，隨著時間的推移，可以觀察到液柱迎風面出現(xiàn)典型RM不穩(wěn)定性現(xiàn)象，出現(xiàn)尖釘和氣泡。在下面給出的算例中，激波馬赫數(shù)M=1.10、初始液柱直徑為2.76 mm。

第一組2D的計算結(jié)果如圖3所示。可以在圖3中觀測到，液柱的具體位移與文獻［5］中的實驗結(jié)果類似。從圖3中不能明顯看清楚微小液滴的剝離情況，但是可以看到液柱表面擾動隨時間不斷的增大，后期形成了RM不穩(wěn)定性中典型的氣泡和尖釘結(jié)構(gòu)。但是，本文的數(shù)值計算結(jié)果與文獻［5］中的實驗結(jié)果有較大的不同。首先，從圖3（4）開始可以看到液柱上下兩端邊角卷起形成了觸角狀。其次，從圖3（8）開始，液柱靠近上下兩端的位置出現(xiàn)中空結(jié)構(gòu)，然而在實驗中并沒有發(fā)現(xiàn)類似現(xiàn)象。液柱的整體呈弓狀結(jié)構(gòu)，中間略窄，上下略粗。液柱迎風面出現(xiàn)擾動增大現(xiàn)象，背風面沒有這一情況。實驗中迎風面也發(fā)現(xiàn)了尖釘氣泡結(jié)構(gòu)，背風面由于液滴剝離形成了尾跡，在數(shù)值模擬結(jié)果中沒有發(fā)現(xiàn)這一情況。最后，在圖3（19）之后，尖釘狀結(jié)構(gòu)不再維持，開始上下擺動、融合、斷裂。在實驗中沒有發(fā)現(xiàn)如此長度的尖釘結(jié)構(gòu)，當尖釘和氣泡結(jié)構(gòu)發(fā)展到一定程度之后，它們會維持一定尺寸，最后隨著整個液柱的裂解霧化一起崩碎。

圖3 第一組2D液柱實驗密度云圖

圖4 第二組2D液柱實驗密度云圖

第二組2D計算的結(jié)果如圖4所示。對比實驗結(jié)果，第一組數(shù)值計算結(jié)果整體液柱的變形速度與實驗接近，尖釘?shù)冉Y(jié)構(gòu)與實驗結(jié)果更為吻合。第二組數(shù)值計算結(jié)果中，由于擾動數(shù)目降低，并且更為規(guī)則，因此引發(fā)的尖釘氣泡結(jié)構(gòu)較為規(guī)則，同時整體變形速度快于實驗。兩組數(shù)值計算結(jié)果的對比說明初始擾動的增加與擾動振幅的降低，會減緩整個液柱變形的速度。在實際過程中，尖釘結(jié)構(gòu)不會拉伸過長，而是維持一定尺寸；在第二組數(shù)值計算結(jié)果中，出現(xiàn)尖釘結(jié)構(gòu)拉伸過長的情況，主要原因是在計算中沒有考慮Kelvin-Helmholtz（KH）不穩(wěn)定性和湍流混合的影響的緣故。

圖5為三維液柱的數(shù)值的計算結(jié)果，左側(cè)為垂直沿軸線方向中央截面密度云圖，右側(cè)為水平沿軸線方向中央截面密度云圖。從圖5（a）圖中可以看到，液柱的模擬結(jié)果中，液柱橫向移動速度不斷變大，其表面有微小液滴剝落，最后發(fā)生液柱斷裂現(xiàn)象。圖5（b）表明，液柱橫截面形狀是沿激波管徑向的直徑增大，沿激波管軸向的直徑降低，呈扁平化發(fā)展，微小液滴從液柱表面不斷剝離，最后發(fā)生崩潰、破碎。對比文獻［5］，整個液柱的失穩(wěn)及破裂過程的模擬結(jié)果與實驗結(jié)果吻合較好。特別地，數(shù)值計算提供了液柱核心破裂的過程：從圖5b（6）-b（7）中可以看到，液柱從橫截面方向上出現(xiàn)破裂，并且軸向也出現(xiàn)斷裂；文獻［5］中的實驗照片中，只能觀察到重疊在一起的影像，而不能分辨出液柱變形的細節(jié)，無法判斷液柱是否斷裂。

圖5 三維液柱實驗數(shù)值計算結(jié)果

圖6是液柱中部在（X，Z）平面里，即在橫向截面上的靜壓云圖。其中，圖6（1）-（3）中，液柱尾部上下各形成了一對稱的低壓區(qū)域。該區(qū)域隨時間增大，有一定的對稱性。對比密度圖，液柱Z方向直徑逐漸增大，其前方形成一塊較規(guī)則的高壓區(qū)域。圖6（5）開始，液柱后方的壓力分布開始出現(xiàn)明顯的不均勻。圖6（6）-（7）中，液柱從軸向和橫截面方向上開始破碎斷裂，流場變得非常紊亂。

圖6 液柱橫向截面上的靜壓云圖變化過程

圖7是模型部（X，Z）平面上的速度云圖。圖7（1）-（4）中，中間面積較大的一塊零速度區(qū)域，即深藍色區(qū)域，其形狀與密度云圖中的液柱橫向截面的形狀完全吻合。液柱左右兩側(cè)的高速區(qū)域不斷擴大，并且速度不斷上升，最高速度達268 m/s，此時，液柱表面存在氣體與剝離的微小水滴的混合流。同時液柱橫向截面尾部后方形成了低速區(qū)域擴張而形成的火焰狀尾翼。圖7（5）開始，液柱橫向截面的火焰狀尾翼擺動明顯，速度場越來越不穩(wěn)定。從圖7（6）-（7）可以看出，液柱周邊的速度場已經(jīng)極度紊亂，此時，可以看到液柱在橫截面方向上已經(jīng)破裂。

圖7 液柱橫向截面上的速度云圖變化過程

利用三維計算結(jié)果，可以從橫向和縱向兩個方向來觀察液柱變形斷裂過程中的細節(jié)?？梢钥闯觯褐紫仍跈M截面方向上發(fā)生破裂，隨后在軸向上也發(fā)生斷裂。對比實驗結(jié)果，在三維數(shù)值計算結(jié)果中，液柱橫向位移發(fā)生的更快，整體變形破裂的速度也要更快；并且在液柱后方?jīng)]有出現(xiàn)實驗中的大規(guī)模微小液滴的剝離、霧化。產(chǎn)生這種現(xiàn)象的原因可能是：第一，計算模型中參數(shù)的設定不能完全符合實驗環(huán)境，以及Fluent軟件本身存在局限性，導致計算結(jié)果與實驗過程存在偏差；第二，液滴剝離時與空氣混合，其密度接近空氣，因此在計算結(jié)果圖中難以分辨。

3 結(jié) 論

本文利用非定常的二維和三維Navier-Stokes（NS）方程、標準k-ε湍流模型和多相流VOF模型，對液柱在激波誘導的高速氣流中的變形斷裂現(xiàn)象進行了數(shù)值計算與分析，得到如下結(jié)論。

a）二維（2D）液柱數(shù)值計算模擬中，由于加入了初始擾動，出現(xiàn)了尖釘與氣泡狀結(jié)構(gòu)，即出現(xiàn)了RM不穩(wěn)定性。在第一組數(shù)值計算中，由于擾動數(shù)目較多，擾動的振幅較小，得到的數(shù)值計算結(jié)果中，液柱整體變形速度與實驗結(jié)果接近。

b）三維（3D）液柱數(shù)值計算中，雖然液柱表面并沒有初始擾動，但由于在垂直于液柱方向上，液柱存在初始邊界形狀，因此在高速氣流的作用下仍然會出現(xiàn)RM不穩(wěn)定性的特征。在時長1.75 ms的數(shù)值模擬過程中，液柱完成了變形斷裂過程，但在實驗中，這一過程歷時4 ms左右，因此數(shù)值計算結(jié)果中液柱的變形速度要快過實驗。從具體的演變過程可以得出結(jié)論，由于液柱在橫截面平面（Z方向）發(fā)生變形失穩(wěn)，導致沿液柱軸向隨之出現(xiàn)變形失穩(wěn)。

c）二維計算接近實驗結(jié)果，但三維計算更符合實驗結(jié)果，在進行數(shù)值研究時，更應采用三維計算方法。在本文二維（2D）數(shù)值計算過程中，存在Fluent軟件系統(tǒng)將液柱辨識為“液墻”的問題；三維（3D）計算存在RM不穩(wěn)定性現(xiàn)象不明顯，即沒有出現(xiàn)文獻［5］中實驗照片上的尖釘氣泡狀結(jié)構(gòu)，以及數(shù)值計算中液柱變形斷裂速度存在過快的問題。

在后續(xù)研究中，為了提高數(shù)值計算的可靠性，需要進一步對計算模型進行優(yōu)化設計，并對Fluent軟件在該算例中的相關參數(shù)設置做進一步的優(yōu)化調(diào)整，以得到更加真實可信的數(shù)值計算結(jié)果。

［1］Igra D，Takayama A.Investigation of aerodynamic breakup of a cylindrical water droplet［J］.Atomization and Sprays，2001，11（2）：167-185.

［2］柏勁松，李 平，陳森華，等.內(nèi)爆加載下果凍內(nèi)外界面不穩(wěn)定性數(shù)值計算［J］.高壓物理學報，2004，18（4）：295-301.

［3］Chen H，Liang S M.Flow visualization of shock/water column interactions［J］.Shock Waves，2008，17（5）：309-321.

［4］王 濤，柏勁松，李 平.二維氣/液界面不穩(wěn)定性數(shù)值模擬［J］.高壓物理學報，2008，22（3）：298-304.

［5］施紅輝，吳 宇，肖 毅，等.激波與液柱相互作用時的氣動特性研究［C］//中國力學學會流體力學專業(yè)委員會實驗流體力學專業(yè)組.第9屆全國實驗流體力學學術會議論文集，2013：36-42.

NumericaI CaIcuIation of RM InstabiIity and Rupture Process of Liquid CoIumn under Impact of Shock Waves

WU Yu，SHI Hong-hui，WANG Chao，YE Bin，ZHANG Ke
（School of Mechanical Engineering&Automation，Zhejiang Sci-Tech University，Hangzhou 310018，China）

This paper calculated the deformation and rupture process of the liquid column under the impact of shock waves and studied specific representations of Richtmyer-Meshkov（RM）instability in the process.Fluent software is used to carry out numerical simulation of evolutionary process of RM instability on gas/liquid interface and flow field around the gas of 2D and 3D liquid columns when shock mach number is 1.10 and initial diameter of shock mach number is 2.76 mm.The results show that the number of initial perturbation imposes significant influences on RM instability；the liquid column deforms and loses stability along the cross section（direction Z），so deformation and instability appear along axial direction of the liquid column.The numerical calculation results well coincide with experimental results.

RM instability；shock wave；liquid column；numerical calculation

O354.5；O359

A

（責任編輯：康 鋒）

1673-3851（2014）05-0502-05

2014-03-03

吳 宇（1983-），男，江西玉山人，碩士研究生，研究方向為湍流與復雜流動以及可壓縮性與瞬態(tài)流動。

施紅輝，E-mail：hhshi@zstu.edu.cn