李旭暉,單肖文,段文洋
(1. 哈爾濱工程大學(xué) 船舶工程學(xué)院,哈爾濱 150001;2. 南方科技大學(xué) 力學(xué)與航空航天工程系,深圳 518055)
在過去三十年,格子玻爾茲曼方法逐漸成為計算流體力學(xué)領(lǐng)域的一個重要分支,并廣泛應(yīng)用到科學(xué)和工程問題當(dāng)中[1-2]。不同于基于Navier-Stokes-Fourier(NSF)方程離散的各種傳統(tǒng)計算流體力學(xué)方法,格子玻爾茲曼方法是一種用統(tǒng)計概率粒子分布函數(shù)的碰撞和遷移來描述流體系統(tǒng)的介觀數(shù)值方法。它與NSF 方程的聯(lián)系可以在小Knudsen 數(shù)(Kn)下通過Chapman-Enskog(CE)展開[3]得到。傳統(tǒng)計算流體力學(xué)方法,如有限差分、有限體積法等,一個難點是處理偏微分方程中的非線性對流項,為了控制數(shù)值耗散和數(shù)值彌散誤差,常采用高階離散格式。在格子玻爾茲曼方法中,對于通常使用的On-Lattice 格子,統(tǒng)計粒子沿特征線進行遷移,即粒子能在一個時間步遷移到達相鄰的直角網(wǎng)格格點上。由于沒有任何插值操作,僅為計算機上的拷貝操作,這一拉格朗日屬性決定了格子玻爾茲曼方法處理對流的簡潔性以及它的低數(shù)值耗散和數(shù)值彌散特性。格子玻爾茲曼方法無需迭代求解壓力泊松方程,時間步顯式推進,非線性的碰撞項為空間完全局部操作,代數(shù)計算為主,數(shù)據(jù)訪存規(guī)則,這些特性與目前在高性能計算領(lǐng)域占越來越重要地位的異構(gòu)并行架構(gòu)—圖形計算處理器(graphic process unit,GPU)架構(gòu)和國產(chǎn)神威架構(gòu)(眾核處理器)特點相符,為超大規(guī)模高精度流體數(shù)值模擬提供了可能性[4-5]。由于格子玻爾茲曼方法的動理學(xué)基礎(chǔ)、執(zhí)行簡單且高效、天然的并行特性,目前它已被廣泛和成功地應(yīng)用于多相流[6-7]、湍流[8-9]、燃燒[10]、流固耦合[11-12]等領(lǐng)域。
然而,格子玻爾茲曼方法很多棘手的問題來自于非線性的碰撞項,如數(shù)值不穩(wěn)定性、輸運系數(shù)不滿足伽利略不變性和不可調(diào)的普朗特數(shù)等。過去幾十年,格子玻爾茲曼方法領(lǐng)域的理論工作很多都集中于處理復(fù)雜的非線性微分積分碰撞項。最簡潔的碰撞模型為錢躍竑等[13]基于Bhatnagar、Gross 和Krook等[14]提出的Boltzmann-BGK(作者三人首字母命名)方程發(fā)展的Lattice BGK(LBGK)模型,該模型中非平衡態(tài)部分的各階成分均以相同的弛豫時間朝平衡態(tài)松弛,該算法因執(zhí)行簡單、易于編程而廣受歡迎。但是,原始的LBGK 模型僅能處理等溫流動,數(shù)值穩(wěn)定性存在問題,壓力場存在很多虛假噪音。隨后很多其他的碰撞模型逐漸被提出,用以實現(xiàn)普朗特數(shù)的可調(diào)、良好的數(shù)值穩(wěn)定性、虛假壓力噪音的壓制和消除,以及輸運系數(shù)的伽利略不變性。對于存在大梯度的流場,幾乎所有的計算流體力學(xué)數(shù)值方法都會遇到數(shù)值穩(wěn)定性問題。對于格子玻爾茲曼方法,在高雷諾數(shù)和高馬赫數(shù)流動中,其數(shù)值穩(wěn)定性問題尤為突出。D'Humières 等[15-16]提 出 了 多 松 弛 時 間(multiplerelaxation-time,MRT)碰撞模型。在該方法中,利用正交化手段,通過離散速度集構(gòu)建了一組正交的基向量,數(shù)量與離散速度個數(shù)相等。通過這組正交的基向量,粒子分布函數(shù)可以計算出不同階的矩,對不同的矩分配不同的松弛時間,朝對應(yīng)矩的平衡態(tài)松弛。該模型雖然沒能通過多松弛因子實現(xiàn)普朗特數(shù)的可調(diào),但獲得了數(shù)值穩(wěn)定性的大幅提高。另一個能顯著提高數(shù)值穩(wěn)定性并改善剪切黏性伽利略不變性的碰撞模型為Geier 等[17-18]提出的級聯(lián)多松弛碰撞模型(Cascaded MRT)。所謂的級聯(lián)多松弛模型即為高階矩(大于等于三階)由低階矩和宏觀速度構(gòu)建,且每階矩獨立松弛。Lycett-Brown等[19-21]則從一維出發(fā),先通過分布函數(shù)與零到二階矩的關(guān)系將分布函數(shù)由這些矩顯式寫出,并用張量積的形式構(gòu)造了二維和三維級聯(lián)多松弛碰撞模型。他們的推導(dǎo)相比Geier 等的推導(dǎo)更加直觀易懂。所謂的級聯(lián)多松弛模型,本質(zhì)就是在中心矩空間里面對分布函數(shù)的中心矩進行松弛,然后通過中心矩和原點矩的數(shù)學(xué)關(guān)系將中心矩松弛轉(zhuǎn)換到原點矩形式,高階矩松弛由低階矩和當(dāng)?shù)睾暧^速度表示,即所謂的級聯(lián)。其他開展級聯(lián)多松弛模型工作的還有Dubois等[22-23]、Fei 和Luo 等[24-25],以及De Rosis 等[26]。然而,由于其采用的離散速度的限制,輸運系數(shù)的伽利略不變性在級聯(lián)多松弛模型中并沒有完全解決。2015 年,Geier 等[27]受Seeger 等[28]的計算氣體動理論累計量方法(Cumulant Method)的啟發(fā),提出了累計量格子玻爾茲曼方法(Cumulant LBM),并在后續(xù)工作中[29-30]進行了參數(shù)優(yōu)化討論。累計量多松弛模型中,以所謂的累計量為獨立松弛對象,二階和三階矩與前述的級聯(lián)碰撞模型一樣,只有在四階及以上的矩才產(chǎn)生偏差。Seeger等[28]在最早的工作中曾指出他的計算氣體動理論累積量方法與Grad[31]等的矩方法在一維空間至少到五階矩是吻合的,從第六階矩開始存在偏差。累計量碰撞模型確實取得了良好的數(shù)值穩(wěn)定性,并被應(yīng)用到湍流邊界層模擬[32]、氣動聲學(xué)模擬[33],以及自由面兩相流模擬[34]。Karlin 等[35-38]提出的熵格子玻爾茲曼模型是另外一個具有良好數(shù)值穩(wěn)定性的碰撞模型。在該模型中,H定理要求在離散分子速度空間得到滿足,平衡態(tài)分布函數(shù)通過迭代求解極值問題得到,高階格子也通過一維格子的張量積形式得到;松弛因子不再是固定值,而是隨流場空間動態(tài)變化,類似于大渦模擬模型。熵格子玻爾茲曼模型雖然能顯著提高數(shù)值穩(wěn)定性,然而由于其極值問題的迭代求解需要在每個網(wǎng)格點每個時間步執(zhí)行,其計算效率并無優(yōu)勢。
正則化松弛模型是另外一種具有良好數(shù)值穩(wěn)定性的碰撞模型。最早的正則化碰撞模型由Ladd[39]于1994 年提出,并成功應(yīng)用到粒子沉降流動中[40-41]。Ladd 認為對于弱可壓縮Navier-Stokes(N-S)方程的求解,格子玻爾茲曼碰撞項中只需保留二階應(yīng)力張量的松弛(包括剪切應(yīng)力的松弛和體積應(yīng)力的松弛),而其他高階非水動力學(xué)動理學(xué)模態(tài)的松弛頻率直接置為1,這樣使得不影響N-S 方程的非水動力學(xué)高階模態(tài)快速松弛。Ladd 模型中二階應(yīng)力張量中的有跡部分實際為速度空間積分精度誤差帶來的數(shù)值誤差項,該項的松弛會調(diào)節(jié)數(shù)值穩(wěn)定性。實際的二階應(yīng)力張量中有跡部分通常只在稠密氣體和考慮分子內(nèi)部自由度激發(fā)的時候才會出現(xiàn),體積黏性松弛頻率應(yīng)置為零,對此Ladd 在論文[41]中做了嚴(yán)格的論述。2005 年,Latt 等[42]也提出了一個正則化碰撞模型(類似于Ladd 的碰撞模型),并發(fā)現(xiàn)正則化執(zhí)行能大幅提高數(shù)值模擬的穩(wěn)定性和精度;但是他們并沒有討論二階應(yīng)力張量中無跡部分和有跡部分的分別松弛,可以視為Ladd 模型的一個特殊版本。幾乎同時,Chen 等[43]也討論了與Latt 等相同的正則化模型在提高高Kn數(shù)流動中數(shù)值格式的旋轉(zhuǎn)不變性的前景。Zhang 等[44]利用Shan 等[45]基于Hermite 多項式展開Boltzmann-BGK 方程的理論體系,將正則化模型延伸到三階,并在該模型中引入了力項的處理,計算了高Kn數(shù)流動。他們認為非平衡態(tài)分布函數(shù)應(yīng)該類似于平衡態(tài)分布函數(shù),在Hermite 截斷階數(shù)空間中進行重構(gòu);然而該模型所使用的是Hermite 原點矩空間,并且在力項的處理中并沒有考慮一階項。2014 年,Chen 等[46]將Shan 和Chen 在2007 年提出的Hermite 展開多松弛模型[47]做了一些改進,將原來的原點矩多松弛改為中心矩多松弛。他們發(fā)現(xiàn)二階矩松弛保持不變,而三階中心矩松弛轉(zhuǎn)換到原點矩空間后多出了一個修正項,該修正項與局部速度和低階矩松弛相關(guān),三階中心矩松弛完全克服了之前模型熱傳導(dǎo)系數(shù)在普朗特數(shù)不為1 時不滿足伽利略不變性的缺陷。2015 年,Malaspinas[48]基于Hermite 展開提出了一個遞歸正則化碰撞模型,證明了非平衡態(tài)Hermite 系數(shù)在等溫層級上存在遞歸關(guān)系。該模型中非平衡態(tài)Hermite 系數(shù)的遞歸關(guān)系主要由以下三個關(guān)系推導(dǎo)而來:1)平衡態(tài)Hermite 系數(shù)的遞歸關(guān)系,該關(guān)系用數(shù)學(xué)歸納法很容易得到;2)零階Chapman-Enskog 展開得到的歐拉方程;3)一階Chapman-Enskog 展開得到的非平衡態(tài)Hermite 系數(shù)與平衡態(tài)Hermite 系數(shù)的時空導(dǎo)數(shù)關(guān)系式。該遞歸正則化碰撞模型被作者證明比原始MRT[15-16]有更好的數(shù)值耗散和色散性質(zhì)。2017 年,Coreixas 等[49]將Malaspinas 的模型推廣到可壓縮傳熱層級上去。Mattila 等[50]利用中心矩空間Hermite基對非平衡態(tài)進行了展開,并在中心矩空間對展開進行了截斷,如:對可壓縮傳熱層級的非平衡態(tài)展開,將四階及以上的展開進行截斷,相當(dāng)于將四階及以上的非平衡態(tài)Hermite 中心矩設(shè)定松弛頻率為1;而對等溫層級的非平衡態(tài)展開,則將三階及以上的展開進行截斷,轉(zhuǎn)換到原點矩空間后發(fā)現(xiàn)各階非平衡態(tài)Hermite 原點矩剛好與Manaspinas[48]的遞歸正則化模型中的非平衡態(tài)Hermite 矩相同。然而,以上三個正則化模型只對非平衡態(tài)Hermite 矩進行了討論,而且只是對各階非平衡態(tài)Hermite 原點矩進行獨立松弛,并沒有在中心矩空間進行松弛,熱傳導(dǎo)系數(shù)的伽利略不變不被保證。2018 年,Jacob 等[51]又在Malaspinas的遞歸正則化模型的基礎(chǔ)上提出了一個混合遞歸正則化模型,對二階應(yīng)力張量的重構(gòu)做了調(diào)整;原來正則化模型中該項由非平衡態(tài)分布函數(shù)的二階矩得到,在Jacob 等的混合正則化模型中,二階應(yīng)力張量也可以用速度場的中心差分得到;他們還對兩種方法得到的二階應(yīng)力張量進行加權(quán),通過調(diào)節(jié)權(quán)值控制額外引入的數(shù)值耗散的大小。后續(xù)該混合正則化模型被Feng 等[52-54]用到了他們的求解器中。在他們的方法中,質(zhì)量和動量方程由格子玻爾茲曼方法在D3Q19和D3Q27 格子上進行求解,能量方程則通過有限體積方法進行求解,并通過修改平衡態(tài)分布函數(shù)和添加源項的方式盡可能地消除標(biāo)準(zhǔn)格子存在的速度離散誤差,實現(xiàn)了跨聲速和超聲速模擬。2019 年,Shan[55]在中心矩空間對非平衡態(tài)和碰撞項分別進行了展開,并在中心矩空間對每階矩進行獨立松弛,然后通過中心矩空間和原點矩空間Hermite 基之間的關(guān)系將碰撞項轉(zhuǎn)換回原點矩空間,得到了遞歸形式的碰撞項。該碰撞模型繼承了2007 年Shan 等[47]提出的原點矩空間Hermite 矩多松弛模型的優(yōu)點,克服了Chen 等[46]提到的熱傳導(dǎo)系數(shù)的伽利略不變性問題,并可以推廣到任意高階松弛。隨后,Li 等[56]在溫度標(biāo)度的中心矩空間中,對非平衡態(tài)和碰撞項進行了展開,并在該空間中對各階矩進行松弛,利用溫度標(biāo)度中心矩空間和原點矩空間Hermite 基之間的關(guān)系將碰撞項轉(zhuǎn)換回原點矩空間,也得到了遞歸形式的碰撞項;與原點矩空間松弛的碰撞項[47]相比,三階及以上的碰撞項包含速度、溫度及低階矩松弛構(gòu)成的修正項。
關(guān)于正則化碰撞模型的邊界條件處理,歷史上也有一些學(xué)者取得了部分進展。2007 年,Niu 等[57]模擬了微尺度稀薄氣體流動,討論了散射邊界條件、松弛時間和正則化碰撞模型。2008 年,Latt 等[58]基于正則化重構(gòu)非平衡態(tài)的思路發(fā)展了適用于平直邊界的正則化邊界條件。2011 年,Malaspinas 等[59]提出了一個通用的處理平直邊界的正則化邊界條件,通過建立邊界點已知分布函數(shù)與宏觀量的超定方程組,解出宏觀矩,再重構(gòu)出邊界點的非平衡態(tài)分布函數(shù)。該邊界處理方法可適用于標(biāo)準(zhǔn)格子(D2Q9、D3Q19 等),也適用于高階格子(D2V37、D3V103 等),但是否適用于任意曲線曲面邊界需進一步研究。2017 年,Wissocq 等[60]發(fā)展了適用于聲學(xué)模擬的正則化特征邊界條件。2019 年,Guo 等[61]在混合正則化模型中討論了可壓縮流動的任意曲面邊界條件和域邊界特征邊界條件。
在將正則化碰撞模型應(yīng)用到多相流方面,也有一些進展。Otomo 等將正則化碰撞模型分別應(yīng)用到多組分偽勢模型[62]和相場多相流模型[63]。Ba 等[64]將正則化碰撞模型應(yīng)用到顏色梯度模型。另外,將正則化碰撞模型用于噪聲模擬也有很多代表性工作,如Zhuo 等[65]、Chen 等[66]、Casalino 等[67-68]、Gendre 等[69]、Astoul 等[70]。本文對正則化碰撞模型在各個領(lǐng)域的應(yīng)用不再做詳細展開。
本文將重點闡述、回顧正則化碰撞模型,通過一個系統(tǒng)的理論分析框架闡述作者與合作者提出的高階正則化碰撞模型,并對以往的正則化碰撞模型的理論關(guān)系進行梳理,同時對正則化多松弛模型的本質(zhì)進行歸納分析。后文將按以下順序進行論述:第1 節(jié)闡述Hermite 展開的基礎(chǔ)理論框架,為后文中的碰撞算子正則化提供分析的準(zhǔn)則;第2 節(jié)則利用第1 節(jié)建立的理論體系,對正則化碰撞模型發(fā)展歷史上有重要貢獻的經(jīng)典模型逐個進行分析、梳理和比較;最后對全文進行歸納總結(jié),得出全文理論分析的結(jié)論。
該部分我們從Boltzmann 方程出發(fā),在Hermite譜空間對分布函數(shù)進行展開,分別在原點矩空間、中心矩空間和溫度標(biāo)度的中心矩空間對Maxwell-Boltzmann 平衡態(tài)、非平衡態(tài)、碰撞項和力項進行展開,為后面的正則化理論模型回顧做準(zhǔn)備。
我們先從Boltzmann 方程進行闡述。Boltzmann方程為基于分子二體碰撞、分子混沌假設(shè)和外力不影響局部碰撞的動力學(xué)行為等三個假設(shè),推導(dǎo)的稀薄氣體系統(tǒng)的控制方程。它描述了統(tǒng)計粒子分布函數(shù)在速度空間和物理空間的演化,其右端的非線性碰撞項極為復(fù)雜,直接求解困難重重。Boltzmann-BGK 方程[14]作為Boltzmann 方程碰撞項的一個極度簡化版本,認為所有動理學(xué)矩都以相同的速度朝平衡態(tài)弛豫,其控制方程如下:
在實際的數(shù)值執(zhí)行中,粒子分布函數(shù)不可能展開到無窮階,必須進行截斷。截斷的階數(shù)直接影響水動力學(xué)宏觀方程的適用范圍及其精度,這意味著實際上需要將分布函數(shù)投影到有限階Hermite 基上去。對于截斷階數(shù)如何影響恢復(fù)的水動力學(xué)方程層級以及其
受益于Hermite 多項式的正交屬性,粒子分布函數(shù)的前若干階矩是可以由截斷形式的Hermite 展開得到的粒子分布函數(shù)獲得。因此,粒子分布函數(shù)可以投影到由前N階Hermite 多項式構(gòu)成的Hilbert 空間,并保證其前N階矩與完整粒子分布函數(shù)的前N階矩嚴(yán)格相等。粒子分布函數(shù)的前N階Hermite 展開截斷為:
為了便于后文對不同空間中Hermite 系數(shù)進行評估,對于不同Hermite 基之間的轉(zhuǎn)換關(guān)系做一個簡單介紹,詳細的推導(dǎo)過程可以參考Shan[55,72]、Li 等[56]、Shan 等[72]的論文附錄部分,他們給出了詳細的數(shù)學(xué)推導(dǎo),這里只給出前四階的數(shù)學(xué)轉(zhuǎn)換關(guān)系。由于在實際的數(shù)值執(zhí)行中,中心矩空間的離散速度和權(quán)值會隨當(dāng)?shù)氐暮暧^速度而變化;而溫度標(biāo)度中心矩空間的離散速度和權(quán)值則會隨當(dāng)?shù)氐暮暧^速度和溫度發(fā)生變化。這會破壞格子玻爾茲曼方法完美的遷移操作,而引入額外的數(shù)值耗散和色散,讓格子玻爾茲曼方法的優(yōu)點喪失。因此,我們的思路是在中心矩空間和溫度標(biāo)度中心矩空間完成碰撞,然后通過它們與原點矩空間的數(shù)學(xué)關(guān)系轉(zhuǎn)換到原點矩空間,最終在原點矩空間實現(xiàn)粒子遷移操作。
首先寫出原點矩Hermite 基的零到四階表達式,如下(具體可參考文獻[55]):
實際上,歷史上發(fā)展的不同的正則化模型的差異主要在以下幾個地方:一個是使用的Hermite 基不同,即在不同的空間(原點矩空間、中心矩空間和溫度標(biāo)度中心矩空間)中展開;另一個是分布函數(shù)Hermite 展開截斷的階數(shù)的差異,尤其是非平衡態(tài)Hermite 展開的截斷;還有一個是碰撞項展開的空間和碰撞項中弛豫的張量矩分量。因此,為了對這些正則化碰撞模型進行系統(tǒng)地回顧,我們將平衡態(tài)、非平衡態(tài)、碰撞項和力項都用三種不同的Hermite 基進行展開,然后對它們進行系統(tǒng)地分析。
1.4.1 平衡態(tài)分布函數(shù)的Hermite 展開
使用權(quán)函數(shù)定義表達式(5),Maxwell-Boltzmann平衡態(tài)分布函數(shù)可寫為下式:
利用中心矩Hermite 基和原點矩Hermite 基的數(shù)學(xué)轉(zhuǎn)換關(guān)系式(22),上式中各階平衡態(tài)中心矩Hermite 系數(shù)可相應(yīng)得到,這里顯式地寫出零到四階:
其中二階和三階非平衡態(tài)系數(shù)可通過Chapman-Enskog 展開[73]得到:
利用溫度標(biāo)度中心矩Hermite 基與原點矩Hermite 基的數(shù)學(xué)關(guān)系,我們可以計算上式,顯式地寫出零到四階:
上面的推導(dǎo)中用到了Maxwell-Boltzmann 平衡態(tài)分布函數(shù),且推導(dǎo)式(63)中用到了等溫假設(shè)。式(62)如果在速度空間進行離散且考慮離散格子效應(yīng),即為Guo 等[76]2002 年推導(dǎo)的力項格式(具體離散在后文中闡述);式(63)即為He 等[75]1998 年設(shè)計的力項格式;式(64)與式(63)等價。
其中多弛豫時間被應(yīng)用,不是原始的Boltzmann-BGK 方程中的單松弛因子,而是在每階矩應(yīng)用了不同的松弛因子,可視為其拓展模型。在我們的工作中,上式的約束條件被拓展到非平衡態(tài)分布函數(shù)的Hermite 矩松弛,這必須與原始的碰撞項的松弛在討論的階數(shù)上分別等價,即:
將上式代入式(67),并將碰撞項和非平衡態(tài)分布函數(shù)的Hermite 展開式也代入,并利用Hermite 多項式的正交特性,可得到下面關(guān)系式:
點評:從式(75)可以看出,不管是Hermite 原點矩、中心矩還是溫度標(biāo)度中心矩,每階都可以獨立松弛。然而從本文1.4.2 節(jié)對非平衡態(tài)Hermite 系數(shù)在不同空間之間的數(shù)學(xué)關(guān)系的分析,我們可以發(fā)現(xiàn),二階矩在不同空間中是等價的,只要離散速度足夠多,能保證二階矩的積分精度,分子速度空間離散形式的非平衡態(tài)分布函數(shù)的二階矩就會等價于連續(xù)形式非平衡態(tài)分布函數(shù)的二階矩。Shan 等[45,71]已證明對于二維至少需要17 個離散速度,三維至少需要39 個離散速度。只要離散速度充分,同時平衡態(tài)展開到三階及以上,不論在哪個空間進行松弛,黏性輸運系數(shù)的伽利略不變性都能得到嚴(yán)格滿足[55]。使用D2Q9 或者D3Q27 離散格子級聯(lián)碰撞模型[17]或者累積量碰撞模型[27]都不能保證嚴(yán)格滿足黏性輸運系數(shù)的伽利略不變性。在不同的平移坐標(biāo)系中,宏觀速度是不一樣的,且是有量綱的,在溫度標(biāo)度中心矩空間進行碰撞,各階矩是與宏觀速度無關(guān)且被局部聲速無量綱化的量。從式(82)~式(84)即可看出,在溫度標(biāo)度中心矩空間進行松弛,松弛的量是真正的非平衡態(tài)量,在原點矩松弛對象中將守恒量(動量、內(nèi)能)剔除,從而保證了每階矩松弛的真正獨立。結(jié)合式(76)、式(82)~式(84)我們可以發(fā)現(xiàn),原點矩空間對每階非平衡態(tài)Hermite 矩進行松弛,實質(zhì)上在三階及以上矩松弛中與低階矩發(fā)生了相互干擾,即所謂的串?dāng)_效應(yīng)。這也說明原點矩空間中的各階矩松弛并不真正的獨立。因此,在碰撞算子中,松弛對象必須滿足坐標(biāo)的平移不變性。Li 等[77]最近的工作討論了松弛對象必須滿足坐標(biāo)的旋轉(zhuǎn)不變性。本綜述中暫不作展開。
包含力項的速度相空間離散形式的格子玻爾茲曼方程為:
將新定義變量代入式(86),可得到如下顯式方程,右端項均為當(dāng)前時間步的量:
特別注意地,當(dāng)不包含力項時,碰撞方程從二階開始,因為質(zhì)量和動量守恒。但是對于新定義的分布函數(shù),如果包含力項,則一階非平衡態(tài)矩不為零,這在多相流[62-63]和浸沒邊界法[11]包含力項時需要注意。先對新定義的粒子分布函數(shù)所求得的宏觀量與實際物理量之間的關(guān)系做一個討論。
對于原點矩空間碰撞,平衡態(tài)部分、非平衡態(tài)部分以及力項均可以直接利用Hermite 展開的截斷階形式重構(gòu)粒子分布函數(shù)。而對于中心矩和溫度標(biāo)度中心矩空間的碰撞形式,由于離散速度依賴于當(dāng)?shù)厮俣群蜏囟?,我們選擇將其轉(zhuǎn)換到原點矩空間,再進行粒子分布函數(shù)重構(gòu)。根據(jù)矩空間碰撞執(zhí)行的空間不同以及截斷階的階數(shù),我們可以對已存在的正則化模型進行系統(tǒng)地分析和討論。為方便起見,后文中非平衡態(tài)矩上面的一拔標(biāo)記全部省略。
2.2.1 Ladd 正則化模型和Latt 正則化模型
Ladd 等在文獻[39](1994 年)和文獻[41](2001 年)中提出了正則化碰撞模型的最早形式。在他們的碰撞模型中,平衡態(tài)分布函數(shù)展開到二階,非平衡態(tài)分布函數(shù)也僅保留到二階。碰撞后的分布函數(shù)寫為下面形式:
Latt 等[42]在2005 年也提出了正則化模型,他們的模型并未分離二階應(yīng)力張量的無跡部分和有跡部分的獨立松弛,僅為
Latt 的模型在Ladd 模型之后很久才提出,而且可視為Ladd 模型的一個特例,即剪切和體積粘性松弛因子相等,本質(zhì)上不能算一個新的模型。Ladd 模型為原點矩空間的正則化碰撞模型,高于三階的超出等溫Navier-Stokes 方程的非平衡態(tài)矩可視為松弛因子都設(shè)置為1,對碰撞后的粒子分布函數(shù)不產(chǎn)生貢獻。
2.2.2 Zhang-Shan-Chen 正則化模型
2006 年,Zhang、Shan 和Chen[43]基 于Hermite 展開提出了高階正則化模型,并用來模擬高Kn數(shù)非連續(xù)流動。他們認為平衡態(tài)分布函數(shù)投影到有限階Hermite 基構(gòu)成的Hilbert 空間,非平衡態(tài)部分也應(yīng)該投影到該Hilbert 空間上,并引入力項。根據(jù)本文2.1 節(jié)的討論,在考慮外力項后,非平衡態(tài)Hermite 系數(shù)的一階項并不為零,因此正則化的碰撞后粒子分布函數(shù)應(yīng)該寫為:
2.2.3 Mattila-Philippi-Hegele 正則化模型
Mattila 等[49]提出了一個基于中心矩空間展開的高階格子正則化模型來處理可壓縮傳熱流動,我們利用本文1.4.4 節(jié)的有關(guān)碰撞項和1.4.2 節(jié)的非平衡態(tài)系數(shù)在中心矩空間和原點矩空間分別展開的數(shù)學(xué)關(guān)系,分析該正則化模型。該模型中不包含力項,因此我們也省略力項。為分析方便,我們把非平衡態(tài)Hermite系數(shù)和碰撞項合并到一起,如式(95)所示,在式(77)中引入新的碰撞形式,包含非平衡的中心矩Hermite 系數(shù),具體如下:
其中四階及以上非平衡態(tài)松弛因子設(shè)為1,處理可壓縮傳熱流動。即認為超出二階應(yīng)力張量和三階熱流矢量的高階動理學(xué)中心矩快速松弛,碰撞后的粒子分布函數(shù)直接為平衡態(tài)。利用Hermite 中心矩和原點矩之間的關(guān)系可整理上式為:
展開基跟平衡態(tài)展開基對應(yīng),二階非平衡態(tài)Hermite 系數(shù)通過非平衡態(tài)粒子分布函數(shù)的二階矩求得,三階和四階非平衡態(tài)Hermite 系數(shù)則遞歸求得。
相應(yīng)地,三維(D3Q27)的平衡態(tài)展開零到二階跟二維一樣,但高階一直展開到六階,一共27 個Hermite 基,見式(129):
最后的正則化碰撞和遷移格式為:
2.2.5 Li-Shi-Shan 正則化模型
Li 等[56]在2019 年提出了一個基于溫度標(biāo)度中心矩空間展開的高階格子正則化模型來處理可壓縮傳熱流動。我們利用本文1.4.4 節(jié)的有關(guān)碰撞項和1.4.2 節(jié)的非平衡態(tài)系數(shù)在溫度標(biāo)度中心矩空間和原點矩空間分別展開的數(shù)學(xué)關(guān)系來分析該正則化模型。為了跟本文2.2.3 節(jié)分析保持一致,我們把溫度標(biāo)度中心矩空間的非平衡態(tài)Hermite 系數(shù)及其松弛進行合并,在式(81)的基礎(chǔ)上定義新的碰撞形式:
因此,Li 等提出的基于溫度標(biāo)度中心矩空間中Hermite 展開的碰撞形式可視為將該空間中四階及以上的非平衡態(tài)Hermite 矩的松弛因子置為1,該空間中的高階矩快速松弛到平衡態(tài)。比較Zhang-Shan-Chen 在原點矩空間的正則化碰撞模型、Mattila-Philippi-Hegele 在中心矩空間的正則化碰撞模型,以及作者等的Li-Shi-Shan 在溫度標(biāo)度中心矩空間的正則化碰撞模型,三者的相同點是將四階及以上的相應(yīng)非平衡態(tài)高階矩的松弛因子置為1,讓這些高階非平衡態(tài)動力學(xué)矩快速松弛到平衡態(tài),而本質(zhì)差異是進行碰撞的空間不同。這些模型的聯(lián)系和差異,在本文的分析框架下變得十分明朗。另外,在溫度標(biāo)度中心矩空間的正則化模型,在包含力項后的形式,我們將在另外的論文中發(fā)表。
2.2.6 Coreixas 等正則化模型
2017 年,Coreixas 等[50]將Malaspinas 提出的等溫層級的遞歸正則化模型推廣到了可壓縮傳熱層級上,發(fā)展了一個高階格子遞歸正則化模型。該正則化模型包含平衡態(tài)Hermite 系數(shù)的遞歸關(guān)系式:
上述遞歸關(guān)系式十分復(fù)雜,其數(shù)學(xué)推導(dǎo)過程也極其復(fù)雜。我們對四階非平衡態(tài)Hermite 系數(shù)進行整理如下:
與Li-Shi-Shan 正則化碰撞模型式(132)~式(134)相比,Coreixas 等的高階遞歸正則化模型即使使用了多松弛模型,也不能保證熱傳導(dǎo)系數(shù)的伽利略不變性;而Li-Shi-Shan 正則化碰撞模型不僅可以得到非平衡態(tài)Hermite 系數(shù)的遞歸形式,還能得到碰撞項的遞歸形式。這充分體現(xiàn)了對溫度標(biāo)度中心矩空間進行展開和正則化在數(shù)學(xué)和物理基本原理上的優(yōu)勢。
本文對格子玻爾茲曼領(lǐng)域的正則化碰撞模型進行了系統(tǒng)的嚴(yán)謹(jǐn)?shù)睦碚摶仡?,通過采用不同自變量的Hermite 基對平衡態(tài)、非平衡態(tài)、力項和碰撞項進行展開,將所有的理論模型在該理論體系下的定位以及它們之間的聯(lián)系進行了分析。分析發(fā)現(xiàn):
1) 每一種正則化碰撞模型都可以在Hermite 展開的理論框架下得到。
2) Malaspinas 和Coreixas 的遞歸正則化模型,非平衡態(tài)可以無限遞歸得到,對于超出所研究的動理學(xué)矩(二階應(yīng)力張量矩和三階熱流矢量矩),其弛豫時間不一定為1。前者由與剪切粘性相關(guān)的弛豫時間統(tǒng)一松弛,后者則對每一階非平衡態(tài)原點矩進行獨立松弛。高階矩的松弛因子如何影響模型的數(shù)值穩(wěn)定性和精度有待進一步探究。
3) 除2)中提到的遞歸正則化模型,其他的正則化模型等價于在原點矩空間、中心矩空間和溫度標(biāo)度中心矩空間,將高于所關(guān)注的非平衡態(tài)動理學(xué)矩的弛豫時間置為1。
4) 可模擬可壓縮傳熱NSF 方程的正則化模型中,只有Chen 等[46]的模型、 Shan 的模型[55]、Li-Shi-Shan 模型[56]是在中心矩或溫度標(biāo)度中心矩空間進行松弛,滿足熱傳導(dǎo)系數(shù)的伽利略不變性。
正則化多松弛碰撞模型的理論本質(zhì)就是將離散速度表達不了的高階動理學(xué)矩進行截斷,或者等價于將高階動理學(xué)矩的松弛因子置為1,以避免虛假模態(tài)對水動力學(xué)方程的負面影響;而高階動理學(xué)矩的松弛因子不為1 時,對數(shù)值穩(wěn)定性的正面影響目前尚不明朗。展望未來,正則化模型相較于其他碰撞模型的優(yōu)勢還有待嚴(yán)格地理論證明,尤其是在湍流模擬、聲學(xué)模擬、多相流以及稀薄氣體流動等領(lǐng)域,還有很大的發(fā)展空間。