激波作用不同橢圓氦氣柱過程中流動混合研究

李冬冬 王革 張斌

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

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

1 引 言

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

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

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

2 數值方法與計算模型

2.1 數值方法

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

其中U為守恒參數組成的向量,F和G為對流通量向量,Fv和Gv為黏性及擴散通量向量,具體表達式為

這里ρ,p和E分別代表混合物密度、壓力以及單位質量的總能量;u和v是混合物速度矢量在x和y方向上的分量;Yi是組分i的質量分數;Ng為物質種類總數;Jx,i和Jy,i分別為組分i的擴散在x和y方向的分量;qx和qy分別為熱擴散量在x和y方向的分量;σ為偏應力張量.多組分流動的理想氣體狀態方程為

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

其中CPi為物質i的定壓比熱,采用美國國家航空航天局(National Aeronautics and Space Administration,NASA)提供的擬合多項式進行計算;h0fi為物質i的標準焓;H為物質的總焓;T0為參考溫度.混合物的定壓比熱CP采用各物質的質量分數Yi和定壓比熱CPi加權平均得到.音速c是氣相混合物的一個重要屬性,定義為

對于理想氣體音速c可以簡化為

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

物質i的質量擴散通量Ji為

其中是物質i的擴散速度,其計算公式為

其中Xi為物質的摩爾分數,Di,mix為物質i在混合物中的混合平均擴散系數.

忽略Dufour效應,熱擴散矢量q計算表達式如下:

其中hi為物質i的焓值,同樣采用由NASA提供的多項式計算;λ為混合物的導熱系數.

流體間正應力和剪應力分別為

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

本文采用簡化的混合方法來計算混合物的輸運性質參數 (ξ=μ,λ,κ)[24]:

對于物質剪切黏性μ、體積黏性κ和熱傳導系數λ,n分別為6,4/3和4.

物質的混合平均擴散系數Di,mix采用下式計算:

其中Dij為物質i在物質j中的二元擴散系數,Xj為物質j的摩爾分數.

單物質的剪切黏性、熱傳導系數、二元擴散系數基于Chapman-Enskog理論[24,25]和Lennard-Jones參數[26]計算獲得,純物質的體積黏性κ依據Ern和Giovangigli[24]的簡化方法獲得.

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

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

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

其半離散格式為

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

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

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

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

圖2 網格界面對流通量Fig.2.Convective flux at cell face.

2.2 計算模型

文中采用如圖3所示的計算模型,計算條件參照Haas和Sturtevant[7]的實驗條件給定,氣柱內輕氣體為氦氣,周圍介質為空氣,誘導激波馬赫數為1.22,位于氣泡右側,波前氣體靜止,溫度和壓力分別為293 K和1 atm(1 atm=1.01325×105Pa),波后氣體參數由Rankine-Hugoniot條件獲得.由于對稱性,選取模型的上半部分進行計算,上邊界采用固體反射邊界,對稱軸采用對稱邊界條件,左右邊界上采用無反射邊界條件,各參數梯度為0.

以圖3所示的計算模型為基礎,并在網格無關性檢驗和算例驗證可靠性的前提下,針對面積相同、幾何構型不同(圓形,兩個激波沿橢圓長軸作用于氣柱,兩個激波沿橢圓短軸作用于氣柱)的氦氣柱與激波的相互作用過程,分析界面變形、波系演化、界面結構參數(界面高度和長度)、氣柱體積壓縮率、總混合率、環量等的變化,探究氣泡內介質和周圍環境氣體的混合機理和不同幾何構型下氣體混合的優劣.以半徑2.5 cm的圓形截面氣柱面積19.63 cm2為標準來設計氣柱截面形態,如圖4所示.為了方便表示,采用統一的標準橢圓方程來約束界面形狀:

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

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

不同幾何構型的橢圓氣柱對應的參數如表1所列,其中Geometry 1和Geometry 2中激波沿橢圓長軸作用氣柱,Geometry 3為圓形氣柱,Geometry 4和Geometry 5中激波沿橢圓短軸作用氣柱.

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

為防止虛假渦量的產生,采用有限厚度擴散層[29]獲得更加光滑的初始氣泡界面,(20)式為氦氣質量分數初始分布函數,

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

3 結果與討論

3.1 網格無關性檢驗與算例驗證

為了考察計算方法的可靠性,選擇與參考文獻[7]中激波作用氦氣柱的實驗結果進行對比.首先,以Geometry 3圓形氣柱為例,采用三種不同尺寸的網格進行網格無關性檢驗,網格數量分別為粗網格(600×83,Grid 1)、中等網格(1200×166,Grid 2)、細網格(2400×330,Grid 3),相應的氣柱內網格數量分別為80,160和320.

