跨流域高超聲速繞流Boltzmann模型方程并行算法

李志輝， 吳俊林， 蔣新宇， 馬強,2

1. 中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所， 綿陽 621000

跨流域高超聲速繞流Boltzmann模型方程并行算法

李志輝1,2,*， 吳俊林1， 蔣新宇1， 馬強1,2

1. 中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所， 綿陽 621000

2. 國家計算流體力學實驗室， 北京 100191

通過對Boltzmann方程碰撞積分進行模型化處理，提出了統(tǒng)一描述各流域復雜高超聲速流動輸運現(xiàn)象的氣體分子速度分布函數(shù)控制方程，使用離散速度坐標法對分布函數(shù)方程所依賴的速度空間離散降維，構(gòu)造出直接求解分子速度分布函數(shù)的氣體動理論耦合迭代數(shù)值格式，研制了復雜飛行器高超聲速繞流氣動熱力學計算模型。基于對氣體動理論數(shù)值計算方法內(nèi)在并行性、變量依賴關(guān)系、數(shù)據(jù)通信與并行可擴展性的分析研究，使用區(qū)域分解并行化方法提出了新型的氣體動理論數(shù)值算法并行方案；研究了數(shù)據(jù)的并行分布與并行執(zhí)行特征，開展了大規(guī)模的并行化程序設(shè)計，構(gòu)造了可穩(wěn)定運行于成千上萬CPU的高性能并行算法，用以模擬各流域復雜飛行器的高超聲速繞流問題。以稀薄流到連續(xù)流環(huán)境下不同Knudsen數(shù)、不同馬赫數(shù)的可重復使用類球錐衛(wèi)星體及翼身組合復雜飛行器等氣動力、熱繞流問題為研究對象展開大規(guī)模并行計算，并進行算法驗證，所得計算結(jié)果與理論分析、直接模擬蒙特卡羅方法(DSMC)的模擬值及有關(guān)實驗數(shù)據(jù)吻合較好，揭示了飛行器跨流域高超聲速下的復雜流動機理與變化規(guī)律，提供了一條能夠可靠模擬高超聲速飛行器跨流域氣動力及熱問題的統(tǒng)一的算法應用研究途徑。

高超聲速飛行器； 跨流域； 氣動熱力學； Boltzmann模型方程； 離散速度坐標法； 統(tǒng)一算法； 并行計算

近空間高超聲速飛行器往返大氣層時，先后經(jīng)歷連續(xù)流、近連續(xù)滑移流、過渡流以及高稀薄流[1]等，不同流域氣體的熱力學性質(zhì)及繞流狀態(tài)互不相同。特別是位于連續(xù)流與高稀薄自由分子流之間的過渡區(qū)流動，無論是從實驗技術(shù)角度還是數(shù)值方法、理論描述與計算分析等層面來說，均是難以處理的流動。近空間高超聲速飛行器跨大氣層飛行，因其高馬赫數(shù)、低環(huán)境大氣壓力等特點，不僅具有嚴重的稀薄氣體非平衡流動熱力學效應，而且由于各部位幾何尺寸差異較大，在飛行器的某些部位甚至會出現(xiàn)連續(xù)流、近連續(xù)滑移流或稀薄過渡流多流區(qū)共存的混合流動現(xiàn)象。如何準確模擬近連續(xù)滑移流區(qū)及過渡流流區(qū)的高超聲速飛行器氣動力/熱環(huán)境，已成為近空間飛行器研制的核心問題[2-3]，亟需相關(guān)基礎(chǔ)理論與數(shù)值計算方法的研究與發(fā)展。

氣體分子運動論(氣體動理論)的基本方程——Boltzmann方程[4]可描述各個流域氣體分子的輸運現(xiàn)象，并能提供任何遠離局部熱平衡態(tài)的氣體流動的詳細特征，使之用于從航天軌道器飛行模擬到微不可見的同位素分離等領(lǐng)域研究。Boltzmann方程是研究包括各種微尺度流動在內(nèi)的氣體動力學問題的一個主要根據(jù)，求解該方程的最大困難在于其碰撞項的高度非線性和高維積分-微分屬性等復雜多尺度剛性問題，因此，精確求解Boltzmann方程一直未成現(xiàn)實。隨著流體力學的發(fā)展，國際上已有眾多學者基于質(zhì)量、動量以及能量的守恒定律，采用數(shù)學上較簡單的統(tǒng)計和碰撞松弛模型來代替Boltzmann方程碰撞項，提出了許多類似于Boltzmann方程各階矩的氣體運動論模型方程[5-7]。根據(jù)文獻[7]～文獻[11]的研究可知，通過求解能夠表征原始Boltzmann方程碰撞松弛和統(tǒng)計物理特性的氣體動理論模型方程，可望得到一個更為經(jīng)濟有效的數(shù)值計算方法，以解決飛行器從高稀薄自由分子流到連續(xù)流不同流域的高超聲速流動問題。

