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

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

0 引 言

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

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

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

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

1 模型方程和數(shù)值方法

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

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

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

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

由于k-ε模型在工程中被廣泛應用,模型中的常數(shù)已經(jīng)過仔細地研究而被確定,即所有常數(shù)通過數(shù)值實驗或與不可壓流體中進行的典型實驗比較來確定。新的常數(shù) λm,λM,Cε0和 σρ根據(jù)激波管實驗標定 。本文中模型常數(shù)的選取為: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 。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

致謝:本工作是在國家自然科學基金委員會-中國工程物理研究院聯(lián)合基金“界面運動不穩(wěn)定性及湍流混合研究”(No.10676005)、中國工程物理研究院科學技術(shù)發(fā)展基金“界面不穩(wěn)定性誘導湍流混合的多相湍流模型研究”(課題批準號: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]周力行.湍流兩相流動與燃燒的數(shù)值模擬[M].清華大學出版社,1991.

[10]張涵信,沈孟育.計算流體力學-差分方法的原理和應用[M].國防工業(yè)出版社,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.