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

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

（浙江理工大學(xué)機(jī)械與自動(dòng)控制學(xué)院，杭州310018）

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

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

（浙江理工大學(xué)機(jī)械與自動(dòng)控制學(xué)院，杭州310018）

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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖3 第一組2D液柱實(shí)驗(yàn)密度云圖

圖4 第二組2D液柱實(shí)驗(yàn)密度云圖

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

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

圖5 三維液柱實(shí)驗(yàn)數(shù)值計(jì)算結(jié)果

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

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

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

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

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

3 結(jié) 論

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

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

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

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

在后續(xù)研究中，為了提高數(shù)值計(jì)算的可靠性，需要進(jìn)一步對(duì)計(jì)算模型進(jìn)行優(yōu)化設(shè)計(jì)，并對(duì)Fluent軟件在該算例中的相關(guān)參數(shù)設(shè)置做進(jìn)一步的優(yōu)化調(diào)整，以得到更加真實(shí)可信的數(shù)值計(jì)算結(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ì)算［J］.高壓物理學(xué)報(bào)，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］.高壓物理學(xué)報(bào)，2008，22（3）：298-304.

［5］施紅輝，吳 宇，肖 毅，等.激波與液柱相互作用時(shí)的氣動(dòng)特性研究［C］//中國(guó)力學(xué)學(xué)會(huì)流體力學(xué)專業(yè)委員會(huì)實(shí)驗(yàn)流體力學(xué)專業(yè)組.第9屆全國(guó)實(shí)驗(yàn)流體力學(xué)學(xué)術(shù)會(huì)議論文集，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

（責(zé)任編輯：康 鋒）

1673-3851（2014）05-0502-05

2014-03-03

吳 宇（1983-），男，江西玉山人，碩士研究生，研究方向?yàn)橥牧髋c復(fù)雜流動(dòng)以及可壓縮性與瞬態(tài)流動(dòng)。

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