楊鯉銘，吳杰，董昊，杜銀杰

南京航空航天大學 航空學院，南京 210016

跨流域流動問題廣泛存在于工程實際中，譬如微電子機械系統(tǒng)[1-2]和載入飛行器[3-5]。對于該類問題，由于分子平均自由程非常接近或遠大于幾何的特征尺度，連續(xù)性假設(shè)不再成立，因而Navier-Stokes方程不再適用。相比于Navier-Stokes方程，Boltzmann方程直接描述氣體分子的分布函數(shù)，沒有引入連續(xù)性假設(shè)，因而理論上可用于從連續(xù)到自由分子流整個流域范圍。為了求解Boltzmann方程，離散速度方法(Discrete Velocity Method，DVM)經(jīng)常被采用[6-7]，該方法在物理空間和分子速度空間同時離散求解Boltzmann方程。

基于不同的通量重構(gòu)方式和分布函數(shù)演化策略，文獻中提出了多種算法。其中，借鑒傳統(tǒng)計算流體力學中求解Navier-Stokes方程的通量重構(gòu)方法，即不考慮分子碰撞效應的迎風格式被廣泛用于重構(gòu)Boltzmann-BGK方程的通量。例如，Yang和Huang[8]、Mieussens[9]采用了二階迎風格式，Wang等[10]采用了三階迎風格式，Kudryavtsev和Shershnev[11]、Diaz等[12]采用了WENO格式。由于通量計算時不考慮分子碰撞的影響，因而無需聯(lián)立所有的離散分布函數(shù)來計算單元界面的平衡態(tài)，所以該方法非常簡單，且易于在離散速度空間并行化計算。然而該方法在連續(xù)和近連續(xù)流區(qū)域也會導致一些問題，主要原因在于這些區(qū)域的分子平均自由程非常小，而物理空間的網(wǎng)格尺度很難小于分子平均自由程，所以分子的碰撞對通量計算的貢獻不可忽略。不考慮分子碰撞的影響時，會使得連續(xù)和近連續(xù)流區(qū)域計算的結(jié)果不準確。另一方面，由于最新時刻的平衡態(tài)是未知的，Boltzmann-BGK方程的碰撞項一般只能采用半隱式來離散，會導致連續(xù)和近連續(xù)流區(qū)域計算收斂速度較慢[13-14]。

為了克服上述傳統(tǒng)DVM的缺陷，一些新的算法被提出用于求解Boltzmann-BGK方程。譬如，氣體動理學統(tǒng)一算法(Gas Kinetic Unified Algorithm，GKUA)[15-17]、統(tǒng)一氣體動理學格式(Unified Gas Kinetic Scheme，UGKS)[18-20]、離散統(tǒng)一氣體動理學格式(Discrete Unified Gas Kinetic Scheme，DUGKS)[21-22]等。其中UGKS在計算通量時采用了含碰撞項Boltzmann-BGK方程的局部積分解，同時同步求解了與Boltzmann-BGK方程對應的宏觀伴隨方程。由于通量重構(gòu)考慮了分子碰撞的影響，UGKS在連續(xù)和近連續(xù)流區(qū)域也可以獲得準確可靠的結(jié)果。另外，由于同步求解了宏觀伴隨方程，其預估的平衡態(tài)可用于實現(xiàn)Boltzmann-BGK方程碰撞項的全隱式離散，從而極大地提高了連續(xù)流區(qū)域的收斂速度。但是，含碰撞項Boltzmann-BGK方程的局部積分解較為復雜，UGKS在過渡流和自由分子流區(qū)域的計算效率通常低于傳統(tǒng)DVM。此外，UGKS需要聯(lián)立所有的離散分布函數(shù)來計算單元界面的宏觀量，進而計算平衡態(tài)分布。

