用湍流模型研究Richtmyer-Meshkov不穩定性誘導的湍流混合

楊 玟,王麗麗 ,張樹道 ,何長江,杭義洪

(北京應用物理與計算數學研究所,北京100094)

0 引 言

界面不穩定性的發展及其引起的湍流混合在各種背景下都可能遇到,如燃燒、旋轉機械、慣性約束聚變(ICF)、超新星爆炸和地球物理學。很多國家都在該領域開展了大量的研究,如美國的三大武器物理實驗室(LANL、LLNL和 SNL)、英國核武器研究院(AWE)、法國原子能委員會(CEA)等[1-3]。目前,隨著計算機技術的迅猛發展,界面不穩定性的高精度高分辨率數值模擬(如直接模擬DNS、大渦模擬LES)[12,14]已成為該領域研究的重要手段之一,它為揭示界面不穩定性發展的內在規律發揮了一定的作用。但是,由于界面不穩定性現象以及與其相關物理過程的復雜性,高精度計算不僅費時,而且,受計算機硬件條件限制,其模擬能力目前尚不能解決許多復雜的實際應用問題。因此,界面不穩定性誘導湍流混合的雷諾平均方法(RANS),即湍流模型研究近年來受到了越來越多的關注[4-6]。國內在這方面與國外的差距較大,至今尚未見相關的報導。

當激波加速兩種不同密度流體間擾動界面時產生Richtmyer-Meshkov不穩定性(RMI)。這種不穩定性理論上由 Richtmyer(1960)發現并描述[7],由Meshkov(1969)從實驗中證實[8]。它是分層流體湍流混合的典型機制,對ICF是相當重要的。ICF中包含氘氚燃料的小靶丸被激光產生的沖擊波壓縮,其目的是在被壓縮的ICF靶丸內部獲得足夠高的壓力和溫度來點燃燃料。外殼和內部燃料間界面上不穩定性產生的混合抑制聚變反應,可能是能量產出的一個限制因素。

RMI通常在激波管中被研究。幾乎在所有激波管實驗中,重流體和輕流體初始采用塑料薄膜分隔,它通常被直接置于薄絲網下面。當激波經過絲網時將薄膜分裂成碎片。激波管實驗大多提供湍流混合區寬度的紋影圖以及采用X射線、紅外吸收或反射和微分干涉測量法得到的平均密度。這些診斷僅給出湍流的間接測量:混合的增強或平均密度脈動的增強。法國原子能委員會的F.Poggi等[13]首次采用激光多普勒測速儀(LDA)測量了激波管中產生于沖擊波的氣相混合物的瞬態速度,給出了混合區與反射激波相互作用前后混合區速度脈動的發展。這些實驗第一次給出了湍流的直接測量,為校核湍流混合模型提供了重要數據。

G.Dimonte[6]發展了描述界面不穩定性誘導湍流混合自相似增長的k-L模型,其中湍動能產生項以簡單的浮阻力模型[4]為基礎來構造。多相流模型[5]具有顯著優于單相模型的特點是能夠描述分離,但這類模型復雜、數值上強度大,且有時不穩定。湍流模型應用于不穩定性時,必須面對如下一些困難:(1)需要給出初始湍流狀態；(2)當模擬湍流與激波相互作用時,為了抹掉激波和消除梯度,強烈的不連續性的數值處理需要數值耗散,使得它們在產生湍動能上人為地效率降低；(3)n階關聯弱依賴于(n-1)階關聯的封閉假設在激波的作用下不再有效；(4)封閉假設無法考慮不均勻性。針對上述難點,本文將一般工程應用中被廣泛采用的k-ε模型運用到不穩定性誘導湍流混合的研究中。

1 模型方程和數值方法

模型以N-S方程為基礎,但是包含依賴于湍動能k和耗散率ε的湍流粘性μT。流動參數,除壓力和密度采用系綜平均之外,其它都采用質量加權平均,即Favre平均。例如,速度可分解為平均分量和脈動分量,即 u=+,將它們代入到流動方程中。對方程取平均,可得到平均和脈動分量各自的演化方程。對于可壓流動,Favre平均以密度加權以致于最終的方程為復雜性漸增的所有階的平均流動的展開式[9]。但是,展開式必須以一組簡單的封閉假設結束。所求解的方程為:

上述方程中雷諾應力τij和湍動能產生項-需要進行封閉。通常雷諾應力采用平均速度來構造,沖擊存在時,尤其是加密網格時,具有梯度乘積的項可能變得不穩定,到目前為止,還未得到激波存在時保持可實現性約束(realizability constraints)的湍流剪切應力的形式。與G.Dimonte等[6]采用的省略偏應力張量的方法不同,本文采用施加通量限制器的方法來排除存在激波時可能會出現的很大的非物理值,即根據量綱分析和 RMI產生的特點(沖擊波作用于不同密度流體間界面),湍動能產生項的計算定義為:

