耦合界面張力的三維流體界面不穩定性的格子Boltzmann 模擬*

馬聰 劉斌 梁宏

(杭州電子科技大學理學院，杭州 310018)

采用介觀格子Boltzmann 方法模擬界面張力作用下三維流體界面的Rayleigh-Taylor (RT)不穩定性的增長過程，主要分析表面張力對流體界面動力學行為及尖釘和氣泡后期增長的影響機制.首先發現三維RT 不穩定性的發生存在臨界表面張力(σc)，其值隨著流體Atwood 數的增大而增大，且數值預測值與理論分析結果 σc (ρh-ρl)g/k2 一致.另外，隨著表面張力的增大，不穩定性演化過程中界面卷吸程度和結構復雜性逐漸減弱，系統中界面破裂形成離散液滴的數目也顯著減少.相界面的后期動力學行為也從非對稱發展轉向始終保持關于中軸線對稱.尖釘與氣泡振幅在表面張力較小時對其變化不顯著，當表面張力增大到一定值后，可以有效地抑制尖釘與氣泡振幅的增長.進一步發現，高雷諾數三維RT 不穩定性在不同表面張力下均經歷4 個不同的發展階段:線性階段、飽和速度階段、重加速和混沌混合階段.尖釘與氣泡在飽和速度階段以近似恒定的速度增長，其漸進速度的值與修正的勢流理論模型結果一致.受非線性Kelvin-Helmholtz 旋渦的剪切作用，尖釘與氣泡隨后的增長被加速，導致在重加速階段的演化速度超過勢流模型的解析解.重加速階段不能持續發展下去，尖釘與氣泡在不穩定性后期的增長速度會隨時間上下波動，這表明不穩定性的演化進入了混沌混合階段.通過數值分析，證實了三維RT 不穩定性在后期的混沌混合階段具有二次增長的規律，并且尖釘與氣泡增長率總體上隨著表面張力的增大而逐漸減少.

1 引言

Rayleigh-Taylor (RT)不穩定性是低密度流體加速高密度流體時兩者密度交界面的不穩定現象，廣泛存在于自然現象(如卷云的形成、超新星爆發、地下鹽丘和火山島的形成)和工程應用(如超燃沖壓發動機、慣性約束核聚變、氣象學、海洋運動學)中，也是多相流體界面不穩定性、湍流混合等基礎問題的理論基礎[1，2].因此，研究RT 不穩定性的發展與演化規律對于理論研究和工程實踐都具有重要的科學價值和現實意義.著名學者Rayleigh[3]和Taylor[4]最早描述了RT 不穩定性現象，提出了著名的線性穩定性理論:交界面處初始擾動隨時間以指數形式增長.后來，Lewis[5]實驗證實了線性穩定性理論可以有效地描述擾動振幅的增長達到 0 .4λ(λ是初始擾動的波長)，并發現不穩定性隨后進入非線性增長階段.2001 年，Waddell等[6]實驗研究了二維單模RT 不穩定性問題，觀察到緊接著線性階段，尖釘與氣泡振幅在非線性階段的平均增長速度接近于常數.Goncharov[7]理論分析單模RT 不穩定性的非線性增長，并給出尖釘與氣泡在非線性階段以恒定速度增長的勢流模型:

其中，ub是氣泡的速度;us是尖釘的速度;g是 重力加速度;At是Atwood 數;k是波數;Cg對于三維流動取 1，而對于二維情形則取3.受Goncharov工作[7]的啟發，Sohn[8]進一步分析了流體黏性和表面張力對氣泡漸進速度的影響，提出了包含這些因素的修正勢流模型，三維情形下可以表述為

