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

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

0 引 言

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

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

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

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

1 GKS簡介

氣體動(dòng)理學(xué)格式基于Boltzmann方程的模型方程，BGK方程。它具有如下形式：

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

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

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

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

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

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

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

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

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

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

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

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

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

2 GKS的拓展及應(yīng)用

2.1 高超聲速湍流模擬

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

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

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

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

2.2 基于CPR框架的高精度GKS

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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