王 濤，陶 鋼，柏勁松，李 平，汪 兵，杜 磊

(1.南京理工大學(xué)能源與動(dòng)力工程學(xué)院，江蘇南京 210094；2.中國工程物理研究院流體物理研究所，四川綿陽 621999)

界面不穩(wěn)定性，包括涉及沖擊驅(qū)動(dòng)的Richtmyer-Meshkov(RM)不穩(wěn)定性[1-2]、慣性力驅(qū)動(dòng)的Rayleigh-Tayloy(RT)不穩(wěn)定性[3-4]以及剪切驅(qū)動(dòng)的Kelvin-Helmholtz(KH)不穩(wěn)定性[5]，是我們在自然界和工程應(yīng)用中經(jīng)常遇到的一類物理現(xiàn)象，比如天體物理中的超新星爆發(fā)、慣性約束聚變、磁約束聚變、核武器爆炸、超音速燃燒系統(tǒng)等。界面不穩(wěn)定性及其誘發(fā)的湍流混合作為流體動(dòng)力學(xué)領(lǐng)域的重要和前沿課題，一直以來受到流體力學(xué)工作者，特別是各國武器實(shí)驗(yàn)室的格外關(guān)注。在實(shí)際問題中，各類界面不穩(wěn)定性經(jīng)常是共存的，使該問題的研究變得更復(fù)雜。

眾所周知，界面不穩(wěn)定性的發(fā)展受初始擾動(dòng)條件影響很大，前人也多次對(duì)界面不穩(wěn)定性發(fā)展受初始擾動(dòng)的影響進(jìn)行過分析和研究。比如，在RM不穩(wěn)定性方面，Thornber等人[6]通過一系列精確控制的數(shù)值實(shí)驗(yàn)，研究了無反射沖擊情況下三維多模態(tài)初始擾動(dòng)條件對(duì)混合區(qū)增長率的影響，結(jié)果顯示，混合區(qū)的增長強(qiáng)烈地依賴于初始條件。Schilling等人[7]證實(shí)，從初始沖擊后至反射沖擊前，混合區(qū)寬度對(duì)初始擾動(dòng)振幅很敏感。Leinov等人[8]研究顯示，反射沖擊波再加載后，混合區(qū)的演化與反射沖擊波到達(dá)時(shí)的界面振幅無關(guān)，但是與反射沖擊波的強(qiáng)度有關(guān)。Balasubramanian等人[9]的實(shí)驗(yàn)顯示，初始擾動(dòng)波長和振幅對(duì)反射沖擊后氣簾不穩(wěn)定性混合有影響。Bai等人[10]關(guān)于初始流場非均勻性對(duì)擾動(dòng)增長影響的研究表明，在反射沖擊波再加載之前的線性和弱非線性階段，初始流場非均勻性對(duì)不穩(wěn)定性的演化具有相當(dāng)大的影響；而在再加載之后的強(qiáng)非線性階段，非均勻性效應(yīng)急劇衰減。王濤等人[11]也分析研究了無反射沖擊波作用下單模態(tài)RM不穩(wěn)定性發(fā)展受初始擾動(dòng)的影響，結(jié)果表明，擾動(dòng)增長對(duì)初始條件有很強(qiáng)的依賴性。但是,反射沖擊波作用后混合區(qū)發(fā)展受初始擾動(dòng)的影響作用還不是很明確，需要進(jìn)一步深入開展相關(guān)研究。本研究通過大渦數(shù)值模擬方法,研究多次沖擊作用下三維多模態(tài)RM不穩(wěn)定性發(fā)展對(duì)初始擾動(dòng)條件的依賴性。

2 計(jì)算方法

數(shù)值模擬采用自研的可壓縮多介質(zhì)黏性流動(dòng)和湍流大渦(Multi-Viscous-Flow and Turbulence,MVFT)模擬程序[12-13]。該大渦模擬的流場控制方程為經(jīng)過Favre濾波后的可壓縮多介質(zhì)黏性流動(dòng)Navier-Stokes方程組

(1)

(2)

