激波作用不同橢圓氦氣柱過(guò)程中流動(dòng)混合研究

李冬冬 王革 張斌

1)(哈爾濱工程大學(xué)航天與建筑工程學(xué)院,哈爾濱 150001)

2)(上海交通大學(xué)航空航天學(xué)院,上海 201100)

1 引 言

當(dāng)激波作用于不同密度的流體界面時(shí)會(huì)發(fā)生折射和反射等現(xiàn)象,流體界面也會(huì)獲得一定的加速,從而導(dǎo)致流場(chǎng)中產(chǎn)生一系列復(fù)雜流動(dòng)現(xiàn)象,這種現(xiàn)象被稱(chēng)為Richtmyer-Meshkov不穩(wěn)定性(Richtmyer-Meshkov instability,RMI)[1,2].許多自然現(xiàn)象和工程問(wèn)題中存在這種現(xiàn)象,例如武器內(nèi)爆、慣性約束核聚變[3,4]、超音速燃燒[5]、超新星爆發(fā)[6]等.對(duì)RMI的研究不僅會(huì)推動(dòng)渦動(dòng)力學(xué)、多相流、湍流等流體力學(xué)基礎(chǔ)難題的研究,而且在工程領(lǐng)域和自然界中有直接應(yīng)用價(jià)值.

近幾十年來(lái),國(guó)內(nèi)外學(xué)者對(duì)激波與不同密度的氣體界面相互作用問(wèn)題中的RMI進(jìn)行了大量的數(shù)值、實(shí)驗(yàn)和理論研究.Haas和Sturtevant[7]采用實(shí)驗(yàn)手段對(duì)弱激波與R22重氣柱、氣泡以及He輕氣柱、氣泡的作用過(guò)程進(jìn)行了研究,對(duì)界面變形和運(yùn)動(dòng)特征進(jìn)行了討論和分析.Jacobs[8]首次采用平面激光誘導(dǎo)熒光(planar laser-induced fluorescence,PLIF)技術(shù)研究了低馬赫數(shù)激波與氦氣柱相互作用的混合問(wèn)題,采用激波壓縮后氣柱面積的減少量來(lái)表征混合.Giordano和Burtschell[9]研究了較低馬赫數(shù)下激波與氣柱及氣泡的相互作用,數(shù)值模擬結(jié)果和一維氣體動(dòng)力學(xué)理論分析表明,不同密度比的氣泡受相同激波壓縮后,體積壓縮率近似是定值,其值與氣泡和周?chē)鷼怏w的比熱之比有關(guān).Ranjan等[10]通過(guò)數(shù)值模擬方法驗(yàn)證了Giordano和Burtschell[9]有關(guān)壓縮率的結(jié)論:高馬赫數(shù)激波與輕氣泡作用時(shí),Giordano等的結(jié)論失效,且隨著馬赫數(shù)以及密度比的增大,體積壓縮率變化存在更大的振蕩.Tomkins等[11]采用PLIF技術(shù)對(duì)馬赫數(shù)1.2的激波與SF6重氣柱的相互作用問(wèn)題進(jìn)行了實(shí)驗(yàn)研究,并以物質(zhì)質(zhì)量分?jǐn)?shù)為基礎(chǔ)的瞬態(tài)擴(kuò)散率和總混合率表征混合效果,對(duì)SF6與空氣的混合機(jī)理進(jìn)行了深入的研究和討論.Shankar等[12]采用局部人工擴(kuò)散方法對(duì)Tomkins等[11]的實(shí)驗(yàn)進(jìn)行了數(shù)值模擬,以總混合率為參考標(biāo)準(zhǔn),研究了SF6初始分布和示蹤物質(zhì)丙酮的含量對(duì)結(jié)果的影響.國(guó)內(nèi)沙莎等[13,14]采用五階加權(quán)基本無(wú)振蕩格式與沉浸邊界法并結(jié)合大渦模擬技術(shù)的數(shù)值方法,對(duì)激波作用R22重氣柱和SF6球形氣泡所致的射流混合進(jìn)行了研究,分析了入射激波以及反射激波在氣泡界面聚焦誘導(dǎo)射流的過(guò)程,詳細(xì)研究了不同反射距離下反射激波與重氣體的作用過(guò)程及流場(chǎng)結(jié)構(gòu),此外對(duì)入射激波與梯形SF6重氣柱相互作用過(guò)程中的復(fù)雜波系演化也進(jìn)行了詳細(xì)的分析[15].Bai等[16]以及廖深飛等[17]采用粒子圖像測(cè)速技術(shù)對(duì)激波誘導(dǎo)單橢圓、雙橢圓SF6重氣柱的不穩(wěn)定性問(wèn)題進(jìn)行了實(shí)驗(yàn)研究,獲得了激波沖擊下的速度場(chǎng)、渦量場(chǎng)和環(huán)量等,定量地表征了激波界面相互作用的RMI現(xiàn)象,并探討了橢圓結(jié)構(gòu)對(duì)界面演化的影響.Zhai等[18]采用VAS2D程序?qū)げㄗ饔肒r/SF6氣泡過(guò)程中射流的形成進(jìn)行了深入研究,結(jié)果表明,壓力擾動(dòng)和斜壓-渦量的沉積是激波與氣柱相互作用(shock bubble interaction,SBI)問(wèn)題中射流形成的兩個(gè)主要因素.隨后采用實(shí)驗(yàn)方法研究了二維V形空氣/SF6界面在入射激波和反射激波作用下的RMI發(fā)展規(guī)律[19].Fan等[20]同樣借助VAS2D程序?qū)げㄗ饔糜陂L(zhǎng)方形、橢圓、三角形等5種不同形狀的SF6氣柱進(jìn)行了數(shù)值模擬,分析了該過(guò)程中氣柱形態(tài)、環(huán)量等的變化.Wang等[21]采用高速攝影技術(shù)和VAS2D數(shù)值模擬方法對(duì)激波作用不同形態(tài)的氦氣柱進(jìn)行了研究,對(duì)界面形態(tài)特征、環(huán)量等進(jìn)行了比較和分析.黃熙龍等[22]采用PLIF技術(shù)對(duì)激波作用于橢圓氣柱過(guò)程中的RMI問(wèn)題進(jìn)行了研究,分析了該過(guò)程中界面氣體聚集、轉(zhuǎn)移、消散等現(xiàn)象.

