可壓流體Rayleigh-Taylor不穩定性的離散Boltzmann模擬?

李德梅賴惠林許愛國張廣財林傳棟4)甘延標5)

1)(福建師范大學數學與信息學院,福建省分析數學及應用重點實驗室,福州 350117)

2)(北京應用物理與計算數學研究所,計算物理國家重點實驗室,北京 100088)

3)(北京大學,應用物理與技術研究中心,高能量密度物理數值模擬教育部重點實驗室,北京 100871)

4)(清華大學能源與動力工程系,燃燒能源中心,北京 100084)

5)(北華航天工業學院,廊坊 065000)

(2017年9月4日收到;2018年1月29日收到修改稿)

1 引 言

當低密度流體支撐或推動較高密度流體時,即重力加速度或慣性加速度由重密度流體指向輕密度流體時,如果流體之間的界面存在擾動,那么界面的擾動幅度將會增長,該物理現象稱為瑞利泰勒(Rayleigh-Taylor,RT)不穩定性.這種不穩定性最早由Rayleigh[1]和Lamb[2]在某種程度上提及,直到1950年,Taylor明確指出不穩定性現象[3].因此,該現象也稱為RT不穩定性或者Rayleigh-Lamb-Taylor不穩定性.由于RT不穩定性現象在慣性約束聚變[4?6]、超新星爆炸[7]、核反應堆[8]等領域中起著重要的作用,因此在過去幾十年里,人們采用各種解析方法和數值方法對其進行研究,包括分子動力學[9]、直接數值模擬[10]、大渦模擬方法[11]等.這些研究對理解RT不穩定性現象的動力學機制提供了許多有用的信息.

作為Boltzmann方程的特殊離散形式,格子Boltzmann方法(lattice Boltzmann method,LBM)在各種復雜流體的研究中取得了巨大的成功[12].LBM在RT不穩定性問題的研究中發展了兩類模型:不可壓LBM[13?15]和可壓LBM[16].這些模型的基本思想上是把LBM看作Navier-Stokes(NS)方程的求解器,能夠模擬得到NS方程一致的結果.近年來,許愛國課題組[17?25]已將LBM發展成為能夠同時描述流動和熱動非平衡效應的離散Boltzmann方法 (discrete Boltzmann method,DBM).在2012年，許愛國等[17]提出構建DBM.DBM與LBM最主要的差異在于:作為偏微分方程解法器的LBM必須忠誠于原始物理模型,而作為流體系統動理學模型的DBM必須具有超越原始物理模型的部分功能;LBM所依賴的演化方程和“矩關系”可以根據算法設計的要求人為構造,即可以沒有物理對應,而DBM所依賴的演化方程和“矩關系”只能是Boltzmann方程及其動理學矩關系,必須與非平衡統計物理學基本理論自洽[18].例如DBM所提供的非平衡行為特征能夠恢復真實分布函數的主要特征[26]、區分不同類型的界面[27]、區分相分離過程的不同階段[18,21],所提供的沖擊波精細物理結果與分子動力學數值模擬結果相互印證,相互補充[28].本文在甘延標等[29]提出的DBM模型的基礎上,進一步驗證了含外力項的DBM模型.通過數值模擬Riemann問題和熱Couette流等問題驗證了DBM的有效性.使用該模型,本文模擬了可壓流體系統多模初始擾動的RT不穩定性現象,能夠得到RT不穩定性的基本物理圖像以及相伴隨的熱動非平衡效應規律,得出一些相關物理解釋.

2 離散Boltzmann模型

考慮含外力項的Bhatnagar-Gross-Krook(BGK)碰撞的Boltzmann方程為

其中fi(r,vi,t)是離散分布函數,r是空間變量,t是時間;vi是離散速度,i=1,2,···,N是離散速度序號;u是宏觀流速;a是加速度;τ是動理學松弛時間;是 Maxwell分布函數的離散化形式;其中Maxwell分布函數的形式如下:

其中,D為空間維數(本文考慮D=2的情形);n是除了平動自由度之外的額外自由度數目;η是自由參數;ρ,T和u分別是密度、溫度和流速.這里考慮的包含外力的方程是不包含外力情形下的拓展,可以處理更加普遍的物理情形,比如重力場存在下的流體不穩定性問題、分子間相互作用下的多相流問題、電場力、磁場力存在下的等離子體輸運問題;其中的外力項使用了f～feq的近似條件,因而該DBM只適用于系統偏離平衡不遠的情形.

前三個方程代表質量守恒、能量守恒和動量守恒.

借助 Chapman-Enskog多尺度分析,可以從離散Boltzmann方程(1)得到NS方程層次的宏觀流體力學方程.首先對密度分布函數、時間導數、空間導數和外力項進行如下多尺度展開:

其中ε?1是一個無量綱小量,正比于克努森數(Knudsen number,Kn)Kn=l/L,l是分子平均自由程或者平均分子間距,L是宏觀上關心的特征尺度.

將方程(10)代入方程(1)中,可以得到一系列關于ε的各階等式:

即分布函數的非平衡部分對宏觀物理量沒有貢獻.

經過一系列代數運算,可以得到可壓NS方程:

分別是壓強和總內能;

為動力黏性系數;

為熱傳導系數.

本文選取如下二維十六速度(D2V16)的離散速度模型: 其中,當i=1,2,···,4 時,ηi=η0,當i=5,6,···,16 時,ηi=0.

DBM擺脫了空間離散化和時間離散化之間的綁定,使得粒子速度可以靈活選擇,并且可以在離散Boltzmann方程的求解中方便地引入多種差分格式.

DBM被認為是Boltzmann方程的特殊離散形式,自然繼承了Boltzmann方程可以用來描述非平衡效應的屬性.在7個動力學矩關系(3)—(9)式中,只有前面3個動力學矩關系(質量、動能和能量的定義),可以被fi取代,而后面的4個動力學矩關系,如果用fi取代則兩側值會產生偏差.這個偏差從物理上來看是描述系統狀態偏離熱力學平衡所引起的宏觀效應,可用于描述系統狀態偏離熱力學平衡的程度[22].本文考慮扣除宏觀流動的微觀粒子熱漲落特征的熱動非平衡效應,對應的中心矩定義如下:

3 數值模擬與驗證

本節通過一維Riemann問題:Sod激波管、沖擊波碰撞和熱Coutte流問題的解析解和數值解的符合程度來驗證DBM的有效性.計算動理學方程(1)時,時間導數采用一階向前差分,空間格式采用無波動無自由參數的耗散(non-oscillatory,containing no free parameters and dissipative,NND)格式[32].事實上,NND格式是二階迎風格式、一階迎風格式、中心差分格式的混合格式,該格式針對激波上下游采用不同的混合格式,其總變差(total variation diminishing,TVD)是減小的,空間上具有實質的二階精度高分辨率,捕捉激波能力較強,可以很好地分辨間斷.

3.1 Sod激波管問題

Sod激波管問題.計算區域 [?1,1],流場的左半部分和右半部分分別給定如下的初始條件:

其中“L”和“R”分別代表遠離間斷界面左右兩側的宏觀量初始值.計算網格為Nx×Ny=2000×2,空間步長為 ?x=?y=0.001,時間步長選取為?t=10?5.其他模型參數選取為τ=10?5,n=3,c=1.0和η=10.0.y方向采用周期邊界條件,對于x方向,左邊界設置為

其中?1和0表示左邊的虛擬點.此類邊界條件指定系統在邊界處一直處于平衡態,即邊界處的宏觀量為

方程(27)和(28)也被稱作微觀和宏觀邊界條件,兩者是互相對應的.

同樣,右邊的微觀邊界設置如下:

則對應的宏觀邊界為