最近，為了保持傳統(tǒng)DVM的簡單性優(yōu)勢，同時克服其在連續(xù)流和近連續(xù)流區(qū)域計算精度差、計算效率低的缺陷，Yang等[23-24]提出了一種改進離散速度方法(Improved Discrete Velocity Method，IDVM)。IDVM的設(shè)計思想是基于傳統(tǒng)DVM在過渡流和自由分子流區(qū)域可以獲得可靠的計算結(jié)果，而在連續(xù)和近連續(xù)流區(qū)域精度差、效率低的表現(xiàn)，通過引入宏觀伴隨方程來完善其在全流域的計算性能。與UGKS類似，該方法也同步求解了Boltzmann-BGK方程和相應的宏觀伴隨方程，且在計算宏觀伴隨方程通量時考慮了分子碰撞的影響，可以視為UGKS的同類算法。但不同于UGKS，IDVM在求解Boltzmann-BGK方程時，單元界面通量的計算沒有耦合碰撞項，保持與傳統(tǒng)DVM一致。沒有耦合碰撞項的原因在于傳統(tǒng)DVM在過渡流和自由分子流區(qū)域可以獲得可靠的解，而在連續(xù)和近連續(xù)流區(qū)域，宏觀伴隨方程將主導流場解，此時由于該方程通量的計算耦合了碰撞項，亦可獲得合理的解。通過這種策略，既不會破壞傳統(tǒng)DVM的簡單性優(yōu)勢，同時還保證了全流域的計算精度。此外，由于同步求解了相應的宏觀伴隨方程，可以利用預估的平衡態(tài)來實現(xiàn)Boltzmann-BGK方程碰撞項的全隱式離散，從而提高算法在連續(xù)流區(qū)域的收斂速度。實際上，在傳統(tǒng)DVM基礎(chǔ)上，通過引入宏觀伴隨方程來提高連續(xù)和近連續(xù)流區(qū)域的計算精度和計算效率的思想在Chen[25]和皮興才[26]等的工作中也已經(jīng)被應用和驗證。但是，由于宏觀伴隨方程的通量和Boltzmann-BGK方程的通量計算不一致，該類方法在理論上存在不自洽的問題，需要在后續(xù)工作中進一步完善。

為了進一步提高IDVM的計算效率，本文將著重研究不同速度空間求積方式和守恒修正對結(jié)果的影響，實現(xiàn)速度空間網(wǎng)格量的盡可能降低。在IDVM中，宏觀物理量由離散分布函數(shù)求積而得到，其實質(zhì)為對基于離散點所構(gòu)造的多項式進行積分，因而會受到多項式誤差的影響。使用均勻節(jié)點構(gòu)造高次多項式時，區(qū)間邊緣的誤差可能很大，從而導致Runge現(xiàn)象[27]。這意味著，高階的求積格式并不一定總能獲得高精度的結(jié)果。所以，在克努森數(shù)(Kn)較大的情形，由于分布函數(shù)嚴重偏離平衡態(tài)，矩形求積規(guī)則反而更容易保證結(jié)果的可靠性；而在宏觀速度較低和Kn數(shù)較小的情形，由于分布函數(shù)接近平衡態(tài)且關(guān)于原點基本對稱，可以采用Gauss-Hermite積分來獲得高精度結(jié)果。此外，由于分布函數(shù)僅在宏觀速度附近時值較大，而遠離該區(qū)域的值逐漸變小，所以采用非結(jié)構(gòu)網(wǎng)格來離散速度空間更能充分利用分布函數(shù)的這一特征，實現(xiàn)速度空間網(wǎng)格量的降低。最后，由于數(shù)值積分獲得的宏觀量與真實積分結(jié)果存在誤差，可能會導致相容性條件不能嚴格滿足，使得誤差累積或計算發(fā)散。為了克服這一問題，本文還引進了守恒修正來強制滿足相容性條件，使得較粗糙的速度空間網(wǎng)格也可以獲得準確的計算結(jié)果，從而進一步提升計算效率。

1 Boltzmann方程和傳統(tǒng)離散速度方法

原始的Boltzmann方程是積分-微分形式，求解較為困難。針對航天領(lǐng)域所關(guān)注的氣動力和氣動熱問題，可以求解采用BGK-Shakhov模型[28]來近似碰撞項的Boltzmann-BGK方程

