多模態瑞利-泰勒不穩定性的離散玻爾茲曼數值研究

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

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

0 引 言

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

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

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

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

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

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

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

1 離散玻爾茲曼模型

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

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

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

2 數值模擬

圖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

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

1）隨機擾動模態不同。在程序中，通過隨機函數在物質界面產生隨機擾動模態，不同軟件平臺的隨機函數產生的結果往往不一致。

2）流體可壓性的影響。與不可壓模型不同，本文使用的可壓縮DBM 可以用于準確預測和研究RT 不穩定性的可壓縮效應。

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

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

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

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

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

圖5 不同隨機擾動下非平衡區域面積占比隨時間演化圖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 給出了系統全局TNE 強度隨時間演化圖。Dˉ的演化情況與宏觀物理量梯度（包括溫度、密度和速度梯度）和非平衡區域面積均密切相關。從圖中可以看到，Dˉ 在 演化早期（大約t＜1.0）首先急劇下降，然后是較長時期的緩慢變化，數值小于0.01。隨著流體界面逐漸拉長，系統中出現了越來越多的小結構，非平衡區域的增大也使得Dˉ呈現出增長的趨勢，在大概t=3.5時到達峰值。隨著系統物理量梯度的減弱，非平衡區域也減小，此時系統逐漸向平衡態發展，Dˉ逐漸減小，直至趨于穩定狀態。

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

3 結 論

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