其中ν是流體黏性，ρh表示重流體的密度.Glimm等[9]基于前追蹤方法模擬了二維單模RT 不穩定性問題，首次發現氣泡與尖釘的演化速度經過恒定速度增長后會出現再加速行為.RT 不穩定性的再加速現象在后續的實驗[10]和數值模擬[11]中被進一步證實，并被歸因于氣泡與尖釘界面處形成的Kelvin-Helmholtz 渦旋的作用.Bian 等[12]采用直接數值模擬方法研究了不同雷諾數和阿特伍德數下渦量對三維單模RT 不穩定性的再加速增長階段的影響，并確定了再加速過程與渦量之間有很強的相關性.2012 年，Ramaprabhu 等[13]數值研究了混相流體的三維單模RT 不穩定性的后期增長規律，首次發現不穩定性在高雷諾數時經歷4 個不同的發展階段，包括線性增長、飽和速度增長、再加速和混沌混合階段.同期，Wei 和Livescu[14]基于直接數值模擬方法高精度模擬了低阿特伍德數流體的二維單模RT 不穩定性的后期增長過程，同樣觀察到高雷諾數的單模RT 不穩定性經歷上述4 個發展階段，并且報道氣泡振幅在后期的混沌混合階段具有平均二次增長的規律.Hu 等[15]模擬了不同雷諾數下的二維單模RT 不穩定性問題，發現在中等雷諾數條件下，氣泡與尖釘在后期階段會出現反復的加速與減速現象.李德梅等[16]近期利用離散Boltzmann 方法研究了可壓縮流體的二維單模RT 不穩定性的非平衡效應.另外，Liang 等[17-19]針對非混相流體的單模RT 不穩定性后期演化過程也開展了一系列介尺度數值研究，分析了廣泛的雷諾數和阿特伍德數對不穩定性的后期增長階段和相界面動力學行為的影響.

上述研究豐富了人們對單模RT 不穩定性演化機制的認識，但均未考慮兩相流體間界面張力對不穩定性增長的影響.已有研究表明，表面張力會顯著影響多相流體界面輸運現象及流體界面不穩定性行為，特別是RT 不穩定性在界面張力作用下可以顯示出毛細波、收縮、破裂等獨特的界面動力學行為.鑒于此，一些學者嘗試探索表面張力對單模RT 不穩定性增長的影響，如Daly[20]采用有限差分方法模擬了界面張力對單模RT 不穩定性早期增長的影響，發現增大表面張力可以減小線性增長率;Zhang 等[21]利用多相流格子Boltzmann 方法模擬了表面張力作用下二維單模RT 不穩定性的演化過程，發現存在表面張力作用可以抑制氣泡與尖釘的增長;Young 和Ham[22]模擬了表面張力作用下單模與多模RT 不穩定性的演化機制，發現增大界面張力可以抑制混合層的增長;Matsuoka[23]采用邊界積分法模擬了耦合表面張力的可壓縮流體的二維單模RT 不穩定性的演化過程，發現存在表面張力有利于誘導界面發生夾斷行為以及抑制不穩定性湍流現象的發生;Sohn[24]與夏同軍等[25]分析了表面張力對單模RT 不穩定性的飽和速度增長階段的影響，給出了包含表面張力效應的氣泡漸進速度的解析表達式;黃皓偉等[26]數值研究了表面張力對高雷諾數的二維單模RT 不穩定性后期增長的影響，發現增大表面張力可以減弱演化過程中相界面結構的復雜程度以及有效抑制后期階段相界面破裂行為，并進一步給出了不同表面張力下氣泡與尖釘的后期增長率.

綜上所述，對耦合界面張力的單模RT 不穩定性的研究已取得了一些進展，但相關研究還不夠深入.絕大多工作[20-25]局限于不穩定性發展的中前期階段，所考慮的雷諾數較小，且所研究的物理工況均為二維情形[20-26]，而對表面張力作用下三維單模RT 不穩定性的演化后期及氣泡與尖釘增長的描述尚未報道.本文將基于介觀格子Boltzmann方法系統研究高雷諾數的三維單模RT 不穩定性的后期演化機制，主要分析界面張力對相界面動力學行為及氣泡與尖釘后期增長的影響.

2 數值模型與問題描述

格子Boltzmann 方法是基于氣體動理學理論的介觀數值方法，具有清晰的物理背景和粒子演化特性，并行計算效率高，可以從底層方便地描述流體內部與環境間的相互作用，因而在模擬多相多組分等復雜流體輸運問題上具有很大的優勢[27].根據流體間作用力的不同物理背景，現有的多相流格子Boltzmann 模型主要可以分為4 類:顏色模型、偽勢模型、自由能模型和相場模型.不同于前3 類模型，相場格子Boltzmann 模型需要求解明確的界面追蹤方程，相界面求解精度高，適用于模擬界面大拓撲變化的多相流問題[28，29].因此，本文采用相場格子Boltzmann 模型研究三維非混相Rayleigh-Taylor 不穩定性的后期增長規律.該模型利用兩套粒子分布函數fi和gi，并采用多松弛碰撞算子以提高模型處理高雷諾數流動的數值穩定性.雙分布函數的多松弛格子Boltzmann 方法的演化方程可以表示為如下的統一形式[30]:

其中，fi(x，t) 為序參數的粒子分布函數，gi(x，t) 為刻畫流場的粒子分布函數，(x，t) 和(x，t) 則為粒子分布函數所對應的平衡態分布函數，Mf和Mg表示碰撞矩陣，Sf和Sg為松弛矩陣，Fi(x，t)和Gi(x，t) 分別表示源項和外力分布函數.為了恢復正確的相場方程和不可壓流體動力學方程，平衡態分布函數(x，t) 和(x，t) 分別構造為[30]

且

其中，η是與遷移率相關的調節參數;cs為聲速;ci和ωi分別為速度空間的離散速度和權系數，取值依賴于選取的格子模型.針對三維的Cahn-Hilliard 方程，格子Boltzmann 方法至少需要7 個離散速度才能滿足恢復宏觀控制方程所約束的矩條件[31]，因此，為了提高計算效率，在序參數分布函數的演化方程(3a)中采用高效的D3Q7 格子模型，其對應的權系數ωi為ω01/4，ω1-61/8，c2/4，離散速度ci定義為

另外，針對三維流體動力學方程，在流場分布函數的演化方程(3b)中采用通用的D3Q15 格子模型，其對應的權系數為ω02/9，ω1-61/9，ω7-141/72，c2/3，且離散速度ci設定為

三維多松弛格子Boltzmann 模型的其他元素中碰撞矩陣Mf和Mg的選取可參考文獻[31，32].為了恢復正確的相界面追蹤的Cahn-Hilliard 方程，Liang 等[30]在演化方程(3a)中引入關于時間導數的局部源項以提高相界面求解的精度:

另一方面，本文考慮了外力引入格子Boltzmann方法所產生的離散效應，并采用了 He 等[33]提出的完全消除離散誤差的外力格式:

其中，G為流體系統所受的外力;Fa是為了恢復正確的動量方程而引入的額外界面力，其定義為表示兩相流體間的表面張力;μ為相場模型中的化學勢，

其中β，k為與界面張力(σ)和界面厚度(D)相關的模型參數，.通過計算粒子分布函數的零階矩和一階矩，可以得到宏觀統計量[30]

另外，流體密度ρ可由序參數φ的線性插值計算，

其中，ρl為輕流體的密度.在實際模擬中，需要采用合適的差分格式離散格子Boltzmann 模型的導數項，本文采用顯式的歐拉格式計算時間導數項:

同時利用各向同性的二階差分格式計算空間梯度和拉普拉斯算子:

其中χ表示任意一個變量.

為了研究流體不穩定性的后期演化規律，本文考慮的物理工況為z軸方向足夠長的長方體通道.如圖 1 所示，Lx，Ly，Lz分別表示管道的長度、寬度與高度，且Lx×Ly×LzW×W×16W.初始時刻，重流體位于管道的上方而輕流體置于管道的下方，并在兩相流體界面z8W處施加一個微小的擾動:

圖1 三維單模RT 不穩定性問題的示意圖Fig.1.Schematic of three-dimensional single-mode RT instability problem.

根據相場理論，序參量的初始分布設定為如下的雙曲正切函數:

其值在體相區取0 和1，而在界面區域從0 連續變化至1.為了表征RT 不穩定性的演化特性，需要引入兩個重要的無量綱參數，即雷諾數(Reynolds number，Re)和阿特伍德數(Atwood number，At)，分別定義為[33]

其中，ν表示流體黏性，重流體的密度ρh給定為 1，而輕流體的密度ρl可以根據阿特伍德數的值來確定.為了描述重力效應，對管道中輕重流體均施加一個豎直方向上的浮力:

本文著重分析表面張力對高雷諾數RT 不穩定性的后期演化特性的影響，除特別聲明，流動的雷諾數和Atwood數分別固定為Re5000，At0.1，其他物理參數給定為W100，D4，松弛矩陣Sf和Sg分別設定為

3 數值結果和討論