為了探索跨流域氣體流動問題的一體化模擬方法，根據(jù)Boltzmann方程的研究現(xiàn)狀與發(fā)展趨勢，文獻[9]、文獻[11]及文獻[12]通過研究確立了統(tǒng)一描述各流域微觀分子輸運現(xiàn)象的Boltzmann模型方程，提出并應用了離散速度坐標法對氣體分子速度分布函數(shù)進行數(shù)值離散，研究發(fā)展了可用于速度空間宏觀流動取矩的離散速度數(shù)值積分法。在此基礎(chǔ)上，將計算流體力學的有限差分方法推廣到了基于時間、位置空間和速度空間的Boltzmann模型方程并進行數(shù)值求解，建立了可模擬各流域一維、二維及三維氣體流動問題的統(tǒng)一算法(GKUA)理論與系列計算技術(shù)。該方法是在描述各流域氣體流動輸運現(xiàn)象的分子速度分布函數(shù)方程的基礎(chǔ)上直接進行求解，因此，首先需要利用氣體動理論離散速度坐標法對氣體分子的速度空間實施離散降維，消除速度分布函數(shù)對于速度空間的連續(xù)依賴性；然后，在各個離散速度坐標點處，數(shù)值求解時間與位置空間具有非線性源項的雙曲型守恒方程，從而使得各個離散速度坐標點之間的計算具有良好的并行獨立性。為了從根本上解決統(tǒng)一算法用于三維復雜繞流計算需要大量計算機內(nèi)存的問題，必須使用現(xiàn)代高性能、大內(nèi)存、多處理器(CPU)的并行計算機資源。

本文從研究如何求解描述各流域復雜高超聲速流動輸運現(xiàn)象統(tǒng)一的Boltzmann模型方程出發(fā)，引入?yún)^(qū)域分解并行化原理，建立氣體動理論大規(guī)模并行算法，并用于求解從稀薄流到連續(xù)流的各流域三維復雜高超聲速流動問題；通過對跨流域衛(wèi)星體飛行器的高超聲速氣動力/熱繞流問題開展并行計算，證實該并行算法的可靠性；在此基礎(chǔ)上，開展翼身組合大尺度復雜飛行器在近空間飛行環(huán)境下，極高超聲速繞流現(xiàn)象與流動機理并行計算分析的應用研究。

1 建立氣體動理論數(shù)值算式

基于相關(guān)文獻的敘述[9,11-12]，引入氣體分子運動論碰撞間隔理論[4]，將氣體分子碰撞松弛參數(shù)ν和當?shù)仄胶鈶B(tài)分布函數(shù)fN與各流域流態(tài)控制參數(shù)、宏觀流動物理量、氣體黏性輸運特性、熱力學效應、氣體分子相互作用規(guī)則以及分子模型聯(lián)系起來，建立Boltzmann方程碰撞積分的簡化式，提出適于表征各流域復雜流動輸運現(xiàn)象統(tǒng)一的Boltzmann模型方程，其無量綱形式為

(1)

式中：

在獲得氣體分子速度分布函數(shù)f與位置空間、速度空間和時間的變化率關(guān)系的基礎(chǔ)上，式(1)從分子運動論的角度給出了氣體的統(tǒng)計描述[4,13]。氣體動力學中感興趣的各種宏觀物理量，如氣體密度、流動速度、溫度、壓力、應力張量、熱流矢量等，可分別通過式(2)～式(6)求得。

(2)

(3)

(4)

(5)

(6)

式中：各物理量下標取值為1、2、3，分別表示x、y、z這3個方向的分量。

應用氣體動理論離散速度坐標法(DVOM)[8-15]可將單一的速度分布函數(shù)方程式(1)轉(zhuǎn)化為在各個離散速度坐標點處基于時間和位置空間(ξ，η,ζ)的具有非線性源項的雙曲型方程組，其守恒形式為

(7)

式中：

令

根據(jù)氣體分子碰撞松馳與對流運動的非定常特點，可將計算流體力學算子分裂有限差分方法應用到基于位置空間、速度空間和時間的七維問題中。運用二階Runge-Kutta方法和基于原始變量的無波動、無自由參數(shù)的耗散差分格式(NND格式)[14]，建立直接求解分子速度分布函數(shù)的氣體動理論耦合迭代數(shù)值算式為

(8)

式中：LS(·)、Lζ(·)、Lη(·)與Lξ(·)分別為式(7)中的源項與沿ζ、η、ξ這3個方向的對流項差分算子。

上述差分格式可分裂為四步進行：

式中：δξ、δη、δζ與δξ2、δη2、δζ2分別為對應的一階與二階差分；Δt為時間步長，由式(9)所示的格式穩(wěn)定條件Δts確定。

(9)

式中：Δξ、Δη和Δζ分別為計算平面內(nèi)沿ξ、η、ζ這3個方向的網(wǎng)格尺度；CFL為時間步長調(diào)節(jié)系數(shù)，可取CFL=0.99。

2 建立跨流域高超聲速氣動力/熱繞流計算模型

由于一般大尺度高超聲速飛行器外形復雜，各部位的幾何尺寸差異很大，在往返大氣層跨流域飛行時，會經(jīng)過近連續(xù)滑移流區(qū)和稀薄過渡流區(qū)，并會在飛行器的不同繞流區(qū)域因分子碰撞頻率的差別巨大而出現(xiàn)連續(xù)流、過渡流或高稀薄流共存的混合流動現(xiàn)象，從而導致物面邊界條件的表述方式對氣動力/熱特性的影響較大。

