氣體動理學格式及其在再入問題中的應用

李詩一, 張 潮, 譚 爽, 李啟兵, 符 松

(清華大學 航天航空學院, 北京 100084)

0 引 言

臨近空間飛行器、跨大氣層飛行器以及返回艙等高超聲速飛行器，在其再入過程中，流動狀態(tài)涵蓋從稀薄流到連續(xù)流，涉及到多尺度、跨流域問題，在近連續(xù)、連續(xù)流區(qū)還會出現(xiàn)黏性干擾、湍流等復雜流動現(xiàn)象，嚴重影響飛行器的氣動力和氣動熱特性。對其進行準確模擬具有非常高的挑戰(zhàn)性。

在高空稀薄流區(qū)，受限于連續(xù)介質(zhì)假設，基于Euler和Navier-Stokes方程的傳統(tǒng)CFD方法不再適用。直接模擬蒙特卡羅方法(DSMC)從分子平均自由程或者分子平均碰撞時間尺度出發(fā)，刻畫宏觀氣體流動的稀薄效應[1]。傳統(tǒng)的離散速度坐標法(DOM或DVM)基于算子分裂[2]，直接在高維空間離散Boltzmann模型方程或其模型方程。對于三維復雜飛行器再入各流域繞流問題，需要在三維位置空間和三維速度空間離散求解分布函數(shù)，計算量太大。為求解高超聲速飛行器面臨的復雜跨流域繞流問題，亟需發(fā)展高效的多尺度計算方法。

香港科技大學徐昆等發(fā)展的氣體動理學格式(GKS)[3-4]基于BGK類方程的解析解來計算網(wǎng)格界面上的數(shù)值通量，耦合了分子間相互碰撞與分子的自由運動，離散尺度僅取決于求解精度的要求，不再受限于平均自由程和平均碰撞時間[3]。對連續(xù)及近連續(xù)流，可以將分布函數(shù)簡化，得到逼近N-S方程的GKS-NS方法，避免了離散速度坐標，極大地減少了存儲量和計算量。GKS-NS在多種流動，特別是高速黏性流動中表現(xiàn)優(yōu)異[5]。

為了將GKS應用于實際再入問題，本文主要開展了以下幾個方面的工作。首先，耦合常用湍流模型對GKS-NS進行拓展，使之能有效模擬高超聲速飛行器飛行高度不太高時面臨的復雜高雷諾數(shù)高速湍流問題。其次，發(fā)展非結(jié)構(gòu)網(wǎng)格上的高精度多維GKS-NS格式，以對高超聲速復雜流動問題進行精細模擬。最后，對結(jié)合速度空間的離散發(fā)展起來的統(tǒng)一氣體動理學格式(UGKS)進行拓展及優(yōu)化，使之能夠模擬跨流域多尺度流動問題。

1 GKS簡介

氣體動理學格式基于Boltzmann方程的模型方程，BGK方程。它具有如下形式：

g=ρ(2πRT)-(K+3)/2e-[(u-U)2+ξ2]/(2RT)(1)

其中，f=f(x,u,ξ,t)為氣體分子速度分布函數(shù)，u為分子平動速度，ξ為內(nèi)自由度，內(nèi)自由度數(shù)為K。g為平衡態(tài)分布函數(shù)，ρ、U、T、E分別為氣體的密度、速度、溫度和總能，τ為分子平均碰撞時間，R為氣體常數(shù)。對f和BGK方程取矩，可以得到宏觀守恒量Q及守恒方程：

Ψ=[1,u,(u2+ξ2)/2]T。(2)

相比于宏觀方程，BGK方程通過分布函數(shù)的變化來刻畫分子的自由運動和相互碰撞的過程，具有更豐富的物理內(nèi)涵[6- 7]。

BGK方程(1)有如下形式的演化解：

e-t/τf0(x-ut,u,ξ)，

x′=x-u(t-t′)(3)

其中，f0是初始分布函數(shù)，g是隨時間演化的平衡態(tài)分布。上述演化解自動耦合了分子的自由運動和相互碰撞，這是確保基于該演化解的GKS成為真正多尺度方法的關鍵[3]。

GKS采用有限體積法，通過演化解(3)來計算單元界面上的數(shù)值通量，從而對宏觀量和分布函數(shù)進行更新：

