李德梅賴惠林許愛國(guó)張廣財(cái)林傳棟4)甘延標(biāo)5)

1)(福建師范大學(xué)數(shù)學(xué)與信息學(xué)院,福建省分析數(shù)學(xué)及應(yīng)用重點(diǎn)實(shí)驗(yàn)室,福州 350117)

2)(北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所,計(jì)算物理國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100088)

3)(北京大學(xué),應(yīng)用物理與技術(shù)研究中心,高能量密度物理數(shù)值模擬教育部重點(diǎn)實(shí)驗(yàn)室,北京 100871)

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

5)(北華航天工業(yè)學(xué)院,廊坊 065000)

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

1 引 言

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

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

2 離散Boltzmann模型

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

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

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

前三個(gè)方程代表質(zhì)量守恒、能量守恒和動(dòng)量守恒.

借助 Chapman-Enskog多尺度分析,可以從離散Boltzmann方程(1)得到NS方程層次的宏觀流體力學(xué)方程.首先對(duì)密度分布函數(shù)、時(shí)間導(dǎo)數(shù)、空間導(dǎo)數(shù)和外力項(xiàng)進(jìn)行如下多尺度展開:

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

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

即分布函數(shù)的非平衡部分對(duì)宏觀物理量沒有貢獻(xiàn).

經(jīng)過一系列代數(shù)運(yùn)算,可以得到可壓NS方程:

分別是壓強(qiáng)和總內(nèi)能;

為動(dòng)力黏性系數(shù);

為熱傳導(dǎo)系數(shù).

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

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

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

3 數(shù)值模擬與驗(yàn)證

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

3.1 Sod激波管問題

Sod激波管問題.計(jì)算區(qū)域 [?1,1],流場(chǎng)的左半部分和右半部分分別給定如下的初始條件:

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

其中?1和0表示左邊的虛擬點(diǎn).此類邊界條件指定系統(tǒng)在邊界處一直處于平衡態(tài),即邊界處的宏觀量為

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

同樣,右邊的微觀邊界設(shè)置如下:

則對(duì)應(yīng)的宏觀邊界為

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

圖1 不同網(wǎng)格數(shù)下t=0.2時(shí)刻一維Sod激波管密度剖面的DBM數(shù)值解與解析解對(duì)比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時(shí)刻一維Sod激波管的密度、壓力、速度和溫度剖面的DBM數(shù)值解與解析解的對(duì)比Fig.2.Comparisons between DBM results and the exact solutions for the one-dimensional Sod problem at t=0.2.

3.2 兩個(gè)強(qiáng)激波碰撞問題

為了充分驗(yàn)證模型,考慮沖擊波碰撞問題,該問題涉及兩個(gè)強(qiáng)激波的碰撞,其初始條件為:

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

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

3.3 熱Coutte流問題

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

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

x方向速度的解析解為

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

當(dāng)系統(tǒng)達(dá)到穩(wěn)態(tài)時(shí),沿y方向溫度場(chǎng)的理論解為

其中cp=γ/(γ?1).圖5展示了不同γ對(duì)應(yīng)的DBM數(shù)值解與解析解在穩(wěn)態(tài)時(shí)的對(duì)比圖.數(shù)值解與解析解符合較好,表明DBM能夠精確模擬不同熱傳導(dǎo)情形下的流體問題.

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

圖4 γ=1.4時(shí)熱 Couette流在不同時(shí)刻速度剖面的DBM數(shù)值解與解析解對(duì)比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流的穩(wěn)態(tài)溫度剖面的DBM數(shù)值解與解析解對(duì)比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不穩(wěn)定性數(shù)值模擬

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

圖6 多模RT不穩(wěn)定性在不同時(shí)刻的密度演化圖: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.

本文考慮二維區(qū)域 [?d/2,d/2]×[?2d,2d],系統(tǒng)處于重力加速度為常數(shù)的重力場(chǎng)下,界面的初始擾動(dòng)滿足

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

所以系統(tǒng)的不穩(wěn)定性初始條件滿足:

其中,p0是上部分流體頂部的初始?jí)簭?qiáng),Tu和Tb代表上下部分流體的初始溫度.在這種條件下,界面處的壓強(qiáng)滿足

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