基于Boltzmann模型方程的氣體運動論統(tǒng)一算法需要求解的是物理空間與速度空間內(nèi)各個離散網(wǎng)格點處的氣體分子速度分布函數(shù)，即求解氣體分子速度分布函數(shù)是建立與實現(xiàn)各類邊界條件與物面氣動熱力學特性數(shù)學模型的必要條件[11-12]，從而避免了應用直接模擬蒙特卡羅方法(DSMC)對純粹微觀粒子進行隨機統(tǒng)計模擬時產(chǎn)生的統(tǒng)計漲落與對低Knudsen數(shù)近連續(xù)滑移過渡流難以進行模擬計算的問題，同時也避開了Navier-Stokes解算器純粹使用宏觀流動量進行邊界表述而導致不同位置的稀薄效應強弱不同、并難于用固定關(guān)系式數(shù)值實現(xiàn)的不足。

根據(jù)第1節(jié)內(nèi)容可知，通過式(8)將各個離散速度坐標點處的氣體分子速度分布函數(shù)進行數(shù)值求解后，須通過離散速度數(shù)值積分，才能及時更新任一時刻物理空間內(nèi)各點的宏觀流動參數(shù)。因此，研究可用于速度空間宏觀流動取矩的離散速度數(shù)值積分方法就顯得格外重要。為了實現(xiàn)高超聲速繞流模擬，提出了基于Legendre多項式的根為積分結(jié)點的Gauss-Legendre復合積分方法，用于解決高超聲速宏觀流動參數(shù)的動態(tài)確定與演化更新問題。Gauss型積分[15]的實質(zhì)是通過選取一套與積分規(guī)則相適應的積分結(jié)點Vσ與積分權(quán)系數(shù)Wσ，使定義在給定區(qū)域D內(nèi)的積分能夠被最佳離散求和，即

(10)

于是，飛行器繞流的各種宏觀流動量均可基于離散速度空間坐標(Vx σ,Vy δ,Vz θ)由所求離散速度分布函數(shù)fσ,δ,θ的三重求和計算確定[11,16]。下面僅列出其中幾個計算公式：

(11)

式中：Wσ、Wδ和Wθ分別為離散速度數(shù)值積分法中離散速度坐標點(Vx σ,Vy δ,Vz θ)所對應的積分權(quán)重因子，分別僅與σ、δ、θ有關(guān)。

根據(jù)上述離散速度數(shù)值積分方法，可獲得任意t時刻物面各網(wǎng)格點(xw,yw,zw)的氣動力/熱特性，由此可求得整個飛行器繞流空氣動力學特征，如當?shù)匚锩鏌崃魇噶縬w在x、y、z這3個方向上的分量值為

(12)

式中：Ma∞為來流馬赫數(shù)；ρ∞為來流氣體密度；γ為比熱比。

將垂直于飛行器物面且沿外法向nw方向的熱流分布記為qnw，則

qnw=qxwnxw+qywnyw+qzwnzw

(13)

式中：nxw、nyw和nzw分別為nw在x、y、z這3個方向上的分量。

3 建立氣體動理論大規(guī)模并行算法

三維繞流氣體動理論數(shù)值算法的計算空間是由離散速度空間和位置空間組成的六維相空間。算法可分成2個主要部分：①利用式(8)求解離散速度分布函數(shù)fσ,δ,θ；②基于已求得的分布函數(shù)，應用氣體動理論離散速度數(shù)值積分法[9,11,12]確定位置空間各點(x，y，z)處的飛行器繞流氣動力/熱環(huán)境宏觀流動參數(shù)。雖然問題的整個求解空間是由位置子空間和速度子空間組成的六維空間，但僅有速度分布函數(shù)fσ,δ,θ是定義在整個六維空間的，而所涉及的其他變量則只定義在位置子空間(x，y，z)或速度子空間(σ,δ,θ)中。因此，算法中的各模塊或基于離散位置空間或基于離散速度空間，彼此之間相互獨立，具有良好的并行獨立性。

從區(qū)域分解并行化原理出發(fā)，根據(jù)求解空間可將算法分為位置空間分解策略、速度空間分解策略以及基于位置空間與速度空間的混合分解策略。根據(jù)文獻[11]和文獻[12]的研究，模擬跨流域的高超聲速高馬赫數(shù)繞流問題，會需要大量的離散速度坐標點數(shù)目，但算法對離散速度空間具有很好的數(shù)據(jù)并行特點。因此，本文在HPF數(shù)據(jù)并行程序設(shè)計環(huán)境基礎(chǔ)上[16]，使用求解偏微分方程常用的區(qū)域分解方法來研究在離散速度空間并行分解的氣體動理論并行算法策略。

假設(shè)Ω為三維繞流氣體動理論數(shù)值模擬的求解空間，Ωr為位置空間(x，y，z)，ΩV為離散速度空間(σ,δ,θ)，則有

Ω=Ωr×ΩV

(14)

一般地，設(shè)處理器數(shù)為Np，則可將Ω分解為Np個子空間Ωi，并滿足

(15)

同時，將子空間Ωi的數(shù)據(jù)映射到相應的處理器PEi中。對離散速度空間ΩV實施并行分解策略，按塊布局或循環(huán)布局將ΩV分解成Np個子空間ΩVi，并滿足

則有

Ωi=Ωr×ΩVi

(16)

