流場非均勻性對非平面激波誘導的Richtmyer-Meshkov不穩定性影響的數值研究*

王 震，王 濤，柏勁松,2，肖佳欣

（1.中國工程物理研究院流體物理研究所，四川 綿陽 621999；2.中國工程物理研究院流體物理研究所沖擊波物理與爆轟物理實驗室，四川 綿陽 621999；3.北京臨近空間飛行器系統工程研究所，北京 100076）

RM不穩定性[1]及其誘導產生的湍流混合現象，由于其廣泛的應用如大特征尺度的超新星爆炸[2]、中等特征尺度的超燃沖壓發動機中超聲速燃燒[3]及小特征尺度的慣性約束聚變[4]等，而備受關注。其產生的關鍵機制為斜壓渦量沉積，即激波加載兩種密度不同的物質界面時，界面處密度梯度與壓力梯度不共線而導致斜壓渦量沉積（即 ? ρ×?p/ρ2≠0），進而界面處的初始擾動會被放大而產生失穩現象。

半個多世紀以來，學者們對RM不穩定性進行了廣泛深入的研究，結果表明界面、流場及激波是影響RM不穩定性發展的重要因素。Yang等[3]通過48組數據數值研究，發現激波強度、密度比、橢圓界面形狀等9種因素將影響界面的不穩定性發展，他的研究為后續RM不穩定性影響因素的研究起了指導作用。Zou等[5]、Bai等[6]分別實驗、數值模擬研究，表明平面激波加載長/短軸比值不同的橢圓形氣柱界面時，界面處沉積的渦量不同，因而橢圓形氣柱界面的形態存在差異。Bai等[7]以高斯分布函數描述初始流場密度的非均勻分布，數值研究表明初始流場非均勻性對平面激波誘導的RM不穩定性演化影響顯著。Xiao等[8]、肖佳欣等[9]在此基礎上更加細致地分析了初始均勻流場與非均勻流場中的渦量及環量分布，揭示了非均勻流場中RM不穩定性產生及演化機理。同時，Bai等[10]通過數值模擬兩種非均勻流場中平面激波誘導的RM不穩定性現象，發現：在線性階段或弱非線性階段，界面擾動增長強烈依賴于初始流場非均勻性；而在非線性階段，該依賴性逐漸減小。

可以發現，Bai等[7,10]、Xiao等[8]、肖佳欣等[9]側重研究流場非均勻性對平面激波誘導產生的RM不穩定性演化的影響。然而，Ishizaki等[11]、Kane等[12]、Zou等[13]和Zhai等[14]分別通過理論分析、實驗研究及數值模擬手段發現，均勻流場中非平面激波沖擊密度不同的平界面時也會產生類RM不穩定性現象。實際上，在工程應用中非平面激波廣泛存在，如在RM不穩定性重要應用領域之一的ICF中，由于非均勻激光的輻照或者靶丸外表面粗糙度的影響，均會不可避免的產生非平面激波。同時，受制造工藝精度的影響，靶丸內部不同物質界面可能存在微小的擾動，并且物質的密度并非絕對均勻分布。因此，在內爆壓縮階段，非均勻流場中非平面激波沖擊擾動界面引起的不穩定性就有可能出現，進而影響靶丸的壓縮效果，甚至造成點火失敗。目前，通過實驗手段實現穩定、可靠的非平面激波生成及其參數測定仍然是一大難點，而且物質內部的密度分布很難精確測量。因此，通過數值模擬手段研究非均勻流場中非平面激波加載擾動界面的RM不穩定性演化規律，就具有重要的學術價值，這對于拓展對界面不穩定性現象的認識具有重要的幫助。

在文獻[6-10, 15-18]的基礎上，本文中，通過大渦模擬方法，數值模擬δ1=0.616 2 m、δ2=0.496 1 m非均勻流場中非平面激波加載擾動Air/SF6界面時的RM不穩定性現象。通過界面形態、擾動振幅、環量及脈動速度相關量，重點分析初始流場非均勻性對非平面激波誘導的RM不穩定性演化的影響，以期為系統地開展非平面激波誘導的RM不穩定性影響因素理論分析及實驗研究提供參考數據。

1 數值方法與計算模型

1.1 數值方法

采用可壓縮多介質黏性流動和湍流大渦模擬程序MVFT數值求解經Favre濾波后的可壓縮多介質黏性流動Navier-Stokes方程組，如下：

