陳 璐，賴惠林，*，林傳棟，李德梅

(1. 福建師范大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，福州 350117；2. 福建師范大學(xué) 福建省分析數(shù)學(xué)及應(yīng)用重點實驗室，福州 350117；3. 福建師范大學(xué) 福建省應(yīng)用數(shù)學(xué)研究中心，福州 350117；4. 中山大學(xué) 中法核工程與技術(shù)學(xué)院，珠海 519082)

0 引 言

在自然界和工程技術(shù)領(lǐng)域中存在三種常見的流體不穩(wěn)定性：瑞利-泰勒（Rayleigh-Taylor，RT）不穩(wěn)定性、瑞奇邁爾-莫西科夫（Richtmyer-Meshkov，RM）不穩(wěn)定性和開爾文-亥姆霍茲（Kelvin-Helmholtz，KH）不穩(wěn)定性。外力場中當(dāng)重流體被輕流體加速或支撐時，兩流體界面處的擾動隨時間發(fā)展起來的一種流體力學(xué)不穩(wěn)定性現(xiàn)象，被稱為RT 不穩(wěn)定性[1-2]。RT 不穩(wěn)定性作為一種基本的流體不穩(wěn)定現(xiàn)象，在科學(xué)工程的許多領(lǐng)域都具有重要意義。例如在慣性約束聚變點火過程中，靶丸在壓縮發(fā)生內(nèi)爆時會產(chǎn)生RT 不穩(wěn)定性，其誘導(dǎo)的湍流混合會直接影響靶丸的能量增益并導(dǎo)致點火失敗[3]。這時發(fā)生的RT 不穩(wěn)定性是具有破壞作用的，我們應(yīng)當(dāng)盡可能減弱它。而對于超燃沖壓發(fā)動機來說，RT 不穩(wěn)定性的存在則會使得燃料充分混合，提高其燃燒效率。RT 不穩(wěn)定性還廣泛存在于天文領(lǐng)域，如旋轉(zhuǎn)恒星[4]、極光斑[5]、行星狀星云[6]、超新星爆炸[7]等。因此，對RT 不穩(wěn)定性的研究不僅具有重要的理論意義，還具有實際的應(yīng)用價值。

一般的，RT 不穩(wěn)定性的研究方法主要有實驗分析[8-9]、理論研究[10-13]和數(shù)值模擬[14-21]三種。然而，自然界的RT 不穩(wěn)定性現(xiàn)象相當(dāng)復(fù)雜，本身的物理系統(tǒng)具有較強的非線性，絕大多數(shù)實際問題的計算均超出了現(xiàn)有解析求解的能力范圍，理論研究需要各種簡化及線性假設(shè)，所得到的結(jié)果往往較為有限。實驗研究則常常會存在儀器昂貴、實驗危險和無法觀察短時間行為等缺點。因而對RT 不穩(wěn)定性的研究在理論求解和實驗方法上都會遇到這樣或那樣的不足。數(shù)值模擬作為一種當(dāng)今新興的研究工具已然成為重要的科學(xué)研究手段。相較于實驗與理論研究，數(shù)值模擬可通過計算機模擬得到復(fù)雜物理系統(tǒng)演化過程中的精細結(jié)構(gòu)和物理機制。對于RT 不穩(wěn)定性現(xiàn)象，研究人員通過數(shù)值模擬已取得了較為豐碩的成果。例如，Gallis 等[14]利用分子氣體動力學(xué)的直接模擬蒙特卡羅方法，數(shù)值模擬再現(xiàn)了RT 不穩(wěn)定性混合層增長的許多定性特征，并且在線性、非線性和自相似狀態(tài)下與理論和經(jīng)驗?zāi)Ｐ途哂卸可系囊恢滦?。Wei 等[15]提出了模擬二維不可壓縮 RT 不穩(wěn)定性的耦合格子玻爾茲曼模型和一種修正的平衡分布函數(shù), 并利用一個簡單的標(biāo)度討論了混合RT 不穩(wěn)定性的區(qū)域。Liang 等[16]利用格子玻爾茲曼方法研究了雷諾數(shù)（Re）對多模、互不相溶的RT 不穩(wěn)定性中的演化界面動力學(xué)和氣泡/尖釘振幅的影響。Lyubimova 等[17]在相場方法的基礎(chǔ)上，研究了限制在水平平面內(nèi)的RT 不穩(wěn)定性中的擴散與對流演化。Scott 等[18]利用基于小波的二維可壓縮直接數(shù)值模擬研究了單模RT 不穩(wěn)定性對渦旋動力學(xué)的影響。Li 等[19]以歐拉（Euler）模擬為基礎(chǔ)對雙界面的RT 不穩(wěn)定性弱非線性區(qū)域內(nèi)的界面耦合效應(yīng)進行了數(shù)值研究。結(jié)果表明，當(dāng)下界面的阿特伍德數(shù)（At）較小時，下界面的微擾增長幅度與上界面的At呈正相關(guān)；但當(dāng)下界面上的At較大時，二者呈負相關(guān)。Luo 等[20]利用二元混合流體模型研究了可壓縮性對RT 不穩(wěn)定性的影響，發(fā)現(xiàn)可壓縮性引起的初始密度分層起穩(wěn)定作用，而膨脹-壓縮效應(yīng)起失穩(wěn)作用。Livescu 等[21]利用直接數(shù)值模擬研究了在重力為零或反轉(zhuǎn)時RT 不穩(wěn)定性的演化。他們發(fā)現(xiàn)，當(dāng)加速度改變時，一些湍流量發(fā)生了顯著變化。以上列舉的是不同學(xué)者利用不同數(shù)值模擬方法對RT 不穩(wěn)定性做出的一系列研究，他們的結(jié)論豐富了我們對RT 不穩(wěn)定性的認知。