變量依賴關(guān)系分析是并行識別的基礎(chǔ)，也是用于指導區(qū)域分解策略的理論依據(jù)之一。通過對各階段算法過程進行變量依賴關(guān)系分析后可得出，在求解離散速度分布函數(shù)fσ,δ,θ的差分格式中，位置空間Ωr的各維方向都存在數(shù)據(jù)相關(guān)性，而離散速度子空間ΩV的各維方向則毫無數(shù)據(jù)相關(guān)。對于使用離散速度數(shù)值積分法對氣體分子離散速度分布函數(shù)進行矩積分并確定飛行器繞流宏觀流動參數(shù)的過程，在位置空間Ωr的各維方向上都不存在數(shù)據(jù)相關(guān)性，而在離散速度空間ΩV的各維方向上則以歸約形式體現(xiàn)了問題的相關(guān)性。因此，對于ΩV分解策略，在計算離散速度分布函數(shù)fσ,δ,θ的階段，如果忽略邊界處理，將會在子空間ΩV內(nèi)出現(xiàn)無數(shù)據(jù)通信的完全并行；而在計算宏觀流動參數(shù)這一階段，則需要在速度空間ΩV內(nèi)進行并行歸約計算，從而會產(chǎn)生數(shù)據(jù)通信。令CV為該策略下整個算法的總數(shù)據(jù)通信量，則依據(jù)二叉樹方式進行并行歸約，可得

CV=14NpσNpδNpθNiNjNk

(17)

其中，假設(shè)位置空間Ωr是Ni×Nj×Nk的網(wǎng)格陣列，處理器陣列為Npσ×Npδ×Npθ。

通常，將ΩV按三維方式等分分解為

(18)

令Ωσ,δ,θ=ΩV,σ,δ,θ×Ωr，并將Ωσ,δ,θ上的變量映射到處理器PEi,j,k中。為了分析離散速度空間的網(wǎng)格陣列與處理器陣列的分布關(guān)系，將式(17)改寫為

(19)

式中：V=NiNjNkNσNδNθ。

從式(19)中可以看出，為獲得較小的CV值，通常設(shè)置Npσ≤Nσ、Npδ≤Nδ、Npθ≤Nθ為宜。因此，基于離散速度空間ΩV的分解策略能夠較為高效地模擬離散速度坐標點數(shù)較多的情況?？缌饔蝻w行器復雜繞流問題通常都是馬赫數(shù)很高的流動，氣體分子速度分布函數(shù)所依賴的速度空間很大，需要覆蓋相應速度空間范圍的離散速度坐標點數(shù)會變得很多，因此，這種情況就特別適合使用ΩV分解策略。對于ΩV分解策略，處理器數(shù)最多可達到Nσ×Nθ×Nδ。對三維復雜高超聲速繞流問題來說，速度空間各離散速度坐標點的數(shù)目會隨流動馬赫數(shù)的增大而迅速增加，因此，基于離散速度空間ΩV并行分解策略能完全實現(xiàn)具有成千上萬CPU核甚至更高并行度的超大規(guī)模并行計算。這為求解Boltzmann模型方程的超大規(guī)模并行數(shù)值模擬奠定了理論基礎(chǔ)。

4 并行計算與結(jié)果分析

4.1 并行算法測試

為了考驗GKUA并行計算策略的可行性與不同規(guī)模CPU下并行計算的加速性能，現(xiàn)將其與文獻[17]中的國際同類研究進行比較，表1為賓西法尼亞大學在美國國家宇航局資助下，依靠高性能并行計算機對不同馬赫數(shù)下的一維激波結(jié)構(gòu)演化問題進行BGK模型方程求解，基于不同的位置空間網(wǎng)格劃分，分別在256、512、1 024個并行處理機進行并行計算的CPU開支情況。

表2為本文并行算法用于計算三維繞流問題所獲得的并行加速比，與文獻[17]中基于256個CPU下的BGK方程計算一維激波結(jié)構(gòu)并行加速比的情況比較。從表中可以看出，二者之間的并行計算加速比相當，表明了本文算法卓越的并行加速性能，不僅率先在國際上創(chuàng)建起求解各流域三維繞流問題的氣體動理論并行算法，而且該并行算法具有良好的并行可擴展性與并行效率。

表1 基于BGK方程256、512、1 024處理器(CPU)的并行計算[17]

Table 1 Parallel computation of the BGK model from using 256， 512 and 1 024 CPU[17]

SpatialnodeNumberofprocessorsMachinetypeCPUtime/sGFlops64K1024CM?23662932K1024CM?21852916K1024CM?21032616K1024CM?2146188K1024CM?2642132K512CM?23461516K512CM?2182158K512CM?21011316K256CM?23410788K256CM?218107416K256 CM?200320085

表2 基于256個CPU的并行計算加速比比較

圖1 統(tǒng)一算法(GKUA)在不同規(guī)模CPU下并行計算的加速比Fig.1 Speedup ratio of parallel computation for the gas-kinetic unified algorithm (GKUA) with different scale of CPU

圖2 統(tǒng)一算法在512～32 768 CPU下的并行計算加速比Fig.2 Speedup ratio of parallel computation for GKUA with 512-32 768 CPU

為了進一步驗證本文求解Boltzmann模型方程的并行算法在不同規(guī)模CPU下并行計算的加速性能，圖1(a)和圖1(b)分別給出了在64～1 024和1 440～7 920 CPU下的并行計算加速比。圖2為統(tǒng)一算法在眾核異構(gòu)并行計算機系統(tǒng)經(jīng)OpenACC眾核程序并行移植優(yōu)化后得到的512～32 768個處理器進程的并行計算加速比，從圖中可以看出，本文算法不僅在中等規(guī)模CPU數(shù)并行計算時，其加速比隨處理機數(shù)增加呈擬線性分布，而且在大規(guī)模CPU并行計算時，加速比仍能夠接近理想加速比。由此表明，算法程序在各CPU之間的通信效率和負載效率很高，證實了本文采用的并行計算策略是高效可行的，并能實現(xiàn)更高并行度的超大規(guī)模并行計算。