(1)

(2)

(3)

式中：Pr為普朗特數(shù)；c=ξ-u為分子的特異速度,u為宏觀速度；q為熱流；ρ為密度；T為溫度；Rg為氣體常數(shù)。另外，文中約定用黑體來表示矢量，用對應的白斜體來表示矢量的模。

離散分布函數(shù)fα與宏觀守恒量W、宏觀通量F、熱流q和應力τ存在如下關(guān)系：

W=[ρρuρE]T=〈ψf〉α

(4)

F=[FρFρuFρE]T=〈ξψf〉α

(5)

(6)

τ=〈ccf〉α

(7)

(8)

式中：上標n表示當前時刻；Δt代表隱式推進的時間步長，由如下CFL條件確定:

Δt=

(9)

式中：σ為CFL數(shù)；V為控制體的體積；ξ1,max、ξ2,max和ξ3,max分別為x、y和z這3個方向的最大離散速度；ΔSx、ΔSy和ΔSz為控制體在y-z平面、x-z平面和x-y平面的投影面積；cs為聲速。

(10)

對方程(10)在控制體Vi積分可得

(11)

(12)

(13)

(14)

(15)

式中：H(nij·ξα)為臺階函數(shù)，當nij·ξα≥0時有H(nij·ξα)=1，而當nij·ξα<0時有H(nij·ξα)=0。

(16)

(17)

(18)

自此，方程(16)可以采用高效的LU-SGS方法來求解，分為如下的前掃和后掃兩步：

(19)

(20)

式中：L(i)和U(i)為N(i)的子集，分別包含單元編號大于i和小于i的單元。

2 改進離散速度方法

2.1 Boltzmann-BGK方程求解

為了克服傳統(tǒng)半隱式DVM在連續(xù)和近連續(xù)流區(qū)域收斂速度較慢的缺陷，IDVM采用了全隱式方法來離散Boltzmann-BGK方程

(21)

(22)

(23)

(24)

2.2 宏觀伴隨方程求解

(25)

(26)

(27)

式中：C(i)為包含i單元和其周圍單元的集合。將方程(27)代入方程(26)可得：

(28)

在應用LU-SGS求解方程(28)時，為了避免矩陣運算，可采用基于Euler方程的通量分裂算法來近似雅克比矩陣?R/?W，即

(29)

將方程(29)代入方程(28)可得：

(30)

式中：

(31)

Fc=[ρuρuu+pI(ρE+p)u]T

(32)

(33)

自此，方程(30)可以通過如下的前掃和后掃求解：

(34)

(35)

(36)

式中：fcol(xij,ξ,Δtp)為單元界面考慮了分子碰撞影響的Boltzmann-BGK方程的局部解，Δtp為通量重構(gòu)的時間步長，其可以通過方程(9)計算，但CFL數(shù)需要設(shè)置為小于1，本文選取0.95。通過相關(guān)推導發(fā)現(xiàn)[14]，fcol(xij,ξ,Δtp)最終可以表示為

fcol(xij,ξ,Δtp)=βfDVM(xij,ξα,Δtp)+

(1-β)fNS(xij,ξα,Δtp)

(37)

β=e-Δtp/τ

(38)

式中：fDVM(xij,ξα,Δtp)為無碰撞Boltzmann方程在單元界面的局部解，即方程(15)。fNS(xij,ξα,Δtp)為對應于連續(xù)流極限的分布函數(shù)。將方程(37)代入方程(36)有

〈ξψfNS(xij,ξ,Δtp)〉α=

(39)

式中：右端第1項可以通過將方程(15)代入進行計算，第2項實際上是Navier-Stokes方程的通量，可以采用常用的連續(xù)流求解器進行計算。采用Yang等[29]提出的LBFS算法進行求解。

2.3 速度空間離散與求解規(guī)則