在RT 不穩(wěn)定性的數(shù)值模擬中，根據(jù)初始擾動界面結(jié)構(gòu)的設(shè)置，可以分為單模RT 不穩(wěn)定性（流體界面被一個具有單個頻率的余弦或者正弦所擾動）和多模RT 不穩(wěn)定性（界面被兩個及以上不同頻率的余弦或者正弦所擾動）。此前，多模RT 不穩(wěn)定性已被許多專家學(xué)者研究，并得到了一些有意義的結(jié)論。例如，Banerjee 等[22]結(jié)合不可壓縮Euler 方程和隱式大渦模擬方法研究初始條件對不可壓縮流體的多模RT 不穩(wěn)定性的影響，發(fā)現(xiàn)RT 混合的整體增長強烈依賴于初始條件。Burton 等[23]使用非線性大渦模擬方法對具有超高At的多模RT 不穩(wěn)定性進行研究，發(fā)現(xiàn)尖釘?shù)母叨群突旌蠈由L速率受到初始密度比的強烈影響。Zhang 等[24]對多模燒蝕可壓縮RT 不穩(wěn)定性的自相似非線性演化進行了二維和三維數(shù)值研究，發(fā)現(xiàn)氣泡發(fā)展速度受到初始條件和燒蝕速度的影響。Liang 等[16]使用改進的相場格子玻爾茲曼方法研究了Re對不可壓縮多模非混相RT 不穩(wěn)定性的影響，研究表明多模RT 不穩(wěn)定性在不同Re下表現(xiàn)出不同的界面動力學(xué)情形。Yilmaz 等[25]利用大渦模擬模擬了高At下的多模三維RT 不穩(wěn)定性，結(jié)果表明高At下RT 不穩(wěn)定性快速發(fā)展，尖釘?shù)纳L速率和速度增加，混合區(qū)域的不對稱性增大。Hamzehloo 等[26]利用相場方法研究了At、Re、表面張力和初始擾動幅值的不同組合對不可壓縮流體的多模非混相RT 不穩(wěn)定性的影響，發(fā)現(xiàn)三維多模RT 不穩(wěn)定性在初始階段表現(xiàn)出與單模RT 不穩(wěn)定性相似的指數(shù)界面增長率。Ding 等[27]利用分子動力學(xué)方法研究了強加速度下的雙模微觀RT 不穩(wěn)定性，發(fā)現(xiàn)微觀RT 不穩(wěn)定性表現(xiàn)出較弱的非線性性質(zhì)，并且雙模RT 不穩(wěn)定性在微觀尺度上的模耦合行為與宏觀尺度上的模耦合行為有明顯差異。

根據(jù)上述研究方法的不同，可以將流體不穩(wěn)定性的模擬工具分為以下三個層次：微觀分子動力學(xué)模型、介觀動理學(xué)模型和宏觀流體力學(xué)模型。宏觀流體力學(xué)模型主要是指基于Euler 方程組或N-S 方程組的各類模型，可以用于描述較大時空尺度的緩變行為。這類模型基于宏觀連續(xù)性假設(shè)，無法描述小尺度上粒子隨機運動對應(yīng)的熱漲落行為，缺乏描述微介觀結(jié)構(gòu)和快模式的能力。為了從更底層的基礎(chǔ)上深入理解流體系統(tǒng)特征機制，需要利用其它數(shù)值模型手段，比如微觀分子動力學(xué)方法[28]。雖然這種方法在納米流等問題中得到了成功應(yīng)用，但由于其對計算資源的高需求，本質(zhì)上局限于小尺度問題。為了解決物理精度與計算效率的矛盾問題，可以借助動理學(xué)理論構(gòu)建介觀層次的數(shù)值模擬方法。

