格子玻爾茲曼正則化碰撞模型的理論進(jìn)展

李旭暉，單肖文，段文洋

(1. 哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001；2. 南方科技大學(xué) 力學(xué)與航空航天工程系，深圳 518055)

0 引 言

在過去三十年，格子玻爾茲曼方法逐漸成為計(jì)算流體力學(xué)領(lǐng)域的一個(gè)重要分支，并廣泛應(yīng)用到科學(xué)和工程問題當(dāng)中[1-2]。不同于基于Navier-Stokes-Fourier（NSF）方程離散的各種傳統(tǒng)計(jì)算流體力學(xué)方法，格子玻爾茲曼方法是一種用統(tǒng)計(jì)概率粒子分布函數(shù)的碰撞和遷移來描述流體系統(tǒng)的介觀數(shù)值方法。它與NSF 方程的聯(lián)系可以在小Knudsen 數(shù)（Kn）下通過Chapman-Enskog（CE）展開[3]得到。傳統(tǒng)計(jì)算流體力學(xué)方法，如有限差分、有限體積法等，一個(gè)難點(diǎn)是處理偏微分方程中的非線性對(duì)流項(xiàng)，為了控制數(shù)值耗散和數(shù)值彌散誤差，常采用高階離散格式。在格子玻爾茲曼方法中，對(duì)于通常使用的On-Lattice 格子，統(tǒng)計(jì)粒子沿特征線進(jìn)行遷移，即粒子能在一個(gè)時(shí)間步遷移到達(dá)相鄰的直角網(wǎng)格格點(diǎn)上。由于沒有任何插值操作，僅為計(jì)算機(jī)上的拷貝操作，這一拉格朗日屬性決定了格子玻爾茲曼方法處理對(duì)流的簡(jiǎn)潔性以及它的低數(shù)值耗散和數(shù)值彌散特性。格子玻爾茲曼方法無需迭代求解壓力泊松方程，時(shí)間步顯式推進(jìn)，非線性的碰撞項(xiàng)為空間完全局部操作，代數(shù)計(jì)算為主，數(shù)據(jù)訪存規(guī)則，這些特性與目前在高性能計(jì)算領(lǐng)域占越來越重要地位的異構(gòu)并行架構(gòu)—圖形計(jì)算處理器（graphic process unit，GPU）架構(gòu)和國(guó)產(chǎn)神威架構(gòu)（眾核處理器）特點(diǎn)相符，為超大規(guī)模高精度流體數(shù)值模擬提供了可能性[4-5]。由于格子玻爾茲曼方法的動(dòng)理學(xué)基礎(chǔ)、執(zhí)行簡(jiǎn)單且高效、天然的并行特性，目前它已被廣泛和成功地應(yīng)用于多相流[6-7]、湍流[8-9]、燃燒[10]、流固耦合[11-12]等領(lǐng)域。

