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

南京航空航天大學(xué) 航空學(xué)院，南京 210016

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

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

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

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

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

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

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

(1)

(2)

(3)

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

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

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

(4)

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

(5)

(6)

τ=〈ccf〉α

(7)

(8)

式中：上標(biāo)n表示當(dāng)前時(shí)刻；Δt代表隱式推進(jìn)的時(shí)間步長，由如下CFL條件確定:

Δt=

(9)

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

(10)

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

(11)

(12)

(13)

(14)

(15)

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

(16)

(17)

(18)

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

(19)

(20)

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

2 改進(jìn)離散速度方法

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)

在應(yīng)用LU-SGS求解方程(28)時(shí)，為了避免矩陣運(yùn)算，可采用基于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)的時(shí)間步長，其可以通過方程(9)計(jì)算，但CFL數(shù)需要設(shè)置為小于1，本文選取0.95。通過相關(guān)推導(dǎo)發(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)為對應(yīng)于連續(xù)流極限的分布函數(shù)。將方程(37)代入方程(36)有

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

(39)

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

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

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

(40)

(41)

式中：HNV-1為NV-1階Hermite多項(xiàng)式。

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

(42)

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

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

(43)

加密區(qū)域取為

(44)

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

2.4 相容性條件與守恒修正

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

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

(45)

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

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

[000 -Prq]T=0

(46)

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

(47)

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

3 算例測試與驗(yàn)證

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

3.1 頂蓋驅(qū)動流

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

μ=μrefTw

(48)

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

(49)

(50)

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

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

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

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

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

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

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

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

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

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

3.2 高超聲速圓柱繞流

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

(51)

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

(52)

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

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

圖6 采用IDVM計(jì)算的速度和溫度云圖

圖7 圓柱表面壓力和應(yīng)力分布

3.3 超聲速繞球流動

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

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

圖9 采用IDVM計(jì)算的速度u云圖

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

4 結(jié) 論

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

1) 當(dāng)前算法在全流域范圍均可獲得準(zhǔn)確的計(jì)算結(jié)果。

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

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