近十年，基于非平衡統(tǒng)計物理的離散玻爾茲曼方法（discrete Boltzmann method，DBM）被成功發(fā)展應(yīng)用于研究復(fù)雜流場中熱力學(xué)非平衡效應(yīng)（thermodynamic non-equilibrium effect，TNE）和流體力學(xué) 非 平 衡 效 應(yīng)（hydrodynamic non-equilibrium effect，HNE）。作為一種粗?；奈锢斫７椒?，DBM 可以從以下兩方面彌補NS 模型的不足：1）適用于模擬研究具有銳利界面的流體流動，如沖擊波；2）可以描述和刻畫更底層的熱力學(xué)非平衡信息。事實上，我們通過查普曼-恩斯柯（Chapman-Enskog，CE）多尺度分析展開，可以得到連續(xù)極限假設(shè)條件下，介觀層次的分布函數(shù)演化與宏觀層次的物理量演化之間的內(nèi)在聯(lián)系。根據(jù)CE 多尺度分析，利用克努森數(shù)（Kn）中的各階項，可以構(gòu)造不同層次的描述熱力學(xué)非平衡信息的DBM。在分布函數(shù)需要滿足的動理學(xué)矩中存在守恒矩和非守恒矩兩種類型，利用守恒矩可得到流體的宏觀物理量信息，同時利用 (f?feq)的動理學(xué)矩可描述系統(tǒng)偏離熱力學(xué)平衡態(tài)的具體信息，其中feq是平衡態(tài)分布函數(shù)[29-33]。通過這些分析方法，一些以前無法提取的信息可被分層、定量地研究。

目前，DBM 已被廣泛用于模擬研究各類復(fù)雜流體流動，包括流體不穩(wěn)定性[34-43]、多相流[44-45]、反應(yīng)流[46-48]、激波[49]和爆轟[50-56]等。這里，我們重點關(guān)注并簡要介紹一下DBM 研究流體不穩(wěn)定性的成果。Lai 等[34]研究了可壓縮性對RT 不穩(wěn)定性的影響，發(fā)現(xiàn)可壓縮性在演化早期抑制了RT 不穩(wěn)定性的發(fā)生，在后期則促進了RT 不穩(wěn)定性。Li 等[35]利用DBM 模擬了可壓縮流體系統(tǒng)的多模RT 不穩(wěn)定性，探討了其演化機理。Chen 等[36]利用多松弛時間DBM 研究了黏性、熱傳導(dǎo)和普朗特數(shù)對RT 不穩(wěn)定性的影響。他們還模擬了二維RM 不穩(wěn)定性和RT 不穩(wěn)定性共存系統(tǒng)，討論了RT 與RM 不穩(wěn)定系統(tǒng)的相似和不同之處；研究了兩種不穩(wěn)定性的協(xié)作和競爭機制，并探討了重力場g和馬赫數(shù)對非平衡的影響，發(fā)現(xiàn)在重力加速度的共同作用下，RT 和RM 不穩(wěn)定共存系統(tǒng)中擾動的增長速度可能會增加[37]。Gan等[38]利用DBM 研究了黏性和熱傳導(dǎo)對KH 不穩(wěn)定性的影響，發(fā)現(xiàn)黏性效應(yīng)穩(wěn)定了KH 不穩(wěn)定性，并提高了局部和全局TNE 強度；而熱傳導(dǎo)效應(yīng)則表現(xiàn)為先抑制后增強KH 不穩(wěn)定性。Lin 等[39]研究了KH 不穩(wěn)定性的動態(tài)非平衡過程，研究表明由于物理梯度和非平衡區(qū)域的共同作用，在KH 不穩(wěn)定性的整個演化周期中，熱力學(xué)非平衡強度先增大后減小，而混合熵的增長速率則呈現(xiàn)先減小后增大最后減小的趨勢，混合自由焓的變化趨勢與混合熵的變化趨勢相反。Chen 等[40]采用多松弛時間DBM 模擬了耦合Rayleigh-Taylor-Kelvin-Helmholtz 不穩(wěn)定系統(tǒng)，并且引入形態(tài)邊界長度和TNE 強度研究了該系統(tǒng)的復(fù)雜流體結(jié)構(gòu)和動力學(xué)過程。Ye 等[41]研究了Kn效應(yīng)對可壓縮RT 不穩(wěn)定性的影響，發(fā)現(xiàn)Kn增大會抑制RT 不穩(wěn)定性的發(fā)展，但會增強全局HNE 和TNE 效應(yīng)。Chen 等[42]研究了比熱比效應(yīng)對可壓縮RT 不穩(wěn)定性的影響，發(fā)現(xiàn)由于宏觀物理量梯度與熱力學(xué)非平衡區(qū)域之間的競爭，TNE 強度先增大后減小，并隨著比熱比的減小而增大。Zhang 等[43]引入示蹤粒子作為可壓縮RT 不穩(wěn)定流動離散玻爾茲曼模擬的補充，并研究了混相雙流體系統(tǒng)界面附近RT 不穩(wěn)定性流動的精細結(jié)構(gòu)和TNE 行為，同時討論了可壓縮性和黏性對RT 混合的影響。DBM 結(jié)果還得到了分子動力學(xué)[28]、直接模擬蒙特卡羅方法[57-58]等的證實和補充。