然而，格子玻爾茲曼方法很多棘手的問題來自于非線性的碰撞項(xiàng)，如數(shù)值不穩(wěn)定性、輸運(yùn)系數(shù)不滿足伽利略不變性和不可調(diào)的普朗特?cái)?shù)等。過去幾十年，格子玻爾茲曼方法領(lǐng)域的理論工作很多都集中于處理復(fù)雜的非線性微分積分碰撞項(xiàng)。最簡(jiǎn)潔的碰撞模型為錢躍竑等[13]基于Bhatnagar、Gross 和Krook等[14]提出的Boltzmann-BGK（作者三人首字母命名）方程發(fā)展的Lattice BGK(LBGK)模型，該模型中非平衡態(tài)部分的各階成分均以相同的弛豫時(shí)間朝平衡態(tài)松弛，該算法因執(zhí)行簡(jiǎn)單、易于編程而廣受歡迎。但是，原始的LBGK 模型僅能處理等溫流動(dòng)，數(shù)值穩(wěn)定性存在問題，壓力場(chǎng)存在很多虛假噪音。隨后很多其他的碰撞模型逐漸被提出，用以實(shí)現(xiàn)普朗特?cái)?shù)的可調(diào)、良好的數(shù)值穩(wěn)定性、虛假壓力噪音的壓制和消除，以及輸運(yùn)系數(shù)的伽利略不變性。對(duì)于存在大梯度的流場(chǎng)，幾乎所有的計(jì)算流體力學(xué)數(shù)值方法都會(huì)遇到數(shù)值穩(wěn)定性問題。對(duì)于格子玻爾茲曼方法，在高雷諾數(shù)和高馬赫數(shù)流動(dòng)中，其數(shù)值穩(wěn)定性問題尤為突出。D'Humières 等[15-16]提 出 了 多 松 弛 時(shí) 間（multiplerelaxation-time，MRT）碰撞模型。在該方法中，利用正交化手段，通過離散速度集構(gòu)建了一組正交的基向量，數(shù)量與離散速度個(gè)數(shù)相等。通過這組正交的基向量，粒子分布函數(shù)可以計(jì)算出不同階的矩，對(duì)不同的矩分配不同的松弛時(shí)間，朝對(duì)應(yīng)矩的平衡態(tài)松弛。該模型雖然沒能通過多松弛因子實(shí)現(xiàn)普朗特?cái)?shù)的可調(diào)，但獲得了數(shù)值穩(wěn)定性的大幅提高。另一個(gè)能顯著提高數(shù)值穩(wěn)定性并改善剪切黏性伽利略不變性的碰撞模型為Geier 等[17-18]提出的級(jí)聯(lián)多松弛碰撞模型（Cascaded MRT）。所謂的級(jí)聯(lián)多松弛模型即為高階矩（大于等于三階）由低階矩和宏觀速度構(gòu)建，且每階矩獨(dú)立松弛。Lycett-Brown等[19-21]則從一維出發(fā)，先通過分布函數(shù)與零到二階矩的關(guān)系將分布函數(shù)由這些矩顯式寫出，并用張量積的形式構(gòu)造了二維和三維級(jí)聯(lián)多松弛碰撞模型。他們的推導(dǎo)相比Geier 等的推導(dǎo)更加直觀易懂。所謂的級(jí)聯(lián)多松弛模型，本質(zhì)就是在中心矩空間里面對(duì)分布函數(shù)的中心矩進(jìn)行松弛，然后通過中心矩和原點(diǎn)矩的數(shù)學(xué)關(guān)系將中心矩松弛轉(zhuǎn)換到原點(diǎn)矩形式，高階矩松弛由低階矩和當(dāng)?shù)睾暧^速度表示，即所謂的級(jí)聯(lián)。其他開展級(jí)聯(lián)多松弛模型工作的還有Dubois等[22-23]、Fei 和Luo 等[24-25]，以及De Rosis 等[26]。然而，由于其采用的離散速度的限制，輸運(yùn)系數(shù)的伽利略不變性在級(jí)聯(lián)多松弛模型中并沒有完全解決。2015 年，Geier 等[27]受Seeger 等[28]的計(jì)算氣體動(dòng)理論累計(jì)量方法（Cumulant Method）的啟發(fā)，提出了累計(jì)量格子玻爾茲曼方法（Cumulant LBM），并在后續(xù)工作中[29-30]進(jìn)行了參數(shù)優(yōu)化討論。累計(jì)量多松弛模型中，以所謂的累計(jì)量為獨(dú)立松弛對(duì)象，二階和三階矩與前述的級(jí)聯(lián)碰撞模型一樣，只有在四階及以上的矩才產(chǎn)生偏差。Seeger等[28]在最早的工作中曾指出他的計(jì)算氣體動(dòng)理論累積量方法與Grad[31]等的矩方法在一維空間至少到五階矩是吻合的，從第六階矩開始存在偏差。累計(jì)量碰撞模型確實(shí)取得了良好的數(shù)值穩(wěn)定性，并被應(yīng)用到湍流邊界層模擬[32]、氣動(dòng)聲學(xué)模擬[33]，以及自由面兩相流模擬[34]。Karlin 等[35-38]提出的熵格子玻爾茲曼模型是另外一個(gè)具有良好數(shù)值穩(wěn)定性的碰撞模型。在該模型中，H定理要求在離散分子速度空間得到滿足，平衡態(tài)分布函數(shù)通過迭代求解極值問題得到，高階格子也通過一維格子的張量積形式得到；松弛因子不再是固定值，而是隨流場(chǎng)空間動(dòng)態(tài)變化，類似于大渦模擬模型。熵格子玻爾茲曼模型雖然能顯著提高數(shù)值穩(wěn)定性，然而由于其極值問題的迭代求解需要在每個(gè)網(wǎng)格點(diǎn)每個(gè)時(shí)間步執(zhí)行，其計(jì)算效率并無優(yōu)勢(shì)。