其中，Sijk為離散物理空間上網(wǎng)格單元(i,j,k)的界面面積，Ωijk為單元體積。演化解中的平衡態(tài)分布g可根據(jù)格式精度需要由單元界面中心關于時、空的Taylor展開逼近。

UGKS同時更新宏觀量(4)和分布函數(shù)(5)，適用于從連續(xù)流到稀薄流區(qū)的全流域。并且，由于通過演化解來計算通量，可以刻畫出不同網(wǎng)格尺度上、不同時間間隔內(nèi)分子的自由運動和相互碰撞的物理過程，網(wǎng)格尺度、時間步長可以不受分子平均自由程和平均碰撞時間的限制，因而是真正的多尺度方法，模擬再入過程中的跨尺度、跨流域流動問題時具有很高的求解效率。

在近連續(xù)流區(qū)，初始分布函數(shù)可以直接由一階Chapmann-Enskog展開逼近，

此時演化解可完全由平衡態(tài)分布及其時、空導數(shù)確定，因此在宏觀量的更新過程中，式(4)，通量F僅與Q及其空間導數(shù)有關[4]。更為重要的是，這種情況下無需直接更新分布函數(shù)，即式(5)，因而避免了速度和內(nèi)自由度空間(u,ξ)的離散，極大地減少了計算量。這樣建立的GKS-NS格式仍然耦合處理了無黏與黏性項，保證了其在高超聲速黏性模擬中的優(yōu)異性能。

2 GKS的拓展及應用

2.1 高超聲速湍流模擬

GKS-NS由于內(nèi)含黏性并且耦合處理無黏與黏性項，在高超聲速黏性流動模擬中表現(xiàn)優(yōu)異。在工程高雷諾數(shù)問題中，為了高效地模擬湍流場，可進一步基于擴展BGK方程[9]發(fā)展耦合湍流模式的拓展GKS[10-11]。擴展BGK方程為：

τeff=τ+τt為流動的有效松弛時間，其中湍流松弛時間通過湍流黏性系數(shù)定義τt=μt/(p+2k/3)。通過μt，GKS可直接與常用的湍流模式耦合，在積分尺度上描述湍流流動。τt定義中k為流場湍動能，通常在高超聲速流動中不能忽略。在該框架下建立的拓展GKS-NS中，擴展BGK方程的多尺度演化過程退化為格式的自適應耗散機制，使格式可在間斷處保持強魯棒性，同時在光滑流場區(qū)保持較高的分辨率。另外，大部分湍流模式中需要求解湍流量的輸運方程。本文中通過Li[8]等發(fā)展的帶有標量輸運的GKS-NS格式，以耦合方式求解湍流量輸運方程，湍流量的數(shù)值精度和耗散機制與守恒量保持一致。

拓展GKS中采用了多種湍流模型，包括零方程模型B-L、一方程模型S-A、兩方程模型k-ε、k-ω、SST以及Menter轉(zhuǎn)捩模式等。這里給出結(jié)合含有可壓縮修正的k-ω-LS模式對高超聲速壓縮拐角流動的計算結(jié)果。來流馬赫數(shù)為9.22，斜坡角度為θ=24°[12]。在斜坡拐角較大時，分離區(qū)的大小以及壁面熱流的準確預測對數(shù)值計算非常具有挑戰(zhàn)性。圖1中為GKS的計算結(jié)果及其與已有實驗測量及數(shù)值結(jié)果的比較。當計算采用相同湍流模式時，拓展GKS的計算結(jié)果比已有數(shù)值結(jié)果[13]與實驗符合得更好。其中，熱流峰值位置基本與實驗一致，熱流峰值大小與實驗的誤差約為41%，精度高于參考計算結(jié)果。拓展GKS在高超聲速模擬中的優(yōu)勢來自于其內(nèi)涵的自適應耗散機制，而湍流量與守恒量耦合的求解方法也有助于提高計算精度。

圖1 高超聲速壓縮拐角壁面壓力和熱流分布(q∞=6.56 W/cm2)Fig.1 Normalized wall pressure and heat flux in hypersonic compression ramp flow

2.2 基于CPR框架的高精度GKS

