大型航天器再入解體氣動力熱特性模擬的直接模擬蒙特卡洛方法研究

梁 杰，李志輝，2，李緒國，杜波強

(1. 中國空氣動力研究與發(fā)展中心超高速空氣動力研究所， 綿陽621000；2. 北京航空航天大學(xué)北京前沿創(chuàng)新中心國家計算流力學(xué)實驗室，北京100191)

1 引言

服役期滿大型航天器如類天宮飛行器離軌再入大氣層跨越自由分子流、稀薄過渡流和連續(xù)流等多個流動區(qū)域，一般在連續(xù)流區(qū)域是通過求解Euler 或N-S 方程的計算流體力學(xué)(CFD)準確模擬，在稀薄流區(qū)域則采用直接模擬蒙特卡洛(Direct Simulation Monte Carlo，DSMC)方法進行計算。 飛行器作為一個整體在連續(xù)流域飛行時，在高超聲速流場中也可能存在局部的稀薄流動區(qū)域或高度非平衡流區(qū)，如激波、邊界層或尾跡區(qū)。 在這些區(qū)域，傳統(tǒng)的CFD 方法失效，而DSMC 方法在計算高密度區(qū)域的粒子碰撞時花費巨大，因此發(fā)展了CFD-DSMC 混合方法來求解這類多尺度復(fù)雜流動問題[1-2]。 這種混合方法的計算精度和計算效率都比較高，但程序代碼復(fù)雜，2 種不同類型方法之間的信息交換困難，DSMC 方法的統(tǒng)計散布或統(tǒng)計波動極大地影響CFD 計算的穩(wěn)定性。同時也發(fā)展了一些基于粒子模擬方法的混合技術(shù)(如平衡取樣技術(shù)[3-4]、基于DSMC 的低擴散粒子方法[5]、eDSMC[6])來模擬近連續(xù)流，這類混合方法大多建立在DSMC 碰撞限制器的概念基礎(chǔ)之上[7-9]。 在一個限制碰撞的DSMC 模擬中，數(shù)值碰撞頻率被限制在一個較低的值，假設(shè)使用碰撞限制器的流場區(qū)域無粘，擴散輸運的影響可以忽略。 如果將DSMC 碰撞限制器應(yīng)用在流場中的強激波和邊界層區(qū)域，明顯的數(shù)值擴散和耗散會導(dǎo)致模擬結(jié)果有較大的誤差[10]。 當采用DSMC 方法模擬高超聲速近連續(xù)流時，可發(fā)現(xiàn)計算區(qū)域中有許多網(wǎng)格，其中的流動處于或接近于熱力學(xué)平衡狀態(tài)。 網(wǎng)格內(nèi)的分子經(jīng)過有限的碰撞次數(shù)就可以使分子速度達到Maxwell 平衡分布，更多的碰撞并不能改變當?shù)氐姆植己瘮?shù)。 如果在這些區(qū)域采用限制網(wǎng)格內(nèi)每個分子碰撞次數(shù)(碰撞限制器)的方法，并與大時間步長和大網(wǎng)格尺度相結(jié)合，可以節(jié)省大量的計算時間。

因此，針對大型航天器再入解體過程氣動力熱問題模擬，本文提出一種基于碰撞限制器技術(shù)的DSMC 混合方法，根據(jù)高超聲速流動中連續(xù)流假設(shè)失效的判斷準則，將流場區(qū)域分解為遵循氣體分子速度函數(shù)平衡態(tài)分布的連續(xù)流區(qū)和強熱非平衡的稀薄流區(qū)。 在連續(xù)流區(qū)域，采用限制碰撞的DSMC 方法，通過采用網(wǎng)格自適應(yīng)和變時間步長技術(shù)來減少數(shù)值擴散的誤差[11]。 同時采用大規(guī)模并行算法[12]，實現(xiàn)大型復(fù)雜飛行器二次解體過程高超聲速近連續(xù)過渡流氣動力熱特性的模擬，以解決類天宮飛行器首次、二次解體過程中的高超聲速氣動熱問題。 通過計算類天宮飛行器兩艙結(jié)構(gòu)體在低密度風(fēng)洞試驗狀態(tài)的氣動力/熱系數(shù)，并與低密度風(fēng)洞試驗數(shù)據(jù)進行對比，以驗證該算法的可靠實用性。 最后計算分析帶太陽電池帆板的類天宮飛行器激波/激波、激波/邊界層干擾及流動分離的復(fù)雜流場結(jié)構(gòu)，為數(shù)值預(yù)報天宮飛行器首次、二次解體的氣動力熱特性可靠求解提供平臺基礎(chǔ)。