式中：i、j代表x、y、z的3個方向分別代表可解尺度流體密度、速度、壓強及單位質量總能量；N為流體種類數；為第s種流體體積分數，滿足為擴散系數，ν為流體運動黏性系數，Sc為Schmidt數；為牛頓流體黏性應力張量，δij為Kronecker符號為亞格子應力張量；分別為單位時空可解尺度、亞格子能量流， λl= μlcp/Prl、 λt= μtcp/Prt分別為可解尺度、亞格子熱傳導系數，μl、 μt為流體動力黏性系數，cp為比定壓熱容，為流體溫度，Prl、Prt為Prandtl數。采用Vreman亞格子模型[19]模化處理，實現控制方程組封閉。

采用理想氣體狀態方程，如下：

式中：γ代表氣體比熱比，對于空氣γ=1.40，e代表比內能。

采用算子分裂技術將方程組（1）描述的物理過程分解成無黏通量、黏性通量及熱通量3個部分進行計算。無黏通量采用多介質高精度分段拋物線方法（PPM）進行計算，兩步Lagrange-Remapping型PPM方法由以下4個步驟組成：（1）物理量分段拋物插值；（2）近似Riemann問題求解；（3）Lagrange方程組推進求解；（4）將物理量映射回靜止的歐拉網格上。在無黏通量的基礎上，采用二階空間中心差分及兩步Runge-Kutta時間推進方法求解黏性通量、熱通量及標量輸運通量，詳細步驟見文獻[20]。

1.2 計算模型

采用圖1所示計算模型，計算域為[-0.02 m,0.25 m]×[0.00 m, 0.20 m]，離散網格數為 540×400。左邊界為自由來流，右、上、下邊界為固壁。擾動air/SF6界面呈現正弦函數分布：

式中：初始平衡位置xi0=0.03 m，波長λi=0.05 m，初始振幅Ai=0.005 m。

馬赫數為1.25的非平面激波初始波陣面呈現正弦函數分布：

式中：波陣面平衡位置xs0=0.01 m，波長λs=0.05 m，初始振幅As=0.005 m，初始相位為φ。

本文中，SF6氣體初始密度呈現高斯函數分布：

式中：yc=0 m；δi(i=1, 2)為流場非均勻系數，反映了流場密度分布的非均勻程度，δi越小，流場的非均勻性越強，密度分布越不均勻；選取δ1=0.616 2 m、δ2=0.496 1 m，即在計算域上部氣體密度為 ρSF6，而下部氣體密度則分別為 0.9ρSF6、 0.85ρSF6，如圖2所示。初始時刻氣體參數見表1。

圖1 計算模型Fig.1 Simplified computational model

圖2 均勻流場及非均勻流場沿y方向的密度分布Fig.2 Density distributions along y direction in uniform and non-uniform flows

表1 Air/SF6初始參數（20 ℃, 101.325 kPa）Table 1 Properties of air and SF6 gases

2 結果與分析

討論的重點在于初始流場密度的非均勻分布對非平面激波導致的RM不穩定性演化的影響，對于非平面激波與平面激波引起的RM不穩定性的差別僅進行簡單對比，其產生機制差異更加深入的分析將在后續工作中開展。

2.1 流場可視化

圖3給出了均勻流場及δ1=0.616 2 m、δ2=0.496 1 m非均勻流場中，非平面激波同相及反相加載擾動Air/SF6界面時的流場密度云圖。同時，還給出了均勻流場中，相同初始位置的平面激波沖擊同一物質界面的演化過程作為對比。當入射激波初始沖擊界面后，流動迅速進入線性階段，界面擾動振幅逐漸增大而形狀變化較小。約0.2 ms后非線性效應逐步增強，流動進入弱非線性階段，界面上形成典型的“氣泡”及“尖釘”結構，如圖3(b)～(d)所示。約1.6 ms反射激波再次沖擊已經變形的界面，非線性效應占據主導作用，流動進入強非線性階段，界面開始出現多尺度渦結構，流動呈現湍流混合狀態，如圖3(e)～(f)所示。可以看到，非平面激波波后流場是非均勻分布的，且其波陣面形狀也隨時間發生變化，而平面激波波后則為均勻流場，如圖3(a)所示。均勻流場中，透射激波波陣面為垂直于激波運動方向的平面，“氣泡”及“尖釘”結構呈軸對稱分布；非均勻流場中，透射激波波陣面則為斜面，“氣泡”及“尖釘”結構非對稱增長，且δi越小，其對稱性越差。同時，無論非平面激波同相加載還是反相加載，非均勻流場中的界面形態在反射激波加載前與均勻流場差異明顯，如圖3(b)～(d)所示；而反射激波加載后這種差異減小，如圖3(e)～(f)所示。此外，反射激波加載前，平面激波加載引起的“氣泡”及“尖釘”結構流向尺寸小于非平面激波同相加載時的尺寸，但大于其反相加載時的數值。