正則化松弛模型是另外一種具有良好數(shù)值穩(wěn)定性的碰撞模型。最早的正則化碰撞模型由Ladd[39]于1994 年提出，并成功應(yīng)用到粒子沉降流動(dòng)中[40-41]。Ladd 認(rèn)為對(duì)于弱可壓縮Navier-Stokes（N-S）方程的求解，格子玻爾茲曼碰撞項(xiàng)中只需保留二階應(yīng)力張量的松弛（包括剪切應(yīng)力的松弛和體積應(yīng)力的松弛），而其他高階非水動(dòng)力學(xué)動(dòng)理學(xué)模態(tài)的松弛頻率直接置為1，這樣使得不影響N-S 方程的非水動(dòng)力學(xué)高階模態(tài)快速松弛。Ladd 模型中二階應(yīng)力張量中的有跡部分實(shí)際為速度空間積分精度誤差帶來的數(shù)值誤差項(xiàng)，該項(xiàng)的松弛會(huì)調(diào)節(jié)數(shù)值穩(wěn)定性。實(shí)際的二階應(yīng)力張量中有跡部分通常只在稠密氣體和考慮分子內(nèi)部自由度激發(fā)的時(shí)候才會(huì)出現(xiàn)，體積黏性松弛頻率應(yīng)置為零，對(duì)此Ladd 在論文[41]中做了嚴(yán)格的論述。2005 年，Latt 等[42]也提出了一個(gè)正則化碰撞模型（類似于Ladd 的碰撞模型），并發(fā)現(xiàn)正則化執(zhí)行能大幅提高數(shù)值模擬的穩(wěn)定性和精度；但是他們并沒有討論二階應(yīng)力張量中無跡部分和有跡部分的分別松弛，可以視為L(zhǎng)add 模型的一個(gè)特殊版本。幾乎同時(shí)，Chen 等[43]也討論了與Latt 等相同的正則化模型在提高高Kn數(shù)流動(dòng)中數(shù)值格式的旋轉(zhuǎn)不變性的前景。Zhang 等[44]利用Shan 等[45]基于Hermite 多項(xiàng)式展開Boltzmann-BGK 方程的理論體系，將正則化模型延伸到三階，并在該模型中引入了力項(xiàng)的處理，計(jì)算了高Kn數(shù)流動(dòng)。他們認(rèn)為非平衡態(tài)分布函數(shù)應(yīng)該類似于平衡態(tài)分布函數(shù)，在Hermite 截?cái)嚯A數(shù)空間中進(jìn)行重構(gòu)；然而該模型所使用的是Hermite 原點(diǎn)矩空間，并且在力項(xiàng)的處理中并沒有考慮一階項(xiàng)。2014 年，Chen 等[46]將Shan 和Chen 在2007 年提出的Hermite 展開多松弛模型[47]做了一些改進(jìn)，將原來的原點(diǎn)矩多松弛改為中心矩多松弛。他們發(fā)現(xiàn)二階矩松弛保持不變，而三階中心矩松弛轉(zhuǎn)換到原點(diǎn)矩空間后多出了一個(gè)修正項(xiàng)，該修正項(xiàng)與局部速度和低階矩松弛相關(guān)，三階中心矩松弛完全克服了之前模型熱傳導(dǎo)系數(shù)在普朗特?cái)?shù)不為1 時(shí)不滿足伽利略不變性的缺陷。2015 年，Malaspinas[48]基于Hermite 展開提出了一個(gè)遞歸正則化碰撞模型，證明了非平衡態(tài)Hermite 系數(shù)在等溫層級(jí)上存在遞歸關(guān)系。該模型中非平衡態(tài)Hermite 系數(shù)的遞歸關(guān)系主要由以下三個(gè)關(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ù)的時(shí)空導(dǎo)數(shù)關(guān)系式。該遞歸正則化碰撞模型被作者證明比原始MRT[15-16]有更好的數(shù)值耗散和色散性質(zhì)。2017 年，Coreixas 等[49]將Malaspinas 的模型推廣到可壓縮傳熱層級(jí)上去。Mattila 等[50]利用中心矩空間Hermite基對(duì)非平衡態(tài)進(jìn)行了展開，并在中心矩空間對(duì)展開進(jìn)行了截?cái)啵纾簩?duì)可壓縮傳熱層級(jí)的非平衡態(tài)展開，將四階及以上的展開進(jìn)行截?cái)啵喈?dāng)于將四階及以上的非平衡態(tài)Hermite 中心矩設(shè)定松弛頻率為1；而對(duì)等溫層級(jí)的非平衡態(tài)展開，則將三階及以上的展開進(jìn)行截?cái)?，轉(zhuǎn)換到原點(diǎn)矩空間后發(fā)現(xiàn)各階非平衡態(tài)Hermite 原點(diǎn)矩剛好與Manaspinas[48]的遞歸正則化模型中的非平衡態(tài)Hermite 矩相同。然而，以上三個(gè)正則化模型只對(duì)非平衡態(tài)Hermite 矩進(jìn)行了討論，而且只是對(duì)各階非平衡態(tài)Hermite 原點(diǎn)矩進(jìn)行獨(dú)立松弛，并沒有在中心矩空間進(jìn)行松弛，熱傳導(dǎo)系數(shù)的伽利略不變不被保證。2018 年，Jacob 等[51]又在Malaspinas的遞歸正則化模型的基礎(chǔ)上提出了一個(gè)混合遞歸正則化模型，對(duì)二階應(yīng)力張量的重構(gòu)做了調(diào)整；原來正則化模型中該項(xiàng)由非平衡態(tài)分布函數(shù)的二階矩得到，在Jacob 等的混合正則化模型中，二階應(yīng)力張量也可以用速度場(chǎng)的中心差分得到；他們還對(duì)兩種方法得到的二階應(yīng)力張量進(jìn)行加權(quán)，通過調(diào)節(jié)權(quán)值控制額外引入的數(shù)值耗散的大小。后續(xù)該混合正則化模型被Feng 等[52-54]用到了他們的求解器中。在他們的方法中，質(zhì)量和動(dòng)量方程由格子玻爾茲曼方法在D3Q19和D3Q27 格子上進(jìn)行求解，能量方程則通過有限體積方法進(jìn)行求解，并通過修改平衡態(tài)分布函數(shù)和添加源項(xiàng)的方式盡可能地消除標(biāo)準(zhǔn)格子存在的速度離散誤差，實(shí)現(xiàn)了跨聲速和超聲速模擬。2019 年，Shan[55]在中心矩空間對(duì)非平衡態(tài)和碰撞項(xiàng)分別進(jìn)行了展開，并在中心矩空間對(duì)每階矩進(jìn)行獨(dú)立松弛，然后通過中心矩空間和原點(diǎn)矩空間Hermite 基之間的關(guān)系將碰撞項(xiàng)轉(zhuǎn)換回原點(diǎn)矩空間，得到了遞歸形式的碰撞項(xiàng)。該碰撞模型繼承了2007 年Shan 等[47]提出的原點(diǎn)矩空間Hermite 矩多松弛模型的優(yōu)點(diǎn)，克服了Chen 等[46]提到的熱傳導(dǎo)系數(shù)的伽利略不變性問題，并可以推廣到任意高階松弛。隨后，Li 等[56]在溫度標(biāo)度的中心矩空間中，對(duì)非平衡態(tài)和碰撞項(xiàng)進(jìn)行了展開，并在該空間中對(duì)各階矩進(jìn)行松弛，利用溫度標(biāo)度中心矩空間和原點(diǎn)矩空間Hermite 基之間的關(guān)系將碰撞項(xiàng)轉(zhuǎn)換回原點(diǎn)矩空間，也得到了遞歸形式的碰撞項(xiàng)；與原點(diǎn)矩空間松弛的碰撞項(xiàng)[47]相比，三階及以上的碰撞項(xiàng)包含速度、溫度及低階矩松弛構(gòu)成的修正項(xiàng)。

