格子玻爾茲曼和氣體動(dòng)理學(xué)通量算法及其應(yīng)用進(jìn)展

舒 昌，楊鯉銘，王 巖，吳 杰

（1.新加坡國(guó)立大學(xué)機(jī)械工程系，新加坡 119260，新加坡；2.南京航空航天大學(xué)航空學(xué)院，南京 210016；3.南京航空航天大學(xué)非定?？諝鈩?dòng)力學(xué)與流動(dòng)控制工信部重點(diǎn)實(shí)驗(yàn)室，南京 210016；4.南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210016）

隨著計(jì)算機(jī)性能的不斷提升，計(jì)算流體力學(xué)（Computational fluid dynamics，CFD）已經(jīng)成為解決流體流動(dòng)問(wèn)題的重要手段之一，在航空、航天、航海等領(lǐng)域均已取得了非常成功的應(yīng)用［1-4］。與理論流體力學(xué)和實(shí)驗(yàn)流體力學(xué)相比，CFD 的應(yīng)用場(chǎng)景更為靈活和廣泛，且大多數(shù)情況下更為經(jīng)濟(jì)。特別是對(duì)于實(shí)驗(yàn)難于復(fù)現(xiàn)的場(chǎng)景以及不便測(cè)量的區(qū)域，例如臨近空間的高速稀薄流動(dòng)和航空發(fā)動(dòng)機(jī)內(nèi)部的精細(xì)流場(chǎng)結(jié)構(gòu)，CFD 更是發(fā)揮著不可替代的作用。CFD 通過(guò)求解流體流動(dòng)控制方程來(lái)獲得流場(chǎng)的解，常用的方法包括有限差分法（Finite difference method，F(xiàn)DM）［5-6］、有 限 體 積 法（Finite volume method，F(xiàn)VM）［7-8］以 及 有 限 元 法（Finite element method，F(xiàn)EM）［9-10］等。在這些方法中，由于離散控制方程的解通常定義在解點(diǎn)上，如何利用解點(diǎn)上的值來(lái)計(jì)算單元界面的通量，是CFD 最為關(guān)心的問(wèn)題之一，該過(guò)程也被稱為通量重構(gòu)。

光滑函數(shù)近似是最為直觀和簡(jiǎn)單的通量重構(gòu)方法，該方法采用光滑函數(shù)將單元界面周圍解點(diǎn)上的物理量連接起來(lái)，由此界面上的物理量便可通過(guò)該光滑函數(shù)插值計(jì)算。由于沒(méi)有考慮流動(dòng)的特征，該類算法可以視為完全的數(shù)學(xué)重構(gòu)。其典型代表為中心格式，即采用單元界面左右兩側(cè)解點(diǎn)的平均值來(lái)重構(gòu)界面通量。雖然該方法極為簡(jiǎn)單，但在求解含有間斷的流動(dòng)問(wèn)題時(shí)，會(huì)導(dǎo)致計(jì)算不穩(wěn)定的現(xiàn)象。克服該問(wèn)題的常用辦法是在計(jì)算通量過(guò)程中增加人工黏性，例如Jameson-Schmidt-Turkel（JST）格式通過(guò)增加2 階和4 階人工黏性項(xiàng)來(lái)改善中心格式的數(shù)值穩(wěn)定性［11-12］。其中，2 階黏性項(xiàng)用于抑制激波附近的振蕩，僅在流場(chǎng)中壓力梯度較大的區(qū)域發(fā)揮作用，而4 階黏性項(xiàng)用于抑制高頻振蕩，使數(shù)值解趨于穩(wěn)定。

不同于數(shù)學(xué)重構(gòu)，物理重構(gòu)在計(jì)算單元界面通量時(shí)會(huì)考慮流動(dòng)的信息或者通過(guò)求解等效的物理方程（或簡(jiǎn)化的物理方程）來(lái)獲得界面的物理量。著名的通量矢量分裂格式，例如van Leer 格式［13-14］和AUSM 格式［15-16］均是依據(jù)波的傳播方向來(lái)計(jì)算無(wú)黏通量。具體地，van Leer 格式依據(jù)當(dāng)?shù)伛R赫數(shù)的符號(hào)直接將無(wú)黏通量分解為左通量和右通量?jī)刹糠?，而AUSM 格式則先將無(wú)黏通量拆分為對(duì)流項(xiàng)和壓力項(xiàng)，然后再依據(jù)當(dāng)?shù)伛R赫數(shù)的符號(hào)對(duì)各自進(jìn)行分解。與通量矢量分裂格式不同，Godunov 格式［17-18］、HLL 格 式［19-20］、HLLC 格 式［21-22］和Roe 格式［23-24］等則利用一維Riemann 問(wèn)題的精確解或近似解來(lái)重構(gòu)單元界面的通量。由于從界面左側(cè)解點(diǎn)和右側(cè)解點(diǎn)插值到界面上的左狀態(tài)和右狀態(tài)通常不相等，會(huì)在界面形成Riemann 問(wèn)題。Godunov格式正是通過(guò)求解一維Riemann 問(wèn)題的解析解來(lái)計(jì)算無(wú)粘通量，HLL 格式和HLLC 格式則是利用一維Riemann 問(wèn)題的近似解來(lái)計(jì)算無(wú)黏通量，而Roe 格式通過(guò)求解近似Riemann 問(wèn)題的解析解來(lái)獲得無(wú)黏通量的表達(dá)式。

無(wú)論是采用數(shù)學(xué)重構(gòu)還是物理重構(gòu)，上述這些通量格式僅考慮了無(wú)黏通量的計(jì)算，黏性通量則是通過(guò)中心差分來(lái)獲得。另外，由于采用一維Riemann 問(wèn)題的精確解或近似解來(lái)重構(gòu)單元界面的通量，Godunov 格式、HLL 格式、HLLC 格式和Roe 格式等只能沿單元界面法向應(yīng)用，切向通量則需要通過(guò)被動(dòng)標(biāo)量的方式來(lái)計(jì)算。由此可見(jiàn)，這些格式在單元界面上并不能精確滿足Navier-Stokes 方程，亦或是多維的Euler 方程。其原因在于，多維Riemann 問(wèn)題的求解非常困難，需要考慮的解的種類非常多，從而導(dǎo)致計(jì)算成本顯著增加。