在DVM中，由于速度空間被離散化，必須通過數(shù)值積分來計算宏觀量以及平衡態(tài)分布函數(shù)。因而，速度空間的離散化策略和求積規(guī)則的選取對計算精度和效率將產(chǎn)生直接的影響。對于低速且接近平衡態(tài)的流動，由于分布函數(shù)基本上呈正態(tài)分布，可以采用精度較高的Gauss-Hermite積分。以φ(ξ)在一維速度空間的積分為例

(40)

(41)

式中：HNV-1為NV-1階Hermite多項式。

對于高速流動和遠離平衡態(tài)的流動問題，由于分布函數(shù)可能會出現(xiàn)“雙峰”等不規(guī)則分布，一般可采用復化Newton-Cotes積分

(42)

目前文獻中應用較多的是四階Newton-Cotes公式，但高階積分可能會出現(xiàn)Runge現(xiàn)象。本文將比較一階到四階Newton-Cotes公式的計算效果。

實際上，當采用一階Newton-Cotes公式時，可以用非結(jié)構(gòu)網(wǎng)格來離散速度空間。此時的求積點即為控制體的中心點，求積權(quán)即為控制體的體積。采用非結(jié)構(gòu)網(wǎng)格來離散速度空間的另一個優(yōu)勢是可以僅在分布函數(shù)值可能較大的區(qū)域布置密網(wǎng)格，而在其他區(qū)域采用粗網(wǎng)格，以減少總的速度空間網(wǎng)格量。非結(jié)構(gòu)網(wǎng)格速度空間離散已在Yuan和Zhong[30]和Chen等[31]的工作中被應用。由于無法預知分布函數(shù)在速度空間中的準確分布，采用非結(jié)構(gòu)網(wǎng)格來離散速度空間時需要人為選擇截斷區(qū)域和加密區(qū)域。文中截斷區(qū)域取為

(43)

加密區(qū)域取為

(44)

式中：Ts為總溫；u0為流場中的最大和最小速度的均值。它們可以通過Navier-Stokes方程解來預估；系數(shù)a取值為4～5之間。

2.4 相容性條件與守恒修正

由于DVM中采用數(shù)值積分來計算宏觀量，相容性條件將不能嚴格滿足，從而可能會導致誤差累積或計算發(fā)散。為了解決這一問題，需要強制滿足相容性條件[32-33]，即

〈ψ1(fS-f)〉α=[000 -Prq]T

(45)

式中：ψ1=[1ξξ2/2cc2/2]T。由于fS是宏觀量的函數(shù)，而f在當前時刻是常數(shù)，方程(45)實際上是變量?=[ρuTq]T的非線性多元函數(shù)。為了便于求解，可以引進如下輔助函數(shù)：

Φ(?)=〈ψ1(fS-f)〉α-

[000 -Prq]T=0

(46)

方程(46)可以采用牛頓迭代求解

(47)

實際計算時，可以將?1取為f求矩的結(jié)果。該算法收斂非?？?，一般迭代2～4步殘差就可以收斂到10-9。

3 算例測試與驗證

本節(jié)通過幾個低速到超聲速、連續(xù)到自由分子流算例來檢驗當前算法的精度和效率。由于所考慮的算例均為定常流動，為了加速計算，宏觀伴隨方程求解時采用了50次內(nèi)迭代，以此來獲得更為準確的預估平衡態(tài)。此外，本文所有算例均只考慮單原子分子情形，對應的比熱比為5/3。所有計算均在個人計算機上完成，其配置為Intel(R) Xeon(R) Gold 6226R CPU@2.90 GHz。

3.1 頂蓋驅(qū)動流

頂蓋驅(qū)動流為經(jīng)典測試算例，可以通過調(diào)節(jié)Kn數(shù)或Re數(shù)來實現(xiàn)流動從連續(xù)流變?yōu)樽杂煞肿恿鳡顟B(tài)。該算例中，頂蓋的無量綱溫度和速度分別為TW=1和uW=0.15，其他壁面保持靜止且溫度固定為TW=1。無量綱動力學黏性系數(shù)計算如下:

μ=μrefTw

(48)