采用算子分裂技術(shù)將方程組(1)式描述的物理過程分解為3個(gè)子過程進(jìn)行計(jì)算，即整個(gè)通量計(jì)算分解為無黏通量、黏性通量和熱通量計(jì)算3部分。無黏通量的計(jì)算采用多介質(zhì)的高精度分段拋物線方法(Piecewise Parabolic Method，PPM)，兩步Lagrange-Remapping型的PPM方法分以下4個(gè)步驟進(jìn)行：(1) 物理量分段拋物插值；(2) 近似Riemann問題求解；(3) Lagrange方程組推進(jìn)求解；(4) 將物理量映射到靜止的歐拉網(wǎng)格上。然后在無黏通量的基礎(chǔ)上，采用二階空間中心差分方法和兩步Runge-Kutta時(shí)間推進(jìn)方法求解黏性通量、熱通量及標(biāo)量輸運(yùn)通量。詳細(xì)的數(shù)值方法參考文獻(xiàn)[15-16]。

3 結(jié)果分析

計(jì)算模型來自Erez等人[17]開展的RM不穩(wěn)定性激波管實(shí)驗(yàn)。激波管測試段長度為180 mm，內(nèi)截面尺寸為80 mm×80 mm。實(shí)驗(yàn)氣體為空氣和六氟化硫(SF6)，初始計(jì)算參數(shù)列于表1。

表1 氣體初始參數(shù)Table 1 Initial properties of air and SF6

計(jì)算模型如圖1所示。馬赫數(shù)為1.2的沖擊波加速空氣/六氟化硫(Air/SF6)界面。計(jì)算域大小為[-60 mm,180 mm]×[0 mm,80 mm]×[0 mm,80 mm]，離散網(wǎng)格數(shù)為768×256×256，網(wǎng)格大小為0.312 5 mm，采用6×2×2個(gè)CPU進(jìn)行并行計(jì)算。多模態(tài)初始擾動(dòng)通過8個(gè)正弦模態(tài)疊加而成

(3)

計(jì)算中，初始擾動(dòng)主模態(tài)的橫向波長為0.80、1.00、1.25、1.60、2.00、2.50、3.20、4.00 mm。為了分析初始擾動(dòng)條件對(duì)RM不穩(wěn)定性及湍流混合發(fā)展的影響，我們通過改變初始擾動(dòng)的振幅來改變初始擾動(dòng)強(qiáng)度。初始擾動(dòng)強(qiáng)度(PS)定義為初始擾動(dòng)振幅和波長之比，對(duì)于多模態(tài)初始擾動(dòng)，PS=η0/λ，其中λ≈2 mm為平均波長。模型的初始參數(shù)列于表2。

圖1 計(jì)算模型Fig.1 Computation model

Caseη0/(mm)PS10.070.03520.140.07030.280.14040.560.28051.120.560

沖擊波加速空氣/SF6擾動(dòng)界面，向SF6氣體中產(chǎn)生一個(gè)透射沖擊波。該透射沖擊波從激波管尾端反射回來再次沖擊界面時(shí)(tres)，向SF6氣體中反射一個(gè)稀疏波。反射稀疏波從激波管尾端反射回來和界面相互作用(texp1)后又向SF6氣體中反射一個(gè)壓縮波，當(dāng)這個(gè)壓縮波再反射回來和界面相互作用時(shí)(tcom1)，又向SF6氣體中反射一個(gè)稀疏波。這樣的一系列壓縮波和稀疏波交替在測試段中運(yùn)動(dòng)，并和界面相互作用，形成對(duì)界面的多次沖擊，但是每次沖擊作用后，能量的轉(zhuǎn)移使波的強(qiáng)度進(jìn)一步降低，長時(shí)間后對(duì)湍流混合區(qū)的演化難以產(chǎn)生影響。

圖2 湍流混合區(qū)定義Fig.2 Definition of turbulent mixing zone

圖3 不同模型湍流混合區(qū)寬度隨時(shí)間的增長歷史Fig.3 Time histories of growth of TMZ width (dTMZ) of different models

圖3(b)給出了不同模型的湍流混合區(qū)寬度的增長歷史。可以看出，初始擾動(dòng)強(qiáng)度大，對(duì)應(yīng)的湍流混合區(qū)寬度較大，相應(yīng)的增長速率也較高，這主要體現(xiàn)在初始沖擊后至反射沖擊前和反射沖擊后至第一次反射稀疏波作用前兩個(gè)階段。而反射稀疏波作用后，初始擾動(dòng)對(duì)湍流混合區(qū)寬度增長的影響逐漸減弱，直至反射壓縮波作用后，湍流混合區(qū)的增長基本不受初始擾動(dòng)的影響，以基本一致的線性規(guī)律增長，增長速度為2.05 m/s，即在第一次反射稀疏波作用后，湍流混合區(qū)的發(fā)展逐漸失去對(duì)初始擾動(dòng)條件的記憶。表3列出了不同模型在初始沖擊后、反射沖擊后以及第一次反射稀疏波作用后湍流混合區(qū)寬度的增長因子。圖4是不同時(shí)刻Case 1模型中，SF6的體積分?jǐn)?shù)為0.1和0.9的等值面所顯示的湍流混合區(qū)圖像，表明混合區(qū)具有復(fù)雜的空間結(jié)構(gòu)。