2 數(shù)值模擬方法

2.1 碰撞限制器技術(shù)的DSMC 方法實現(xiàn)

對于大型航天器再入解體過程中的繞流，因其大尺度非規(guī)則，前后不同區(qū)域的局部高超聲速流場中，混合有氣體分子速度分布函數(shù)呈現(xiàn)平衡與非平衡的流動，根據(jù)基于DSMC 碰撞限制器的混合方法，利用連續(xù)流失效參數(shù)來指定采用限制碰撞的DSMC 流場區(qū)域和當?shù)刈赃m應(yīng)DSMC 計算區(qū)域。 連續(xù)流失效參數(shù)根據(jù)Boyd 等的研究[13]，采用式(1)所示基于當?shù)亓鲃犹荻乳L度(GLL)的克努森數(shù)(Kn)來判斷:

式中，梯度取流動參數(shù)的最大梯度方向， Q是某些量如密度、壓力、溫度或速度。 本文采用基于密度的KnGLL來進行區(qū)域分解，KnGLL＞0.02 的區(qū)域定義為非平衡區(qū)域，其它區(qū)域則為平衡區(qū)域。采用碰撞限制器的DSMC 方法計算。 由于整個流場都統(tǒng)一采用基于DSMC 的方法進行計算，因此在穿越連續(xù)流失效邊界時并不存在額外的信息交換。

對于精確模擬，標準的DSMC 方法要求足夠小的計算網(wǎng)格、足夠小的時間步長和每個網(wǎng)格內(nèi)足夠多的模擬粒子數(shù)。 其網(wǎng)格尺寸應(yīng)該小于當?shù)胤肿悠骄杂沙?，時間步長間隔應(yīng)小于當?shù)胤肿悠骄鲎矔r間，才能得到有效的結(jié)果。 但在高超聲速流場的平衡區(qū)域，由于穿過每個網(wǎng)格的流動都為常量，所以時間步長和網(wǎng)格尺寸可以大許多倍。 網(wǎng)格中的每個粒子經(jīng)過一定次數(shù)的碰撞后，當?shù)氐乃俣确植己瘮?shù)就能接近Maxwell 平衡態(tài)分布。 碰撞限制器就是把網(wǎng)格中每個粒子在一個時間步長之內(nèi)的碰撞次數(shù)限制為一個常數(shù)值χ，即在一個時間步長內(nèi)把一個網(wǎng)格內(nèi)粒子總的碰撞次數(shù)限制為χ·N 次，其中N 是這個網(wǎng)格內(nèi)模擬的粒子數(shù)。 不少文獻對碰撞限制器DSMC 中χ 值的影響進行了檢驗[14-16]，本文算例中選取χ ＝2。與在連續(xù)流平衡區(qū)域采用真實的高碰撞頻率相比，每個時間步長內(nèi)每個網(wǎng)格總的粒子碰撞次數(shù)限制為2N 次，可以減少粒子碰撞處理過程的計算時間。

2.2 計算網(wǎng)格及網(wǎng)格自適應(yīng)