該算例采用變徑硬球模型，取w=0.81。另外，對于給定Kn數(shù)或Re數(shù)的情形，分別采用方程(49)或方程(50)來計算參考黏性系數(shù)μref

(49)

(50)

式中：ρ0=1和T0=TW分別為參考的無量綱密度和溫度。計算時，采用相鄰兩步原始變量的1-范數(shù)誤差最大值小于10-7作為收斂判據(jù)。

首先，以Kn=10的算例來測試不同階數(shù)的Newton-Cotes積分對計算結(jié)果的影響。計算時，將速度空間截斷為[-4, 4]×[-4, 4]，并用24×24的均勻結(jié)構(gòu)網(wǎng)格對其進行離散。計算得到的中截面速度分布如圖1所示?？梢钥闯觯斔俣瓤臻g網(wǎng)格非常稀疏時，較小階數(shù)的Newton-Cotes積分(n=1)的結(jié)果反而優(yōu)于高階Newton-Cotes積分(n=4)。由此表明，高Kn數(shù)時，高階Newton-Cotes積分可能會導致Runge現(xiàn)象。所以在后續(xù)計算中，如果Kn數(shù)較大或者流動速度較高時，均采用一階Newton-Cotes積分(即矩形律)，并利用非結(jié)構(gòu)網(wǎng)格來離散速度空間。

圖1 采用不同階數(shù)Newton-Cotes積分獲得的速度分布比較

其次，驗證當前算法在全流域計算時的精度。測試狀態(tài)包含Re=100,1 000的連續(xù)流，Kn=0.075的滑移流，Kn=1的過渡流和Kn=10的自由分子流。連續(xù)流和滑移流計算時分別采用8×8和28×28個離散點的Gauss-Hermite積分，而過渡流和自由分子流采用包含1 648個三角形單元的非結(jié)構(gòu)網(wǎng)格對速度空間進行離散，如圖2所示。另外，連續(xù)流計算時采用包含10 000個四邊形單元的均勻非結(jié)構(gòu)網(wǎng)格對物理空間進行離散，而其他工況則采用包含1 600個四邊形單元的均勻非結(jié)構(gòu)網(wǎng)格。圖3展示了Re=1 000，Kn=0.075和Kn=10時，采用IDVM計算的結(jié)果。可以看出，無論在連續(xù)流區(qū)域(Ghia等的結(jié)果[34])還是稀薄流區(qū)域(DSMC的計算結(jié)果[35])，IDVM都可以獲得準確合理的結(jié)果。Re=100和Kn=1的計算結(jié)果也與參考解吻合，但限于篇幅，圖中未展示。

最后，驗證速度空間非結(jié)構(gòu)網(wǎng)格和守恒修正對IDVM計算效率和內(nèi)存開銷的影響，同時比較IDVM和UGKS的性能。頂蓋驅(qū)動流問題中，原始IDVM僅在Kn=1的過渡流和Kn=10的自由分子流計算時需要采用Newton-Cotes積分，因此選用這兩個工況來比較不同算法的性能。為了獲得光滑的流場分布，原始IDVM一般需要采用101×101的均勻結(jié)構(gòu)網(wǎng)格來對速度空間進行離散。而采用非結(jié)構(gòu)網(wǎng)格離散速度空間時僅需要1 648個網(wǎng)格單元，內(nèi)存開銷約為原始IDVM的16%。

圖4比較了采用和不采用非結(jié)構(gòu)網(wǎng)格速度空間和守恒修正IDVM以及采用非結(jié)構(gòu)網(wǎng)格速度空間UGKS計算得到的溫度云圖。圖中云圖背景和黑線為UGKS的計算結(jié)果，紅線和白線分別為采用和不采用非結(jié)構(gòu)網(wǎng)格速度空間和守恒修正IDVM的計算結(jié)果?？梢钥闯?種格式計算的結(jié)果基本一致。另外，表1比較了采用和不采用非結(jié)構(gòu)網(wǎng)格速度空間和守恒修正IDVM的迭代次數(shù)和計算時間。顯然，引入非結(jié)構(gòu)網(wǎng)格速度空間和守恒修正后，計算效率得到了明顯提升。