表3 不同模型湍流混合區(qū)寬度的增長因子Table 3 Growth factors of TMZ width of different models

圖4 SF6的體積分?jǐn)?shù)為0.1和0.9的等值面所顯示的不同時(shí)刻湍流混合區(qū)圖像Fig.4 Images of TMZ visualized by the volume fraction isosurface of 0.1 and 0.9 for SF6 at different times

下面我們通過比較湍流混合區(qū)的統(tǒng)計(jì)量進(jìn)一步分析初始擾動(dòng)條件的影響。統(tǒng)計(jì)量包括湍動(dòng)能k、湍動(dòng)能耗散率ε和擬渦能Ω

(4)

(5)

(6)

圖5～圖7分別給出了不同模型湍流混合區(qū)的湍動(dòng)能峰值、湍動(dòng)能耗散率峰值及擬渦能峰值隨時(shí)間變化的歷史?？梢钥闯觯寒?dāng)界面受到波的沖擊作用時(shí)，沖擊波攜帶的能量通過斜壓作用沉積到界面上，促使界面快速發(fā)展，此刻混合區(qū)統(tǒng)計(jì)量也隨之急速增大；之后由于波的沖擊作用消失，統(tǒng)計(jì)量在耗散和擴(kuò)散作用下逐漸衰減。初始沖擊后，湍動(dòng)能和擬渦能以冪次律衰減，而湍動(dòng)能耗散率以指數(shù)規(guī)律衰減；反射沖擊后，湍動(dòng)能、湍動(dòng)能耗散率和擬渦能均以指數(shù)規(guī)律衰減；反射稀疏波作用后，這些統(tǒng)計(jì)量也以指數(shù)規(guī)率衰減，但是衰減因子不同；其后，由于波的強(qiáng)度太小，沒有足夠的能量沉積，所以統(tǒng)計(jì)量以漸近方式衰減。

湍動(dòng)能、湍動(dòng)能耗散率及擬渦能都隨著初始擾動(dòng)強(qiáng)度的增大而增大，也就是說，在相同強(qiáng)度沖擊波作用下，較大的初始擾動(dòng)會(huì)導(dǎo)致更強(qiáng)的湍流脈動(dòng)，湍動(dòng)能的傳輸更快，湍動(dòng)能耗散也更快，產(chǎn)生的斜壓渦量場也更強(qiáng)。在初始沖擊后至反射沖擊前和反射沖擊后至反射稀疏波作用這兩個(gè)階段，受初始擾動(dòng)的影響較大。第一次反射稀疏波作用后，湍流混合區(qū)的發(fā)展基本不受初始擾動(dòng)條件的影響。

圖5 不同模型湍流混合區(qū)湍動(dòng)能峰值時(shí)間歷史Fig.5 Time history of turbulent kinetic energy of different models in TMZ

圖6 不同模型湍流混合區(qū)湍動(dòng)能耗散率峰值時(shí)間歷史Fig.6 Time histories of turbulent kinetic energy dissipation rate of different models in TMZ

圖7 不同模型湍流混合區(qū)擬渦能峰值時(shí)間歷史Fig.7 Time histories of enstrophy of different models in TMZ

4 結(jié) 論

利用發(fā)展的適用于可壓縮多介質(zhì)黏性流動(dòng)和湍流大渦模擬的程序MVFT，對(duì)輕/重氣體(空氣/SF6)界面構(gòu)型下多次沖擊作用的三維多模態(tài)RM不穩(wěn)定性發(fā)展及其對(duì)初始擾動(dòng)條件的依賴性進(jìn)行了數(shù)值模擬。研究結(jié)果顯示，初始沖擊后，湍流混合區(qū)寬度以冪次律增長；反射沖擊后和第一次反射稀疏波作用后，湍流混合區(qū)寬度以指數(shù)規(guī)律增長，但其增長因子不同；第一次反射壓縮波作用后，近似以線性規(guī)律增長。而湍流混合區(qū)統(tǒng)計(jì)量則以相似的規(guī)律衰減。多模態(tài)RM不穩(wěn)定性發(fā)展對(duì)初始擾動(dòng)條件有很強(qiáng)的依賴性，主要體現(xiàn)在初始沖擊后至反射沖擊前和反射沖擊后至第一次反射稀疏波作用前兩個(gè)階段，而在第一次反射稀疏波與界面作用后，湍流混合區(qū)的發(fā)展逐漸失去對(duì)初始擾動(dòng)條件的記憶。