關(guān)于正則化碰撞模型的邊界條件處理，歷史上也有一些學(xué)者取得了部分進(jìn)展。2007 年，Niu 等[57]模擬了微尺度稀薄氣體流動(dòng)，討論了散射邊界條件、松弛時(shí)間和正則化碰撞模型。2008 年，Latt 等[58]基于正則化重構(gòu)非平衡態(tài)的思路發(fā)展了適用于平直邊界的正則化邊界條件。2011 年，Malaspinas 等[59]提出了一個(gè)通用的處理平直邊界的正則化邊界條件，通過建立邊界點(diǎn)已知分布函數(shù)與宏觀量的超定方程組，解出宏觀矩，再重構(gòu)出邊界點(diǎn)的非平衡態(tài)分布函數(shù)。該邊界處理方法可適用于標(biāo)準(zhǔn)格子（D2Q9、D3Q19 等），也適用于高階格子（D2V37、D3V103 等），但是否適用于任意曲線曲面邊界需進(jìn)一步研究。2017 年，Wissocq 等[60]發(fā)展了適用于聲學(xué)模擬的正則化特征邊界條件。2019 年，Guo 等[61]在混合正則化模型中討論了可壓縮流動(dòng)的任意曲面邊界條件和域邊界特征邊界條件。

在將正則化碰撞模型應(yīng)用到多相流方面，也有一些進(jìn)展。Otomo 等將正則化碰撞模型分別應(yīng)用到多組分偽勢(shì)模型[62]和相場(chǎng)多相流模型[63]。Ba 等[64]將正則化碰撞模型應(yīng)用到顏色梯度模型。另外，將正則化碰撞模型用于噪聲模擬也有很多代表性工作，如Zhuo 等[65]、Chen 等[66]、Casalino 等[67-68]、Gendre 等[69]、Astoul 等[70]。本文對(duì)正則化碰撞模型在各個(gè)領(lǐng)域的應(yīng)用不再做詳細(xì)展開。

本文將重點(diǎn)闡述、回顧正則化碰撞模型，通過一個(gè)系統(tǒng)的理論分析框架闡述作者與合作者提出的高階正則化碰撞模型，并對(duì)以往的正則化碰撞模型的理論關(guān)系進(jìn)行梳理，同時(shí)對(duì)正則化多松弛模型的本質(zhì)進(jìn)行歸納分析。后文將按以下順序進(jìn)行論述：第1 節(jié)闡述Hermite 展開的基礎(chǔ)理論框架，為后文中的碰撞算子正則化提供分析的準(zhǔn)則；第2 節(jié)則利用第1 節(jié)建立的理論體系，對(duì)正則化碰撞模型發(fā)展歷史上有重要貢獻(xiàn)的經(jīng)典模型逐個(gè)進(jìn)行分析、梳理和比較；最后對(duì)全文進(jìn)行歸納總結(jié)，得出全文理論分析的結(jié)論。

1 理論基礎(chǔ)

該部分我們從Boltzmann 方程出發(fā)，在Hermite譜空間對(duì)分布函數(shù)進(jìn)行展開，分別在原點(diǎn)矩空間、中心矩空間和溫度標(biāo)度的中心矩空間對(duì)Maxwell-Boltzmann 平衡態(tài)、非平衡態(tài)、碰撞項(xiàng)和力項(xiàng)進(jìn)行展開，為后面的正則化理論模型回顧做準(zhǔn)備。