本文采用直角與表面非結(jié)構(gòu)混合網(wǎng)格的DSMC 數(shù)值方法。 在描述物面幾何形狀的非結(jié)構(gòu)網(wǎng)格建立以后，直接將其嵌入直角網(wǎng)格的流場中，使DSMC 計算對流場網(wǎng)格的依賴程度大大降低。同時通過判斷分子運動軌跡方程和物面三角形面元上任一點的位置方程，唯一確定出分子與物面的碰撞點坐標[17]，解決了這種混合網(wǎng)格流場分子運動與物面碰撞的難題。 對分子在物體三角形面元上碰撞、反射前后的流場參數(shù)進行統(tǒng)計取樣，就可以獲得飛行器的整體氣動力特性以及表面力、熱載荷分布。

在流動梯度較大的區(qū)域，碰撞分子應(yīng)該通過有效的最相鄰搜索方法選取。 本文在高度非平衡區(qū)采用基于密度梯度的動態(tài)自適應(yīng)碰撞網(wǎng)格和取樣網(wǎng)格，避免網(wǎng)格內(nèi)隨機選取的碰撞分子之間的距離遠大于當?shù)氐姆肿悠骄杂沙獭?即在背景網(wǎng)格的基礎(chǔ)上，根據(jù)流場中密度梯度的變化分別對碰撞網(wǎng)格和取樣網(wǎng)格進行細化，碰撞分子則是在自適應(yīng)后最小的亞網(wǎng)格內(nèi)選取，以保證計算的空間精度。 但是為保證高密度區(qū)域精細的空間分辨率而增加總的模擬分子數(shù)，會導(dǎo)致低密度區(qū)域的粒子過度增多，嚴重降低計算效率。 為了使計算區(qū)域每個網(wǎng)格更均勻地分布模擬粒子，解決方法之一就是采用變時間步長格式。

2.3 變時間步長格式

當計算區(qū)域中流場密度有較大的變化時(如強烈的激波壓縮或激波干擾流動)，按照DSMC方法的模擬要求，時間步長要小于當?shù)氐姆肿悠骄鲎矔r間，如果整個計算區(qū)域統(tǒng)一采用較小的時間步長，那么整個流場要達到穩(wěn)定狀態(tài)需要花費非常多的計算時間，計算效率很低。 如果流場中每個模擬粒子的權(quán)重相同，則在高密度區(qū)域模擬的粒子數(shù)非常少，而在低密度區(qū)域粒子數(shù)又大大超過求解的需要[18]。 過多的時間浪費在低密度區(qū)域的計算上，而在高密度區(qū)域取樣嚴重不足。由于每個網(wǎng)格內(nèi)的粒子數(shù)與氣體密度呈反比，為了使計算區(qū)域中每個網(wǎng)格的模擬分子數(shù)分布更加均勻，并且不會被分子“復(fù)制”產(chǎn)生有害的影響，必須將當?shù)氐臅r間步長和分子的權(quán)重相匹配。 為了保證模擬粒子穿越網(wǎng)格交界面時通量守恒(包括質(zhì)量、動量和動力學(xué)能量守恒)，Wi/Δti在整個計算區(qū)域都保持相同。 模擬分子穿過網(wǎng)格交界面時的剩余時間，需要按照分子運動的初始網(wǎng)格和目的網(wǎng)格的時間步長的比值重新調(diào)整。 經(jīng)過時間步長的自適應(yīng)過程后，網(wǎng)格i 的時間步長可用式(2)計算:

式中，Δt∞和n∞是時間步長和數(shù)密度的參考值， ·[ ] 表示截斷。 如果某個網(wǎng)格的數(shù)密度大于參考值，則其時間步長就等于參考時間步長。 網(wǎng)格內(nèi)的模擬粒子數(shù)僅根據(jù)其自身的權(quán)重自適應(yīng)，而粒子權(quán)重又與時間步長成正比，因此在高密度區(qū)域粒子數(shù)會增加，密度低的區(qū)域粒子數(shù)則減少。這就有可能使自適應(yīng)后的碰撞網(wǎng)格尺度接近當?shù)氐姆肿悠骄杂沙獭?/p>