良好的并行加速性能可以保證在較高并行效率的情況下，通過增加處理機數(shù)來大大增加計算規(guī)模，使依靠傳統(tǒng)計算條件難以解決的各流域三維復雜高超聲速氣動力熱繞流問題計算求解成為可能，也足以證明發(fā)展國家高性能大規(guī)模并行數(shù)值模擬應用環(huán)境對流體力學數(shù)值方法研究與應用發(fā)展具有巨大的推動作用。

4.2 跨流域復雜飛行器高超聲速繞流并行計算

為了證明本文給出的跨流域復雜氣動力/熱繞流問題統(tǒng)一算法對近連續(xù)流到高稀薄流區(qū)高超聲速流動的模擬能力與并行計算的準確性，使用來自文獻[18]中的可重復使用衛(wèi)星體類球錐飛行器 (底部半徑R=503.5 mm，飛行器總長L=1 410mm，飛行器錐身段半錐角θ=11.4°)進行仿真驗證，并選取飛行器的底部半徑為特征尺度。擬定Ma∞=5時的2種繞流狀態(tài)：H=70.8km、Kn∞=0.002、雷諾數(shù)Re∞=4 074.37；H=105km、Kn∞=0.74、Re∞=10.19。在位置空間網(wǎng)格51×25×31和速度空間網(wǎng)格60×40×40的設(shè)置下，使用具有160MB/CPU的512個處理機并行計算。圖3和圖4分別為上述2種飛行高度下繞流流場軸對稱面內(nèi)的馬赫數(shù)等值線與繞流流場和物面流線結(jié)構(gòu)?？梢钥闯觯S著飛行器飛行高度從70.8km上升到105km，來流Knudsen數(shù)增大，氣體繞流稀薄效應迅速加重，飛行器繞流流場的受擾動區(qū)域大大增加；對于H=70.8km、Kn∞=0.002條件下的近連續(xù)滑移流，飛行器繞流在物體前面產(chǎn)生明晰清楚的脫體激波罩在飛行器前面，能夠很好地捕捉到包括駐點域、附著激波以及在飛行器底部出現(xiàn)真空區(qū)、尾渦流場結(jié)構(gòu)和飛行器后部流場再壓縮尾跡流動現(xiàn)象；對于H=105km、Kn∞=0.74條件下的高稀薄流，飛行器繞流并不存在激波等強間斷現(xiàn)象與背風漩渦回流結(jié)構(gòu)，而在飛行器周圍形成厚厚的強擾動區(qū)域，氣流附著物面流動，反映了稀薄氣體繞流流場完全不同于近連續(xù)流區(qū)繞流的面貌。

圖3 不同高度下的繞流流場馬赫數(shù)等值線分布Fig.3 Mach number contours distribution in flow field under different heights

圖4 不同高度下的繞流流場與物面流線結(jié)構(gòu)Fig.4 Flow field and stream structure under different heights

圖5 飛行器繞流迎風物面熱流分布的GKUA計算結(jié)果與DSMC模擬值比較Fig.5 Comparison of GKUA computational results and DSMC simulation of heat flux distribution along vehicles windward surface

圖5為GKUA計算該飛行器繞流迎風物面熱流分布與DSMC模擬值[18]定量化比較，計算狀態(tài)為：H=105km、Ma∞=5、Tw/T∞=1(Tw為物面溫度，T∞為來流溫度)。從圖中可看出，從球頭前駐點沿物面向后，由GKUA得到的不同物面角熱流計算值均與DSMC結(jié)果具有良好的吻合度，2種結(jié)果偏差范圍僅為0.27%～8.56%，而且GKUA得到的迎風物面熱流分布的非線性效應更明顯。

圖6為在Ma∞=10時，由GKUA計算上述衛(wèi)星體飛行器繞流的駐點線密度及溫度分布并與文獻[18]中結(jié)果的定量化比較，計算狀態(tài)為：H=75.9 km、Kn∞=0.004、Tw/T∞=1.64。從圖中可以看出，GKUA的計算值與文獻結(jié)果的變化規(guī)律吻合度很高，2種結(jié)果的計算偏差范圍僅為0.39%～3.26%。另外，圖中還顯示出由于設(shè)置了Tw/T∞=1.64的冷壁面，在飛行器前面一定區(qū)域會出現(xiàn)流動高溫區(qū)，但從整體上看，GKUA對駐點線流動參數(shù)的計算分辨率要優(yōu)于DSMC模擬值。

圖6 飛行器繞流駐點線密度、溫度分布的GKUA計算值與文獻[18]結(jié)果的比較Fig.6 Comparison of GKUA calculation and results of Ref.[18] of stagnation line density and temperature distribution for vehicle shape

圖7為該飛行器繞流的物面壓力系數(shù)Cp分布與文獻[18]中的結(jié)果比較，計算狀態(tài)為：H=75.9km、Ma∞=10、Re∞=4 074.5、Tw/T∞=1。從圖中可以看出，GKUA的計算結(jié)果與文獻[18]中結(jié)果的變化趨勢較為一致。其中，在離球頭前駐點較遠的物面，GKUA的計算結(jié)果較文獻結(jié)果稍大些，文獻結(jié)果所表現(xiàn)非線性效應要差些。