我們先從Boltzmann 方程進(jìn)行闡述。Boltzmann方程為基于分子二體碰撞、分子混沌假設(shè)和外力不影響局部碰撞的動(dòng)力學(xué)行為等三個(gè)假設(shè)，推導(dǎo)的稀薄氣體系統(tǒng)的控制方程。它描述了統(tǒng)計(jì)粒子分布函數(shù)在速度空間和物理空間的演化，其右端的非線性碰撞項(xiàng)極為復(fù)雜，直接求解困難重重。Boltzmann-BGK 方程[14]作為Boltzmann 方程碰撞項(xiàng)的一個(gè)極度簡(jiǎn)化版本，認(rèn)為所有動(dòng)理學(xué)矩都以相同的速度朝平衡態(tài)弛豫，其控制方程如下：

1.1 粒子分布函數(shù)的Hermite 展開

1.2 Hermite 展開的截?cái)嗪退俣瓤臻g離散

在實(shí)際的數(shù)值執(zhí)行中，粒子分布函數(shù)不可能展開到無窮階，必須進(jìn)行截?cái)?。截?cái)嗟碾A數(shù)直接影響水動(dòng)力學(xué)宏觀方程的適用范圍及其精度，這意味著實(shí)際上需要將分布函數(shù)投影到有限階Hermite 基上去。對(duì)于截?cái)嚯A數(shù)如何影響恢復(fù)的水動(dòng)力學(xué)方程層級(jí)以及其

受益于Hermite 多項(xiàng)式的正交屬性，粒子分布函數(shù)的前若干階矩是可以由截?cái)嘈问降腍ermite 展開得到的粒子分布函數(shù)獲得。因此，粒子分布函數(shù)可以投影到由前N階Hermite 多項(xiàng)式構(gòu)成的Hilbert 空間，并保證其前N階矩與完整粒子分布函數(shù)的前N階矩嚴(yán)格相等。粒子分布函數(shù)的前N階Hermite 展開截?cái)酁椋?/p>

1.3 不同Hermite 基之間的數(shù)學(xué)聯(lián)系

為了便于后文對(duì)不同空間中Hermite 系數(shù)進(jìn)行評(píng)估，對(duì)于不同Hermite 基之間的轉(zhuǎn)換關(guān)系做一個(gè)簡(jiǎn)單介紹，詳細(xì)的推導(dǎo)過程可以參考Shan[55,72]、Li 等[56]、Shan 等[72]的論文附錄部分，他們給出了詳細(xì)的數(shù)學(xué)推導(dǎo)，這里只給出前四階的數(shù)學(xué)轉(zhuǎn)換關(guān)系。由于在實(shí)際的數(shù)值執(zhí)行中，中心矩空間的離散速度和權(quán)值會(huì)隨當(dāng)?shù)氐暮暧^速度而變化；而溫度標(biāo)度中心矩空間的離散速度和權(quán)值則會(huì)隨當(dāng)?shù)氐暮暧^速度和溫度發(fā)生變化。這會(huì)破壞格子玻爾茲曼方法完美的遷移操作，而引入額外的數(shù)值耗散和色散，讓格子玻爾茲曼方法的優(yōu)點(diǎn)喪失。因此，我們的思路是在中心矩空間和溫度標(biāo)度中心矩空間完成碰撞，然后通過它們與原點(diǎn)矩空間的數(shù)學(xué)關(guān)系轉(zhuǎn)換到原點(diǎn)矩空間，最終在原點(diǎn)矩空間實(shí)現(xiàn)粒子遷移操作。

首先寫出原點(diǎn)矩Hermite 基的零到四階表達(dá)式，如下（具體可參考文獻(xiàn)[55]）：

1.4 不同Hermite 基級(jí)數(shù)展開及其系數(shù)評(píng)估

實(shí)際上，歷史上發(fā)展的不同的正則化模型的差異主要在以下幾個(gè)地方：一個(gè)是使用的Hermite 基不同，即在不同的空間（原點(diǎn)矩空間、中心矩空間和溫度標(biāo)度中心矩空間）中展開；另一個(gè)是分布函數(shù)Hermite 展開截?cái)嗟碾A數(shù)的差異，尤其是非平衡態(tài)Hermite 展開的截?cái)啵贿€有一個(gè)是碰撞項(xiàng)展開的空間和碰撞項(xiàng)中弛豫的張量矩分量。因此，為了對(duì)這些正則化碰撞模型進(jìn)行系統(tǒng)地回顧，我們將平衡態(tài)、非平衡態(tài)、碰撞項(xiàng)和力項(xiàng)都用三種不同的Hermite 基進(jìn)行展開，然后對(duì)它們進(jìn)行系統(tǒng)地分析。

1.4.1 平衡態(tài)分布函數(shù)的Hermite 展開

使用權(quán)函數(shù)定義表達(dá)式（5），Maxwell-Boltzmann平衡態(tài)分布函數(shù)可寫為下式：

利用中心矩Hermite 基和原點(diǎn)矩Hermite 基的數(shù)學(xué)轉(zhuǎn)換關(guān)系式（22），上式中各階平衡態(tài)中心矩Hermite 系數(shù)可相應(yīng)得到，這里顯式地寫出零到四階：

其中二階和三階非平衡態(tài)系數(shù)可通過Chapman-Enskog 展開[73]得到：

利用溫度標(biāo)度中心矩Hermite 基與原點(diǎn)矩Hermite 基的數(shù)學(xué)關(guān)系，我們可以計(jì)算上式，顯式地寫出零到四階：

