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

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

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

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

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

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

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

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

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

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

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

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

本文從研究如何求解描述各流域復雜高超聲速流動輸運現象統一的Boltzmann模型方程出發，引入區域分解并行化原理，建立氣體動理論大規模并行算法，并用于求解從稀薄流到連續流的各流域三維復雜高超聲速流動問題；通過對跨流域衛星體飛行器的高超聲速氣動力/熱繞流問題開展并行計算，證實該并行算法的可靠性；在此基礎上，開展翼身組合大尺度復雜飛行器在近空間飛行環境下，極高超聲速繞流現象與流動機理并行計算分析的應用研究。

1 建立氣體動理論數值算式

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

(1)

式中：

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

(2)

(3)

(4)

(5)

(6)

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

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

(7)

式中：

令

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

(8)

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

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

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

(9)

式中：Δξ、Δη和Δζ分別為計算平面內沿ξ、η、ζ這3個方向的網格尺度；CFL為時間步長調節系數，可取CFL=0.99。

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

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

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

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

(10)

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

(11)

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

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

(12)

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

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

qnw=qxwnxw+qywnyw+qzwnzw

(13)

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

3 建立氣體動理論大規模并行算法

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

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

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

Ω=Ωr×ΩV

(14)

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

(15)

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

則有

Ωi=Ωr×ΩVi

(16)

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

CV=14NpσNpδNpθNiNjNk

(17)

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

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

(18)

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

(19)

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

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

4 并行計算與結果分析

4.1 并行算法測試

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

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

表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 統一算法(GKUA)在不同規模CPU下并行計算的加速比Fig.1 Speedup ratio of parallel computation for the gas-kinetic unified algorithm (GKUA) with different scale of CPU

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

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

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

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

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

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

圖4 不同高度下的繞流流場與物面流線結構Fig.4 Flow field and stream structure under different heights

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結 論

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

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

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

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

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

致 謝

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

[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). 莊逢甘， 崔爾杰， 張涵信. 未來空間飛行器的某些發展和空氣動力學的任務[C]∥中國第一屆近代空氣動力學與氣動熱力學會議論文集. 北京： 國防工業出版社， 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). 李志輝， 張涵信． 稀薄流到連續流的氣體運動論統一數值算法初步研究[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模型方程的氣體運動論統一算法研究[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]. 北京： 國防工業出版社， 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模型方程不同流區復雜三維繞流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

馬強 男，博士后。主要研究方向：氣動熱力學繞流環境結構耦合計算。

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; 網絡出版時間： 2014-10-20 10：09

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

國家“973”計劃 (2014CB744100)； 國家自然科學基金 (91016027， 91130018， 11325212)； 國防基礎科研項目 (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