[1] RICHTMYER R D.Taylor instability in shock acceleration of compressible fluids [J].Commun Pure Appl Math,1960,13(2):297-319.

[2] MESHKOV E E.Instability of the interface of two gases accelerated by a shock wave [J].Sov Fluid Dyn,1969,4(5):101-104.

[3] RAYLEIGH LORD.Scientific papers Ⅱ [M].Cambridge:Cambridge University Press,1900:200.

[4] TAYLOR G I.The instability of liquid surfaces when accelerated in a direction perpendicular to their plane [J].Proc R Soc London (Ser A),1950,201:192.

[5] CHANDRASEKHAR S.Hydrodynamic and hydromagnetic stability [M].London:Oxford University Press,1961.

[6] THORNBER B,DRIKAKIS D,YOUNGS D L,et al.The influence of initial conditions on turbulent mixing due to Richtmyer-Meshkov instability [J].J Fluid Mech,2010,654:99-139.

[7] SCHILLING O,LATINI M.High-order WENO simulations of three-dimensional reshocked Richtmyer-Meshkov instability to late times:dynamics,dependence on initial conditions,and comparisons to experimental data [J].Acta Math Sci,2010,30(2):595-620.

[8] LEINOV E,MALAMUD G,ELBAZ Y,et al.Experimental and numerical investigation of the Richtmyer-Meshkov instability under re-shock conditions [J].J Fluid Mech,2009,626:449-475.

[9] BALASUBRAMANIAN S,ORLICZ G C,PRESTRIDGE K P,et al.Experimental study of initial condition dependence on Richtmyer-Meshkov instability in the presence of reshock [J].Phys Fluids,2012,24(3):034103.

[10] BAI J S,WANG B,WANG T,et al.Numerical simulation of the Richtmyer-Meshkov instability in initially nonuniform flows and mixing with reshock [J].Phys Rev E,2012,86(6):066319.

[11] 王 濤，柏勁松，李 平,等.單模態(tài)RM不穩(wěn)定性的二維和三維數(shù)值計(jì)算研究 [J].高壓物理學(xué)報(bào)，2013,27(2):277-286.(in English)

WANG T,BAI J S,LI P,et al.Two and three dimensional numerical investigations of the single-mode Richtmyer-Meshkov instability [J].Chinese Journal of High Pressure Physics,2013,27(2):277-286.

[12] WANG T,BAI J S,LI P,et al.Large-eddy simulations of the Richtmyer-Meshkov instability of rectangular interface accelerated by shock waves [J].Sci China Ser G,2010,53(5):905-914.

[13] BAI J S,LIU J H,WANG T,et al.Investigation of the Richtmyer-Meshkov instability with double perturbation interface in nonuniform flows [J].Phys Rev E,2010,81(2):056302.

[14] VREMAN A W.An eddy-viscosity subgrid-scale model for turbulent shear flow:algebraic theory and applications [J].Phys Fluids,2004,16(10):3670-3681.

[15] WANG T,BAI J S,LI P,et al.The numerical study of shock-induced hydrodynamic instability and mixing [J].Chin Phys B,2009,18(3):1127-1135.

[16] BAI J S,WANG T,LI P,et al.Numerical simulation of the hydrodynamic instability experiments and flow mixing [J].Sci China Ser G,2009,52(12):2027-2040.

[17] EREZ L,SADOT O,ORON D,et al.Study of the membrane effect on turbulent mixing measurements in shock tubes [J].Shock Waves,2000,10:241-251.

[18] SREBRO Y,ELBAZ Y,SADOT O,et al.A general buoyancy-drag model for the evolution of the Rayleigh-Taylor and Richtmyer-Meshkov instabilities [J].Laser Part Beams,2003,21:347-353.

[19] WANG T,BAI J S,LI P,et al.Large-eddy simulation of the multi-mode Richtmyer-Meshkov instability and turbulent mixing under reshock [J].High Energy Density Physics,2016,19:65-75.