求解的程序中采用具有總變差遞減(Total Variation Diminishing,簡稱 TVD)保持性質的三階Runge-Kutta方法進行時間積分[10]:

對流項采用高階WENO方法重構網格邊界的通量,擴散項采用交替方向隱式迭代(ADI)法求解。

由于k-ε模型在工程中被廣泛應用,模型中的常數已經過仔細地研究而被確定,即所有常數通過數值實驗或與不可壓流體中進行的典型實驗比較來確定。新的常數 λm,λM,Cε0和 σρ根據激波管實驗標定 。本文中模型常數的選取為:CD=0.09,λm=0.1,λM=1.25,σe=0.9,σc=1.0,σk=0.87,σε=1.3,σρ=0.3,Cε0=1.06,Cε1=1.47,Cε2=1.90,Cε3=0,αρ=0.67 。

充分發展湍流的統計模型大多用于定常流動,其中初始條件僅承受數值約束。實際上,為了達到收斂和最小化計算花費,初始條件選得與想要得到的解盡可能接近。在非定常問題中情況很不一樣。因為封閉模型不計算不穩定性的初始階段,湍流狀態必須采用現象學原則和可得到的實驗數據來初始化。換句話說,為了采用相應平均量的第一梯度定理來封閉湍流通量,已假設充分湍流狀態。混合區湍動能的初始分布取為高度為kinit(=0.0028U)的三角形,耗散率的初值分布為εinit=0.164k/Linit,其中UI為激波經過后界面的運動速度,Linit為模型開啟時刻混合區寬度。

2 結果分析與討論

采用上述模型對Andronov等的激波管實驗[11]進行了模擬。在Andronov實驗中,分隔空氣和氦氣的有機薄膜位于距離管右端16.9cm處,Atwood數(A=(ρ2-ρ1)/(ρ2+ρ1))為0.75。馬赫數為1.3 的激波產生于管左側的空氣中,激波向右運動。計算域取為53.6cm×5cm。界面初始擾動由 ζ(y)=arcos給出。系數ar從一高斯分布中選取,W為y方向寬度,初始擾動中存在的波長范圍為0.25cm～1cm,擾動振幅被歸一化來使擾動標準差/2等于0.02cm。計算中界面在t=0.76ms時刻被加速,k-ε模型在t=0.8ms時打開(數值實驗發現模型打開時間的變化(t=0.9ms)對結果影響不大)。

圖1 空氣/氦氣湍流混合區寬度和軌跡隨時間的變化Fig.1 Widths and trajectory of Air/He turbulent mixing zone at different times

圖1給出了混合區寬度以及混合區邊界隨時間的變化曲線。從圖中可見,這些宏觀量的計算值與實驗值基本吻合,誤差在10%之內。雖然后期混合區寬度的計算值與實驗值吻合得較好,但是混合區邊界都被略微低估了。這可能是由于第一次反射激波與界面相互作用后界面運動速度的誤差造成的。但是,計算得到的入射激波賦予界面的運動速度為220m/s,與從沖擊波關系式中得到的理論值(218.72m/s)一致。

隨著計算機技術的突飛猛進,被實驗驗證過的直接模擬為模型的校核提供了大量有用的信息。圖2給出了不同時刻混合區內空氣摩爾份額沿橫向的分布。由圖可知,空氣與氦氣的相互滲透隨時間逐漸增強,模型計算結果與D.L.Youngs直接模擬結果[12]吻合得很好。混合區內湍動能的分布(如圖3所示)定性上是合理的,它很好地捕捉了第一次反射波到達界面時湍動能的增強。但是,由于湍動能衰減的弛豫時間遠大于激波與界面的作用時間,而計算的時間還不足夠長,因此圖3中未呈現明顯的湍流衰減。

圖2 不同時刻空氣/氦氣湍流混合區濃度沿橫向的分布Fig.2 Concentration distributions of Air/He turbulent mixing zone at different times

圖3 空氣/氦氣湍流混合區湍動能沿橫向的分布Fig.3 Distributions of turbulent kinetic energy of Air/Heturbulent mixing zone

正如上文中提到過的,Poggi等的激波管實驗[13]首次對湍流進行了直接測量(圖5(a))。該實驗在橫截面為8cm寬的正方形、長30cm的測試段進行。馬赫數為1.45的激波入射到六氟化硫中,下游的空氣由薄膜和絲網分隔(Atwood數為0.67)。計算域取為70cm×8cm。計算中初始條件的給定以實驗為基礎。界面初始擾動的波長取為與實驗中的絲網尺寸同一量級,λ=0.05cm～0.25cm；擾動振幅沒有任何實驗信息,本文對所有振幅取同一個值,該值(0.02cm)由數值實驗確定。

與Andronov實驗的計算不同的是該模擬中t=0時刻對應于激波沖擊擾動左端的時刻。圖4給出了混合區寬度隨時間的變化,從圖中可見,不穩定性發展早期計算值與實驗值吻合得較好,而后期則略低于實驗值,這與直接模擬結果定性上一致[14]。