圖7 飛行器繞流物面壓力系數(shù)分布的GKUA計算值與文獻[18]結(jié)果的比較Fig.7 Comparison of GKUA computation and results of Ref.[18] of pressure coefficient distribution along vehicle surface

作為翼身組合大尺度復雜飛行器在近空間飛行環(huán)境極高超聲速下的繞流算例，本文開展了H=60～90 km、Ma∞=20～25、不同攻角下的繞流狀態(tài)大規(guī)模并行計算與流動分析。圖8為該飛行器在稀薄過渡流區(qū)繞流狀態(tài)(H=90km、Kn∞=0.012、Ma∞=25、α=20°、Re∞=3 153、Tw/T0=0.05、Pr=0.72)下的軸對稱面內(nèi)流場馬赫數(shù)等值線分布，該狀態(tài)設(shè)置位置空間網(wǎng)格為101×101×81、速度空間網(wǎng)格為110×100×80，總計算內(nèi)存3 998GB，使用具有168MB/CPU的23 800個處理機并行計算。從圖中可以看出，對此稀薄效應很強的高高度高超聲速繞流，在飛行器前部流場存在劇烈的流動壓縮強擾動，離飛行器前緣駐點一定距離處產(chǎn)生厚厚的脫體激波層，并罩在飛行器上下一定距離的空間位置，形成較密集的馬赫數(shù)等值線過渡帶，而且有趣的是，該強擾動激波過渡帶輪廓與飛行器外形相似，反映了大尺度復雜飛行器稀薄過渡區(qū)特殊的繞流面貌。

圖8 翼身組合復雜飛行器繞流的馬赫數(shù)等值線Fig.8 Mach number contours of flow around complex wing-body combination shape vehicle

為了進一步揭示該飛行器繞流變化規(guī)律，圖9給出了飛行器的表面流線與繞流流場結(jié)構(gòu)。經(jīng)分析表明，雖然其飛行高度較高，已進入稀薄過渡流區(qū)，但由于飛行器特征尺寸很大，達數(shù)米量級，使得此種狀態(tài)的來流Knudsen數(shù)較小(Kn∞=0.012)，氣體繞流整體上仍呈現(xiàn)為近連續(xù)過渡流，在飛行器腹部下面的氣流形成一撮厚厚的壓縮激波流線。該壓縮匯聚流線帶起源于毫米量級的飛行器前部端頭強擾動斜激波，附著在飛行器細長翼身一定位置，并延續(xù)至飛行器遠后方，同時氣流遇到飛行器后端面會迅速膨脹，使一部分氣流改變方向，與來自飛行器肩頂部氣流壓縮匯聚，向飛行器遠后方流去。

圖9 翼身組合復雜飛行器表面流線與繞流流場結(jié)構(gòu)Fig.9 Surface streamlines and flow structure around complex wing-body combination shape vehicle

由于本文發(fā)展的氣體動理論統(tǒng)一算法通過直接求解跟蹤氣體分子速度分布函數(shù)的時間演化，來更新物理空間各點處的宏觀流動參數(shù)，使得通常所說的非連續(xù)流區(qū)繞流物面所出現(xiàn)的速度滑移、溫度跳躍等流動非平衡稀薄效應及復雜飛行器不同部位出現(xiàn)多流區(qū)混合繞流物面力熱流動信息能被準確捕捉并得到自然模擬，避免了DSMC方法針對微觀粒子隨機模擬易產(chǎn)生統(tǒng)計漲落與Navier-Stokes解算器使用宏觀流動量難于用固定關(guān)系式表示不同部位稀薄效應強弱的局限性。

為了分析并揭示該翼身組合大尺度復雜飛行器近空間飛行時的繞流物面氣動力和熱變化規(guī)律，在擬定近連續(xù)流區(qū)繞流狀態(tài)下(H=60km、Kn∞=0.000 14、Ma∞=20、α=12°、Re∞=2.164×105、Tw/T0=0.25、Pr=0.72、γ=Cp/CV=1.4)進行大規(guī)模并行計算。圖10為物面熱流系數(shù)與邊緣線弧長位置的變化關(guān)系。首次計算發(fā)現(xiàn)，該類飛行器后端面邊緣線不同位置的物面熱流在不同流區(qū)和不同飛行高度時，呈現(xiàn)出不同的變化規(guī)律，尾部兩端水平舵邊緣存在較高熱流，熱流最大值發(fā)生在左右水平舵外側(cè)拐角處，為駐點熱流的1/6量級。該發(fā)現(xiàn)對此類飛行器的工程設(shè)計具有重要意義與應用價值。

圖10 翼身組合復雜飛行器繞流后端邊緣線物面熱流分布Fig.10 Surface heat flux distribution of flow along afterbody edge of complex wing-body combination shape vehicle

5 結(jié) 論

1) 本文在確立描述各流域復雜高超聲速流動輸運現(xiàn)象統(tǒng)一的Boltzmann模型方程基礎(chǔ)上，借助大規(guī)模并行計算機系統(tǒng)，應用氣體動理論離散速度坐標法與顯式算子分裂規(guī)則，構(gòu)造直接求解氣體分子速度分布函數(shù)演化更新的氣體動理論耦合迭代數(shù)值格式；使用離散速度空間區(qū)域分解并行化策略研究新型的氣體動理論數(shù)值算法并行方案，建立了求解Boltzmann模型方程可靠模擬稀薄流到連續(xù)流跨流域復雜高超聲速氣動力/熱繞流問題統(tǒng)一算法與穩(wěn)定運行于數(shù)千數(shù)萬以至更大規(guī)模CPU高性能并行計算應用研究平臺。