為了進一步從動理學(xué)角度深入研究可壓縮流體RT 不穩(wěn)定性的演化規(guī)律和物理機制，本文使用Bhatnagar-Gross-Krook（BGK）DBM 模擬分析多模RT系統(tǒng)中的HNE 和TNE 的動態(tài)演化過程和宏觀表征，并通過分析流體系統(tǒng)中溫度梯度與非平衡面積占比解釋非平衡分量和全局非平衡效應(yīng)的演化規(guī)律，進而從介觀與宏觀兩方面了解多模RT 不穩(wěn)定性。

1 離散玻爾茲曼模型

圖1 離散速度示意圖Fig. 1 Schematic diagram of discrete velocities

基于以上幾種非平衡量，我們可進一步定義一種可以描述流體系統(tǒng)整體偏離平衡態(tài)程度的全局熱力學(xué)非平衡量：

另外，對于離散玻爾茲曼演化方程（2）中的時間導(dǎo)數(shù)采用一階精度的向前Euler 格式進行處理，空間導(dǎo)數(shù)采用二階精度的無波動、無自由參數(shù)的耗散有限 差 分 格 式（ nonoscillatory and nonfree-parameter dissipation finite difference scheme，NND 格式）[59]。

2 數(shù)值模擬

圖2 t 分別為0 、0.5、1.0、1.5、2.0、2.5、3.0和 3.5時的溫度輪廓圖Fig. 2 Contours of the temperature at t =0,0.5,1.0,1.5,2.0,2.5,3.0 and 3 .5 respectively

本文的模擬結(jié)果與經(jīng)典多模RT 不穩(wěn)定性的主要特征基本一致[60-61]，但同時也存在些許不同，造成這種差異的原因有以下幾點：

1）隨機擾動模態(tài)不同。在程序中，通過隨機函數(shù)在物質(zhì)界面產(chǎn)生隨機擾動模態(tài)，不同軟件平臺的隨機函數(shù)產(chǎn)生的結(jié)果往往不一致。

2）流體可壓性的影響。與不可壓模型不同，本文使用的可壓縮DBM 可以用于準(zhǔn)確預(yù)測和研究RT 不穩(wěn)定性的可壓縮效應(yīng)。

3）與Euler、N-S 等傳統(tǒng)流體力學(xué)模型不同，本文使用的DBM 不僅包含了黏性和熱傳導(dǎo)的影響，還同時包含了其他重要的非平衡效應(yīng)。

由于多模隨機初始擾動中擾動種子的不確定性，我們模擬了10 組不同隨機擾動情況下的RT 不穩(wěn)定性，并根據(jù)其統(tǒng)計特性進行分析。首先，我們分析了RT 不穩(wěn)定系統(tǒng)中系統(tǒng)溫度梯度的變化趨勢。圖3（a）

圖3 在10 組不同隨機擾動情況下全局平均溫度梯度隨時間的演化圖Fig. 3 Time evolution of the global average temperature gradients under ten initial conditions with different random perturbations

圖4 不同隨機擾動下與熱傳導(dǎo)相關(guān)的非平衡量的全局平均強度演化圖Fig. 4 Time evolution of the global non-organized energy flux with different initial random perturbations