圖4 六氟化硫/空氣湍流混合區寬度隨時間的變化Fig.4 Widths of SF6/Air turbulent mixing zone at different times

圖5(c)給出了不同時刻湍動能的分布,與圖5(b)的直接模擬結果[14]相比較:一方面,兩者都很好地再現了第一次反射激波與湍流混合區相互作用之前湍流的衰減和它被反射激波增強的特征；另一方面,模型計算無法捕捉入射激波與界面作用時湍動能的增強,這是由于模型在入射激波經過后不穩定性發展至非線性階段才打開。此外,模型中的湍流耗散機制似乎還需要進一步完善,因為模型計算的湍流衰減程度明顯低于直接模擬的結果。

圖5 六氟化硫/空氣混合區湍流量的實驗與模擬結果的比較Fig.5 Comparison of turbulence variants of SF6/Air mixing zonebetween experimental results and simulation results

實驗結果清楚地顯示了湍流的各向異性(如圖5(a)所示),但是k-ε模型無法描述該現象,因此我們進一步推導了考慮密度脈動和各向異性湍流的二階矩模型,并將它應用于不穩定性誘導湍流混合的研究中,該工作正在進行中。

3 結 論

本文將k-ε模型應用到 RMI誘導混合的研究中,經典的封閉關系采用代數關系式來補充,給出了合理描述RMI特征的湍動能產生項的表達式。采用上述模型對Andronov(1976)和Poggi(1998)的激波管實驗進行了模擬,結果與實驗和細觀模擬(DNS)吻合得較好,并且較好地再現了第一次反射激波與湍流混合區相互作用之前湍流的衰減和它被反射激波增強的特征。這些結果表明本文采用的模型封閉、模型常數、數值算法和程序實現是合適的。此外,該工作也為不穩定性誘導湍流混合模型研究的深入打下了良好的基礎。

致謝:本工作是在國家自然科學基金委員會-中國工程物理研究院聯合基金“界面運動不穩定性及湍流混合研究”(No.10676005)、中國工程物理研究院科學技術發展基金“界面不穩定性誘導湍流混合的多相湍流模型研究”(課題批準號:2008B0101006)和國防科技重點實驗室基金“含強間斷三維多尺度復雜流場的高精度歐拉方法研究”(編號:9140C6901010702)的資助下完成的。

[1]LLOR A.Statistical hydrodynamic models for developed mixing instability flows[A].Lecture Notes in Physics[M].Springer,Heidelberg,2005.

[2]ALON U,HECHT J,OFER D,SHVART SD.Power laws and similarity of Rayleigh-Taylor and Richtmyer-Meshkov mixing fronts at all density ratios[J].Phys.Rev.Lett.,1995,74(4):534-537.

[3]CRANFILL CW.A new multifluid turbulent-mix model[R].LANL,Report LA-UR-92-2484,1992.

[4]DIMONTE G.Spanwise homogeneous buoyancy-drag model for Rayleigh-Taylor mixing and experimental evaluation[J].Phys.Plasmas,2000,7(6):2255-2269.

[5]CHEN Y,GLIMM J,SHARP DA,ZHANG Q.A twophase flow model of the Rayleigh-Taylor mixing zone[J].Phys.Fluids,1996,8(3):816-825.

[6]DIMONTEG,TIPTON R.K-L turbulence model for the self-similar growth of the Rayleigh-Taylor and Richtmyer-Meshkov instabilities[J].Phys.Fluids,2006,18:085101.

[7]RICHTMYER R D.Taylor instability in shock acceleration of compressible fluids[J].Commun.Pure Ap pl.Math.,1960,13:297-319.

[8]MESHKOV E E.Interface of two gases accelerated by a shock wave[J].Fluid Dynamics,1969,4:101.

[9]周力行.湍流兩相流動與燃燒的數值模擬[M].清華大學出版社,1991.

[10]張涵信,沈孟育.計算流體力學-差分方法的原理和應用[M].國防工業出版社,2003.

[11]ANDRONOV V,BAKHRAKH SM,MESHKOV E E,MOKHOV V N,NIKIFOROV V V,PEVNITSKII A V,TOLSHMYAKOV A I.Turbulent mixing at contact surface accelerated by shock waves[J].Sov.Phys.JETP,1976,44:424.

[12]YOUNGS D L.Numerical simulation of mixing by Rayleigh-Taylor and Richtmyer-Meshkov instabilities[J].Laser and Particle Beams,1994,12(4):725-750.

[13]POGGIF,THOREMBEY M H,RODRIGUEZ G.Velocity measurements in turbulent gaseous mixtures induced by Richtmyer-Meshkov instability[J].Phys.Fluids,1998,10(11):2698-2700.

[14]M UGLERC,GAUTHIER S.Two-dimensional Navier-Stokes simulations of gaseous mixtures induced by Richtmyer-Meshkov instability[J].Phys.Fluids,2000,12(7):1783-1798.