圖3 密度云圖Fig.3 Simulated density contours

2.2 擾動振幅分析

圖4 界面擾動振幅演化歷史Fig.4 Perturbation amplitude evolution histories of interface

為定量描述流場非均勻性對非平面激波誘導的RM不穩定性的影響，定義沿激波運動方向、SF6氣體體積分數為1%～99%區域寬度的一半為界面擾動振幅。界面擾動振幅的演化過程可以宏觀地表征RM不穩定性中物質界面的變形及后期湍流混合的程度的變化。入射激波初始加載及反射激波加載后界面擾動振幅的演化過程如圖4所示，圖4(b)中tre=0 ms即對應反射激波加載后界面擾動振幅再次增長的時刻。從圖4(a)可以看出，由于非平面激波初始加載，Air/SF6界面擾動振幅開始增長，兩種非均勻流場中界面擾動振幅均大于均勻流場中的數值，并且δi越小，其界面擾動振幅偏離均勻流場條件下的數值越遠。約1.6 ms時刻反射激波到達并再次沖擊變形的物質界面，1.9 ms時刻（即圖4(b)中的0 ms時刻）反射激波與界面作用完成，界面擾動振幅再次增長，而且增長速度更快。反射激波加載后的湍流混合階段，均勻流場與兩種非均勻流場中界面擾動振幅增長規律性不強，但二者的差異總體上呈減小趨勢。前述規律說明，非平面激波無論是同相加載還是反相加載，在反射激波加載前階段，初始流場的非均勻性會顯著影響界面失穩的發展；而反射激波加載后階段，雖然界面擾動增長加劇了，但是流場逐步趨向于均勻分布，因而初始流場的非均勻性對界面擾動振幅增長的影響有所減小。此外，從圖中還可以看到，平面激波加載引起的界面擾動振幅大小總體上介于非平面激波同相加載與反相加載時的數值之間。

2.3 環量分析

在RM不穩定性中，渦量的生成與分布決定了界面的演化，在二維流動中，忽略質量力、體積力、黏性效應及可壓縮效應的渦量輸運方程為：

這說明流場中密度梯度與壓力梯度方向的不共線是流體微團運動過程中渦量變化的主要原因。二維流動中，渦量表達式為：

式中：u、v分別為x、y方向速度分量。渦量在整個計算域中的分布決定了環量的分布，環量是解釋界面擾動演化快慢的統計量之一，其定義為速度矢量沿封閉曲線的積分，由斯托克斯定理可變換為封閉曲線所圍成面積上渦量的積分。即：

式中：L、A分別代表計算域[-0.02 m, 0.25 m]×[0.00 m, 0.20 m]所圍成的封閉曲線及面積。

圖5給出了均勻流場及δ1=0.616 2 m、δ2=0.496 1 m非均勻流場中，非平面激波同相及反相加載時流場中正、負及總環量的變化規律。總環量反映了激波與界面相互作用過程中整個流場生成的渦量的多少，總環量為負表示整個流場中產生的正渦量小于負渦量，反之則表示正渦量大于負渦量。從圖中可以看出，均勻流場中正、負環量高度對稱增長，總環量為0；非均勻流場中，負環量主導了總環量的大小，總環量并不為0。同時，圖5還展示了均勻流場中平面激波誘導下流場中的環量變化。可以看到，由于非平面激波x方向速度分量u存在沿y方向的速度梯度，因而在入射激波沖擊界面前，非平面激波波陣面處存在渦量，從而流場中正、負環量不為0，而入射激波為平面激波時，流場中正、負環量則從0開始增長，這個特點是非平面激波與平面激波誘導的RM不穩定性最為明顯的區別。此外，平面激波加載時，流場中的正、負環量明顯小于非平面激波引起的數值。

圖5 環量演化歷史Fig.5 Evolution histories of circulations

為了進一步分析流場非均勻性對非平面激波誘導產生的RM不穩定性演化的影響，引入了非均勻流場與均勻流場中正、負環量的相對差值，表示非均勻流場中環量偏離均勻流場的程度。定義為：

圖6 正、負環量相對差值Fig.6 Relative differences of positive and negative circulations

表2 正、負環量相對差值極大值Table 2 Maximum values of relative difference for positive and negative circulations

圖7 沿x方向的脈動速度一點二階相關量分布Fig.7 Correlation distribution of fluctuating velocities along x direction