與上述通量重構(gòu)的思路不同，氣體動(dòng)理學(xué)格式（Gas kinetic scheme，GKS）通過(guò)在單元界面應(yīng)用Boltzmann 模型方程的局部解來(lái)計(jì)算宏觀Navier-Stokes 方程的通量［25-30］，實(shí)現(xiàn)了通量的完全物理重構(gòu)。在連續(xù)介質(zhì)假設(shè)的前提下，Boltzmann模型方程可以還原回Navier-Stokes 方程，該關(guān)系為采用Boltzmann 模型方程的局部解來(lái)計(jì)算Navier-Stokes 方程的通量奠定了基礎(chǔ)。由于Boltzmann 模型方程僅為分布函數(shù)的單變量方程，形式比Navier-Stokes 方程更簡(jiǎn)單，因而可以方便地獲得該方程的局部解。 GKS 正是采用Boltzmann 模型方程在單元界面的局部積分解結(jié)合對(duì)應(yīng)于Navier-Stokes 方程截?cái)嘈问降某跏挤植己瘮?shù)來(lái)計(jì)算Navier-Stokes 方程的通量，從而使得單元界面的物理量也滿足Navier-Stokes 方程的等價(jià)形式，并且無(wú)黏通量和黏性通量可以采用統(tǒng)一的方式計(jì)算。另一方面，GKS 由于在求解過(guò)程中引入了許多參數(shù)，其計(jì)算效率要低于傳統(tǒng)的Navier-Stokes 方程求解器。

本文所要介紹的格子玻爾茲曼通量算法（Lattice Boltzmann flux solver，LBFS）和氣體動(dòng)理學(xué)通量算法（Gas kinetic flux solver，GKFS）與GKS 類似，也是借助于Boltzmann 模型方程的局部解來(lái)計(jì)算宏觀Navier-Stokes 方程的通量。但與GKS 不同，LBFS 和GKFS 采 用的是Boltzmann 模型方程的漸進(jìn)展開(kāi)解［30-35］，其形式比GKS 采用的局部積分解更為簡(jiǎn)單而且其計(jì)算量和傳統(tǒng)的Navier-Stokes 方程求解器相當(dāng)。該漸進(jìn)展開(kāi)解的不同階截?cái)鄬?duì)應(yīng)于不同層次的宏觀控制方程，其中采用一階截?cái)嗉纯蛇€原Navier-Stokes 方程。利用漸進(jìn)展開(kāi)解的一階截?cái)嘟Y(jié)合不同的平衡態(tài)分布函數(shù)，即可發(fā)展相應(yīng)的LBFS 和GKFS 用于宏觀控制方程的求解。除了應(yīng)用于連續(xù)流動(dòng)問(wèn)題，GKFS 也可拓展到稀薄流動(dòng)問(wèn)題求解［36-38］。但在稀薄流動(dòng)問(wèn)題求解時(shí)，單元界面的初始分布函數(shù)不能再采用平衡態(tài)來(lái)近似，需要替換為考慮非平衡效應(yīng)的分布函數(shù)，例如取為Boltzmann 模型方程的解或Grad 分布函數(shù)。將單元界面的初始分布函數(shù)直接取為Boltzmann 模型方程的解可以實(shí)現(xiàn)全流域流動(dòng)問(wèn)題的求解，但該方式需要在分子速度空間離散求解Boltzmann 模型方程。通過(guò)Grad 分布函數(shù)來(lái)近似單元界面的初始分布函數(shù)可以避開(kāi)Boltzmann 模型方程的求解，提高計(jì)算效率并降低內(nèi)存需求，但該方式適用的克努森數(shù)范圍有限，一般與矩方法相當(dāng)。鑒于GKS、LBFS 和GKFS 均采用完全物理重構(gòu)的方式來(lái)計(jì)算單元界面的通量，因而通量也可視為滿足宏觀控制方程的等價(jià)形式。大量數(shù)值實(shí)驗(yàn)表明，該類算法從低速到高超聲速流動(dòng)問(wèn)題求解均具有非常好的穩(wěn)定性［39-40］，因而已被擴(kuò)展應(yīng)用于多相流、動(dòng)邊界、流固共軛傳熱以及化學(xué)反應(yīng)流等問(wèn)題的求解。

1 Boltzmann 模型方程和宏觀守恒律方程

1.1 Boltzmann 模型方程及其漸進(jìn)展開(kāi)解

Boltzmann 方程通過(guò)引入氣體分子速度分布函數(shù)f來(lái)描述微觀系統(tǒng)中分子的遷移和碰撞過(guò)程。原始Boltzmann 方程的碰撞項(xiàng)為分布函數(shù)的五重積分，而且與分子的模型和碰撞截面積相關(guān)，因而表達(dá)式極為復(fù)雜，難以在實(shí)際工程中直接應(yīng)用。為了簡(jiǎn)化該方程的碰撞項(xiàng)，相繼提出了各種簡(jiǎn)化模型方程［41-43］，其中BGK 模型由于其形式最為簡(jiǎn)單而獲得了廣泛的應(yīng)用。BGK 模型由Bhatnagar、Gross和Krook 提出，它可以滿足質(zhì)量、動(dòng)量和能量守恒，滿足熵增條件，并且導(dǎo)出的平衡態(tài)即為Maxwell 分布。但是，該模型對(duì)應(yīng)的普朗特?cái)?shù)恒為1，與通過(guò)原始Boltzmann 方程推導(dǎo)得到的正確值2/3 不一致，因此在實(shí)際應(yīng)用中需要對(duì)其進(jìn)行修正。采用BGK 模型作為碰撞項(xiàng)的Boltzmann 方程可以改寫為

式中：t為時(shí)間；ξ為分子速度；?為空間導(dǎo)數(shù)；τ為分子碰撞時(shí)間，對(duì)于硬球模型分子可以表示為