上面的推導(dǎo)中用到了Maxwell-Boltzmann 平衡態(tài)分布函數(shù)，且推導(dǎo)式（63）中用到了等溫假設(shè)。式（62）如果在速度空間進(jìn)行離散且考慮離散格子效應(yīng)，即為Guo 等[76]2002 年推導(dǎo)的力項(xiàng)格式（具體離散在后文中闡述）；式（63）即為He 等[75]1998 年設(shè)計(jì)的力項(xiàng)格式；式（64）與式（63）等價(jià)。

其中多弛豫時(shí)間被應(yīng)用，不是原始的Boltzmann-BGK 方程中的單松弛因子，而是在每階矩應(yīng)用了不同的松弛因子，可視為其拓展模型。在我們的工作中，上式的約束條件被拓展到非平衡態(tài)分布函數(shù)的Hermite 矩松弛，這必須與原始的碰撞項(xiàng)的松弛在討論的階數(shù)上分別等價(jià)，即：

將上式代入式（67），并將碰撞項(xiàng)和非平衡態(tài)分布函數(shù)的Hermite 展開式也代入，并利用Hermite 多項(xiàng)式的正交特性，可得到下面關(guān)系式：

點(diǎn)評(píng)：從式（75）可以看出，不管是Hermite 原點(diǎn)矩、中心矩還是溫度標(biāo)度中心矩，每階都可以獨(dú)立松弛。然而從本文1.4.2 節(jié)對(duì)非平衡態(tài)Hermite 系數(shù)在不同空間之間的數(shù)學(xué)關(guān)系的分析，我們可以發(fā)現(xiàn)，二階矩在不同空間中是等價(jià)的，只要離散速度足夠多，能保證二階矩的積分精度，分子速度空間離散形式的非平衡態(tài)分布函數(shù)的二階矩就會(huì)等價(jià)于連續(xù)形式非平衡態(tài)分布函數(shù)的二階矩。Shan 等[45,71]已證明對(duì)于二維至少需要17 個(gè)離散速度，三維至少需要39 個(gè)離散速度。只要離散速度充分，同時(shí)平衡態(tài)展開到三階及以上，不論在哪個(gè)空間進(jìn)行松弛，黏性輸運(yùn)系數(shù)的伽利略不變性都能得到嚴(yán)格滿足[55]。使用D2Q9 或者D3Q27 離散格子級(jí)聯(lián)碰撞模型[17]或者累積量碰撞模型[27]都不能保證嚴(yán)格滿足黏性輸運(yùn)系數(shù)的伽利略不變性。在不同的平移坐標(biāo)系中，宏觀速度是不一樣的，且是有量綱的，在溫度標(biāo)度中心矩空間進(jìn)行碰撞，各階矩是與宏觀速度無關(guān)且被局部聲速無量綱化的量。從式（82）～式（84）即可看出，在溫度標(biāo)度中心矩空間進(jìn)行松弛，松弛的量是真正的非平衡態(tài)量，在原點(diǎn)矩松弛對(duì)象中將守恒量（動(dòng)量、內(nèi)能）剔除，從而保證了每階矩松弛的真正獨(dú)立。結(jié)合式（76）、式（82）～式（84）我們可以發(fā)現(xiàn)，原點(diǎn)矩空間對(duì)每階非平衡態(tài)Hermite 矩進(jìn)行松弛，實(shí)質(zhì)上在三階及以上矩松弛中與低階矩發(fā)生了相互干擾，即所謂的串?dāng)_效應(yīng)。這也說明原點(diǎn)矩空間中的各階矩松弛并不真正的獨(dú)立。因此，在碰撞算子中，松弛對(duì)象必須滿足坐標(biāo)的平移不變性。Li 等[77]最近的工作討論了松弛對(duì)象必須滿足坐標(biāo)的旋轉(zhuǎn)不變性。本綜述中暫不作展開。

2 碰撞算子的正則化

2.1 格子玻爾茲曼方程的時(shí)空離散

包含力項(xiàng)的速度相空間離散形式的格子玻爾茲曼方程為：

將新定義變量代入式（86），可得到如下顯式方程，右端項(xiàng)均為當(dāng)前時(shí)間步的量：

特別注意地，當(dāng)不包含力項(xiàng)時(shí)，碰撞方程從二階開始，因?yàn)橘|(zhì)量和動(dòng)量守恒。但是對(duì)于新定義的分布函數(shù)，如果包含力項(xiàng)，則一階非平衡態(tài)矩不為零，這在多相流[62-63]和浸沒邊界法[11]包含力項(xiàng)時(shí)需要注意。先對(duì)新定義的粒子分布函數(shù)所求得的宏觀量與實(shí)際物理量之間的關(guān)系做一個(gè)討論。