2) 基于離散速度空間并行分布策略在整個算法過程具有較好的靜態(tài)穩(wěn)定性，不存在臨界與邊界離散通信，算法通信效率與并行計算效率較高，不僅負載平衡好，基本達到了線性加速的并行效果，且并行化代價較低，具有良好的并行可擴展性，顯示出求解Boltzmann模型方程的氣體運動論統(tǒng)一算法具有相當高性能大規(guī)模并行計算優(yōu)勢。

3) 對跨流域不同Knudsen數(shù)、高低不同馬赫數(shù)、不同攻角繞流問題，如可重復使用衛(wèi)星球錐體、翼身組合大尺度復雜飛行器跨流域高超聲速氣動力/熱繞流問題并行計算與算法驗證，計算結(jié)果與典型文獻、DSMC模擬值及理論分析吻合較好，揭示了稀薄流到連續(xù)流不同流域復雜高超聲速繞流現(xiàn)象與物面氣動力熱變化規(guī)律。

4) 為了將統(tǒng)一算法應用于大尺度復雜飛行器跨越飛行各流域高超聲速氣動力/熱繞流問題工程實際，特別需要大規(guī)模并行計算機資源，顯示出國家高性能并行計算環(huán)境對發(fā)展流體力學數(shù)值方法，解決跨流域空氣動力學問題的意義和作用。

5) 通過對尖前緣翼身組合復雜飛行器近空間飛行環(huán)境高超聲速繞流問題大規(guī)模并行計算，揭示了近連續(xù)過渡流區(qū)特殊繞流現(xiàn)象與氣動熱環(huán)境變化規(guī)律。為進一步發(fā)展求解尖前緣大尺度復雜飛行器、大型航天器從外層空間再入跨流域高超聲速繞流氣動力、熱問題高效數(shù)值模擬方法指明了方向。

致 謝

本文工作得到中國空氣動力研究與發(fā)展中心張涵信院士、無錫江南計算技術(shù)研究所徐金秀高級工程師、第一作者課題組彭傲平、方明等的支持幫助，部分計算在國家超級計算濟南中心、總參三部北方計算中心進行，特此表示感謝。

[1] Tsien H S. Superaerodynamics，mechanics of rarefied gases[J]. Journal of Aeronautics Science， 1946， 13(12)： 653-664.

[2] Whitehead A H, Jr. NASP aerodynamics， AIAA-1989-5013[R]. Reston： AIAA， 1989.

[3] Zhuang F G， Cui E J， Zhang H X. Some development of future spacecrafts and aerodynamics tasks[C]∥Proceedings of First Aerodynamics and Aerothermodynamics. Beijing： National Defense Industry Press， 2006： 1-12 (in Chinese). 莊逢甘， 崔爾杰， 張涵信. 未來空間飛行器的某些發(fā)展和空氣動力學的任務(wù)[C]∥中國第一屆近代空氣動力學與氣動熱力學會議論文集. 北京： 國防工業(yè)出版社， 2006： 1-12.

[4] Chapmann S， Cowling T G. The mathematical theory of non-uniform gases[M]. 3rd ed. Cambridge： Cambridge University Press， 1990.

[5] Bhatnagar P L， Gross E P， Krook M． A model for collision processes in gases： I．small amplitude processes is charged and neutral one-component system[J]．Physics of Review， 1954， 94： 511-525.

[6] Abe T， Oguchi H． A hierarchy kinetic model and its applications[C]∥Potter J I． Rarefied Gas Dynamics. Reston： AIAA， 1977： 781-793.

[7] Shakhov E M． Kinetic model equations and numerical results[C]∥Oguchi H． Proceedings of 14th International Symposium on Rarefied Gas Dynamics. Tokyo： University of Tokyo Press， 1984: 137-148.

[8] Yang J Y， Huang J C. Rarefied flow computations using nonlinear model Boltzmann equtions[J]. Journal of Computational Physics， 1995， 120(2)： 323-339.

[9] Li Z H， Zhang H X. Study on gas kinetic algorithm for flows from rarefied transition to continuum[J]. Acta Aerodynamica Sinica， 2000， 18(3)： 255-263 (in Chinese). 李志輝， 張涵信． 稀薄流到連續(xù)流的氣體運動論統(tǒng)一數(shù)值算法初步研究[J]． 空氣動力學學報， 2000， 18(3)： 255-263.

[10] Mieussens L. Discrete-velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries[J]． Journal of Computational Physics， 2000， 162(2)： 429-466.

[11] Li Z H， Zhang H X. Study on gas kinetic unified algorithm for flows from rarefied transition to continuum[J]. Journal of Computational Physics， 2004， 193(2)： 708-738.

[12] Li Z H， Zhang H X. Study on gas kinetic numerical algorithm using Boltzmann model equation[J]．Advances in Mechanics， 2005， 35(4)： 559-576 (in Chinese). 李志輝， 張涵信. 基于Boltzmann模型方程的氣體運動論統(tǒng)一算法研究[J]． 力學進展， 2005， 35(4)： 559-576.

[13] Cercignani C. Kinetic theories and the Boltzmann equation[M]. Berlin： Springer-Verlag， 1984.