在數(shù)值模擬中,計(jì)算區(qū)域?yàn)?12×512的均為網(wǎng)格,空間步長(zhǎng)為?x=?y=0.001,頂部初始?jí)簭?qiáng)為p0=1.0,時(shí)間步長(zhǎng)為 ?t=1×10?5,松弛因子為τ=1×10?5,上部分溫度為Tu=1.0,下部分溫度為Tb=4.0,因此,初始At=0.6.其他參數(shù)為c=1.3,η=15,n=3,ax=0.0,ay=?g=?1.0.

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

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

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

5 結(jié) 論

應(yīng)用含外力項(xiàng)的DBM數(shù)值模擬研究可壓流體多模初始擾動(dòng)的RT不穩(wěn)定性問題.Chapman-Enskog多尺度分析表明該模型在連續(xù)極限可恢復(fù)到Navier-Stokes方程.模型通過了熱Coutte流問題和三個(gè)一維Riemann問題的檢測(cè),表明模型能夠精確模擬黏性耗散和熱傳導(dǎo)以及復(fù)雜激波之間的相互作用.采用DBM對(duì)多模、可壓、具有間斷界面的多模初始擾動(dòng)RT不穩(wěn)定性進(jìn)行數(shù)值模擬.結(jié)果表明,在RT不穩(wěn)定性發(fā)展的初期由于多模的設(shè)置,界面處的黏性和熱傳導(dǎo)效應(yīng)突出,這些耗散效應(yīng)會(huì)“抹平”界面,多模之間相互競(jìng)爭(zhēng)和吸收,形成較少的主導(dǎo)模式;在這一階段系統(tǒng)內(nèi)沒有形成明顯的“氣泡”和“尖釘”結(jié)構(gòu).在RT不穩(wěn)定性的中后期,由于模式的合并導(dǎo)致界處的耗散效應(yīng)減弱,重力占主導(dǎo)地位,擾動(dòng)界面逐漸變形、長(zhǎng)大,形成典型的“氣泡-尖釘”結(jié)構(gòu),即出現(xiàn)典型的“蘑菇頭”形狀,而后進(jìn)入湍流混合階段.這些現(xiàn)象與經(jīng)典的實(shí)驗(yàn)結(jié)果一致.同時(shí)給出系統(tǒng)整體非平衡程度隨時(shí)間發(fā)展的演化情況,一開始系統(tǒng)先趨于平衡態(tài),這是由于系統(tǒng)處于調(diào)整階段,從多模初始界面擾動(dòng)調(diào)整到本征模階段;而后系統(tǒng)以線性形式偏離平衡態(tài),這是由于系統(tǒng)界面被抹平,壓縮能部分轉(zhuǎn)化為內(nèi)能;然后系統(tǒng)又趨于平衡態(tài),這是由于模式的耦合與擾動(dòng)界面進(jìn)一步被“抹平”,系統(tǒng)處于相對(duì)穩(wěn)定狀態(tài);最后系統(tǒng)越來越遠(yuǎn)離平衡態(tài),此時(shí)是由于系統(tǒng)輕重流體的重力勢(shì)能相互轉(zhuǎn)換,系統(tǒng)的壓縮能進(jìn)一步被釋放出來,系統(tǒng)動(dòng)能進(jìn)一步增加所致.在最近的一系列學(xué)術(shù)報(bào)告中,許愛國(guó)等[34?37]進(jìn)一步給出了非平衡程度更深、超越Navier-Stokes描述能力的復(fù)雜流動(dòng)系統(tǒng)的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)[許愛國(guó),張廣財(cái),甘延標(biāo)2016力學(xué)與實(shí)踐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)[許愛國(guó),張廣財(cái),李英駿,李華2014物理學(xué)進(jìn)展34 136]

[24]Xu A G,Zhang G C,Ying Y J 2015Acta Phys.Sin.64 184701(in Chinese)[許愛國(guó),張廣財(cái),應(yīng)陽(yáng)君2015物理學(xué)報(bào)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空氣動(dòng)力學(xué)學(xué)報(bào)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)[許愛國(guó),張廣財(cái)2016第九屆全國(guó)流體力學(xué)學(xué)術(shù)會(huì)議,南京,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)[許愛國(guó),張廣財(cái)2016南京林業(yè)大學(xué)機(jī)電學(xué)院專題學(xué)術(shù)報(bào)告,中國(guó)南京,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)[許愛國(guó),張廣財(cái) 2016中國(guó)人民大學(xué)物理系專題學(xué)術(shù)報(bào)告,中國(guó)北京,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)[許愛國(guó),張廣財(cái) 2016第四屆LBM及其應(yīng)用學(xué)術(shù)研討會(huì),中國(guó)北京,2016年11月26日]