通過數值計算發現，當表面張力較大時，RT 不穩定性的初始擾動會發生振蕩直至相界面變為穩定的水平平面，即RT 不穩定性現象將不會發生.這表明RT 不穩定性的發生存在一個臨界表面張力，當表面張力大于該臨界值時，RT 不穩定性行為將不會出現，而當表面張力小于該臨界值時，RT 不穩定性現象將會發生.本文統計了不同流體Atwood 數下三維RT 不穩定性發生的臨界表面張力，結果如圖 2 所示.可以發現，臨界表面張力隨著Atwood 數的增大而呈現遞增規律.進一步理論分析了RT 不穩定性現象發生的臨界表面張力.根據線性穩定性理論[6，10]，RT 不穩定性的擾動振幅在初始階段的增長具有指數形式，即ha1eγt+a2e-γt，其中h表示擾動振幅，t是演化時間，a1和a2為擬合系數，γ是線性增長因子且可表示為[6]

要使RT 不穩定性能夠發生，線性增長因子中被開方數的值必須大于 0，從而可以推導出臨界表面張力的解析表達式為

其中波數k2π/W.圖2 進一步給出了不同Atwood 數下臨界表面張力的理論值，可以發現理論值基本落在數值模擬所預測的范圍內，即數值模擬結果與理論分析結果一致.在接下來的研究中，只考慮RT 不穩定性能夠發生的情形，即σ＜σc.

圖2 不同Atwood 數下三維RT 不穩定性的臨界表面張力Fig.2.Critical surface tensions of three-dimensional RT instability at various Atwood numbers.

圖3 給出了不同界面張力作用下三維單模RT 不穩定性的相界面的演化過程.可以看出，三維RT 不穩定性在演化初期顯示出相似的界面動力學特征，重流體向下運動而輕流體向上升起，即輕流體和重流體之間相互滲透，且滲透深度隨著時間的演化而逐漸增加，從而輕流體上升形成了氣泡，而重流體下落形成了尖釘.緊接著，尖釘繼續往下運動，氣泡繼續上升.與此同時，兩流體剪切作用產生的Kelvin-Helmholtz 不穩定性開始發展，其強度也隨時間逐漸增強，并開始影響相界面的動力學行為.在Kelvin-Helmholtz 不穩定性的作用下，尖釘向上卷起，形成了經典的蘑菇結構(見圖3(a)—(c)中t8 時刻與圖3(d)的t10 時刻).同時也注意到尖釘的卷起幅度在表面張力較小時差異不大，但繼續增大表面張力會減弱尖釘卷起程度，并且卷起發生的時刻也會隨之推遲.接下來觀察到三維單模RT 不穩定性在不同的表面張力作用下，表現出的顯著不同的界面動力學行為.當表面張力較小(σ10-6或 1 0-5)時，尖釘隨著時間繼續往下運動，在流體間剪切力的作用下，卷起表面變得不光滑，出現了凹凸不平的結構，并且在卷起尾端形成了4 個尺寸相同的類似于“卷發”的界面結構(見t10 時刻).隨后，不穩定性系統的非線性效應越來越劇烈，在高強度的Kelvin-Helmholtz 不穩定性的剪切作用下，產生了許多不同尺度的旋渦，進而使得流體界面變得異常不穩定，演化后期相界面發生了多層次卷起、劇烈變形、混沌破裂，形成了非常復雜的拓撲結構，并且使得流體系統中產生了大量的游離的離散液滴.另外可以發現，在低界面張力的三維RT 不穩定性的增長后期，相界面的演化圖案關于中心軸失去了對稱性.相界面的非對稱性發展在高雷諾數下混相流體的三維單模RT 不穩定性現象[13]中同樣被觀察到.當表面張力增至σ2×10-4，流體間的剪切作用減弱，流體界面卷起表面相對比較光滑，未出現凹凸不平的結構，形成的“卷發”長度也隨之減小.此后，非線性Kelvin-Helmholtz 不穩定性的強度隨時間逐漸增強，致使相界面在多個位置發生卷起甚至出現破裂行為，最終在演化后期也形成了非常復雜的拓撲結構.但不同于上述表面張力較低時的情形，高雷諾數三維單模RT 不穩定性在較大表面張力的作用下，相界面演化圖案始終保持著關于中軸線對稱的特性，并且整個演化過程中界面破裂的現象減少，系統中形成離散液滴的數量也相應地減少.特別地，將表面張力增大至σ5×10-4，在卷起尾端未觀察到明顯的“卷發”結構，雖然在演化后期仍然觀察到界面的多層次卷起行為，但界面卷起程度和界面結構的復雜性顯著減弱，系統中生成的離散液滴數目也顯著減少.另外，在整個不穩定性的演化過程中，相界面始終維持關于中軸線的對稱性.