[14] Zhang H X， Shen M Y. Computational fluid dynamics—principle and application of difference method[M]. Beijing： National Defense Industry Press， 2003 (in Chinese). 張涵信， 沈孟育．計算流體力學——差分方法的原理和應用[M]. 北京： 國防工業(yè)出版社， 2003.

[15] Golub G H，Welsch J. Calculation of Gauss quadrature rules [J]. Mathematics of Computation， 1969， 23(106)： 221-230.

[16] Li Z H， Zhang H X. HPF parallel computation based on Boltzmann model equation for flows past complex body from various flow regimes[J]. Acta Aeronautica et Astronautica Sinica， 2006， 27(2)： 175-181 (in Chinese). 李志輝， 張涵信. 基于Boltzmann模型方程不同流區(qū)復雜三維繞流HPF并行計算[J]. 航空學報， 2006， 27(2)： 175-181.

[17] Long L N， Kamon M， Myczkowski J. A massively parallel algorithm to solve the Boltzmann (BGK) equation，AIAA-1992-0563[R]. Reston： AIAA， 1992.

[18] Sharipov F. Hypersonic flow of rarefied gas near the Brazilian satellite during its reentry into atmosphere[J]. Brazilian Journal of Physics， 2003， 33(2)： 398-405.

Tel: 010-82330957

E-mail: zhli0097@x263.net

吳俊林 男，碩士，助理研究員。主要研究方向：稀薄氣體動力學。

Tel: 0816-2465261

E-mail: wujunlin130@aliyun.com

蔣新宇 男，碩士，助理工程師。主要研究方向：稀薄氣體動力學。

Tel: 0816-2465261

E-mail: janxy1987@163.com

馬強 男，博士后。主要研究方向：氣動熱力學繞流環(huán)境結(jié)構(gòu)耦合計算。

Tel: 010-82330957

E-mail: maqiang@lsec.cc.ac.cn

*Corresponding author. Tel.： 010-82330957 E-mail: zhli0097@x263.net

A massively parallel algorithm for hypersonic covering various flow regimes to solve Boltzmann model equation

LI Zhihui1,2,*, WU Junlin1, JIANG Xinyu1, MA Qiang1,2

1.HypervelocityAerodynamicsInstitute,ChinaAerodynamicsResearch&DevelopmentCenter,Mianyang621000,China2.NationalLaboratoryforComputationalFluidDynamics,Beijing100191,China

The unified equation on the molecular velocity distribution function is presented for describing complex hypersonic flow transport phenomena covering various flow regimes by the computable model of Boltzmann collision integral. The discrete velocity ordinate method is used to discretize and reduce velocity space dimensionality of the velocity distribution function, and the gas-kinetic numerical schemes of coupling iteration are constructed directly to solve the molecular velocity distribution function. The computing models of hypersonic aerothermodynamics for the complex vehicles are developed by the evolution and updating based on the molecular velocity distribution function. The new parallel strategy of the gas-kinetic numerical algorithm is established by using the parallelizing technique of domain decomposition with the analysis from variable dependency relations, data communication and parallel expansibility. The data parallel distribution and parallel implementation are researched, the large-scale parallel program design is carried out and then the high-performance parallel algorithm has been established to simulate the hypersonic flow problems around complex vehicles covering various flow regimes, which can run stably in the tens of thousands of CPU or more scale. The hypersonic aerothermodynamics problems from high rarefied transition to continuum flow regimes around three-dimensional sphere-cone satellite body and complex wing-body combination shape with various Knudsen numbers, different Mach numbers, and diverse flying of angles have been computed and verified in high-performance computer with massive scale parallel. The computed results are found in high resolution of the flow fields and good agreement with the related reference experimental data, direct simulation Monte Carlo (DSMC) and theoretical predictions, and the hypersonic complex flow mechanism and changing laws are revealed. It is probably practical that the applying research approaches of the gas-kinetic unified algorithm can be provided to simulate complex hypersonic flow problems covering the whole of flow regimes.

hypersonic vehicles； covering various flow regimes； aerothermodynamics； Bolztamnn model equation； discrete velocity oridinate method； gas-kinetic unified algorithm； parallel computation

2014-07-02； Revised： 2014-09-16； Accepted： 2014-10-20； Published online： 2014-10-20 10:09

s： National Basic Research Program of China(2014CB744100)； National Natural Science Foundation of China (91016027，91130018，11325212)； National Defense Basic Research Program (51313030104)

2014-07-02; 退修日期： 2014-09-16; 錄用日期： 2014-10-20; 網(wǎng)絡(luò)出版時間： 2014-10-20 10：09

www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0219.html

國家“973”計劃 (2014CB744100)； 國家自然科學基金 (91016027， 91130018， 11325212)； 國防基礎(chǔ)科研項目 (51313030104)

Li Z H， Wu J L， Jiang X Y， et al. A massively parallel algorithm for hypersonic covering various flow regimes to solve Boltzmann model equation[J]. Acta Aeronautica et Astronautica Sinica， 2015， 36(1)： 201-212. 李志輝， 吳俊林， 蔣新宇， 等． 跨流域高超聲速繞流Boltzmann模型方程并行算法[J]． 航空學報， 2015， 36(1)： 201-212.

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2014.0219

V411.3

A

1000-6893(2015)01-0201-12

李志輝 男，博士，研究員，博士生導師。主要研究方向：稀薄氣體動力學與計算流體力學。

*通訊作者.Tel.： 010-82330957 E-mail: zhli0097@x263.net

URL： www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0219.html