高精度格式能夠以相對較少的計算量得到精度更高的數(shù)值結(jié)果，從而提高計算效率。近年來，Huynh和Wang等人所發(fā)展的通量重構(gòu)方法(FR/CPR)[14-15]為DG、SV/SD等經(jīng)典的高精度方法提供了一個統(tǒng)一框架。相比于傳統(tǒng)的DG方法，F(xiàn)R/CPR方法擺脫了數(shù)值積分的計算量，并且通過特定求解點和通量點的選取，能夠有效減少重構(gòu)的計算量，因而FR/CPR方法具有很高的效率。而考慮到通量演化，基于氣體動理學格式所構(gòu)造的高精度GKS-NS通量求解器耦合了無黏通量和黏性通量，并且包含了通量隨時間的演化，因而只需要單步推進即可同時達到時空高階精度。此外，GKS-NS具有真正的多維特性，可以直接構(gòu)造出通量在單元內(nèi)的空間分布。因此，結(jié)合高效的FR/CPR框架和高精度GKS-NS通量求解器，發(fā)展適用于飛行器復雜外形的非結(jié)構(gòu)網(wǎng)格上的新型高精度格式具有很好的應用前景。

其利用單元內(nèi)設置的若干求解點可以插值得到守恒量和通量在單元內(nèi)的分布?？紤]到單元界面兩側(cè)可能存在間斷，還需在每個單元界面上設置若干通量點，計算耦合相鄰單元影響的通量值，并由此對單元內(nèi)的通量分布進行修正，得到修正項δij，進而更新每個求解點上的守恒變量Qij。

為了計算每個求解點和通量點上的通量值，需要分別在三角形單元內(nèi)和單元界面上構(gòu)造出通量分布。首先對氣體分子速度分布函數(shù)進行高階泰勒展開，通過重構(gòu)得到的守恒量分布及其空間導數(shù)，以及相容條件可以計算出展開式中的各項系數(shù)。高精度GKS通量求解器所具有的多維特性使得所有求解點以及通量點上的通量值均可以同時獲得，而傳統(tǒng)的通量求解器則需要逐點計算通量。

為了計算含間斷的可壓縮流動問題，格式中引入了緊致型WENO限制器[16]。該限制器能夠較好地抑制間斷附近的振蕩，同時保證光滑區(qū)域的高精度。多種典型算例表明，新發(fā)展的時-空三階精度CPR-GKS在可壓縮流動中兼具高精度和較好的激波捕捉性能，并且計算效率能和原CPR格式相當。本文采用CPR-GKS計算了Re=200的黏性激波管問題。它包含了強激波和邊界層的相互作用，不僅要求格式具有良好的穩(wěn)健性，同時還對格式的精度和分辨率有較高的要求，因而是驗證格式在非定常黏性流動中捕捉激波和復雜流動結(jié)構(gòu)的經(jīng)典算例。圖2為計算得到的密度等值線圖，其中網(wǎng)格尺度為1/300，計算時間為t=1。由圖中可以看出，格式較好地捕捉到了激波和邊界層干擾產(chǎn)生的復雜流場，尤其是對渦結(jié)構(gòu)的捕捉展現(xiàn)了格式的高精度。計算出的中間渦結(jié)構(gòu)的高度約為0.172，與有限體積GKS得到的0.171吻合較好[18]。該有限體積GKS采用了五階重構(gòu)，時間精度為四階，且網(wǎng)格尺度為1/500。由此也進一步表明了基于內(nèi)點重構(gòu)的CPR框架的高精度。

圖2 黏性激波管：密度分布(t=1)Fig.2 Viscous shock tube: density distribution at t=1

2.3 適合大規(guī)模并行化計算的UGKS

UGKS能夠有效模擬跨流域稀薄流動，但需要同時求解式(4)和式(5)。除了需要離散物理空間，還需要對速度空間進行離散。一般地，模擬三維稀薄流動問題，相比于直接求解宏觀方程，如Navier-Stokes方程所適用的連續(xù)流問題，計算規(guī)模將有量級上的增加。對于高超聲速流動，由于分布函數(shù)更寬，速度空間的離散規(guī)模將進一步增加。因此，為了模擬高超聲速跨流域稀薄流動，迫切需要發(fā)展適合UGKS的大規(guī)模高效并行算法。