激波與氣柱以及液滴的相互作用可加速其與外部氣體的混合,從而提高特定條件下的燃燒性能,相關(guān)研究對(duì)提高超燃發(fā)動(dòng)機(jī)的性能具有極其重要的作用,因而對(duì)激波加速氣柱與外部氣體混合機(jī)理的研究具有重要意義.本文基于雙通量模型[23],結(jié)合五階WENO方法求解二維多組分Navier-Stokes(N-S)方程,對(duì)激波與不同結(jié)構(gòu)的橢圓氦氣柱作用過(guò)程進(jìn)行數(shù)值模擬,研究了橢圓氦氣柱與周?chē)橘|(zhì)之間的混合情況及幾何構(gòu)型對(duì)流動(dòng)和混合的影響,以期為進(jìn)一步研究混合燃燒效率提供一定的理論基礎(chǔ).

2 數(shù)值方法與計(jì)算模型

2.1 數(shù)值方法

考慮可壓縮多組分的二維N-S方程組,其形式為

其中U為守恒參數(shù)組成的向量,F和G為對(duì)流通量向量,Fv和Gv為黏性及擴(kuò)散通量向量,具體表達(dá)式為

這里ρ,p和E分別代表混合物密度、壓力以及單位質(zhì)量的總能量;u和v是混合物速度矢量在x和y方向上的分量;Yi是組分i的質(zhì)量分?jǐn)?shù);Ng為物質(zhì)種類(lèi)總數(shù);Jx,i和Jy,i分別為組分i的擴(kuò)散在x和y方向的分量;qx和qy分別為熱擴(kuò)散量在x和y方向的分量;σ為偏應(yīng)力張量.多組分流動(dòng)的理想氣體狀態(tài)方程為

其中Mwi為物質(zhì)i的摩爾分子量,Ru為通用氣體常數(shù),T為混合物溫度.比總能量的表達(dá)式為

其中CPi為物質(zhì)i的定壓比熱,采用美國(guó)國(guó)家航空航天局(National Aeronautics and Space Administration,NASA)提供的擬合多項(xiàng)式進(jìn)行計(jì)算;h0fi為物質(zhì)i的標(biāo)準(zhǔn)焓;H為物質(zhì)的總焓;T0為參考溫度.混合物的定壓比熱CP采用各物質(zhì)的質(zhì)量分?jǐn)?shù)Yi和定壓比熱CPi加權(quán)平均得到.音速c是氣相混合物的一個(gè)重要屬性,定義為

對(duì)于理想氣體音速c可以簡(jiǎn)化為

其中γ為混合物的比熱,定義為

物質(zhì)i的質(zhì)量擴(kuò)散通量Ji為

其中是物質(zhì)i的擴(kuò)散速度,其計(jì)算公式為

其中Xi為物質(zhì)的摩爾分?jǐn)?shù),Di,mix為物質(zhì)i在混合物中的混合平均擴(kuò)散系數(shù).

忽略Dufour效應(yīng),熱擴(kuò)散矢量q計(jì)算表達(dá)式如下:

其中hi為物質(zhì)i的焓值,同樣采用由NASA提供的多項(xiàng)式計(jì)算;λ為混合物的導(dǎo)熱系數(shù).

流體間正應(yīng)力和剪應(yīng)力分別為

其中μ為混合物的剪切黏性,κ為體積黏性.

本文采用簡(jiǎn)化的混合方法來(lái)計(jì)算混合物的輸運(yùn)性質(zhì)參數(shù) (ξ=μ,λ,κ)[24]:

對(duì)于物質(zhì)剪切黏性μ、體積黏性κ和熱傳導(dǎo)系數(shù)λ,n分別為6,4/3和4.

物質(zhì)的混合平均擴(kuò)散系數(shù)Di,mix采用下式計(jì)算:

其中Dij為物質(zhì)i在物質(zhì)j中的二元擴(kuò)散系數(shù),Xj為物質(zhì)j的摩爾分?jǐn)?shù).

單物質(zhì)的剪切黏性、熱傳導(dǎo)系數(shù)、二元擴(kuò)散系數(shù)基于Chapman-Enskog理論[24,25]和Lennard-Jones參數(shù)[26]計(jì)算獲得,純物質(zhì)的體積黏性κ依據(jù)Ern和Giovangigli[24]的簡(jiǎn)化方法獲得.

考慮到控制方程(1)各部分的物理性質(zhì),采用時(shí)間分裂的方法,將控制方程分裂成對(duì)流項(xiàng)和黏性項(xiàng)依次求解,分裂方式如下:

對(duì)流項(xiàng)的求解在空間上采用五階WENO格式并結(jié)合雙通量模型消減多物質(zhì)比熱不同引起的界面參數(shù)振蕩的影響,時(shí)間上采用三階Runge-Kutta格式[27];黏性項(xiàng)的求解中空間項(xiàng)采用四階中心差分,時(shí)間上采用二階Runge-Kutta-Chebyshev方法[28]進(jìn)行積分完成.