3 計算結(jié)果分析

3.1 類天宮飛行器風(fēng)洞試驗狀態(tài)氣動力特性模擬驗證

類天宮飛行器兩艙結(jié)構(gòu)模型的氣動力試驗是在中國空氣動力研究與發(fā)展中心的高超聲速低密度風(fēng)洞完成的，試驗來流氣體是M12 型面噴管噴出的氮氣，采用2%的模型縮比，模型總長為0.202 89 m，模型參考面積為3.525 66×10-3m2，俯仰力矩的參考點位于頭部頂點。 表1 給出了試驗狀態(tài)參數(shù)，試驗的迎角范圍為-2 ° ～25 °。DSMC 計算中采用VHS 分子模型、L-B 內(nèi)能交換模型，壁面采用固定壁溫300 K 的完全漫反射模型。

表1 類天宮飛行器氣動力試驗狀態(tài)參數(shù)Table 1 Test conditions for Tiangong-type vehicle

雖然試驗狀態(tài)的克努森數(shù)Kn 非常低，按傳統(tǒng)流區(qū)劃分，其模擬高度處于連續(xù)流區(qū)域，但是采用本文構(gòu)建的碰撞限制器技術(shù)和當?shù)貐?shù)自適應(yīng)的混合DSMC 方法進行計算，既能保持計算的高精度又能大幅度提高計算效率。 圖1 分別給出了計算的軸向力、法向力和對頭部頂點的俯仰力矩系數(shù)與試驗結(jié)果的對比情況，可看出計算與試驗得到的氣動力系數(shù)隨迎角的變化規(guī)律完全一致。本文DSMC 模擬得到的軸向力系數(shù)與試驗結(jié)果在0°迎角的最大偏差為5%，法向力系數(shù)在20°迎角的最大偏差7%左右，俯仰力矩系數(shù)在20°迎角的最大偏差在5%附近。 對于相對細長的類天宮飛行器模型，計算的氣動力系數(shù)與低密度風(fēng)洞試驗數(shù)據(jù)的最大偏差均在7%以內(nèi)，表明本文算法可以初步實現(xiàn)大型復(fù)雜航天器跨流域氣動力特性的模擬能力。

圖1 類天宮飛行器氣動力系數(shù)DSMC 計算與試驗比較Fig.1 Comparison of Tiangong-type aerodynamic coefficients between DSMC results and experimental data

圖2～3 分別給出了0°和20°迎角下試驗狀態(tài)流場的密度、壓力和馬赫數(shù)等值線分布云圖，可以看出復(fù)雜繞流現(xiàn)象及頭部區(qū)域的激波/邊界層流動干擾特征，表明基于當?shù)亓鲃訁?shù)自適應(yīng)的DSMC 方法對于捕捉復(fù)雜的流動干擾細節(jié)是有效的。

圖2 類天宮飛行器0°迎角流場等值線分布Fig.2 Flowfield contours distribution of Tiangongtype vehicle at 0 degree angle of attack

3.2 帶太陽電池帆板的類天宮飛行器氣動力/熱復(fù)雜流動模擬分析

圖3 類天宮飛行器20°迎角流場等值線分布Fig.3 Flowfield contours distribution of Tiangongtype vehicle at 20 degree angle of attack