圖3 表面張力對高雷諾數三維單模RT 不穩定性相界面演化的影響，R e=5000，A t=0.1 (a) σ=1×10-6 ;(b)σ=1×10-5 ;(c) σ=2×10-4 ;(d)σ=5×10-4Fig.3.Effect of the surface tension on the evolution of phase interface in three-dimensional single-mode RT instability with high Reynolds number，R e=5000，A t=0.1 :(a) σ=1×10-6 ;(b) σ=1×10-5 ;(c) σ=2×10-4 ;(d) σ=5×10-4 .

為了更細致地顯示表面張力對界面動力學行為的影響，給出了上述界面張力作用下三維RT 不穩定性中對角平面 (xy) 的界面演化圖案，結果如圖4 所示.可以看出，在最初階段，流體界面在不同界面張力作用下均表現出相似的行為，重流體在重力作用下向下運動形成尖釘，而密度較小的流體向上升起形成了氣泡，并受Kelvin-Helmholtz不穩定性的影響，界面在兩層位置發生卷吸現象，從而形成了兩對旋轉方向相反的旋渦.界面卷起的程度會隨著表面張力的增大而逐漸減弱，卷起發生的時間也會隨之而推遲.界面兩層卷起行為是三維單模RT 不穩定性獨有的特征，在二維單模RT 不穩定性的演化過程中并未出現.接下來，對于低表面張力(σ10-6或 1 0-5)情形，形成的兩對旋渦隨著時間演化而逐漸增大，導致在各自卷起尾端處分別形成了一對二級旋渦.在多個旋渦相互作用下，界面在多個位置發生卷起行為，形成了一些不同尺度的旋渦結構.最終，在高流體界面剪切作用下，中軸線上界面的多處位置發生卷吸、變形、混沌的破裂，形成了許多離散的小液滴.另外，進一步可以發現，對于低表面張力情形，流體界面在演化后期難以維持關于中軸線的對稱性.當表面張力較大(σ2×10-4)時，兩對一級旋渦隨時間繼續增長，伴隨著卷起部分的長度越來越長，并逐漸開始收縮，與中軸線附近的流體界面發生接觸，相比低表面張力情形，在界面卷起的尾端未出現二級旋渦的界面結構.隨后，非線性Kelvin-Helmholtz 不穩定性的影響逐漸加強，導致界面多個位置發生卷起變形，形成較復雜的結構，但相界面在整個演化過程中始終保持中軸線的對稱性.當界面張力繼續增大時，流體界面變得相對簡單與光滑，界面卷起的長度和程度也隨之減小，未觀察到劇烈拓撲變化的復雜界面結構.

圖4 表面張力對高雷諾數三維單模RT 不穩定性的對角平面(x=y 平面)的相界面演化的影響，R e=5000，A t=0.1 (a)σ=1×10-6 ;(b) σ=1×10-5 ;(c) σ=2×10-4 ;(d)σ=5×10-4Fig.4.Effect of the surface tension on the evolution of phase interface in the diagonal plane of three-dimensional single-mode RT instability with high Reynolds number，R e=5000，A t=0.1 :(a) σ=1×10-6 ;(b) σ=1×10-5 ;(c) σ=2×10-4 ;(d)σ=5×10-4.