以一維歐拉方程為例對(duì)對(duì)流項(xiàng)的求解方法(五階WENO格式結(jié)合雙通量模型)進(jìn)行必要的描述.一維歐拉方程形式如下:

其半離散格式為

其中Ui為網(wǎng)格節(jié)點(diǎn)上的守恒變量,Fi±1/2是網(wǎng)格左右界面上的對(duì)流通量.定義離散算子~(U):

方程(15)可以寫(xiě)成如下形式:

文中采用了Abgrall和 Karni[23]提出的雙通量模型消除由于變比熱帶來(lái)的壓力和速度振蕩.在雙通量模型中,網(wǎng)格界面上的對(duì)流通量計(jì)算兩次,分別標(biāo)記為如圖1所示,其中采用五階WENO格式獲得.當(dāng)所有網(wǎng)格界面的兩次通量計(jì)算完成后,如圖2所示,空間算子~(U)可以改寫(xiě)為

圖1 雙通量模型在高階格式中的應(yīng)用Fig.1.Application of double flux model in high order scheme.

圖2 網(wǎng)格界面對(duì)流通量Fig.2.Convective flux at cell face.

2.2 計(jì)算模型

文中采用如圖3所示的計(jì)算模型,計(jì)算條件參照Haas和Sturtevant[7]的實(shí)驗(yàn)條件給定,氣柱內(nèi)輕氣體為氦氣,周?chē)橘|(zhì)為空氣,誘導(dǎo)激波馬赫數(shù)為1.22,位于氣泡右側(cè),波前氣體靜止,溫度和壓力分別為293 K和1 atm(1 atm=1.01325×105Pa),波后氣體參數(shù)由Rankine-Hugoniot條件獲得.由于對(duì)稱(chēng)性,選取模型的上半部分進(jìn)行計(jì)算,上邊界采用固體反射邊界,對(duì)稱(chēng)軸采用對(duì)稱(chēng)邊界條件,左右邊界上采用無(wú)反射邊界條件,各參數(shù)梯度為0.

以圖3所示的計(jì)算模型為基礎(chǔ),并在網(wǎng)格無(wú)關(guān)性檢驗(yàn)和算例驗(yàn)證可靠性的前提下,針對(duì)面積相同、幾何構(gòu)型不同(圓形,兩個(gè)激波沿橢圓長(zhǎng)軸作用于氣柱,兩個(gè)激波沿橢圓短軸作用于氣柱)的氦氣柱與激波的相互作用過(guò)程,分析界面變形、波系演化、界面結(jié)構(gòu)參數(shù)(界面高度和長(zhǎng)度)、氣柱體積壓縮率、總混合率、環(huán)量等的變化,探究氣泡內(nèi)介質(zhì)和周?chē)h(huán)境氣體的混合機(jī)理和不同幾何構(gòu)型下氣體混合的優(yōu)劣.以半徑2.5 cm的圓形截面氣柱面積19.63 cm2為標(biāo)準(zhǔn)來(lái)設(shè)計(jì)氣柱截面形態(tài),如圖4所示.為了方便表示,采用統(tǒng)一的標(biāo)準(zhǔn)橢圓方程來(lái)約束界面形狀:

圖3 SBI計(jì)算模型示意圖Fig.3.A schematic of computational domain for SBI.

圖4 氣柱幾何構(gòu)型 (a)圓形;(b)激波沿橢圓長(zhǎng)軸作用于氣柱;(c)激波沿橢圓短軸作用于氣柱Fig.4.Helium cylinders geometry:(a)Circle;(b)shock hitting along major axis of ellipse;(c)shock hitting along minor axis of ellipse.

不同幾何構(gòu)型的橢圓氣柱對(duì)應(yīng)的參數(shù)如表1所列,其中Geometry 1和Geometry 2中激波沿橢圓長(zhǎng)軸作用氣柱,Geometry 3為圓形氣柱,Geometry 4和Geometry 5中激波沿橢圓短軸作用氣柱.

表1 氣柱邊界方程參數(shù)Table 1.Parameters of gas cylinder boundary equation.

為防止虛假渦量的產(chǎn)生,采用有限厚度擴(kuò)散層[29]獲得更加光滑的初始?xì)馀萁缑?(20)式為氦氣質(zhì)量分?jǐn)?shù)初始分布函數(shù),

其中fin,fout分別表示氣柱內(nèi)外氦氣質(zhì)量分?jǐn)?shù);d為橢圓上一點(diǎn)到橢圓中心的距離;r為計(jì)算域內(nèi)任意一點(diǎn)到橢圓中心的距離;r0為橢圓中心位置;Cr為常數(shù),值越大界面越光滑,所有算例中Cr=10000.在笛卡爾坐標(biāo)系中,假定橢圓長(zhǎng)軸或者短軸位于計(jì)算域的對(duì)稱(chēng)軸上,幾何中心位于坐標(biāo)原點(diǎn),計(jì)算域內(nèi)有一點(diǎn)坐標(biāo)為(x0,y0),則r和d的表達(dá)式為

3 結(jié)果與討論

3.1 網(wǎng)格無(wú)關(guān)性檢驗(yàn)與算例驗(yàn)證

為了考察計(jì)算方法的可靠性,選擇與參考文獻(xiàn)[7]中激波作用氦氣柱的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比.首先,以Geometry 3圓形氣柱為例,采用三種不同尺寸的網(wǎng)格進(jìn)行網(wǎng)格無(wú)關(guān)性檢驗(yàn),網(wǎng)格數(shù)量分別為粗網(wǎng)格(600×83,Grid 1)、中等網(wǎng)格(1200×166,Grid 2)、細(xì)網(wǎng)格(2400×330,Grid 3),相應(yīng)的氣柱內(nèi)網(wǎng)格數(shù)量分別為80,160和320.