式中分別表示非均勻流場、均勻流場中的正、負環量。圖6給出了非平面激波同相和反相加載時，δ1=0.616 2 m、δ2=0.496 1 m非均勻流場中正、負環量與均勻流場中相對差值的變化規律。同時，表2列出了反射激波加載前、過渡階段及反射激波加載后兩種非均勻流場中正、負環量與均勻流場相對差值的極大值。可以看出，反射激波加載前，負環量的相對差值大于正環量的相對差值，并且δi越小，差別越明顯。這種現象產生的原因為，由于初始流場的非均勻分布，增加了物質界面處的密度梯度，導致界面區域各方向速度發生改變，使得流場中生成的負渦量大于正渦量，從而流場中負環量值大于正環量值[8]；而且，界面處密度梯度隨著δi減小而增大（見圖2），因而流場中負渦量值與正渦量值的差異相應增大，環量也呈現相同的規律。此外，δ2=0.496 1 m非均勻流場中的正、負環量相對差值大于δ1=0.616 2 m非均勻流場中的數值。在過渡階段，兩種非均勻流場與均勻流場的正、負環量相對差值急劇增大，仍然是負環量的相對差值大于正環量的相對差值。反射激波加載后，兩種非均勻流場與均勻流場中正、負環量相對差值明顯減小。這些說明，流場非均勻性引起的環量差別主要存在于反射激波加載前及過渡階段，而反射激波加載后差異則顯著減小。由于流場中渦量主要集中在變形的界面區域，環量的分布則由渦量分布決定。因此，反射激波加載后，兩種非均勻流場與均勻流場正、負環量差別的縮小就可以解釋反射激波加載后初始流場非均勻性對界面擾動振幅影響減小的現象。

2.4 脈動速度相關量分析

脈動速度相關量可以表征激波加載后流場脈動的劇烈程度，其實質是通過速度的脈動統計量來反映界面的演化。在二維流動中，脈動速度一點二階相關量表示脈動運動的平均動量輸運：

式中分別代表速度分量u、v的脈動量。圖7給出了反射激波加載前、后流場中沿x方向的脈動速度一點二階相關量分布，其中圖7(b)、(d)分別是圖7(a)、(c)中均勻流場條件下非平面激波及平面激波加載引起的脈動速度相關量的局部放大圖。由圖可以看出，反射激波加載前，非均勻流場中非平面激波導致的脈動速度一點二階相關量比均勻流場中的數值高出約兩個數量級，并且δi越小，非均勻流場中脈動速度一點二階相關量越大，如圖7(a)、(c)。反射激波加載后，非均勻流場與均勻流場中脈動速度一點二階相關量差別縮小至一個數量級，而且兩種非均勻流場中的脈動速度一點二階相關量相對大小差異也較反射激波加載前有所減小，如圖7(e)、(f)。這與前述界面擾動振幅及環量的演化規律相吻合，這從高階統計量的角度分析了初始流場非均勻性對非平面激波誘導的RM不穩定性演化發展的影響。此外，平面激波引起的流場脈動速度一點二階相關量總體上小于非平面激波誘導的速度脈動量數值，圖7(b)、(d)。

3 結 論

運用MVFT程序，數值模擬了均勻流場及兩種不同系數的非均勻流場中非平面激波同相及反相加載擾動界面時的RM不穩定性現象，并與平面激波誘導的RM不穩定性進行了簡要對比分析。主要考察了界面擾動振幅、環量及流場脈動速度的演化規律，具體結論如下。

（1）無論非平面激波同相加載還是反相加載，入射激波加載后至反射激波加載前，流動進入線性與弱非線性階段，界面的不穩定性增長對初始流場的非均勻性十分敏感；流場非均勻系數δi越小，初始流場密度分布越不均勻，非平面激波加載后，流場中脈動速度相關量及負環量數值越大，界面擾動振幅越大。反射激波加載后，非線性效應占據主導作用，流場逐步趨向于均勻分布，初始流場非均勻性對界面演化的影響逐漸減弱。

（2）非平面激波與平面激波誘導的RM不穩定性有著明顯的區別：入射激波沖擊界面前，非平面激波波后呈現非均勻流場，且波陣面形狀隨時間發生變化。此外，由于x方向速度分量u存在y方向的速度梯度，從而非平面激波波陣面區域存在渦量，該渦量與非平面激波加載物質界面時生成的渦量共同作用，使得流場中環量及速度脈動與平面激波加載時存在顯著的差異，因而呈現出與平面激波誘導的RM不穩定性不一樣的界面演化過程。