圖2 頂蓋驅(qū)動流問題的速度空間非結(jié)構(gòu)網(wǎng)格

圖3 采用IDVM計算的頂蓋驅(qū)動流速度分布

圖4 采用不同方法計算的頂蓋驅(qū)動流溫度云圖比較

表1 采用和不采用非結(jié)構(gòu)網(wǎng)格速度空間和守恒修正IDVM迭代次數(shù)和計算時間的比較

3.2 高超聲速圓柱繞流

該算例為Ma=10和Kn=0.003 03的高超聲速圓柱繞流問題。其來流溫度為T0=200 K，壁面溫度固定為TW=500 K。動力學黏性系數(shù)采用方程(48)計算，并將參考黏性取為

(51)

式中：ρ0和R0分別為來流密度和圓柱半徑。由此可得Ma、Re和Kn之間的關(guān)系如下：

(52)

計算時，物理空間采用100×80的結(jié)構(gòu)網(wǎng)格進行離散。該算例如果速度空間采用四階Newton-Cotes積分將需要約121×121個網(wǎng)格單元[23]，而采用非結(jié)構(gòu)網(wǎng)格時僅需要6 338個網(wǎng)格單元，內(nèi)存開銷小于原來的1/2。圖5為所采用的速度空間非結(jié)構(gòu)網(wǎng)格，中間的圓形范圍為速度空間的加密區(qū)域。采用和不采用非結(jié)構(gòu)網(wǎng)格速度空間和守恒修正IDVM計算該算例時所需的CPU時間分別為1.00 h和2.30 h。圖6為采用IDVM計算的速度和溫度云圖。由于Kn數(shù)較小，在圓柱尾部可見兩個反向的渦。圖7比較了物面壓力和應力分布，IDVM與DS2V和UGKS的結(jié)果[36]均吻合較好。另外，采用IDVM計算得到的阻力系數(shù)為1.259，其與DS2V的計算結(jié)果1.252和UGKS的計算結(jié)果1.258的相對誤差均小于1%。

圖5 高超聲速圓柱繞流問題的速度空間非結(jié)構(gòu)網(wǎng)格

圖6 采用IDVM計算的速度和溫度云圖

圖7 圓柱表面壓力和應力分布

3.3 超聲速繞球流動

該算例為Ma=3.834和Kn=0.03的超聲速繞球流動問題。來流無量綱溫度取為T0=1，物面溫度固定為TW=[1+(γ-1)Ma2/2]T0。動力學黏性系數(shù)采用w=0.75的變徑硬球模型計算，物理空間采用包含20 736個六面體單元的非結(jié)構(gòu)網(wǎng)格進行離散，速度空間采用包含36 053個四面體單元的非結(jié)構(gòu)網(wǎng)格進行離散，如圖8所示，中間的球形范圍為速度空間的加密區(qū)域。計算的速度云圖如圖9所示。駐點線速度和溫度分布如圖10所示，IDVM與DSMC結(jié)果基本吻合。

圖8 超聲速繞球流動問題的速度空間非結(jié)構(gòu)網(wǎng)格

圖9 采用IDVM計算的速度u云圖

同時，采用IDVM計算的阻力系數(shù)為1.4295，與Li和Zhang[37]結(jié)果1.374 9以及DSMC結(jié)果1.412 2均吻合較好。此外，該算例采用了24線程的OpenMP并行計算時，總耗時為13.59 h。

4 結(jié) 論

本文提出了一種基于非結(jié)構(gòu)網(wǎng)格速度空間和守恒修正的改進離散速度方法。

1) 當前算法在全流域范圍均可獲得準確的計算結(jié)果。

2) 采用非結(jié)構(gòu)網(wǎng)格速度空間和守恒修正可以進一步提高IDVM的計算效率。

3) 采用非結(jié)構(gòu)網(wǎng)格速度空間相較于基于四階Newton-Cotes求積的內(nèi)存開銷更小。