圖5是以氦氣質量分數顯示的不同網格尺度下的界面變形情況.結果表明在三種網格尺度下界面變形過程基本相似,較細的網格體現出更多的界面細節特征.圖6為不同網格尺度下界面三個特征點(特征點位置如圖6中右下角所示)在激波作用He氣泡過程中的位置變化(界面以氦氣質量分數0.01作為標準),三種不同網格下所獲得的結果基本相同.在后續的計算分析中為了獲得較為清晰的流動細節均采用細網格(2400×330).

圖5 網格無關性檢驗(左,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 不同網格尺寸下特征點位置的變化Fig.6.Movement of characteristic points for different grids.

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

為進一步檢驗方法的準確性,將在細網格下獲得的數值紋影結果與文獻[7]實驗采集到的紋影結果進行定性的對比,如圖8所示,數值紋影圖清晰地表明了作用過程中界面變形和波系演化過程,進一步說明所采用數值方法的可靠性和精度.

圖7 激波作用于氦氣泡前期,特征點位置的變化及與其他計算結果[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 數值紋影與實驗結果[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 不同幾何構型下氦氣混合過程分析

3.2.1 界面形態和特征

1)界面變形和波系演化

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

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

圖9 幾何構型1數值紋影圖(激波沿橢圓長軸作用氣柱;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分別為激波沿橢圓長軸作用于半長軸a=3.25 cm、半短軸b=1.9231 cm的氦氣柱與半徑為2.5 cm的圓形氦氣柱過程中界面和波系演化情況.其中波系演化模式及界面變形情況與幾何構型1基本相似,但不難發現由于長軸的變短和短軸長度的增加(離心率變小),激波入射角度、壓力梯度和密度梯度之間的方向發生了改變,一方面導致了界面前端空氣射流形成時間增加,射流初期的形態也在增大,由其主導的流動不穩定性的發展有所滯后.不規則反射形成的激波結構中三叉點的位置及自由前體激波更貼近界面(圖14(b)和圖14(c)).

圖12和圖13分別為激波沿橢圓短軸作用于半長軸b=3.25 cm、半短軸a=1.9231 cm的氦氣柱與半長軸b=4.0 cm、半短軸a=1.5625 cm的氦氣柱過程中界面和波系演化情況.由于橢圓界面曲率的變化,誘導激波角度較小,初始壓力梯度和密度梯度的夾角也較小,因而此時界面變形和波系演化與上述的三種幾何構型表現不同.在波系演化上,激波作用后產生的是規則反射現象(圖12(a)和圖13(a)),圖14(d)和圖14(e)給出了詳細的波系情況,激波作用于重/輕界面后,在氦氣柱內產生了一道向下游傳播的透射激波,同時反射出一道在空氣中向上游傳播的稀疏波.在界面變形方面,激波作用后,界面的前端并未迅速產生空氣射流,相反地在激波壓縮作用下,界面上游較大區域界面曲率變小,產生了類似平面結構的狀態(圖12(b)和圖12(c),圖13(b)—圖13(e)),其后隨著渦在這一平面結構的末端產生,界面彎曲并向下游滲透,流動繼續發展,最終形成與上述三種情況類似的兩個獨立渦團結構(圖12(d)—(h)、圖13(f)—(h)).

圖10 幾何構型2數值紋影圖(激波沿橢圓長軸作用于氣柱;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 幾何構型3數值紋影圖(激波作用于半徑為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 幾何構型4數值紋影圖(激波沿橢圓短軸作用于氣柱;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 幾何構型5數值紋影圖(激波沿橢圓短軸作用于氣柱;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后的典型波系結構 (a)幾何構型1;(b)幾何構型2;(c)幾何構型3;(d)幾何構型4;(e)幾何構型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)界面特征點位置變化

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

圖15 界面特征點位置及界面長度和高度的變化 (a)幾何構型1;(b)幾何構型2;(c)幾何構型3;(d)幾何構型4;(e)幾何構型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.

對于所有五種幾何構型界面,上游界面位置upstream隨時間的變化可以分為兩個階段:流動的初期,上游界面在誘導激波的作用下獲得一個較大的速度ui;其后,隨著透射激波在氣柱下游界面的反射激波作用于上游界面后,上游界面重新獲得一個新的速度uf. 在兩個階段上,幾何構型1、幾何構型2、幾何構型3、幾何構型4和幾何構型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.下游界面在透射激波作用前基本保持靜止,隨后同樣以一個基本恒定的速度沿流動方向發展,其平均運動速度分別為150.5,148.4,138.8,128.6,121.8 m/s,后期渦環(射流處)的速度分別為184.4,159.1,143.2,143.3,144.0 m/s.界面長度(length)定義為界面上下游之間的距離,在流動發展初期,由于激波的壓縮作用,界面長度首先呈現變小的趨勢,流動過程中存在一個長度基本不變的平臺期,主要是因為此時界面上下游速度相近,其后由于界面下游速度高于上游速度,界面長度又有所上升,平臺期的長度隨著橢圓參數a值的減小而減小.界面高度(height)在透射激波作用到氣柱上下緣之前基本不變,而后在透射激波的作用下獲得一個向外的速度,使得氣柱高度有所增加,增長速率隨著橢圓參數a的減小而增加,其后界面高度變化主要取決于該過程中空氣射流或者渦旋的發展,特別是對于激波沿橢圓短軸作用的幾何構型3和4,后期由于渦旋在內側產生,發展過程中不斷彎曲內側的界面,在獨立渦團形成之間基本不對界面高度產生影響(圖15(d)和圖15(e)).

3)氣柱體積壓縮率

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

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

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

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

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

對界面變形、特征點位置和特征點速度的分析表明,激波強度和氣柱面積一定的情況下,激波沿橢圓長軸作用于氣柱時,離心率越大,流動發展得越迅速,界面的變形及氦氣與周圍氣體的混合越快;激波沿橢圓短軸的作用規律與此相反,并且整體流動混合也慢于前者.所有的幾何構型下,體積壓縮率變化趨勢基本相同,最終穩定值的差異并不大,這與Giordano和Burtschell[9]關于體積壓縮率的結論也基本相符(體積壓縮率只與激波強度和氣體比熱比有關),認為氣體體積壓縮率并不能很好地反映氣柱幾何結構對流動和混合的影響.

3.2.2 環 量

激波作用于不連續氣柱的過程中,由于斜壓機制的作用會在物質界面處產生大量的渦旋.在這些沉積于界面的渦量的驅動下,氣柱形態產生大的變形.由壓力梯度和密度梯度不一致而產生的斜壓渦w可以表示為

進一步對氦氣柱與激波作用中產生的環量進行計算,結果如圖18所示.數值模擬中環量Γ采用如下方式計算:

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

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

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

結合對界面及波系的分析,正、負環量的變化主要經歷三個時期:1)作用初期,環量主要由誘導激波與界面的壓縮作用產生,此時正環量值緩慢增加,負環量值迅速增加;2)作用中期,環量主要受RMI產生的主渦控制,此時正、負環量值均迅速增加;3)作用后期,環量主要由RMI產生的主渦和小的速度剪切引起的次級渦共同控制,主渦傾向于增加環量值,次渦傾向于減小環量值.由于正負環量在中后期變化的一致性,總環量在整個發展過程中快速穩定在一個基本恒定的值上,其中幾何構型1達到穩定期經歷的時間最長,且最終穩定值大.環量的變化同樣可以體現SBI過程中的激波壓縮效應、RMI引起的主渦的產生發展及其對流動和混合的作用,定量地反映流動和混合機制.

3.2.3 總混合率

物質擴散率(混合率)在時間和空間上的分布可以很好地反映物質之間的混合特征.Tomkins等[11]以物質質量分數為基礎,定義物質瞬時標量擴散率或混合率(instantaneous scalar dissipation rate or mixing rate)χ(x,t)為

總混合率由瞬時混合率進行空間積分得到,定義如下:

其中DHe為氦氣在混合氣體中的擴散系數,為初始氦氣最大質量分數,YHe為氦氣質量分數.

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

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

4 結 論

本文采用數值模擬方法研究了激波作用于面積相同,結構不同的橢圓型氣柱過程中的流動混合情況,通過對該過程中氣柱界面變形、波系演化、體積壓縮率、環量和總混合率的分析,得到以下結論.

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

2)當激波沿橢圓短軸作用于橢圓氣柱時,此時短軸越短,激波誘導角度較小,密度梯度與壓力梯度方向的不一致性也較小,誘導激波作用后界面發生規則反射現象,前緣界面被壓縮成近平面結構,隨后渦在平面結構的邊緣產生,主導流動和混合的發展,這種作用情況下RM不穩定性發展較為緩慢.此時離心率越大,流動發展越慢.

3)對環量和總混合率的定量對比分析表明,SBI過程主要分為三個階段:激波壓縮期、RMI引起的主渦發展期、速度剪切引起的次級渦發展期.并且橢圓結構參數a越大,環量值和總混合率也越大.結合對界面變形的分析可以推斷激波沿橢圓長軸作用于氣柱,流動過程中的RMI產生越早發展越快,越有利于氦氣與周圍環境氣體的混合,并且離心率越大,混合發展得越快越充分.