上述研究分析了表面張力對三維RT 不穩定性的相界面動力學的影響，而尖釘與氣泡振幅及增長速度是描述RT 不穩定性演化特征的幾個重要物理量.為了進一步表征界面張力效應，定量統計了廣泛表面張力條件下尖釘與氣泡振幅、增長速度隨時間的演化規律.尖釘與氣泡的振幅定義為尖釘與氣泡的實時位置與初始位置在z方向的間距，因此在初始時刻，尖釘與氣泡振幅均為零.圖 5 給出了尖釘與氣泡振幅在不同表面張力下隨時間變化的演化曲線.可以看出，尖釘與氣泡振幅在不同表面張力下隨時間演化不斷增長，并在較小表面張力時，同一時刻可以獲得更大的尖釘與氣泡的振幅.具體而言，當表面張力在較小范圍內，尖釘與氣泡振幅的演化曲線幾乎相互重合，即界面張力對不穩定性增長的影響不再顯著.在非混相不可壓流體的RT 不穩定性的演化中，尖釘與氣泡的動力學特征理論上由單位質量的重力、黏性耗散力、界面張力之間的競爭決定[34].針對表面張力極小的情形，表面張力遠遠小于系統中重力與耗散力，因此對不穩定性中尖釘與氣泡發展的影響可以忽略.當表面張力增大到一定值后，它可以有效地抑制尖釘與氣泡振幅的增長，即尖釘與氣泡振幅隨著表面張力的增長而顯著減小.進一步模擬了不同Atwood 數下界面張力對三維單模RT 不穩定性增長的影響，同樣發現表面張力較小時對三維不穩定性中尖釘與氣泡的增長影響不大，繼續增大表面張力可以有效抑制不穩定性的發展.圖 6 描述了不同界面張力下尖釘與氣泡隨時間演化的增長速度.根據速度演化曲線，發現高雷諾數的三維單模RT 不穩定性在不同界面張力下的增長仍然可以劃分為4 個階段:線性增長階段、飽和速度階段、重加速階段和混沌混合階段.當表面張力充分小時，氣泡和尖釘的增長速度曲線在線性階段相互重合，繼續增大表面張力可以抑制尖釘與氣泡的線性增長，這符合線性穩定性理論的分析結果:線性增長因子如(23)式所示，在表面張力較小時差別不大，隨著表面張力的增大而逐漸減小.緊接著，尖釘與氣泡均以近似恒定的速度增長，這表明三維RT 不穩定性的演化進入了飽和速度增長階段.經典的勢流理論模型[7]預測了理想流體的單模RT 不穩定性中尖釘與氣泡在飽和速度階段的漸進速度.隨后，Sohn[24]基于上述勢流模型分析了流體黏性和界面張力對氣泡飽和速度增長的影響，給出了氣泡漸進速度的理論解，如(2)式所示.受Goncharov 勢能模型[7]的啟發，尖釘在飽和速度增長階段的漸進速度可以表示為

圖5 表面張力對隨時間演化的尖釘(左)和氣泡(右)振幅的影響Fig.5.Influence of surface tension on the time evolutions of spike (left) and bubble (right) amplitudes.

圖6 表面張力對隨時間演化的尖釘(左)和氣泡(右)增長速度的影響.黑色和藍色虛線分別表示修正的勢能模型所預測尖釘與氣泡漸進速度在 σ=1×10-5 σ=5×10-4 時的解析解Fig.6.Influence of surface tension on the time evolutions of spike (left) and bubble (right) growth velocities.The black and blue dotted lines represent the analytical solutions of the spike and bubble asymptotic velocities from the modified potential flow model at σ=1×10-5 and σ=5×10-4 .