為驗證網格無關性,先固定其他模型參數,x方向采用三種不同的網格數:Nx=1000,2000,4000,模擬結果見圖1.可見,三種不同空間分辨率都能夠清晰捕捉激波、接觸間斷和稀疏波.采用Nx=2000的模擬結果與采用Nx=4000的模擬結果區別不大.為了更好地展示該物理問題不同物理量的非線性間斷結構,圖2給出DBM數值解與解析解在t=0.2的對比圖,圖中圓圈為DBM 數值解,直線為精確解.結果顯示,DBM數值解與解析解符合較好,驗證了模型的準確性和健壯性.

圖1 不同網格數下t=0.2時刻一維Sod激波管密度剖面的DBM數值解與解析解對比Fig.1. Comparisons between DBM results with dif f erent grids and the exact solution for the onedimensional Sod problem,at t=0.2.

圖2 t=0.2時刻一維Sod激波管的密度、壓力、速度和溫度剖面的DBM數值解與解析解的對比Fig.2.Comparisons between DBM results and the exact solutions for the one-dimensional Sod problem at t=0.2.

3.2 兩個強激波碰撞問題

為了充分驗證模型,考慮沖擊波碰撞問題,該問題涉及兩個強激波的碰撞,其初始條件為:

該問題的精確解包含了一個緩慢向右傳播的左激波、向右的接觸界面和一個左行激波.其中,左激波向右傳播很慢給數值方法帶來額外的困難,對模型的穩定性和魯棒性要求較高.

數值模擬中,選取參數為：網格參數為Nx×Ny=2000×2,?x= ?y=0.003,時間步長為?t=10?5.其他參數選取為τ=2×10?5,n=3,c=8.0和η=40.0.圖3給出了t=0.08時刻γ=1.4的密度、壓力、速度和溫度剖面的DBM數值解與解析解的對比.對比結果表明,DBM數值解與解析解符合較好,進一步說明DBM模型具有較好的穩定性和魯棒性.

3.3 熱Coutte流問題

作為經典熱傳導問題,熱Coutte流能夠用來檢測DBM模擬流體黏性熱傳導問題.該問題描述如下:考慮介于兩個無限長平行板之間的黏性流體,平板之間距離為H.初始條件為(ρ,u,v,T)|t=0=(1.0,0,0,1.0).當t>0 時,溫度為T0的上板以速度u0=0.8移動,溫度為T0的下板保持靜止不動.

網格參數選取為Nx×Ny=1×200,空間步長為?x=?y=2×10?3,其他參數選取為:n=3,τ=10?3,c=1.0,η=10.0,?t=10?5.x方向采用周期邊界條件,y方向采用非平衡外推格式[33].

x方向速度的解析解為

圖4給出了DBM數值解與解析解在不同時刻的對比圖,兩者十分符合,表明DBM能夠精確計算黏性耗散下的流體問題.計算結果與NS模型得到的結果一致.

當系統達到穩態時,沿y方向溫度場的理論解為

其中cp=γ/(γ?1).圖5展示了不同γ對應的DBM數值解與解析解在穩態時的對比圖.數值解與解析解符合較好,表明DBM能夠精確模擬不同熱傳導情形下的流體問題.

圖3 t=0.08時刻兩個強激波碰撞問題的密度、壓力、速度和溫度剖面的DBM數值解與解析解對比Fig.3.Comparisons between DBM results and the exact solutions for collision of two strong shocks problem at t=0.08.

圖4 γ=1.4時熱 Couette流在不同時刻速度剖面的DBM數值解與解析解對比Fig.4.Comparisons between DBM results and the exact solutions for the velocity profiles in thermal Couette flow for the case with γ=1.4 at various times.

圖5 不同γ值下熱Couette流的穩態溫度剖面的DBM數值解與解析解對比Fig.5.Comparisons between DBM results and the exact solutions for the temperature profiles in steady thermal Couette flow for various values of γ.

4 可壓流體RT不穩定性數值模擬

