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

（1.新加坡國立大學機械工程系，新加坡 119260，新加坡；2.南京航空航天大學航空學院，南京 210016；3.南京航空航天大學非定常空氣動力學與流動控制工信部重點實驗室，南京 210016；4.南京航空航天大學機械結構力學及控制國家重點實驗室，南京 210016）

隨著計算機性能的不斷提升，計算流體力學（Computational fluid dynamics，CFD）已經成為解決流體流動問題的重要手段之一，在航空、航天、航海等領域均已取得了非常成功的應用［1-4］。與理論流體力學和實驗流體力學相比，CFD 的應用場景更為靈活和廣泛，且大多數(shù)情況下更為經濟。特別是對于實驗難于復現(xiàn)的場景以及不便測量的區(qū)域，例如臨近空間的高速稀薄流動和航空發(fā)動機內部的精細流場結構，CFD 更是發(fā)揮著不可替代的作用。CFD 通過求解流體流動控制方程來獲得流場的解，常用的方法包括有限差分法（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］等。在這些方法中，由于離散控制方程的解通常定義在解點上，如何利用解點上的值來計算單元界面的通量，是CFD 最為關心的問題之一，該過程也被稱為通量重構。

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

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

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

與上述通量重構的思路不同，氣體動理學格式（Gas kinetic scheme，GKS）通過在單元界面應用Boltzmann 模型方程的局部解來計算宏觀Navier-Stokes 方程的通量［25-30］，實現(xiàn)了通量的完全物理重構。在連續(xù)介質假設的前提下，Boltzmann模型方程可以還原回Navier-Stokes 方程，該關系為采用Boltzmann 模型方程的局部解來計算Navier-Stokes 方程的通量奠定了基礎。由于Boltzmann 模型方程僅為分布函數(shù)的單變量方程，形式比Navier-Stokes 方程更簡單，因而可以方便地獲得該方程的局部解。 GKS 正是采用Boltzmann 模型方程在單元界面的局部積分解結合對應于Navier-Stokes 方程截斷形式的初始分布函數(shù)來計算Navier-Stokes 方程的通量，從而使得單元界面的物理量也滿足Navier-Stokes 方程的等價形式，并且無黏通量和黏性通量可以采用統(tǒng)一的方式計算。另一方面，GKS 由于在求解過程中引入了許多參數(shù)，其計算效率要低于傳統(tǒng)的Navier-Stokes 方程求解器。

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

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

1.1 Boltzmann 模型方程及其漸進展開解

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

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

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

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

為了便于推導Boltzmann-BGK 模型方程的漸進展開解，可將式（1）改寫為

式中Df=?t f+ξ·?f為分布函數(shù)的物質導數(shù)。將上述f的表達式不斷地代入方程最后1 項，即可獲得Boltzmann-BGK 模型方程的漸進展開解，即

式中已將τ近似為局部常數(shù)，因而可以提到物質導數(shù)算子前面。實際上，式（5）的不同階截斷即代表不同的宏觀控制方程。

1.2 宏觀守恒律方程

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

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

如果將分布函數(shù)f近似為其平衡態(tài)feq，則對應的宏觀控制方程為Euler 方程。此時相當于將τ取為0，即克努森數(shù)Kn設置為0，等效于引入了無粘假設。為了還原Navier-Stokes 方程，可以將分布函數(shù)f取為其一階截斷形式［44］

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

并且：μb為體積黏性系數(shù)；cV為等容比熱容；cp為等壓比熱容；κ為熱傳導系數(shù)；Pr為普朗特數(shù)。上述推導表明，采用Boltzmann-BGK 模型方程漸進展開解的一階截斷即可還原得到可壓縮Navier-Stokes 方程，并且截斷誤差為O(τ2)。實際上分布函數(shù)f的一階截斷相當于引入了連續(xù)性假設，即克努森數(shù)遠小于1 的情形，因此O(τ2)屬于高階小量，可以忽略。另外，由于采用了BGK 碰撞模型，還原得到的Navier-Stokes 方程的普朗特數(shù)恒為1。

值得指出，在還原上述可壓縮Navier-Stokes方程的過程中，僅需要用到平衡態(tài)分布函數(shù)滿足的7 個矩關系，分別對應質量、動量、能量、動量方程的無黏和黏性通量、能量方程的無黏和黏性通量。換言之，采用其他平衡態(tài)分布函數(shù)也可還原宏觀Navier-Stokes 方程，其要求僅為滿足所需的矩關系。另外，通過該推導過程可知，將分布函數(shù)取為Boltzmann-BGK 模型方程漸進展開解的一階截斷即可還原Navier-Stokes 方程。故可以將該漸進展開解的一階截斷直接代入式（8）來設計Navier-Stokes 方程的通量計算格式。

2 格子玻爾茲曼通量算法

2.1 算法設計