對(duì)于原點(diǎn)矩空間碰撞，平衡態(tài)部分、非平衡態(tài)部分以及力項(xiàng)均可以直接利用Hermite 展開的截?cái)嚯A形式重構(gòu)粒子分布函數(shù)。而對(duì)于中心矩和溫度標(biāo)度中心矩空間的碰撞形式，由于離散速度依賴于當(dāng)?shù)厮俣群蜏囟?，我們選擇將其轉(zhuǎn)換到原點(diǎn)矩空間，再進(jìn)行粒子分布函數(shù)重構(gòu)。根據(jù)矩空間碰撞執(zhí)行的空間不同以及截?cái)嚯A的階數(shù)，我們可以對(duì)已存在的正則化模型進(jìn)行系統(tǒng)地分析和討論。為方便起見，后文中非平衡態(tài)矩上面的一拔標(biāo)記全部省略。

2.2 正則化碰撞模型分析與討論

2.2.1 Ladd 正則化模型和Latt 正則化模型

Ladd 等在文獻(xiàn)[39]（1994 年）和文獻(xiàn)[41]（2001 年）中提出了正則化碰撞模型的最早形式。在他們的碰撞模型中，平衡態(tài)分布函數(shù)展開到二階，非平衡態(tài)分布函數(shù)也僅保留到二階。碰撞后的分布函數(shù)寫為下面形式：

Latt 等[42]在2005 年也提出了正則化模型，他們的模型并未分離二階應(yīng)力張量的無跡部分和有跡部分的獨(dú)立松弛，僅為

Latt 的模型在Ladd 模型之后很久才提出，而且可視為L(zhǎng)add 模型的一個(gè)特例，即剪切和體積粘性松弛因子相等，本質(zhì)上不能算一個(gè)新的模型。Ladd 模型為原點(diǎn)矩空間的正則化碰撞模型，高于三階的超出等溫Navier-Stokes 方程的非平衡態(tài)矩可視為松弛因子都設(shè)置為1，對(duì)碰撞后的粒子分布函數(shù)不產(chǎn)生貢獻(xiàn)。

2.2.2 Zhang-Shan-Chen 正則化模型

2006 年，Zhang、Shan 和Chen[43]基 于Hermite 展開提出了高階正則化模型，并用來模擬高Kn數(shù)非連續(xù)流動(dòng)。他們認(rèn)為平衡態(tài)分布函數(shù)投影到有限階Hermite 基構(gòu)成的Hilbert 空間，非平衡態(tài)部分也應(yīng)該投影到該Hilbert 空間上，并引入力項(xiàng)。根據(jù)本文2.1 節(jié)的討論，在考慮外力項(xiàng)后，非平衡態(tài)Hermite 系數(shù)的一階項(xiàng)并不為零，因此正則化的碰撞后粒子分布函數(shù)應(yīng)該寫為：

2.2.3 Mattila-Philippi-Hegele 正則化模型

Mattila 等[49]提出了一個(gè)基于中心矩空間展開的高階格子正則化模型來處理可壓縮傳熱流動(dòng)，我們利用本文1.4.4 節(jié)的有關(guān)碰撞項(xiàng)和1.4.2 節(jié)的非平衡態(tài)系數(shù)在中心矩空間和原點(diǎn)矩空間分別展開的數(shù)學(xué)關(guān)系，分析該正則化模型。該模型中不包含力項(xiàng)，因此我們也省略力項(xiàng)。為分析方便，我們把非平衡態(tài)Hermite系數(shù)和碰撞項(xiàng)合并到一起，如式（95）所示，在式（77）中引入新的碰撞形式，包含非平衡的中心矩Hermite 系數(shù)，具體如下：

其中四階及以上非平衡態(tài)松弛因子設(shè)為1，處理可壓縮傳熱流動(dòng)。即認(rèn)為超出二階應(yīng)力張量和三階熱流矢量的高階動(dòng)理學(xué)中心矩快速松弛，碰撞后的粒子分布函數(shù)直接為平衡態(tài)。利用Hermite 中心矩和原點(diǎn)矩之間的關(guān)系可整理上式為：

展開基跟平衡態(tài)展開基對(duì)應(yīng)，二階非平衡態(tài)Hermite 系數(shù)通過非平衡態(tài)粒子分布函數(shù)的二階矩求得，三階和四階非平衡態(tài)Hermite 系數(shù)則遞歸求得。

相應(yīng)地，三維（D3Q27）的平衡態(tài)展開零到二階跟二維一樣，但高階一直展開到六階，一共27 個(gè)Hermite 基，見式（129）：

最后的正則化碰撞和遷移格式為：

2.2.5 Li-Shi-Shan 正則化模型

Li 等[56]在2019 年提出了一個(gè)基于溫度標(biāo)度中心矩空間展開的高階格子正則化模型來處理可壓縮傳熱流動(dòng)。我們利用本文1.4.4 節(jié)的有關(guān)碰撞項(xiàng)和1.4.2 節(jié)的非平衡態(tài)系數(shù)在溫度標(biāo)度中心矩空間和原點(diǎn)矩空間分別展開的數(shù)學(xué)關(guān)系來分析該正則化模型。為了跟本文2.2.3 節(jié)分析保持一致，我們把溫度標(biāo)度中心矩空間的非平衡態(tài)Hermite 系數(shù)及其松弛進(jìn)行合并，在式（81）的基礎(chǔ)上定義新的碰撞形式：