對于RT不穩定性的數值模擬,以往模型主要采用等溫不可壓模型,即上下密度是常數而溫度始終不變的情形,而實際系統往往是可壓的且溫度是變化的.本文考慮單介質流體的可壓非等溫情形,即溫度自適應情形.該流體系統由上下兩部分組成,上下溫度不同,系統密度滿足力學平衡條件呈指數分布[16?27].例如,考慮上流體是冷空氣下流體是熱空氣.當中間界面處沒有發生擾動,則系統只有熱擴散作用,界面始終處在中間位置.當中間界面出現小擾動之后,由于重力的作用,擾動會隨著時間的演變而慢慢放大,形成“氣泡-尖釘”結構,而后出現典型的“蘑菇頭”形狀,即RT不穩定性發生.在數值模擬過程中,邊界影響比較大,本文采用如下邊界條件:上下邊界采用絕熱、無滑移邊界條件;左右采用周期邊界條件.模型從最簡單的理想氣體狀態方程出發,暫時忽略表面張力的影響.

圖6 多模RT不穩定性在不同時刻的密度演化圖:t=0,0.5,1.0,1.5,1.8,2.0,2.5,3.0Fig.6.Density evolution of Rayleigh-Taylor instability from a multiple mode perturbation at dif f erent times:t=0,0.5,1.0,1.5,1.8,2.0,2.5,3.0.

本文考慮二維區域 [?d/2,d/2]×[?2d,2d],系統處于重力加速度為常數的重力場下,界面的初始擾動滿足

其中kn=2nπ/Lx,an,bn是 0—1之間均勻分布的隨機數.上下部分流體的溫度不同,每部分流體的密度分布滿足如下靜力學平衡條件:

所以系統的不穩定性初始條件滿足:

其中,p0是上部分流體頂部的初始壓強,Tu和Tb代表上下部分流體的初始溫度.在這種條件下,界面處的壓強滿足

其中ρu和ρb是上下部分流體臨近界面兩側網格處的密度,則界面處初始Atwood數可以定義為[16]

在數值模擬中,計算區域為512×512的均為網格,空間步長為?x=?y=0.001,頂部初始壓強為p0=1.0,時間步長為 ?t=1×10?5,松弛因子為τ=1×10?5,上部分溫度為Tu=1.0,下部分溫度為Tb=4.0,因此,初始At=0.6.其他參數為c=1.3,η=15,n=3,ax=0.0,ay=?g=?1.0.

圖6展示了RT不穩定性的密度分布隨時間變化的時空演化圖,可以看出,初始階段,熱擴散作用迅速抹平了間斷界面,產生有限寬度的過渡層,降低了界面處局部At數.經過短暫的線性階段,RT不穩定性進入了非線性階段.在重力場的作用下,隨著時間的發展,重流體下降,輕流體上升,又由于重流體相對較“硬”,輕流體相對較“軟”,因而呈現典型的“氣泡”和“尖釘”的界面結構.之所以形成這種結構,是因為當密度較大時,慣性力較大,較難改變速度,從而向上的擾動形成較平的“氣泡”結構,向下的擾動形成較尖銳的“尖釘”結構.后期由于界面切向速度差變大(即KH不穩定性慢慢起作用),“尖頂”尾部翻轉起來,形成“蘑菇頭”形狀.由于熱擴散和黏性作用,“蘑菇頭”尾部漸漸模糊且變狹長.事實上,一開始(t=0.5之前)演化較慢,且界面整體下移,這是由于一開始熱傳導起主導作用,在界面附近的上下流體交換內能,上流體吸收熱量,體積膨脹,界面附近的上流體下移,而下流體釋放熱量,體積縮小,界面附近的下流體下移.同時,初始多模互相競爭合并,模式慢慢變少,界面被“抹平”;而后(t=0.5)之后演化加速,界面演化變成重力主導,上下流體開始以交換重力勢能為主,呈現非線性演化階段.后期兩流體在界面附近相互滲透,相互混合,進入湍流混合階段.

