基于速度空間非結構網格和守恒修正的改進離散速度方法

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

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

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

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

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

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

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

1 Boltzmann方程和傳統離散速度方法

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

(1)

(2)

(3)

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

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

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

(4)

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

(5)

(6)

τ=〈ccf〉α

(7)

(8)

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

Δt=

(9)

式中：σ為CFL數；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·ξα)為臺階函數，當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方程求解

為了克服傳統半隱式DVM在連續和近連續流區域收斂速度較慢的缺陷，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為通量重構的時間步長，其可以通過方程(9)計算，但CFL數需要設置為小于1，本文選取0.95。通過相關推導發現[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)為對應于連續流極限的分布函數。將方程(37)代入方程(36)有

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

(39)

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

2.3 速度空間離散與求解規則

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

(40)

(41)

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

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

(42)

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

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

(43)

加密區域取為

(44)

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

2.4 相容性條件與守恒修正

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

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

(45)

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

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

[000 -Prq]T=0

(46)

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

(47)

實際計算時，可以將?1取為f求矩的結果。該算法收斂非?？欤话愕?～4步殘差就可以收斂到10-9。

3 算例測試與驗證

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

3.1 頂蓋驅動流

頂蓋驅動流為經典測試算例，可以通過調節Kn數或Re數來實現流動從連續流變為自由分子流狀態。該算例中，頂蓋的無量綱溫度和速度分別為TW=1和uW=0.15，其他壁面保持靜止且溫度固定為TW=1。無量綱動力學黏性系數計算如下:

μ=μrefTw

(48)

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

(49)

(50)

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

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

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

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

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

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

圖2 頂蓋驅動流問題的速度空間非結構網格

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

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

表1 采用和不采用非結構網格速度空間和守恒修正IDVM迭代次數和計算時間的比較

3.2 高超聲速圓柱繞流

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

(51)

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

(52)

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

圖5 高超聲速圓柱繞流問題的速度空間非結構網格

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

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

3.3 超聲速繞球流動

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

圖8 超聲速繞球流動問題的速度空間非結構網格

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

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

4 結 論

本文提出了一種基于非結構網格速度空間和守恒修正的改進離散速度方法。

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

2) 采用非結構網格速度空間和守恒修正可以進一步提高IDVM的計算效率。

3) 采用非結構網格速度空間相較于基于四階Newton-Cotes求積的內存開銷更小。