式中：μ和p分別為動(dòng)力學(xué)黏性和壓力；Ma、Re和Kn分別為馬赫數(shù)、雷諾數(shù)和克努森數(shù)；γ為比熱容比；u∞為參考速度；L∞為參考長(zhǎng)度；c∞為聲速。另外，式（1）中feq為Maxwell 平衡態(tài)分布函數(shù)，即

式中：c=ξ-u為分子的最可幾熱運(yùn)動(dòng)速度，u為宏觀速度；ρ為密度；T為溫度；R為氣體常數(shù)。除特殊說(shuō)明外，本文中約定用黑體來(lái)表示矢量，用對(duì)應(yīng)的白斜體來(lái)表示矢量的長(zhǎng)度。

為了便于推導(dǎo)Boltzmann-BGK 模型方程的漸進(jìn)展開(kāi)解，可將式（1）改寫為

式中Df=?t f+ξ·?f為分布函數(shù)的物質(zhì)導(dǎo)數(shù)。將上述f的表達(dá)式不斷地代入方程最后1 項(xiàng)，即可獲得Boltzmann-BGK 模型方程的漸進(jìn)展開(kāi)解，即

式中已將τ近似為局部常數(shù)，因而可以提到物質(zhì)導(dǎo)數(shù)算子前面。實(shí)際上，式（5）的不同階截?cái)嗉创聿煌暮暧^控制方程。

1.2 宏觀守恒律方程

不同于Boltzmann 方程，宏觀控制方程直接在宏觀層次上利用質(zhì)量，動(dòng)量和能量守恒規(guī)律來(lái)建立流體流動(dòng)的控制方程。由于同為描述流體問(wèn)題的控制方程，Boltzmann 方程和宏觀控制方程間存在必然的聯(lián)系。在式（1）左右兩邊同時(shí)乘以微觀守恒矩，并在分子速度空間積分，可得

式中：E為總能；ψ=(1，ξ，ξ22)T為微觀守恒矩矢量。符號(hào) · 表示在分子速度空間的積分，即

如果將分布函數(shù)f近似為其平衡態(tài)feq，則對(duì)應(yīng)的宏觀控制方程為Euler 方程。此時(shí)相當(dāng)于將τ取為0，即克努森數(shù)Kn設(shè)置為0，等效于引入了無(wú)粘假設(shè)。為了還原Navier-Stokes 方程，可以將分布函數(shù)f取為其一階截?cái)嘈问剑?4］

將平衡態(tài)分布函數(shù)式（3）代入上述積分，即可得到如下形式的Navier-Stokes 方程