圖7展示了在不同初始多模擾動下總平均熱動非平衡效應的演化情形.由于初始條件處于熱動非平衡,系統有趨于熱動平衡態趨勢,D?有下降的趨勢.而后,隨著模式的耦合以及混合層厚度不斷增加,界面越來越復雜,系統偏離熱動平衡態的演化以線性形式增長.而后,在t?=0.7后系統趨向平衡態,t?=1.2后系統又慢慢遠離平衡態,這是因為系統重力勢能和壓縮能得到釋放,部分轉化為動能,促進了RT不穩定性的發展,界面越來越復雜,非平衡模式越來越豐富.

圖7 不同初始多模擾動下RT不穩定性演化引起的總平均熱動非平衡效應隨時間的演化Fig.7.The time evolution of the global average TNE strength due to Rayleigh-Taylor instability with different multi-mode initial conditions.

5 結 論

應用含外力項的DBM數值模擬研究可壓流體多模初始擾動的RT不穩定性問題.Chapman-Enskog多尺度分析表明該模型在連續極限可恢復到Navier-Stokes方程.模型通過了熱Coutte流問題和三個一維Riemann問題的檢測,表明模型能夠精確模擬黏性耗散和熱傳導以及復雜激波之間的相互作用.采用DBM對多模、可壓、具有間斷界面的多模初始擾動RT不穩定性進行數值模擬.結果表明,在RT不穩定性發展的初期由于多模的設置,界面處的黏性和熱傳導效應突出,這些耗散效應會“抹平”界面,多模之間相互競爭和吸收,形成較少的主導模式;在這一階段系統內沒有形成明顯的“氣泡”和“尖釘”結構.在RT不穩定性的中后期,由于模式的合并導致界處的耗散效應減弱,重力占主導地位,擾動界面逐漸變形、長大,形成典型的“氣泡-尖釘”結構,即出現典型的“蘑菇頭”形狀,而后進入湍流混合階段.這些現象與經典的實驗結果一致.同時給出系統整體非平衡程度隨時間發展的演化情況,一開始系統先趨于平衡態,這是由于系統處于調整階段,從多模初始界面擾動調整到本征模階段;而后系統以線性形式偏離平衡態,這是由于系統界面被抹平,壓縮能部分轉化為內能;然后系統又趨于平衡態,這是由于模式的耦合與擾動界面進一步被“抹平”,系統處于相對穩定狀態;最后系統越來越遠離平衡態,此時是由于系統輕重流體的重力勢能相互轉換,系統的壓縮能進一步被釋放出來,系統動能進一步增加所致.在最近的一系列學術報告中,許愛國等[34?37]進一步給出了非平衡程度更深、超越Navier-Stokes描述能力的復雜流動系統的DBM建模思路.

[1]Rayleigh L 1882Proc.London Math.Soc.s1-14 170

[2]Lamb H 1932Hydrodynamics(6th Ed.)(London:Cambridge University press)p501

[3]Taylor G 1950Proc.R.Soc.London A201 192

[4]Betti R,Goncharov V,McCrory R,Verdon C 1998Phys.Plasmas(1994–present)5 1446

[5]Wang L F,Ye W H,Wu J F,Liu J,Zhang W Y,He X T 2016Phys.Plasmas23 052713

[6]Wang L F,Ye W H,He X T,Wu J F,Fan Z F,Xue C,Guo H Y,Miao W Y,Yuan Y T,Dong J Q,Jia G,Zhang J,Li Y J,Liu J,Wang L M,Ding Y K,Zhang W Y 2017Sci.China:Phys.Mech.Astron.60 055201

[7]Cabot W,Cook A 2006Nat.Phys.2 562

[8]Berthoud G 2000Annu.Rev.Fluid Mech.32 573

[9]Barber J L,Kadau K,Germann T C,Alder B J 2008Eur.Phys.J.B64 271

[10]Celani A,Mazzino A,Vozella L 2006Phys.Rev.L.96 134504