進一步可以發(fā)現(xiàn)，全局平均TNE 強度Dˉ與宏觀物理量梯度（包括溫度、密度和速度梯度）和非平衡區(qū)域面積密切相關(guān)。數(shù)值研究表明，在RT 不穩(wěn)定性發(fā)展的早期階段，局部物理量梯度減小，局部TNE 強度減小。隨著RT 不穩(wěn)定性的發(fā)展，非平衡區(qū)域面積增加，使得全局平均TNE 強度增強，即宏觀物理梯度效應(yīng)和非平衡區(qū)域效應(yīng)是相互競爭的。進一步，我們給出了非平衡區(qū)域面積占比（Sr）的演化圖，見圖5。為了直觀了解非平衡強度的演化，圖6 展示了系統(tǒng)局部非平衡強度在不同時刻的演化云圖（從藍色到紅色表示 非 平 衡 強 度 增 加）。在 前 期（大 約 0 ＜t＜1.2），非平衡區(qū)域先增加后減小。增加的原因是：接觸界面過渡層變寬，非平衡區(qū)域增加；減小的原因是：界面處的非平衡強度會隨著擴散作用逐漸減小，這導(dǎo)致了大于閾值的非平衡區(qū)域占比的減小。之后（大約t＞1.2），非平衡區(qū)域先增加、后減小，其原因是：流體開始形成尖釘與氣泡結(jié)構(gòu)，并且隨著尖釘（氣泡）下降（上升），非平衡區(qū)域不斷伸展，Sr不斷增大，并在t=3.5左右到達峰值。隨著尖釘氣泡到達上下壁面，非平衡區(qū)域增加空間受限，流體的融合也更加充分，此時系統(tǒng)中的物理量梯度變得光滑，非平衡強度逐漸減小，使得Sr下降，最后保持在一個定值附近。同時，系統(tǒng)中出現(xiàn)的渦結(jié)構(gòu)使得流場更加復(fù)雜，并在后期呈現(xiàn)出不規(guī)則的振蕩。

圖5 不同隨機擾動下非平衡區(qū)域面積占比隨時間演化圖Fig. 5 Time evolution of the proportion of the non-equilibrium region with different initial random perturbations

圖6 t 分別為0 .02、0.5、1.0、1.5、2.0、2.5、3.0 和 3.5時的局部非平衡強度輪廓圖Fig. 6 Contours of the local non-equilibrium effects at t =0.02, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0 and 3 .5 respectively

圖7 給出了系統(tǒng)全局TNE 強度隨時間演化圖。Dˉ的演化情況與宏觀物理量梯度（包括溫度、密度和速度梯度）和非平衡區(qū)域面積均密切相關(guān)。從圖中可以看到，Dˉ 在 演化早期（大約t＜1.0）首先急劇下降，然后是較長時期的緩慢變化，數(shù)值小于0.01。隨著流體界面逐漸拉長，系統(tǒng)中出現(xiàn)了越來越多的小結(jié)構(gòu)，非平衡區(qū)域的增大也使得Dˉ呈現(xiàn)出增長的趨勢，在大概t=3.5時到達峰值。隨著系統(tǒng)物理量梯度的減弱，非平衡區(qū)域也減小，此時系統(tǒng)逐漸向平衡態(tài)發(fā)展，Dˉ逐漸減小，直至趨于穩(wěn)定狀態(tài)。

圖7 不同隨機擾動下全局平均TNE 強度隨時間演化圖Fig. 7 Evolution of the global average TNE intensity with different initial random perturbations

3 結(jié) 論

本文利用DBM 研究了可壓縮流體中的多模初始擾動的RT 不穩(wěn)定性。首先，分析了與熱通量相關(guān)的熱力學(xué)非平衡分量的演化趨勢，有兩個物理機制起主要作用：一是溫度梯度的增大使非平衡強度增大；二是兩種流體接觸面積的增大促進熱交換，從而使非平衡區(qū)域增大。兩個因素相互競爭共同影響全局非平衡強度的發(fā)展趨勢。其次，分析了非平衡面積占比的變化趨勢。界面熱擴散作用使得界面上的局部非平衡強度有減小的趨勢，同時界面的拉伸增加非平衡的區(qū)域。兩種競爭機制使得非平衡區(qū)域呈現(xiàn)出復(fù)雜的變化趨勢。最后，分析了全局平均TNE 強度，全局平均TNE 強度先增后減最后趨于穩(wěn)定。此過程也存在兩個競爭機制：非平衡區(qū)域的面積增大（減小）會增強（減弱）非平衡強度，物質(zhì)界面物理梯度的增大（減小）對全局平均熱力學(xué)非平衡強度有相同的影響，二者相互競爭，使其呈現(xiàn)出先減、后增、再減的趨勢。這些現(xiàn)象和規(guī)律有助于我們更好地理解和分析可壓縮RT 不穩(wěn)定性背后的物理機制和機理，為今后相關(guān)的流體不穩(wěn)定性研究提供介觀尺度的物理參考。