并且：μb為體積黏性系數(shù)；cV為等容比熱容；cp為等壓比熱容；κ為熱傳導(dǎo)系數(shù)；Pr為普朗特?cái)?shù)。上述推導(dǎo)表明，采用Boltzmann-BGK 模型方程漸進(jìn)展開(kāi)解的一階截?cái)嗉纯蛇€原得到可壓縮Navier-Stokes 方程，并且截?cái)嗾`差為O(τ2)。實(shí)際上分布函數(shù)f的一階截?cái)嘞喈?dāng)于引入了連續(xù)性假設(shè)，即克努森數(shù)遠(yuǎn)小于1 的情形，因此O(τ2)屬于高階小量，可以忽略。另外，由于采用了BGK 碰撞模型，還原得到的Navier-Stokes 方程的普朗特?cái)?shù)恒為1。

值得指出，在還原上述可壓縮Navier-Stokes方程的過(guò)程中，僅需要用到平衡態(tài)分布函數(shù)滿足的7 個(gè)矩關(guān)系，分別對(duì)應(yīng)質(zhì)量、動(dòng)量、能量、動(dòng)量方程的無(wú)黏和黏性通量、能量方程的無(wú)黏和黏性通量。換言之，采用其他平衡態(tài)分布函數(shù)也可還原宏觀Navier-Stokes 方程，其要求僅為滿足所需的矩關(guān)系。另外，通過(guò)該推導(dǎo)過(guò)程可知，將分布函數(shù)取為Boltzmann-BGK 模型方程漸進(jìn)展開(kāi)解的一階截?cái)嗉纯蛇€原Navier-Stokes 方程。故可以將該漸進(jìn)展開(kāi)解的一階截?cái)嘀苯哟胧剑?）來(lái)設(shè)計(jì)Navier-Stokes 方程的通量計(jì)算格式。

2 格子玻爾茲曼通量算法

2.1 算法設(shè)計(jì)

建立Boltzmann-BGK模型方程與Navier-Stokes方程之間的聯(lián)系之后，即可采用該關(guān)系來(lái)設(shè)計(jì)Navier-Stokes 方程的通量計(jì)算格式。本小節(jié)采用離散的格子Boltzmann 模型來(lái)替代Maxwell 平衡態(tài)分布函數(shù)，發(fā)展基于離散模型的格子Boltzmann 通量算法。以常用的不可壓縮D2Q9 模型為例（圖1（a）），平衡態(tài)分布函數(shù)可以寫為

式 中：Vi為 控 制 體i的 體 積；N(i)為 控 制 體i所 包含鄰近單元的集合；nij為控制體i和控制體j交界面上的單位外法向矢量；Sij為該面的面積。依據(jù)Boltzmann-BGK 模型方程與上述弱可壓縮Navier-Stokes 方 程 之 間 的 聯(lián) 系，式（18）中 的 通量Fij為

式（20）中包含了tn時(shí)刻單元界面周圍的平衡態(tài) 分 布 函 數(shù)feq(xij-ξαδt，ξα，tn)和tn+δt界 面 上的 平 衡 態(tài) 分 布 函 數(shù)feq(xij，ξα，tn+δt)。 其 中，feq(xij-ξαδt，ξα，tn)可 由 相 應(yīng) 位 置 的 宏 觀 量 來(lái) 確定，而這些宏觀量可通過(guò)解點(diǎn)上的值插值得到，如圖1（b）所示。值得注意的是，單元界面周圍點(diǎn)（1～8）的位置可通過(guò)局部坐標(biāo)系來(lái)確定，該方式適用于任意網(wǎng)格。

圖1 D2Q9 模型和基于該模型的格子玻爾茲曼通量算法示意圖Fig.1 D2Q9 model and schematic diagram of D2Q9 model-based lattice Boltzmann flux solver

利用式（21）計(jì)算出(xij-ξαδt，tn)位置處的ρ和u，即可通過(guò)式（14）獲得對(duì)應(yīng)位置的平衡態(tài)分布函數(shù)。

式（22）表明，W(xij，tn+δt)可由tn時(shí)刻界面周 圍 平 衡 態(tài) 分 布 函 數(shù)feq(xij-ξαδt，ξα，tn)求 得。將W(xij，tn+δt) 代 入 式（14）即 可 獲 得feq(xij，ξα，tn+δt)。由此，單元界面上的一階截?cái)喾植己瘮?shù)f(xij，ξα，tn+δt)可完全確定，將其代入式（19）即可獲得宏觀方程通量Fij。在LBFS 中，δt與宏觀式（18）的推進(jìn)時(shí)間步長(zhǎng)Δt無(wú)關(guān)，只要xijξαδt滿足無(wú)外插即可。一種可選的方案為δt=δx=0.4×min{ΔrL，ΔrR，0.5lmin}，式中ΔrL、ΔrR和lmin分別為左側(cè)單元中心到界面的距離、右側(cè)單元中心到界面的距離和左右兩側(cè)單元的最小邊長(zhǎng)。

式（18）中的通量確定之后，便可應(yīng)用Runge-Kutta 方法或者LU-SGS 方法沿時(shí)間方向推進(jìn)求解。整體而言，LBFS仍然求解的是宏觀控制方程，只是在計(jì)算單元界面通量時(shí)采用了Boltzmann-BGK 模型方程漸進(jìn)展開(kāi)解的一階截?cái)唷Ｔ撘浑A截?cái)嗾脤?duì)應(yīng)于所需求解的宏觀Navier-Stokes方程，所以利用LBFS所構(gòu)造的通量在單元界面處也滿足所求解的控制方程，而且無(wú)粘和粘性通量可以采用一致的方式計(jì)算。相比于傳統(tǒng)CFD方法，由于僅通量計(jì)算格式發(fā)生了變化，傳統(tǒng)方法中的各種加速收斂技巧均可直接應(yīng)用于當(dāng)前格式。

2.2 不可壓縮流應(yīng)用

上述基于離散模型的LBFS 已被廣泛應(yīng)用于不 可 壓 等 溫 流 動(dòng)［45］、不 可 壓 傳 熱 流 動(dòng)［46］、多 相流［47-49］、不可壓共軛傳熱問(wèn)題［50］以及不可壓動(dòng)邊界問(wèn)題［51］等。相比于傳統(tǒng)的不可壓算法，當(dāng)前方法避免了迭代求解壓力泊松方程；相比于人工壓縮方法，當(dāng)前方法的無(wú)黏和黏性通量可以統(tǒng)一計(jì)算；相比于通過(guò)求解Boltzmann 方程來(lái)獲得連續(xù)流解的格子Boltzmann 方法，當(dāng)前方法更易用于非結(jié)構(gòu)和非均勻網(wǎng)格，并且穩(wěn)定性更好。

（1）該方法在大密度比多相流問(wèn)題中的應(yīng)用［47-49］。由于相界面附近包含密度、黏性和速度在內(nèi)的多個(gè)物理量在幾個(gè)網(wǎng)格內(nèi)急劇變化，傳統(tǒng)基于格子Boltzmann 算法模擬大密度比和高雷諾數(shù)多相流問(wèn)題一直是重要的挑戰(zhàn)，并對(duì)方法的數(shù)值穩(wěn)定性提出了苛刻的要求。與傳統(tǒng)多相流格子Boltzmann 算法不同，多相流LBFS 采用有限體積法直接求解宏觀控制方程，界面上的通量由標(biāo)準(zhǔn)的格子Boltzmann 解進(jìn)行重構(gòu)。另外，相界面的捕捉利用WENO 格式求解Cahn-Hilliard 方程來(lái)實(shí)現(xiàn)。該方法成功模擬了多個(gè)具有挑戰(zhàn)性的大密度比多相流問(wèn)題，如Rayleigh-Taylor 不穩(wěn)定性、液滴撞擊固壁以及液滴對(duì)撞和液滴撞擊液膜等，驗(yàn)證了算法的正確性。圖2 展示了雷諾數(shù)（Re）為2 000，韋 伯 數(shù)（We）為800，密 度 比（ρH/ρL）為1 000，直 徑 為55 μm 的 液 滴 以32 m/s 的 速 度 撞 擊液膜的界面演化過(guò)程［48］（T為無(wú)量綱時(shí)間）。該計(jì)算采用了501×501×261 的網(wǎng)格來(lái)精確捕捉液滴撞擊過(guò)程中”皇冠”的形成、及其上部邊緣由于不穩(wěn)定而逐漸析出多個(gè)微小液滴的狀態(tài)。

圖2 Re=2 000 時(shí)液滴撞擊液膜界面演化過(guò)程[48]Fig.2 Evolution of interface morphology for droplet splashing on a thin film at Re=2 000[48]

（2）該方法在不可壓共軛傳熱問(wèn)題中的應(yīng)用。共軛傳熱問(wèn)題相比于不可壓等溫流動(dòng)問(wèn)題，增加了用于流場(chǎng)和結(jié)構(gòu)的傳熱方程，同時(shí)動(dòng)量方程中增加了浮力項(xiàng)。因此，需要引入額外的溫度分布函數(shù)來(lái)計(jì)算傳熱方程的通量。具體細(xì)節(jié)參見(jiàn)文獻(xiàn)［52］。應(yīng)用該方法，模擬了帶肋板的三維圓環(huán)自然對(duì)流問(wèn)題，如圖3 所示。圖4 展示了該問(wèn)題在瑞利數(shù)Ra=105時(shí)的溫度云圖和速度剖面圖，該方法的計(jì)算結(jié)果與CFD 商業(yè)軟件FLUENT 的計(jì)算結(jié)果吻合良好。計(jì)算效率方面，該方法耗時(shí)5.325 h，而FLUENT 耗時(shí)5.332 h，表明其效率與成熟商業(yè)軟件相當(dāng)。

圖3 帶肋板的三維圓環(huán)自然對(duì)流問(wèn)題示意圖[52]Fig.3 Schematic diagram of natural convection in a finned 3D annulus[52]

圖4 帶肋板的三維圓環(huán)自然對(duì)流問(wèn)題在Ra = 105時(shí)的溫度云圖和速度剖面圖[52]Fig.4 Temperature contours and velocity profile for natural convection in a finned 3D annulus at Ra = 105[52]

3 氣體動(dòng)理學(xué)通量算法

3.1 算法設(shè)計(jì)

圖5 圓盤在G=100 和m*=0.1 時(shí)的擺動(dòng)模式瞬時(shí)渦系結(jié)構(gòu)[53]Fig.5 Instantaneous vortical structures for the fluttering disk at G=100 and m*=0.1[53]

圖6 圓盤在G=200 和m*=0.75 時(shí)的翻騰模式瞬時(shí)渦系結(jié)構(gòu)[53]Fig.6 Instantaneous vortical structures for the tumbling disk at G=200 and m*=0.75[53]

除了離散的平衡態(tài)分布函數(shù)，連續(xù)的平衡態(tài)分布函數(shù)也可用于構(gòu)造Navier-Stokes 方程的通量計(jì)算格式，即氣體動(dòng)理學(xué)通量算法。正如1.2 小節(jié)中的介紹，采用Boltzmann-BGK 模型方程漸進(jìn)展開(kāi)解的一階截?cái)嘟Y(jié)合Maxwell 分布函數(shù)可以還原可壓縮Navier-Stokes 方程。實(shí)際上，推導(dǎo)過(guò)程中僅使用了滿足Navier-Stokes 方程的7 個(gè)矩關(guān)系，因此Maxwell 分布函數(shù)也可以替換為其他的平衡態(tài)分布函數(shù)。基于此，作者還發(fā)展了圓函數(shù)和球函數(shù)等簡(jiǎn)化平衡態(tài)分布函數(shù)，并據(jù)此構(gòu)造了相應(yīng)的GKFS［54-55］。 不 失 一 般 性 ，該 小 節(jié) 采 用Boltzmann-BGK 模型方程漸進(jìn)展開(kāi)解的一階截?cái)嘟Y(jié)合 Maxwell 分布函數(shù)來(lái)構(gòu)造可壓縮Navier-Stokes 方程的通量計(jì)算格式。首先，將式（11～13）改寫為

與基于離散模型的LBFS 類似，采用有限體積法對(duì)式（23）在空間網(wǎng)格Vi進(jìn)行離散可得

其關(guān)鍵在于計(jì)算單元界面的一階截?cái)喾植己瘮?shù)。

在推導(dǎo)具體的通量表達(dá)式時(shí)，由于需要在速度空間積分，可以在單元界面處引入局部坐標(biāo)系以方便推導(dǎo)和應(yīng)用。以三維問(wèn)題為例，可將局部坐標(biāo)系的1 方向定義為單元界面的法向，2 方向和3 方向均為界面的切向，并且3 個(gè)坐標(biāo)方向構(gòu)成正交右手坐標(biāo)系。局部坐標(biāo)系和全局坐標(biāo)系中的通量滿足如下變換關(guān)系

式中：ψˉ=(1，ξ1，ξ2，ξ3，ξ22)T為局部坐標(biāo)系下的微觀守恒矩矢量；ξ=(ξ1，ξ2，ξ3)為局部坐標(biāo)系下的分子速度分量。

與基于離散模型的LBFS 類似，可將單元界面分布函數(shù)的一階截?cái)啾硎緸?/p>

其中界面周圍的平衡態(tài)分布函數(shù)feq(xijξδt，ξ，tn)可采用泰勒級(jí)數(shù)展開(kāi)為

式（32）方程右端項(xiàng)為左右兩側(cè)單元中心守恒量在局部坐標(biāo)系下的導(dǎo)數(shù)。由此，式（29）可以改寫為

式中H(ξ1)為臺(tái)階函數(shù)。式（33）中唯一未確定的物理量為tn+δt時(shí)刻界面上的平衡態(tài)分布函 數(shù)feq(xij，ξ，tn+δt)，其 可 由 相 容 性 條 件 來(lái)計(jì)算

式 中 ·≥0和 ·＜0分 別 表 示 在 速 度 空 間 中ξ1≥0和ξ1＜0 的半空間積分。通過(guò)式（34）計(jì)算得到單元界面守恒量W(xij，tn+δt)，并由此計(jì)算界面上的平衡態(tài)分布函數(shù)feq(xij，ξ，tn+δt)，則界面上的一階截?cái)喾植己瘮?shù)f(xij，ξ，tn+δt)可以完全確定。將f(xij，ξ，tn+δt)代入式（28）便可求得局部坐標(biāo)系下的Navier-Stokes 方程通量，應(yīng)用方程（27）可變換得到全局坐標(biāo)系下的通量。在該過(guò)程中可以同時(shí)計(jì)算無(wú)黏和黏性通量。

由于采用Boltzmann-BGK 模型方程漸進(jìn)展開(kāi)解的一階截?cái)嘟Y(jié)合Maxwell 分布函數(shù)還原得到的Navier-Stokes 方程普朗特?cái)?shù)恒為1，需要對(duì)能量方程的通量進(jìn)行修正

式中：pL和pR分別為單元界面左右兩側(cè)的壓力；Δt為式（25）顯示推進(jìn)的最大時(shí)間步長(zhǎng)。

3.2 可壓縮流應(yīng)用

由于Boltzmann-BGK 模型方程漸進(jìn)展開(kāi)解的一階截?cái)嘟Y(jié)合Maxwell 分布函數(shù)可以還原回Navier-Stokes 方程，所以采用該漸進(jìn)展開(kāi)解來(lái)計(jì)算Navier-Stokes 方程的通量等效于在單元界面也求解了等價(jià)的流動(dòng)控制方程。該格式已被成功應(yīng)用于低速［56］、跨聲速［57］、超/高超聲速［58］、可壓縮流固耦 合 傳 熱［59］、可 壓 縮 多 組 分 流［60］以 及 化 學(xué) 反 應(yīng)流［61］等流動(dòng)問(wèn)題的求解。數(shù)值結(jié)果表明，該格式具有較好的穩(wěn)定性，不會(huì)出現(xiàn)“紅寶石”等非物理解。

（1）該算法在航空氣動(dòng)力計(jì)算方面的應(yīng)用。算例1 為M6 機(jī) 翼 跨 聲 速 繞 流 問(wèn) 題［62］，馬 赫 數(shù) 為0.839 5，來(lái) 流 迎 角 為3.06°，側(cè) 滑 角 為0°。采 用NASA 網(wǎng)站上提供的標(biāo)準(zhǔn)網(wǎng)格，網(wǎng)格單元數(shù)為294 912。圖7 給出了機(jī)翼的壓力云圖和截面壓力系數(shù)Cp分布。圖中可以清晰觀察到機(jī)翼上表面“λ形狀”的激波，而且截面壓力系數(shù)分布與實(shí)驗(yàn)值也吻合較好。算例2 為DLR-F6 翼身組合體含有和不含F(xiàn)X2B 導(dǎo)流罩兩種情形的跨聲速繞流問(wèn)題［63］。馬赫數(shù)為0.75，雷諾數(shù)為3×106，來(lái)流迎角為0.49°，側(cè)滑角為0°。采用NASA 網(wǎng)站上提供的標(biāo)準(zhǔn)網(wǎng)格，包含26 個(gè)塊和2 298 880 個(gè)網(wǎng)格節(jié)點(diǎn)。圖8 展示了DLR-F6 翼身組合體的壓力云圖和截面壓力系數(shù)分布。同樣，當(dāng)前方法的計(jì)算結(jié)果與實(shí)驗(yàn)值以及參考數(shù)值解均吻合較好。計(jì)算效率方面，當(dāng)前算法與采用Roe 格式計(jì)算無(wú)黏通量結(jié)合中心格式計(jì)算黏性通量基本相當(dāng)。

圖7 M6 機(jī)翼的壓力云圖和截面壓力系數(shù)分布[62]Fig.7 Pressure contours and pressure coefficient distribution at a selected position for M6 wing[62]

圖8 DLR-F6 翼身組合體的壓力云圖和截面壓力系數(shù)分布[63]Fig.8 Pressure contours and pressure coefficient distribution at a selected position for DLR-F6 wing-body[63]

（2）該方法在考慮化學(xué)反應(yīng)的導(dǎo)彈噴流方面的應(yīng)用。噴流控制技術(shù)在航空航天領(lǐng)域有重要的應(yīng)用。高速飛行器在需要進(jìn)行快速機(jī)動(dòng)時(shí)，通過(guò)氣動(dòng)控制舵面控制可能會(huì)需要一個(gè)較長(zhǎng)的響應(yīng)時(shí)間，而噴流控制技術(shù)可以得到快速響應(yīng)的直接力，能夠滿足高機(jī)動(dòng)性的要求。由于噴流與主流之間相互干擾的流場(chǎng)結(jié)構(gòu)十分復(fù)雜，研究多化學(xué)反應(yīng)效應(yīng)下的側(cè)向噴流流場(chǎng)對(duì)高速飛行器直接力控制手段的工程應(yīng)用有很大的參考價(jià)值。為了考慮化學(xué)反應(yīng)中的不同組分，當(dāng)前算法需要增加組分控制方程，具體細(xì)節(jié)可以參見(jiàn)文獻(xiàn)［64］。圖9 給出了應(yīng)用該方法計(jì)算得到的錐-筒-裙結(jié)構(gòu)的普通導(dǎo)彈外形側(cè)噴流場(chǎng)壓力云圖和使用不同化學(xué)反應(yīng)模型時(shí)物面第1 層網(wǎng)格的切向速度。圖9（b）包含了無(wú)化學(xué)反應(yīng)（Nochem）、空氣化學(xué)反應(yīng)（Air）、燃燒化學(xué)反應(yīng)（Combustion）和混合化學(xué)反應(yīng)（Mix）4 種情形的計(jì)算結(jié)果?？梢钥吹剑煌瘜W(xué)反應(yīng)模型作用下的物面流場(chǎng)拓?fù)浣Y(jié)構(gòu)基本一致，在噴口上游，燃燒化學(xué)反應(yīng)模型算例的分離最靠前，而在噴口下游，燃燒化學(xué)反應(yīng)模型算例的分離最靠后。圖10 給出了導(dǎo)彈外形側(cè)噴流混合化學(xué)反應(yīng)模型反應(yīng)焓變功率密度及最大焓變反應(yīng)分布。對(duì)于混合化學(xué)反應(yīng)模型，由于該算例噴口附近的燃燒反應(yīng)的化學(xué)反應(yīng)焓變功率密度比空氣化學(xué)反應(yīng)高出約4 個(gè)數(shù)量級(jí)，因此其化學(xué)反應(yīng)熱量變化與燃燒化學(xué)反應(yīng)模型算例基本一致。

圖9 導(dǎo)彈外形側(cè)噴流壓力云圖和使用不同化學(xué)反應(yīng)模型時(shí)物面第一層網(wǎng)格的切向速度[64]Fig.9 Pressure contours and tangential velocity at the first layer grid of wall using different chemical reaction models for missile jet flow[64]

圖10 導(dǎo)彈外形側(cè)噴流混合化學(xué)反應(yīng)模型反應(yīng)焓變功率密度及最大焓變反應(yīng)分布[64]Fig.10 Enthalpy change power density and maximum enthalpy change reaction distribution for missile jet flow using the mixed chemical reaction model[64]

4 基于Grad 分布函數(shù)的氣體動(dòng)理學(xué)通量算法

4.1 算法設(shè)計(jì)

由于引入了連續(xù)性假設(shè)，Navier-Stokes 方程僅適合于連續(xù)流動(dòng)問(wèn)題求解。 相比于Navier-Stokes 方程，Boltzmann 方程不受連續(xù)性假設(shè)的限制，因而理論上適用于連續(xù)到稀薄整個(gè)流域范圍。但是，采用確定性方法求解Boltzmann方程時(shí)不僅需要在物理空間離散，還需要在分子速度空間離散，因而三維復(fù)雜流動(dòng)問(wèn)題求解時(shí)會(huì)導(dǎo)致維度災(zāi)難。實(shí)際上Boltzmann 方程存在對(duì)應(yīng)的宏觀守恒律方程，只是宏觀控制方程的通量依賴于分布函數(shù)的具體形式。顯然，對(duì)于稀薄程度非常高的流動(dòng)問(wèn)題，分布函數(shù)的具體形式只能通過(guò)求解Boltzmann 方程來(lái)獲得。但對(duì)于稀薄程度不太高的流動(dòng)問(wèn)題，借助于矩方法的思想，可以采用某些具有顯式表達(dá)式的非平衡分布函數(shù)來(lái)計(jì)算宏觀控制方程的通量?；谠摬呗?，本文發(fā)展了基于Grad 分布函數(shù)的GKFS 用于適度稀薄流動(dòng)問(wèn)題的求解。該方法僅求解宏觀方程式（6），采用有限體積法對(duì)該方程在空間網(wǎng)格Vi進(jìn)行離散可得

求解式（37）的關(guān)鍵為計(jì)算單元界面通量Fij。由于Fij依賴于分布函數(shù)，需要通過(guò)式（8）進(jìn)行計(jì)算。因此，需要先確定單元界面的分布函數(shù)f(xij，ξ，tn+δt)。

為了適用于稀薄流動(dòng)問(wèn)題，f(xij，ξ，tn+δt)可以由Boltzmann-BGK 模型方程的特征解來(lái)計(jì)算。對(duì)式（1）沿特征線積分有

由此可見(jiàn)，feq(xij，ξ，tn+δt)也取決于f(xijξδt，ξ，tn)。故f(xij-ξδt，ξ，tn)是計(jì)算宏觀守恒律式（37）通量的關(guān)鍵。

為了避免在分子速度空間離散，GKFS 采用Grad 分布函數(shù)來(lái)近似單元界面周圍的初始分布函數(shù)f(xij-ξδt，ξ，tn)。Grad 13 分 布 函 數(shù) 可 以 表示為

式中：σij為應(yīng)力張量分量；ci為分子最可幾熱運(yùn)動(dòng)速度分量；qi為熱流分量。由式（41）可知，fG13可以顯式表示為宏觀量的函數(shù)。通過(guò)插值可以求得界面周圍的宏觀量，進(jìn)而應(yīng)用Grad 分布函數(shù)計(jì)算f(xij-ξδt，ξ，tn)，然后計(jì)算單元界面守恒量W(xij，tn+δt) 和 平 衡 態(tài) 分 布 函 數(shù)feq(xij，ξ，tn+δt)，由 此 便 可 獲 得 界 面 分 布 函 數(shù)f(xij，ξ，tn+δt)。最后，式（37）的通量可通過(guò)將f(xij，ξ，tn+δt)代入式（8）求得。具體細(xì)節(jié)可以參見(jiàn)文獻(xiàn)［65-66］。

通過(guò)沿時(shí)間方向推進(jìn)求解宏觀守恒律式（37），可以實(shí)現(xiàn)單元中心守恒量的更新。但是在基于Grad 分布函數(shù)的GKFS 中，還需要更新單元中心的應(yīng)力和熱流，以便計(jì)算下一時(shí)刻的Grad 分布函數(shù)。受離散速度方法（Discrete velocity method，DVM）的啟發(fā)，應(yīng)力和熱流也可借助于f(xij，ξ，tn+δt)來(lái)更新。具體地，單元界面上的應(yīng)力和熱流可以表示為

式中角括號(hào)中的指標(biāo)表示張量的對(duì)稱和無(wú)跡部分。由此，單元中心的應(yīng)力和熱流便可取為單元界面的平均值。應(yīng)用該方法可以避免直接求解高階非守恒量的控制方程。

基于Grad 分布函數(shù)的GKFS 所能求解稀薄流動(dòng)的克努森數(shù)范圍依賴于所采用的近似分布函數(shù)。為了提高該方法適用的克努森數(shù)上限，作者團(tuán)隊(duì)還發(fā)展了基于Grad 45 分布函數(shù)的GKFS［67］。相比于離散速度方法，基于Grad 分布函數(shù)的GKFS 無(wú)需在分子速度空間離散，因而計(jì)算量和內(nèi)存開(kāi)銷更小，基本上與Navier-Stokes 方程求解算法保持同一數(shù)量級(jí)。另外相比于矩方法，由于當(dāng)前方法僅求解基于守恒律的控制方程，因而避開(kāi)了復(fù)雜高階非守恒量控制方程的求解，也不需要實(shí)施高階非守恒量的邊界條件。但是，由于引入了Grad 分布函數(shù)來(lái)近似真實(shí)的分布函數(shù)，當(dāng)前方法適用的克努森數(shù)范圍與基于Grad 分布函數(shù)的矩方法基本一致，不能直接用于全流域流動(dòng)問(wèn)題的求解。

4.2 稀薄流應(yīng)用

本節(jié)通過(guò)兩個(gè)算例來(lái)驗(yàn)證基于Grad 分布函數(shù)的GKFS 在求解稀薄流問(wèn)題中的性能。首先考慮不同克努森數(shù)下的三維頂蓋驅(qū)動(dòng)流問(wèn)題［66］。該問(wèn)題中，方腔頂蓋以u(píng)W=0.15 2RT0的速度沿x方向運(yùn)動(dòng)，其余物面靜止不動(dòng)，T0為參考溫度。所有物面采用等溫邊界條件，物面溫度取為T0。圖11給出了采用當(dāng)前方法和DVM 計(jì)算得到的不同克努森數(shù)時(shí)的速度剖面。采用DVM 計(jì)算時(shí)，需要在分子速度空間進(jìn)行離散，離散網(wǎng)格選為18×18×18，并使用Gauss-Hermite 求積來(lái)計(jì)算宏觀量。結(jié)果表明克努森數(shù)在0.025～0.1 范圍內(nèi)，當(dāng)前方法的計(jì)算結(jié)果均與直接求解Boltzmann-BGK 模型方程的DVM 結(jié)果一致。由于無(wú)需在分子速度空間離散，當(dāng)前方法的計(jì)算效率比DVM 高近兩個(gè)數(shù)量級(jí)，結(jié)果如表1 所示。

表1 計(jì)算時(shí)間比較[66]Table 1 Comparison of computational time[66] s

圖11 不同克努森數(shù)下三維頂蓋驅(qū)動(dòng)流的速度剖面圖[66]Fig.11 Velocity profiles for 3D lid-driven cavity flow at different Knudsen numbers[66]

算例2 為熱流逸問(wèn)題，即流動(dòng)由溫度梯度驅(qū)動(dòng)［67］。如圖12 所示，封閉方腔左右壁面的溫度固定為263 K，而上下壁面溫度呈先升后降的折線段分布，中點(diǎn)處最高溫為283 K。圖13 給出了采用當(dāng)前方法和DVM 計(jì)算得到的不同克努森數(shù)時(shí)的溫度剖面。采用DVM 計(jì)算時(shí)，需要在分子速度空間進(jìn)行離散，離散網(wǎng)格選為28×28，并使用Gauss-Hermite 求積來(lái)計(jì)算宏觀量。其他設(shè)置方面，當(dāng)前方法和DVM 保持一致。結(jié)果表明，在較低克努森數(shù)時(shí)（Kn= 0.01），無(wú)論是采用Grad 13分布函數(shù)還是Grad 45 分布函數(shù)，當(dāng)前方法的結(jié)果均與DVM 吻合較好。但當(dāng)克努森數(shù)逐漸增大時(shí)，Grad 13 分布函數(shù)的結(jié)果偏離DVM 結(jié)果，但Grad 45 分布函數(shù)的結(jié)果仍與DVM 結(jié)果基本吻合。由此表明，隨著克努森數(shù)的增加，稀薄效應(yīng)變強(qiáng)，需要選用能合理描述更強(qiáng)非平衡效應(yīng)的分布函數(shù)才能獲得準(zhǔn)確的計(jì)算結(jié)果。

圖12 熱驅(qū)動(dòng)方腔流示意圖[67]Fig.12 Schematic diagram of thermally driven cavity flow[67]

圖13 不同克努森數(shù)時(shí)熱驅(qū)動(dòng)方腔流的垂直中心線（上）和水平中心線（下）溫度分布圖[67]Fig.13 Temperature profiles along the vertical (upper) and horizontal (lower) center lines for thermally driven cavity flow at different Knudsen numbers[67]

5 結(jié) 論

本文首先簡(jiǎn)單回顧了傳統(tǒng)CFD 中的幾種典型的通量重構(gòu)算法，它們?cè)趩卧缑娌捎脭?shù)學(xué)重構(gòu)或部分物理重構(gòu)的方式來(lái)計(jì)算宏觀方程的通量。因而在單元界面處并不能保證完全滿足所求解的宏觀控制方程，而且無(wú)粘和粘性通量通常需要分開(kāi)計(jì)算。隨后，本文重點(diǎn)介紹了LBFS 和GKFS 及其應(yīng)用，展示了該算法在連續(xù)流以及適度稀薄流動(dòng)問(wèn)題中的應(yīng)用。

（1）Boltzmann-BGK 模型方程漸進(jìn)展開(kāi)解的一階截?cái)嗾脤?duì)應(yīng)于Navier-Stokes 方程，這給重構(gòu)Navier-Stokes 方程通量提供了另一種思路。在連續(xù)流動(dòng)問(wèn)題求解時(shí)，LBFS/GKFS 正是采用了該漸進(jìn)展開(kāi)解的一階截?cái)鄟?lái)計(jì)算Navier-Stokes 方程通量，因而單元界面也等效于求解了相應(yīng)的宏觀控制方程。應(yīng)用該方法，可以同時(shí)考慮法向和切向速度對(duì)通量的貢獻(xiàn)，無(wú)需采用被動(dòng)標(biāo)量的方式來(lái)計(jì)算切向通量，并且也可以采用一致的方式計(jì)算無(wú)粘和粘性通量。無(wú)論是離散平衡態(tài)分布函數(shù)或連續(xù)平衡態(tài)分布函數(shù)，只要它們滿足還原Navier-Stokes 方程所需的矩關(guān)系，均可用于構(gòu)造相應(yīng)的通量計(jì)算格式。目前該類算法在低速到高超聲速流動(dòng)、多相流、動(dòng)邊界等方面均已獲得了成功的應(yīng)用，并且展現(xiàn)出良好的穩(wěn)定性。但是，該方法基本上還局限于二階精度，如何將其推廣到具有緊致屬性的高精度格式值得進(jìn)一步研究。

（2）采用Boltzmann-BGK 模型方程離散特征解結(jié)合Grad 分布函數(shù)來(lái)計(jì)算宏觀守恒律方程的通量，可以將GKFS 推廣應(yīng)用于適度稀薄流動(dòng)問(wèn)題的求解。該方法相比于DVM，可以避開(kāi)在分子速度空間離散，其計(jì)算量和存儲(chǔ)量與Navier-Stokes 方程算法基本上保持在同一量級(jí)。另外相比于矩方法，由于當(dāng)前方法無(wú)需求解高階非守恒量的控制方程，因而更為簡(jiǎn)單。但是，該方法適用的克努森數(shù)上限與所采用的非平衡分布函數(shù)相關(guān)，為了將其推廣到更高的克努森數(shù)范圍，需要發(fā)展適用于更強(qiáng)非平衡效應(yīng)的分布函數(shù)。借助于機(jī)器學(xué)習(xí)手段，從大量DVM 模擬結(jié)果中訓(xùn)練出適用于更強(qiáng)非平衡效應(yīng)的分布函數(shù)是一個(gè)非常值得探索的方向。另外，將當(dāng)前算法與DVM 搭接，在低克努森數(shù)流場(chǎng)區(qū)域采用GKFS 求解，而其他區(qū)域采用DVM 求解，以實(shí)現(xiàn)全流域范圍的高效計(jì)算，具有極大的實(shí)用價(jià)值。