一方面，UGKS算法需要同時對物理空間和速度空間進行離散，核心部分為氣體分子速度分布函數(shù)在物理空間上的重構(gòu)、演化和投影，在演化和投影過程中涉及到分布函數(shù)在速度空間上的積分求矩?？紤]到UGKS算法本身在離散物理空間和速度空間上計算過程的差異，在并行算法的設計中可以對二者進行區(qū)分。另一方面，考慮到在實際三維問題中，需要離散三維速度空間，特別是在非定常高超聲速稀薄流模擬中，其存儲量尤其巨大，因此有必要對速度空間進行并行分塊。

本文發(fā)展的UGKS并行程序適合物理空間分塊搭接的結(jié)構(gòu)網(wǎng)格，且物理空間和速度空間同時并行分塊。為了加速定常問題收斂，還引入了當?shù)貢r間步長、隱式格式[24]等措施。在數(shù)據(jù)交換中，采用了減帶寬技術[17]優(yōu)化了分塊界面兩側(cè)不同物理分塊之間的宏觀物理量和分布函數(shù)的數(shù)據(jù)交換，有效減少了數(shù)據(jù)交換的搜索和等待時間。針對速度空間并行分塊策略，采用邏輯網(wǎng)格MPI分組算法[19-20]，劃分不同子通信域，將分塊編號排列如圖3，物理空間分Nx塊，數(shù)據(jù)傳遞在子域(Col_common)中進行，速度空間分Nv塊，在子通信域(Row_common)中歸約，分塊編號矩陣保持了僅對物理空間進行分塊時的大小，速度空間同時進行分塊并沒有對數(shù)據(jù)傳遞中分塊編號的搜索和等待造成更大的負擔。該并行策略能充分考慮到UGKS算法中物理空間和速度空間上的算法差異，并能顯著提高大規(guī)模并行計算的效率。對其進行并行效率測試，其加速比呈線性增長[21]，顯示出了UGKS并行程序的高效并行特性。

圖3 UGKS并行分塊策略Fig.3 Parallel strategy for UGKS

基于本文發(fā)展的大規(guī)模高效并行UGKS程序，對多個典型稀薄流動，如三維方腔流動、圓球繞流、衛(wèi)星體飛行器繞流進行了驗證[21-22]。在此基礎上，對典型返回艙定常繞流問題進行了模擬。模擬真實飛行器復雜外形下的高超聲速稀薄流動，尤其是捕捉多尺度精細流動結(jié)構(gòu)，對格式的穩(wěn)健性和計算效率提出了很高的要求。在此算例中，來流迎角20°，馬赫數(shù)23.5，努森數(shù)0.27?？傆嬎阋?guī)模為140億，其中速度空間離散點數(shù)為2萬。物理空間和速度空間均分塊并行，使用了5040個進程進行計算，總計算時間為25萬CPU小時。圖4和圖5所示為返回艙周圍的溫度和壓力分布云圖。

圖4 返回艙表面及對稱面上的溫度分布Fig.4 Temperature contours of re-entry vehicle

圖5 返回艙壓力分布Fig.5 Pressure contours of re-entry vehicle

3 結(jié) 論

本文介紹了氣體動理學格式GKS在飛行器再入過程中面臨的高超聲速跨流域問題中的拓展和應用。

研究工作主要包括三個方面：1) 基于拓展BGK方程發(fā)展了耦合湍流模式的拓展GKS，格式在典型高超聲速流動模擬中具有良好的表現(xiàn)； 2) 基于CPR框架，結(jié)合緊致型WENO限制器發(fā)展了非結(jié)構(gòu)網(wǎng)格上的高精度動理學格式CPR-GKS，通過黏性激波管等典型算例驗證了其在可壓縮黏性流動中的高精度和較好的激波捕捉性能；3) 發(fā)展了適合復雜分塊結(jié)構(gòu)網(wǎng)格上跨流域稀薄流動大規(guī)模高效并行計算的UGKS方法，初步對典型飛行器的再入進行了精細模擬。

本文研究揭示了GKS在再入飛行器面臨的高超聲速湍流及跨流域稀薄流動模擬中的優(yōu)異性能，具有廣闊的應用前景。

致謝:感謝江南計算技術研究所的徐金秀高級工程師和浙江省水利河口研究院的于普兵博士的幫助。清華探索100高性能計算平臺、并行科技高性能計算支持服務為本文工作提供部分計算機時。