圖5是以氦氣質(zhì)量分?jǐn)?shù)顯示的不同網(wǎng)格尺度下的界面變形情況.結(jié)果表明在三種網(wǎng)格尺度下界面變形過(guò)程基本相似,較細(xì)的網(wǎng)格體現(xiàn)出更多的界面細(xì)節(jié)特征.圖6為不同網(wǎng)格尺度下界面三個(gè)特征點(diǎn)(特征點(diǎn)位置如圖6中右下角所示)在激波作用He氣泡過(guò)程中的位置變化(界面以氦氣質(zhì)量分?jǐn)?shù)0.01作為標(biāo)準(zhǔn)),三種不同網(wǎng)格下所獲得的結(jié)果基本相同.在后續(xù)的計(jì)算分析中為了獲得較為清晰的流動(dòng)細(xì)節(jié)均采用細(xì)網(wǎng)格(2400×330).

圖5 網(wǎng)格無(wú)關(guān)性檢驗(yàn)(左,600×83;中,1200×166;右,2400×330) (a)62μs;(b)240μs;(c)427μs;(d)674μsFig.5. Grid re finement test(left,600×83;middle,1200×166;right,2400×330):(a)62μs;(b)240μs;(c)427μs;(d)674μs.

圖6 不同網(wǎng)格尺寸下特征點(diǎn)位置的變化Fig.6.Movement of characteristic points for different grids.

圖7給出了激波作用于氦氣泡前期,特征點(diǎn)位置的變化及與其他計(jì)算結(jié)果[30,31]的對(duì)比,圖中曲線(xiàn)數(shù)據(jù)符合較好,表明本文所采用的數(shù)值方法具有較高的精度,可以用于激波作用氦氣柱的研究中.

為進(jìn)一步檢驗(yàn)方法的準(zhǔn)確性,將在細(xì)網(wǎng)格下獲得的數(shù)值紋影結(jié)果與文獻(xiàn)[7]實(shí)驗(yàn)采集到的紋影結(jié)果進(jìn)行定性的對(duì)比,如圖8所示,數(shù)值紋影圖清晰地表明了作用過(guò)程中界面變形和波系演化過(guò)程,進(jìn)一步說(shuō)明所采用數(shù)值方法的可靠性和精度.

圖7 激波作用于氦氣泡前期,特征點(diǎn)位置的變化及與其他計(jì)算結(jié)果[30,31]的比較Fig.7.Space-time diagram for the interaction of a shock wave with a helium bubble at initial stage;comparisons with the results of Ref.[30,31].

圖8 數(shù)值紋影與實(shí)驗(yàn)結(jié)果[7]的比較 (a)62μs;(b)102μs;(c)467μs;(d)674μsFig.8.Comparison between numerical schlieren and experimental results[7]:(a)62 μs;(b)102 μs;(c)467 μs;(d)674μs.

3.2 不同幾何構(gòu)型下氦氣混合過(guò)程分析

3.2.1 界面形態(tài)和特征

1)界面變形和波系演化

激波作用于不同幾何構(gòu)型的橢圓氦氣柱中界面形態(tài)及波系演化分別如圖9—圖13所示.在所有的算例中,激波由右向左作用于氣柱,將激波接觸界面最右端的時(shí)刻定義為0μs時(shí)刻.

圖9為激波沿橢圓長(zhǎng)軸作用于半長(zhǎng)軸a=4 cm、半短軸b=1.5625 cm的氦氣柱過(guò)程中界面和波系演化情況.當(dāng)誘導(dǎo)激波與氣柱最右端接觸并沿橢圓氣柱邊界向下游運(yùn)動(dòng)的過(guò)程中,在界面處產(chǎn)生了不規(guī)則反射現(xiàn)象,如圖9(a)和圖14(a)所示.由于氣柱內(nèi)氦氣的聲阻抗較空氣低,進(jìn)入氣柱內(nèi)的透射激波(transmitted shock wave)傳播速度要快于在空氣中傳播的誘導(dǎo)激波(incident shock wave),透射激波在界面處發(fā)生折射,在空氣中形成一道新的激波(free-precursor shock wave),與誘導(dǎo)激波相遇后形成三叉激波結(jié)構(gòu).當(dāng)透射激波運(yùn)動(dòng)到下游界面后,在空氣中形成一道向下游傳播的二次透射激波,同時(shí)在氣柱內(nèi)形成一道向上游界面?zhèn)鞑サ姆瓷浼げ?如圖9(b)所示.其后波系經(jīng)過(guò)復(fù)雜的反射、折射及相互作用,變得極其復(fù)雜,從圖9(c)和圖9(d)可以清晰地觀(guān)察到這一現(xiàn)象.當(dāng)激波接觸界面后,界面頂端會(huì)形成一個(gè)小的空氣射流(air jet)(圖9(b)),隨著流動(dòng)的發(fā)展,空氣射流不斷增長(zhǎng)并向下游滲透,該過(guò)程中伴隨著由RMI所致的主渦及小的速度剪切引起的次級(jí)渦的出現(xiàn)和發(fā)展(圖9(c)—(f)),初始的小空氣射流最終導(dǎo)致流動(dòng)分離成兩個(gè)獨(dú)立發(fā)展的渦團(tuán)(圖9(g)和圖9(h)).

圖9 幾何構(gòu)型1數(shù)值紋影圖(激波沿橢圓長(zhǎng)軸作用氣柱;a=4.0 cm,b=1.5625 cm)Fig.9.Sequence of numerical schlieren images of shock-accelerated elliptic helium cylinder for Geometry 1(shock hitting along major axis of ellipse;a=4.0 cm,b=1.5625 cm).