根據Sohn[24]改進的勢流模型，可以發現氣泡與尖釘在界面張力σ＜1×10-5時飽和速度的解析解差別極小，因此在圖6 中僅給出σ1×10-5和σ5×10-4時尖釘與氣泡漸進速度的理論值.可以發現，本文通過格子Boltzmann 方法所預測的尖釘與氣泡漸進速度與勢流理論模型的結果一致，進一步可以發現，尖釘與氣泡在表面張力較小時更早地進入飽和速度階段，并且該階段的持續時間隨著表面張力的增大而增長.接下來，受系統中非線性Kelvin-Helmholtz 旋渦的剪切作用，使得尖釘和氣泡的增長速度超過勢流模型的理論值，這預示著不穩定性的演化進入了重加速階段.重加速階段是由遠離初始中心平面的輕流體與尖釘的界面處產生的渦旋向氣泡尖端傳播，導致氣泡重新加速，其發生必須滿足兩個條件:(a)渦旋需要在豎直方向上比氣泡更快地移動;(b)渦旋的結構必須要保持足夠長時間.結合相界面動力學圖案，可以發現對所有表面張力情形，高雷諾數三維RT 不穩定性均可以產生保留時間足夠長的旋渦，導致出現重加速階段，但由于旋渦的強度隨著界面張力的增大而減小，表面張力較小的RT 不穩定性可以更早地進入重加速階段.重加速階段不能持續地發展下去，從圖 6 可以發現，對所有表面張力情形，尖釘與氣泡在演化后期的增長速度隨時間上下波動，其值反復地加速與減速，這表明RT 不穩定性的發展最終進入了混沌混合階段.已有的高精度直接數值模擬研究顯示單模RT 不穩定性在高雷諾數條件下會出現后期的混沌混合增長階段，在該階段尖釘與氣泡的振幅呈現平均二次增長的規律，即hs，bαs，bgAtt2，其中αs，b表示尖釘與氣泡的后期增長率.αs，b反映了不穩定性在混沌混合階段的發展規律，因此是RT 不穩定性后期研究中最為重要的物理統計量.為了計算不穩定性湍流混合階段的增長率，一些學者基于振幅與演化時間的二次關系提出多種不同的測量方法，Liang 等[19]近期在二維單模RT 不穩定性的后期研究中詳細地比較了這些測量方法，指出 Olson 等[35]提出的計算增長率的方法具有良好的收斂性與一致性，并被用于本文測量不同表面張力下三維單模RT 不穩定性的后期增長率.Olson 等[35]對演化后期的振幅與時間的二次關系式兩邊開方得到，因此通過與演化時間t的線性擬合可以獲得擬合直線的斜率，進而通過計算可以得到尖釘與氣泡的后期增長率αs，b.圖 7 給出了不同表面張力條件下尖釘與氣泡振幅的開方與的關系曲線以及在后期相應的線性擬合直線.從圖 7 可以發現，兩者在不穩定性的演化后期確實成線性關系，這也驗證了三維單模RT 不穩定性在后期混沌混合階段具有二次增長的規律.進一步通過線性擬合，可以獲得尖釘增長率在表面張力σ1×10-7，1×10-6，1×10-5，2×10-4，5×10-4時分別為0.1293，0.1304，0.1348，0.1106，0.1069，而對應的氣泡增長率依次為0.1286，0.1268，0.1275，0.0976，0.0691.可以發現，對于相同的界面張力，尖釘的后期增長率大于氣泡的增長率，即尖釘比氣泡運動得更快.另外，當界面張力較小時，界面張力對尖釘與氣泡增長率的影響不再顯著，而當界面張力增加至一定值時，尖釘與氣泡增長率則隨著界面張力的繼續增大而減小，即界面張力會抑制尖釘與氣泡的后期增長.

圖7 不同界面張力下三維RT 不穩定性中尖釘(左)與氣泡(右)振幅的開方 與演化時間 的關系曲線，實線表示在后期混沌混合階段的線性擬合結果Fig.7.Relationsbetweenthesquare roots of the spike and bubble amplitudes and the time in the three-dimension-al RT instability with different surface tensions.The solid lines represent the linear fitting results at the late-time chaotic mixing stage.

4 結論

本文采用多相流的介觀格子Boltzmann 方法模擬了三維非混相、不可壓流體的RT 不穩定性問題，系統分析表面張力對不穩定性中界面動態過程和尖釘與氣泡定量增長的影響.通過數值模擬發現，三維RT 不穩定性的發生存在著一個臨界表面張力，當表面張力小于該臨界值時，RT 不穩定性現象才會發生.根據線性穩定性理論建立了臨界表面張力的解析表達式，發現臨界表面張力隨著流體Atwood 數的增大而增大，且數值預測的臨界值與理論分析結果一致.另外發現，界面張力對流體界面動力學起著穩定因子的作用，即增大表面張力可以抑制界面的卷吸與破裂行為，相界面的后期演化也從混沌的非對稱發展轉化為關于中軸線的對稱發展.還定量統計了表面張力對尖釘與氣泡振幅、演化速度和無量綱加速度的影響.在表面張力較小時，界面張力對尖釘與氣泡增長的影響不再顯著，而繼續增大界面張力則可以抑制尖釘與氣泡的增長.根據速度演化曲線，高雷諾的三維RT 不穩定性在所有界面張力情形下均呈現線性增長、飽和速度增長、重加速、混沌混合4 個不同的發展階段.尖釘與氣泡在飽和速度階段的漸進速度與包含界面張力效應的勢流理論模型相符合.進一步發現，隨著界面張力的增大，由于系統中形成旋渦的數目與強度逐漸減弱，導致不穩定性進入飽和速度增長和重加速階段的時間隨之推遲.最后，通過尖釘與氣泡振幅的開方與時間的線性擬合，從數值上證明了三維單模RT 不穩定性在后期的混沌混合階段符合二次增長的規律，并且發現尖釘與氣泡后期增長率與表面張力呈現出遞減關系.