為了研究壽命末期的天宮1 號目標飛行器無控再入過程中的熔融燒蝕和解體情況，需要準確模擬跨流域不同高度的氣動力、氣動熱詳細分布及流場干擾特征。 圖4 繪出了帶太陽電池帆板的類天宮飛行器在克努森數(shù)Kn ＝6.7×10-4(特征長度10.4 m)、0°迎角下的壓力等值線分布，計算的來流條件是:速度7500 m/s、高度85 km。 圖中可以清晰地看到復(fù)雜激波干擾下的流動特征:頭部簡化的對接機構(gòu)產(chǎn)生的脫體激波，在錐段壓縮激波以及與實驗艙連接處膨脹波的干擾下向外延伸，并與太陽電池帆板前方的脫體激波相交，產(chǎn)生的反射激波與帆板邊界層相互干擾，使帆板局部表面壓力急劇上升。 另有多道反射激波從帆板脫體激波內(nèi)部與之相交，使帆板脫體激波形成波浪形狀。 在艙體表面附近區(qū)域，由于氣體的膨脹，壓力迅速下降，因此在頭部脫體激波、艙體和太陽帆板表面之間形成較大的逆壓梯度，引起較強的流動分離，產(chǎn)生一個較大的分離主渦和兩個較小的分離附著渦，如圖5 所示。 還可看出，從反射激波作用在帆板表面邊界層的位置開始，沿帆板左右兩側(cè)各形成一股射流，朝向帆板外側(cè)的射流在帆板表面附近僅引起微弱的流動分離(參見圖5)；而朝向艙體的射流在較大的壓力梯度下加速至超音速并沖向資源艙的表面，形成局部的脫體激波并在資源艙表面引起局部高達120 kW/m2左右的熱流分布，如圖6 所示。 太陽電池帆板與艙體間構(gòu)成的這種特殊非規(guī)則繞流沖擊形成的強氣動熱，導(dǎo)致的結(jié)構(gòu)響應(yīng)變形最終撕裂了太陽電池帆板，使之發(fā)生首次解體。

圖4 帶帆板的類天宮飛行器0°迎角壓力等值線分布Fig.4 Pressure contours of Tiangong-type vehicle with solar panels at 0 degree angle of attack

圖5 帶帆板的類天宮飛行器0°迎角流線分布Fig.5 Streamline distribution of Tiangong-type vehicle with solar panels at 0 degree angle of attack

圖6 激波干擾引起的局部表面高熱流Fig.6 High heat flux on local surface induced by shock wave interactions

4 結(jié)論

1)本文針對大尺度復(fù)雜結(jié)構(gòu)航天器再入繞流不同位置呈現(xiàn)多流區(qū)局部流場特點，構(gòu)建了基于碰撞限制器技術(shù)的全粒子混合方法，并進一步發(fā)展了基于密度梯度的動態(tài)自適應(yīng)混合網(wǎng)格處理技術(shù)與變時間步長方案。 根據(jù)連續(xù)流失效參數(shù)，將流場區(qū)域分解為平衡的連續(xù)流區(qū)域和非平衡的稀薄流區(qū)域，在平衡區(qū)域采用大尺度網(wǎng)格和大時間步長的碰撞限制器DSMC 方法，而在非平衡區(qū)域則采用網(wǎng)格自適應(yīng)和變時間步長的DSMC 方法，在提高計算精度的同時，節(jié)省了大量計算時間；

2)使用本文算法對類天宮飛行器兩艙體風(fēng)洞試驗狀態(tài)氣動力特性進行不同迎角的氣動特性計算，解得的氣動力/力矩系數(shù)與低密度風(fēng)洞試驗數(shù)據(jù)對比一致性較好，偏差7%以內(nèi)，驗證了本文提出的碰撞限制器和當?shù)貐?shù)自適應(yīng)的全粒子混合算法對大型復(fù)雜航天器再入氣動特性高精度模擬能力；

3)模擬了帶太陽電池帆板的類天宮飛行器85 km 高度再入時的高超聲速復(fù)雜氣動熱及流場激波干擾結(jié)構(gòu)，揭示了頭部脫體激波、艙體和太陽電池帆板非規(guī)則繞流激波引起流動分離、產(chǎn)生強氣動熱和結(jié)構(gòu)響應(yīng)變形導(dǎo)致太陽電池帆板毀壞發(fā)生首次解體的機制。