圖10和圖11分別為激波沿橢圓長(zhǎng)軸作用于半長(zhǎng)軸a=3.25 cm、半短軸b=1.9231 cm的氦氣柱與半徑為2.5 cm的圓形氦氣柱過(guò)程中界面和波系演化情況.其中波系演化模式及界面變形情況與幾何構(gòu)型1基本相似,但不難發(fā)現(xiàn)由于長(zhǎng)軸的變短和短軸長(zhǎng)度的增加(離心率變小),激波入射角度、壓力梯度和密度梯度之間的方向發(fā)生了改變,一方面導(dǎo)致了界面前端空氣射流形成時(shí)間增加,射流初期的形態(tài)也在增大,由其主導(dǎo)的流動(dòng)不穩(wěn)定性的發(fā)展有所滯后.不規(guī)則反射形成的激波結(jié)構(gòu)中三叉點(diǎn)的位置及自由前體激波更貼近界面(圖14(b)和圖14(c)).

圖12和圖13分別為激波沿橢圓短軸作用于半長(zhǎng)軸b=3.25 cm、半短軸a=1.9231 cm的氦氣柱與半長(zhǎng)軸b=4.0 cm、半短軸a=1.5625 cm的氦氣柱過(guò)程中界面和波系演化情況.由于橢圓界面曲率的變化,誘導(dǎo)激波角度較小,初始?jí)毫μ荻群兔芏忍荻鹊膴A角也較小,因而此時(shí)界面變形和波系演化與上述的三種幾何構(gòu)型表現(xiàn)不同.在波系演化上,激波作用后產(chǎn)生的是規(guī)則反射現(xiàn)象(圖12(a)和圖13(a)),圖14(d)和圖14(e)給出了詳細(xì)的波系情況,激波作用于重/輕界面后,在氦氣柱內(nèi)產(chǎn)生了一道向下游傳播的透射激波,同時(shí)反射出一道在空氣中向上游傳播的稀疏波.在界面變形方面,激波作用后,界面的前端并未迅速產(chǎn)生空氣射流,相反地在激波壓縮作用下,界面上游較大區(qū)域界面曲率變小,產(chǎn)生了類(lèi)似平面結(jié)構(gòu)的狀態(tài)(圖12(b)和圖12(c),圖13(b)—圖13(e)),其后隨著渦在這一平面結(jié)構(gòu)的末端產(chǎn)生,界面彎曲并向下游滲透,流動(dòng)繼續(xù)發(fā)展,最終形成與上述三種情況類(lèi)似的兩個(gè)獨(dú)立渦團(tuán)結(jié)構(gòu)(圖12(d)—(h)、圖13(f)—(h)).

圖10 幾何構(gòu)型2數(shù)值紋影圖(激波沿橢圓長(zhǎng)軸作用于氣柱;a=3.25 cm,b=1.9231 cm)Fig.10.Sequence of numerical schlieren images of shock-accelerated elliptic helium cylinder for Geometry 2(shock hitting along major axis of ellipse;a=3.25 cm,b=1.9231 cm).

圖11 幾何構(gòu)型3數(shù)值紋影圖(激波作用于半徑為2.5 cm的圓形氣柱)Fig.11.Sequence of numerical schlieren images of shock-accelerated elliptic helium cylinder for Geometry 3(shock hitting circular helium cylinder;a=b=2.5 cm).

圖12 幾何構(gòu)型4數(shù)值紋影圖(激波沿橢圓短軸作用于氣柱;a=1.9231 cm,b=3.25 cm)Fig.12.Sequence of numerical schlieren images of shock-accelerated elliptic helium cylinder for Geometry 4(shock hitting along minor axis of ellipse;a=1.9231 cm,b=3.25 cm).

圖13 幾何構(gòu)型5數(shù)值紋影圖(激波沿橢圓短軸作用于氣柱;a=1.5625 cm,b=4.0 cm)Fig.13.Sequence of numerical schlieren images of shock-accelerated elliptic helium cylinder for Geometry 5(shock hitting along minor axis of ellipse;a=1.5625 cm,b=4.0 cm).

圖14 不同形狀下激波作用氣柱20μs后的典型波系結(jié)構(gòu) (a)幾何構(gòu)型1;(b)幾何構(gòu)型2;(c)幾何構(gòu)型3;(d)幾何構(gòu)型4;(e)幾何構(gòu)型5Fig.14.Typical wave systems derived from numerical simulation for shock-accelerated elliptic cylinder at 20μs:(a)Geometry 1;(b)Geometry 2;(c)Geometry 3;(d)Geometry 4;(e)Geometry 5.

2)界面特征點(diǎn)位置變化

為了定量化地研究界面形態(tài)的變化,依據(jù)數(shù)值模擬得到的參數(shù)分布,以氦氣質(zhì)量分?jǐn)?shù)0.01作為物質(zhì)界面,得到界面三個(gè)特征點(diǎn)的位置(界面上游位置upstream、界面下游位置downstream、射流位置jet)變化如圖15所示,圖中同樣給出了特征點(diǎn)的位置示意圖.

圖15 界面特征點(diǎn)位置及界面長(zhǎng)度和高度的變化 (a)幾何構(gòu)型1;(b)幾何構(gòu)型2;(c)幾何構(gòu)型3;(d)幾何構(gòu)型4;(e)幾何構(gòu)型5Fig.15.Movement of distorted upstream,jet and downstream interfaces,and the interface length and height:(a)Geometry 1;(b)Geometry 2;(c)Geometry 3;(d)Geometry 4;(e)Geometry 5.