因此，Li 等提出的基于溫度標(biāo)度中心矩空間中Hermite 展開的碰撞形式可視為將該空間中四階及以上的非平衡態(tài)Hermite 矩的松弛因子置為1，該空間中的高階矩快速松弛到平衡態(tài)。比較Zhang-Shan-Chen 在原點(diǎn)矩空間的正則化碰撞模型、Mattila-Philippi-Hegele 在中心矩空間的正則化碰撞模型，以及作者等的Li-Shi-Shan 在溫度標(biāo)度中心矩空間的正則化碰撞模型，三者的相同點(diǎn)是將四階及以上的相應(yīng)非平衡態(tài)高階矩的松弛因子置為1，讓這些高階非平衡態(tài)動(dòng)力學(xué)矩快速松弛到平衡態(tài)，而本質(zhì)差異是進(jìn)行碰撞的空間不同。這些模型的聯(lián)系和差異，在本文的分析框架下變得十分明朗。另外，在溫度標(biāo)度中心矩空間的正則化模型，在包含力項(xiàng)后的形式，我們將在另外的論文中發(fā)表。

2.2.6 Coreixas 等正則化模型

2017 年，Coreixas 等[50]將Malaspinas 提出的等溫層級(jí)的遞歸正則化模型推廣到了可壓縮傳熱層級(jí)上，發(fā)展了一個(gè)高階格子遞歸正則化模型。該正則化模型包含平衡態(tài)Hermite 系數(shù)的遞歸關(guān)系式：

上述遞歸關(guān)系式十分復(fù)雜，其數(shù)學(xué)推導(dǎo)過程也極其復(fù)雜。我們對(duì)四階非平衡態(tài)Hermite 系數(shù)進(jìn)行整理如下：

與Li-Shi-Shan 正則化碰撞模型式（132）～式（134）相比，Coreixas 等的高階遞歸正則化模型即使使用了多松弛模型，也不能保證熱傳導(dǎo)系數(shù)的伽利略不變性；而Li-Shi-Shan 正則化碰撞模型不僅可以得到非平衡態(tài)Hermite 系數(shù)的遞歸形式，還能得到碰撞項(xiàng)的遞歸形式。這充分體現(xiàn)了對(duì)溫度標(biāo)度中心矩空間進(jìn)行展開和正則化在數(shù)學(xué)和物理基本原理上的優(yōu)勢(shì)。

3 結(jié) 論

本文對(duì)格子玻爾茲曼領(lǐng)域的正則化碰撞模型進(jìn)行了系統(tǒng)的嚴(yán)謹(jǐn)?shù)睦碚摶仡?，通過采用不同自變量的Hermite 基對(duì)平衡態(tài)、非平衡態(tài)、力項(xiàng)和碰撞項(xiàng)進(jìn)行展開，將所有的理論模型在該理論體系下的定位以及它們之間的聯(lián)系進(jìn)行了分析。分析發(fā)現(xiàn)：

1） 每一種正則化碰撞模型都可以在Hermite 展開的理論框架下得到。

2） Malaspinas 和Coreixas 的遞歸正則化模型，非平衡態(tài)可以無限遞歸得到，對(duì)于超出所研究的動(dòng)理學(xué)矩（二階應(yīng)力張量矩和三階熱流矢量矩），其弛豫時(shí)間不一定為1。前者由與剪切粘性相關(guān)的弛豫時(shí)間統(tǒng)一松弛，后者則對(duì)每一階非平衡態(tài)原點(diǎn)矩進(jìn)行獨(dú)立松弛。高階矩的松弛因子如何影響模型的數(shù)值穩(wěn)定性和精度有待進(jìn)一步探究。

3） 除2）中提到的遞歸正則化模型，其他的正則化模型等價(jià)于在原點(diǎn)矩空間、中心矩空間和溫度標(biāo)度中心矩空間，將高于所關(guān)注的非平衡態(tài)動(dòng)理學(xué)矩的弛豫時(shí)間置為1。

4） 可模擬可壓縮傳熱NSF 方程的正則化模型中，只有Chen 等[46]的模型、 Shan 的模型[55]、Li-Shi-Shan 模型[56]是在中心矩或溫度標(biāo)度中心矩空間進(jìn)行松弛，滿足熱傳導(dǎo)系數(shù)的伽利略不變性。

正則化多松弛碰撞模型的理論本質(zhì)就是將離散速度表達(dá)不了的高階動(dòng)理學(xué)矩進(jìn)行截?cái)?，或者等價(jià)于將高階動(dòng)理學(xué)矩的松弛因子置為1，以避免虛假模態(tài)對(duì)水動(dòng)力學(xué)方程的負(fù)面影響；而高階動(dòng)理學(xué)矩的松弛因子不為1 時(shí)，對(duì)數(shù)值穩(wěn)定性的正面影響目前尚不明朗。展望未來，正則化模型相較于其他碰撞模型的優(yōu)勢(shì)還有待嚴(yán)格地理論證明，尤其是在湍流模擬、聲學(xué)模擬、多相流以及稀薄氣體流動(dòng)等領(lǐng)域，還有很大的發(fā)展空間。