建立Boltzmann-BGK模型方程與Navier-Stokes方程之間的聯(lián)系之后，即可采用該關系來設計Navier-Stokes 方程的通量計算格式。本小節(jié)采用離散的格子Boltzmann 模型來替代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時刻單元界面周圍的平衡態(tài) 分 布 函 數(shù)feq(xij-ξαδt，ξα，tn)和tn+δt界 面 上的 平 衡 態(tài) 分 布 函 數(shù)feq(xij，ξα，tn+δt)。 其 中，feq(xij-ξαδt，ξα，tn)可 由 相 應 位 置 的 宏 觀 量 來 確定，而這些宏觀量可通過解點上的值插值得到，如圖1（b）所示。值得注意的是，單元界面周圍點（1～8）的位置可通過局部坐標系來確定，該方式適用于任意網格。

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

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

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

式（18）中的通量確定之后，便可應用Runge-Kutta 方法或者LU-SGS 方法沿時間方向推進求解。整體而言，LBFS仍然求解的是宏觀控制方程，只是在計算單元界面通量時采用了Boltzmann-BGK 模型方程漸進展開解的一階截斷。該一階截斷正好對應于所需求解的宏觀Navier-Stokes方程，所以利用LBFS所構造的通量在單元界面處也滿足所求解的控制方程，而且無粘和粘性通量可以采用一致的方式計算。相比于傳統(tǒng)CFD方法，由于僅通量計算格式發(fā)生了變化，傳統(tǒng)方法中的各種加速收斂技巧均可直接應用于當前格式。

2.2 不可壓縮流應用

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

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

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

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

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

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

3 氣體動理學通量算法

3.1 算法設計

圖5 圓盤在G=100 和m*=0.1 時的擺動模式瞬時渦系結構[53]Fig.5 Instantaneous vortical structures for the fluttering disk at G=100 and m*=0.1[53]

圖6 圓盤在G=200 和m*=0.75 時的翻騰模式瞬時渦系結構[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ù)也可用于構造Navier-Stokes 方程的通量計算格式，即氣體動理學通量算法。正如1.2 小節(jié)中的介紹，采用Boltzmann-BGK 模型方程漸進展開解的一階截斷結合Maxwell 分布函數(shù)可以還原可壓縮Navier-Stokes 方程。實際上，推導過程中僅使用了滿足Navier-Stokes 方程的7 個矩關系，因此Maxwell 分布函數(shù)也可以替換為其他的平衡態(tài)分布函數(shù)?；诖?，作者還發(fā)展了圓函數(shù)和球函數(shù)等簡化平衡態(tài)分布函數(shù)，并據(jù)此構造了相應的GKFS［54-55］。 不 失 一 般 性 ，該 小 節(jié) 采 用Boltzmann-BGK 模型方程漸進展開解的一階截斷結合 Maxwell 分布函數(shù)來構造可壓縮Navier-Stokes 方程的通量計算格式。首先，將式（11～13）改寫為

與基于離散模型的LBFS 類似，采用有限體積法對式（23）在空間網格Vi進行離散可得

其關鍵在于計算單元界面的一階截斷分布函數(shù)。

在推導具體的通量表達式時，由于需要在速度空間積分，可以在單元界面處引入局部坐標系以方便推導和應用。以三維問題為例，可將局部坐標系的1 方向定義為單元界面的法向，2 方向和3 方向均為界面的切向，并且3 個坐標方向構成正交右手坐標系。局部坐標系和全局坐標系中的通量滿足如下變換關系

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

與基于離散模型的LBFS 類似，可將單元界面分布函數(shù)的一階截斷表示為

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

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

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

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

由于采用Boltzmann-BGK 模型方程漸進展開解的一階截斷結合Maxwell 分布函數(shù)還原得到的Navier-Stokes 方程普朗特數(shù)恒為1，需要對能量方程的通量進行修正

式中：pL和pR分別為單元界面左右兩側的壓力；Δt為式（25）顯示推進的最大時間步長。

3.2 可壓縮流應用

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

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

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

圖9 導彈外形側噴流壓力云圖和使用不同化學反應模型時物面第一層網格的切向速度[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 導彈外形側噴流混合化學反應模型反應焓變功率密度及最大焓變反應分布[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ù)的氣體動理學通量算法

4.1 算法設計

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

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

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

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

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

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

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

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

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

4.2 稀薄流應用

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

表1 計算時間比較[66]Table 1 Comparison of computational time[66] s

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

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

圖12 熱驅動方腔流示意圖[67]Fig.12 Schematic diagram of thermally driven cavity flow[67]

圖13 不同克努森數(shù)時熱驅動方腔流的垂直中心線（上）和水平中心線（下）溫度分布圖[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 結 論

本文首先簡單回顧了傳統(tǒng)CFD 中的幾種典型的通量重構算法，它們在單元界面采用數(shù)學重構或部分物理重構的方式來計算宏觀方程的通量。因而在單元界面處并不能保證完全滿足所求解的宏觀控制方程，而且無粘和粘性通量通常需要分開計算。隨后，本文重點介紹了LBFS 和GKFS 及其應用，展示了該算法在連續(xù)流以及適度稀薄流動問題中的應用。

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

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