對(duì)于所有五種幾何構(gòu)型界面,上游界面位置upstream隨時(shí)間的變化可以分為兩個(gè)階段:流動(dòng)的初期,上游界面在誘導(dǎo)激波的作用下獲得一個(gè)較大的速度ui;其后,隨著透射激波在氣柱下游界面的反射激波作用于上游界面后,上游界面重新獲得一個(gè)新的速度uf. 在兩個(gè)階段上,幾何構(gòu)型1、幾何構(gòu)型2、幾何構(gòu)型3、幾何構(gòu)型4和幾何構(gòu)型5的平均速度ui/uf分別為196.3/134.2 m/s,193.4/120.2 m/s,181.3/105.2 m/s,153.1/93.5 m/s,136.9/83.3 m/s.下游界面在透射激波作用前基本保持靜止,隨后同樣以一個(gè)基本恒定的速度沿流動(dòng)方向發(fā)展,其平均運(yùn)動(dòng)速度分別為150.5,148.4,138.8,128.6,121.8 m/s,后期渦環(huán)(射流處)的速度分別為184.4,159.1,143.2,143.3,144.0 m/s.界面長(zhǎng)度(length)定義為界面上下游之間的距離,在流動(dòng)發(fā)展初期,由于激波的壓縮作用,界面長(zhǎng)度首先呈現(xiàn)變小的趨勢(shì),流動(dòng)過(guò)程中存在一個(gè)長(zhǎng)度基本不變的平臺(tái)期,主要是因?yàn)榇藭r(shí)界面上下游速度相近,其后由于界面下游速度高于上游速度,界面長(zhǎng)度又有所上升,平臺(tái)期的長(zhǎng)度隨著橢圓參數(shù)a值的減小而減小.界面高度(height)在透射激波作用到氣柱上下緣之前基本不變,而后在透射激波的作用下獲得一個(gè)向外的速度,使得氣柱高度有所增加,增長(zhǎng)速率隨著橢圓參數(shù)a的減小而增加,其后界面高度變化主要取決于該過(guò)程中空氣射流或者渦旋的發(fā)展,特別是對(duì)于激波沿橢圓短軸作用的幾何構(gòu)型3和4,后期由于渦旋在內(nèi)側(cè)產(chǎn)生,發(fā)展過(guò)程中不斷彎曲內(nèi)側(cè)的界面,在獨(dú)立渦團(tuán)形成之間基本不對(duì)界面高度產(chǎn)生影響(圖15(d)和圖15(e)).

3)氣柱體積壓縮率

氣泡體積變化可以很好地反映流動(dòng)過(guò)程中不同介質(zhì)的混合情況,因而常被作為SBI問(wèn)題中的重要分析參數(shù).Giordano和Burtschell[9]定義了以體積壓縮率γc表征氣泡體積隨時(shí)間的變化,

圖16 不同橢圓界面下氣柱體積壓縮率變化Fig.16.Volume compressibility for different elliptic interfaces.

上述的數(shù)值結(jié)果與分析表明,界面初始形狀對(duì)SBI過(guò)程中界面變形和波系演化有很大的影響.界面變形和波系演化特點(diǎn)主要由橢圓結(jié)構(gòu)自身曲率的變化和激波作用方向共同決定.激波強(qiáng)度一定時(shí),激波在界面處的反射類(lèi)型依賴(lài)于激波誘導(dǎo)角度θ(incident angleθ)的大小.由于橢圓界面各處曲率不同,針對(duì)上述五種橢圓構(gòu)型及激波作用方向,采用一種簡(jiǎn)單的近似方法估計(jì)激波誘導(dǎo)角度,如圖17所示.在這種近似方法下得到的誘導(dǎo)角度θ分別為68.6?,59.4?,45.0?,30.6?和21.3?.分析結(jié)果顯示在激波誘導(dǎo)角度θ為68.6?,59.4?和45.0?時(shí)界面發(fā)生的是不規(guī)則反射現(xiàn)象,如圖14(a)—(c)所示,而30.6?和21.3?的誘導(dǎo)角度下發(fā)生簡(jiǎn)單的規(guī)則反射現(xiàn)象.盡管研究的界面結(jié)構(gòu)形式不同,Zhai等[32]在對(duì)正方形、等腰三角形和菱形氣柱的研究中,同樣指出在誘導(dǎo)角度為60.0?,60.3?和90.0?時(shí),激波在重/輕界面處的反射類(lèi)型為不規(guī)則反射.

同樣地,若將該角度用于衡量壓力梯度和密度梯度方向的不一致性,分析界面形態(tài)的變化.激波沿長(zhǎng)軸作用于氣柱(此時(shí)不一致性較大,平均夾角分別為68.6?,59.4?),氣柱前端產(chǎn)生空氣射流,向氦氣柱內(nèi)浸入,進(jìn)而主導(dǎo)流動(dòng)和混合的發(fā)展;反之,當(dāng)激波沿短軸作用于氣柱(此時(shí)不一致性較小,平均夾角分別為30.6?和21.3?),由于激波的壓縮作用,在上游界面形成平行于誘導(dǎo)激波的近平面結(jié)構(gòu),隨后渦旋在該平面結(jié)構(gòu)的末端產(chǎn)生并主導(dǎo)流動(dòng)和混合的發(fā)展過(guò)程.界面特征點(diǎn)速度及結(jié)構(gòu)參數(shù)也因此有所差異.

圖17 激波誘導(dǎo)角度θ的近似方法Fig.17.Approximate method of shock incident angle θ.