[11]Moin P 1991Comput.Meth.Appl.Mech.Eng.87 329

[12]Succi S 2001The Lattice Boltzmann Equation for Fluid Dynamics and Beyond(New York:Oxford University Press)pp179–255

[13]He X Y,Chen S Y,Zhang R Y 1999J.Comput.Phys.152 642

[14]Li Q,Luo K H,Gao Y J,He Y L 2012Phys.Rev.E85 026704

[15]Liu G J,Guo Z L 2013Int.J.Numer.Method H.23 176

[16]Scagliarini A,Biferale L,Sbragaglia M,Sugiyama K,Toschi F 2010Phys.Fluids22 055101

[17]Xu A G,Zhang G C,Gan Y B,Chen F,Yu X J 2012Front.Phys.7 582

[18]Xu A G,Zhang G C,Gan Y B 2016Mech.Eng.38 361(in Chinese)[許愛國,張廣財,甘延標2016力學與實踐38 361]

[19]Gan Y B,Xu A G,Zhang G C,Yu X J,Li Y J 2008Physica A387 1721

[20]Gan Y B,Xu A G,Zhang G C,Li Y J 2011Phys.Rev.E83 056704

[21]Gan Y B,Xu A G,Zhang G C,Li Y J,Li H 2011Phys.Rev.E84 046715

[22]Yan B,Xu A G,Zhang G C,Ying Y J,Li H 2013Front.Phys.8 94

[23]Xu A G,Zhang G C,Li Y J,Li H 2014Prog.Phys.34 136(in Chinese)[許愛國,張廣財,李英駿,李華2014物理學進展34 136]

[24]Xu A G,Zhang G C,Ying Y J 2015Acta Phys.Sin.64 184701(in Chinese)[許愛國,張廣財,應陽君2015物理學報64 184701]

[25]Xu A G,Zhang G C,Ying Y J,Wang C 2016Sci.China:Phys.Mech.Astron.59 650501

[26]Lin C D,Xu A G,Zhang G C,Li Y J,Succi S 2014Phys.Rev.E89 013307

[27]Lai H L,Xu A G,Zhang G C,Gan Y B,Ying Y J,Succi S 2016Phys.Rev.E94 023106

[28]Liu H,Kang W,Zhang Q,Zhang Y,Duan H L,He X T 2016Front.Phys.11 115206

[29]Gan Y B,Xu A G,Zhang G C,Yang Y 2013Europhys.Lett.103 24003

[30]Gan Y B,Xu A G,Zhang G C,Succi S 2015Soft Matter11 5336

[31]Watari M,Tsutahara M 2004Phys.Rev.E70 016703

[32]Zhang H X 1988Acta Aerodyn.Sin.6 43(in Chinese)[張涵信1988空氣動力學學報6 43]

[33]Guo Z L,Zheng C G,Shi B C 2002Phys.Fluids14 2007

[34]Xu A G,Zhang G C 2016The 9th National Conference on Fluid MechanicsNanjing,China Oct.20–23,2016(in Chinese)[許愛國,張廣財2016第九屆全國流體力學學術會議,南京,2016年10月20—23日]

[35]Xu A G,Zhang G C 2016Special Academic Report of Electromechanical College of Nanjing Forestry UniversityNanjing,China,Oct.25,2016(in Chinese)[許愛國,張廣財2016南京林業大學機電學院專題學術報告,中國南京,2016年10月25日]

[36]Xu A G,Zhang G C 2016Academic Report on Physics Department of Renmin University of ChinaBeijing,China,Nov.23,2016(in Chinese)[許愛國,張廣財 2016中國人民大學物理系專題學術報告,中國北京,2016年11月23日]

[37]Xu A G,Zhang G C 2016The 4th Academic Seminar of LBM and Its ApplicationsBeijing,China,Nov.26,2016(in Chinese)[許愛國,張廣財 2016第四屆LBM及其應用學術研討會,中國北京,2016年11月26日]