對(duì)界面變形、特征點(diǎn)位置和特征點(diǎn)速度的分析表明,激波強(qiáng)度和氣柱面積一定的情況下,激波沿橢圓長(zhǎng)軸作用于氣柱時(shí),離心率越大,流動(dòng)發(fā)展得越迅速,界面的變形及氦氣與周?chē)鷼怏w的混合越快;激波沿橢圓短軸的作用規(guī)律與此相反,并且整體流動(dòng)混合也慢于前者.所有的幾何構(gòu)型下,體積壓縮率變化趨勢(shì)基本相同,最終穩(wěn)定值的差異并不大,這與Giordano和Burtschell[9]關(guān)于體積壓縮率的結(jié)論也基本相符(體積壓縮率只與激波強(qiáng)度和氣體比熱比有關(guān)),認(rèn)為氣體體積壓縮率并不能很好地反映氣柱幾何結(jié)構(gòu)對(duì)流動(dòng)和混合的影響.

3.2.2 環(huán) 量

激波作用于不連續(xù)氣柱的過(guò)程中,由于斜壓機(jī)制的作用會(huì)在物質(zhì)界面處產(chǎn)生大量的渦旋.在這些沉積于界面的渦量的驅(qū)動(dòng)下,氣柱形態(tài)產(chǎn)生大的變形.由壓力梯度和密度梯度不一致而產(chǎn)生的斜壓渦w可以表示為

進(jìn)一步對(duì)氦氣柱與激波作用中產(chǎn)生的環(huán)量進(jìn)行計(jì)算,結(jié)果如圖18所示.數(shù)值模擬中環(huán)量Γ采用如下方式計(jì)算:

其中D表示計(jì)算域,w(x,y,t)表示垂直于區(qū)域D的渦量.

圖18 不同橢圓界面下上半部分環(huán)量的變化 (a)正環(huán)量;(b)負(fù)環(huán)量;(c)總環(huán)量Fig.18.Circulation versus time for different elliptic interfaces:(a)Positive circulation;(b)negative circulation;(c)total circulation.

環(huán)量的計(jì)算中只考慮模型的上半部分,負(fù)環(huán)量在其中起主導(dǎo)作用.激波沿橢圓長(zhǎng)軸作用于橢圓氣柱和圓形氣柱的過(guò)程中,正、負(fù)環(huán)量的變化過(guò)程基本相似.在作用的初期,由于激波和界面的作用,負(fù)環(huán)量迅速增加;其后隨著氣柱前緣空氣射流的產(chǎn)生和發(fā)展,正、負(fù)環(huán)量以基本相似的趨勢(shì)在增加,直到空氣射流基本發(fā)展完全,在小的速度剪切引起的次級(jí)渦的共同作用下,此后環(huán)量的變化相對(duì)復(fù)雜.激波沿橢圓短軸作用于氣柱的過(guò)程中環(huán)量的變化基本與上述過(guò)程相似,同樣的負(fù)環(huán)量在激波的作用下迅速增加,但是由于橢圓自身的幾何結(jié)構(gòu),此種作用方式下,激波與界面接觸時(shí)間較短,因而導(dǎo)致負(fù)環(huán)量增加時(shí)間較短,增加值也較小,其后在其他波系的作用下緩慢增加,并且隨著渦的產(chǎn)生,增長(zhǎng)速率有所增加,對(duì)于幾何構(gòu)型4,小的速度剪切造成的次級(jí)渦在后期同樣導(dǎo)致了環(huán)量值的下降.總體來(lái)看,在激波與氦氣柱相互作用的過(guò)程中,盡管橢圓氣柱結(jié)構(gòu)不同,但環(huán)量的變化趨勢(shì)基本相同,并且隨著橢圓結(jié)構(gòu)參數(shù)a的減小(從幾何1到幾何5),該過(guò)程中正、負(fù)環(huán)量所能達(dá)到的最大值(大小)在不斷減小,總環(huán)量的穩(wěn)定值也在不斷減小.

結(jié)合對(duì)界面及波系的分析,正、負(fù)環(huán)量的變化主要經(jīng)歷三個(gè)時(shí)期:1)作用初期,環(huán)量主要由誘導(dǎo)激波與界面的壓縮作用產(chǎn)生,此時(shí)正環(huán)量值緩慢增加,負(fù)環(huán)量值迅速增加;2)作用中期,環(huán)量主要受RMI產(chǎn)生的主渦控制,此時(shí)正、負(fù)環(huán)量值均迅速增加;3)作用后期,環(huán)量主要由RMI產(chǎn)生的主渦和小的速度剪切引起的次級(jí)渦共同控制,主渦傾向于增加環(huán)量值,次渦傾向于減小環(huán)量值.由于正負(fù)環(huán)量在中后期變化的一致性,總環(huán)量在整個(gè)發(fā)展過(guò)程中快速穩(wěn)定在一個(gè)基本恒定的值上,其中幾何構(gòu)型1達(dá)到穩(wěn)定期經(jīng)歷的時(shí)間最長(zhǎng),且最終穩(wěn)定值大.環(huán)量的變化同樣可以體現(xiàn)SBI過(guò)程中的激波壓縮效應(yīng)、RMI引起的主渦的產(chǎn)生發(fā)展及其對(duì)流動(dòng)和混合的作用,定量地反映流動(dòng)和混合機(jī)制.

3.2.3 總混合率

物質(zhì)擴(kuò)散率(混合率)在時(shí)間和空間上的分布可以很好地反映物質(zhì)之間的混合特征.Tomkins等[11]以物質(zhì)質(zhì)量分?jǐn)?shù)為基礎(chǔ),定義物質(zhì)瞬時(shí)標(biāo)量擴(kuò)散率或混合率(instantaneous scalar dissipation rate or mixing rate)χ(x,t)為

總混合率由瞬時(shí)混合率進(jìn)行空間積分得到,定義如下:

其中DHe為氦氣在混合氣體中的擴(kuò)散系數(shù),為初始氦氣最大質(zhì)量分?jǐn)?shù),YHe為氦氣質(zhì)量分?jǐn)?shù).

圖19為不同橢圓界面下氦氣總混合率隨時(shí)間的變化.根據(jù)總混合率的定義可知其大小主要取決于氦氣濃度梯度和界面面積的大小(二維中表現(xiàn)為接觸界面長(zhǎng)度).在作用初期,誘導(dǎo)激波一方面導(dǎo)致界面面積迅速減小,另一方面導(dǎo)致部分激波掠過(guò)的區(qū)域氦氣濃度梯度增加,從總混合率迅速下降的結(jié)果可以得知在這一階段,界面面積的減小占據(jù)著主導(dǎo)地位.其后隨著空氣射流和渦的產(chǎn)生,使得氦氣與空氣的接觸界面面積增加,進(jìn)而導(dǎo)致混合率又迅速變大.后期總混合率存在一個(gè)極大值,分析認(rèn)為此后雖然界面面積在主渦的作用下依然增長(zhǎng),但流動(dòng)中小的次級(jí)渦的出現(xiàn),傾向于將氦氣的濃度在空間上呈現(xiàn)均勻化分布的趨勢(shì),因此此時(shí)的濃度梯度有所減小,可以認(rèn)為作用的后期氦氣與空氣處在一個(gè)充分混合的狀態(tài)中,混合率的變化也因此較為復(fù)雜.總而言之,從總混合率的角度出發(fā),同樣可以將混合分為三個(gè)時(shí)期:1)激波壓縮期,此階段波系的作用占據(jù)主導(dǎo)地位,一方面導(dǎo)致界面面積急劇減小,另一方面使得部分激波掠過(guò)的區(qū)域氦氣濃度梯度增加,總混合率呈下降趨勢(shì);2)RMI階段,在此時(shí)期,RMI引起的主渦在流動(dòng)過(guò)程中不斷彎曲和伸展進(jìn)入氦氣的空氣,使得接觸界面面積增加,總混合率也隨之迅速增加;3)充分混合階段,此階段以總混合率達(dá)到最大值為開(kāi)始的標(biāo)志,流動(dòng)中出現(xiàn)的次級(jí)渦一方面使得氦氣在核心渦區(qū)分布更加均勻,另一方面和主渦共同作用使得界面面積進(jìn)一步有所增加,這一階段混合率變化比較復(fù)雜,該階段氣體處于充分混和狀態(tài).不同幾何構(gòu)型的對(duì)比表明,總混合率的時(shí)間和空間分布雖然在整體上基本相似,但具體的流動(dòng)和混合細(xì)節(jié)存在一定的差異,隨著橢圓幾何參數(shù)a的減小,流動(dòng)和混合發(fā)展的越來(lái)越慢,進(jìn)入充分混合階段所需的時(shí)間越來(lái)越長(zhǎng),尤是激波沿橢圓短軸作用氣柱的情況.

圖19 不同橢圓界面下氦氣的總混合率Fig.19.Total mixing rate versus time for different elliptic interfaces.

4 結(jié) 論

本文采用數(shù)值模擬方法研究了激波作用于面積相同,結(jié)構(gòu)不同的橢圓型氣柱過(guò)程中的流動(dòng)混合情況,通過(guò)對(duì)該過(guò)程中氣柱界面變形、波系演化、體積壓縮率、環(huán)量和總混合率的分析,得到以下結(jié)論.

1)當(dāng)激波沿橢圓長(zhǎng)軸作用于橢圓氣柱時(shí),由于橢圓前緣附近界面的曲率特征,激波入射角度較大,密度梯度與壓力梯度方向的不一致性也較大,此時(shí)誘導(dǎo)激波作用后界面發(fā)生不規(guī)則發(fā)射現(xiàn)象,界面最前端形成小的空氣射流,主導(dǎo)其后的流動(dòng)和混合過(guò)程,流動(dòng)的RMI出現(xiàn)的也較早,小的速度剪切引起的Kelvin-Helmholtz不穩(wěn)定性影響也越大,有利于氦氣更大程度地與環(huán)境介質(zhì)混合.并且此時(shí)離心率越大,流動(dòng)發(fā)展越快.

2)當(dāng)激波沿橢圓短軸作用于橢圓氣柱時(shí),此時(shí)短軸越短,激波誘導(dǎo)角度較小,密度梯度與壓力梯度方向的不一致性也較小,誘導(dǎo)激波作用后界面發(fā)生規(guī)則反射現(xiàn)象,前緣界面被壓縮成近平面結(jié)構(gòu),隨后渦在平面結(jié)構(gòu)的邊緣產(chǎn)生,主導(dǎo)流動(dòng)和混合的發(fā)展,這種作用情況下RM不穩(wěn)定性發(fā)展較為緩慢.此時(shí)離心率越大,流動(dòng)發(fā)展越慢.

3)對(duì)環(huán)量和總混合率的定量對(duì)比分析表明,SBI過(guò)程主要分為三個(gè)階段:激波壓縮期、RMI引起的主渦發(fā)展期、速度剪切引起的次級(jí)渦發(fā)展期.并且橢圓結(jié)構(gòu)參數(shù)a越大,環(huán)量值和總混合率也越大.結(jié)合對(duì)界面變形的分析可以推斷激波沿橢圓長(zhǎng)軸作用于氣柱,流動(dòng)過(guò)程中的RMI產(chǎn)生越早發(fā)展越快,越有利于氦氣與周?chē)h(huán)境氣體的混合,并且離心率越大,混合發(fā)展得越快越充分.