大型航天器離軌再入氣動融合結構變形失效解體落區數值預報與應用

李志輝，彭傲平，馬 強，石衛波，黨雷寧，梁 杰，蔣新宇，唐小偉

（1.中國空氣動力研究與發展中心超高速空氣動力研究所，綿陽621000；2.國家計算流力學實驗室，北京100191；3.四川大學數學學院，成都610041）

1 引言

近地軌道高度300～1000 km運行的板艙桁架結構大型航天器（空間實驗室、貨運飛船、空間站等），服役期滿或失效時，如已處于失聯無控飛行狀態的歐洲最大對地觀測衛星ENVISAT，將面臨發生可能危及其他航天航空器的碰撞、再入墜毀航跡預示、人財物安全評估、危害性分析與處置等問題［1-2］。正常情況下，這類航天器在進行返回制動后，控制系統就不再介入返回過程，在過渡段、再入段都不進行控制，與之前的再入航天器（如載人飛船返回艙）有很大區別，稱無控隕落航天器，從離軌到再入過程是一個軌道高度和能量逐漸降低、連續變化過程［3］。這個過程可以劃分為2個階段：軌道衰降階段和再入損毀階段。在軌道衰降階段，航天器仍然能夠環繞地球、以螺旋形橢圓軌道飛行，并在稀薄空氣動力、地球重力、磁場力等外力持續作用以及太陽活動的間歇性影響下，航天器的飛行軌道高度逐漸降低，軌道衰降速度越來越快；當軌道高度低于120 km時，進入再入損毀階段［4］，如圖1所示，航天器在大氣劇烈摩擦氣動力/熱作用下，摩擦生熱使航天器金屬（合金）桁架結構變形失效熔融，在氣動熱、力、減速過載等綜合作用下，發生劇烈破壞解體，面臨再入墜毀處置問題。

圖1 大型航天器離軌再入氣動/結構響應熔融解體過程示意Fig.1 Sketch of disintegration process of aerothermodynam ics/structural responsemelting for large-scale spacecraft during deorbit reentry

圍繞這類航天器壽命末期離軌飛行軌道衰降與隕落再入飛行航跡預報試驗［5-9］，開展離軌再入空氣動力學計算模型、空間環境、強氣動力熱環境致金屬（合金）桁架結構響應變形/軟化/熔融/嵌層復合材料熱解燒蝕/解體相關研究，對建設與檢驗國家航天監測、系統跟蹤、定軌預報、風險評估具有極重要意義。如中國天宮一號目標飛行器超期服役兩年半后突發功能失效，無控飛行軌道衰降與再入過程中如何進行飛行航跡落區數值預報，直接依賴如何準確可靠求解再入各流域復雜多尺度非平衡條件下空氣流動與傳熱問題；如何開展航天器金屬（合金）桁架結構在超高速再入強氣動力熱環境變形軟化、熔融毀壞非線性力學行為可計算建模，一直是流體力學、材料結構力學、計算數學與應用物理研究工作者面臨的挑戰［10-15］。

航天器軌道衰降與再入飛行過程，周圍流場出現強烈的間斷粒子稀薄非平衡流動效應，從航天器離軌飛行軌道衰降到再入解體飛行最后落區，整個隕落解體過程將從外層空間自由分子流過渡至近地面非規則解體物連續介質繞流。雖然連續流體力學和稀薄氣體動力學均各有研究方法［16-18］，但近連續滑移過渡區真實氣體高超聲速非平衡流動介于連續流與稀薄氣體自由分子流之間，無論實驗模擬還是理論計算均很難處理，目前尚沒有統一的、準確可靠的模擬方法，急待發展相關基礎理論、計算與實驗模擬方法。為此，中國國家杰出青年基金項目、國家重點基礎研究發展計劃（973計劃）項目“航天飛行器跨流域空氣動力學與飛行控制關鍵基礎問題研究（2014CB744100）”頂層設計制定了核心前沿基礎研究方向“服役期滿大型航天器離軌再入跨流域氣動/結構響應熔融/燒蝕/解體飛行航跡預報模擬”，確立了項目核心子課題群圍繞自制類天舟一號貨運飛船再入大氣層過程，建立跨流域氣動力熱一體化模擬融合彈（軌）道/結構響應計算分析方法，初步形成大型航天器跨流域氣動力熱/變形失效/熱解燒蝕/解體飛行航跡數值預報平臺框架，對接實施載人航天工程辦公室針對天宮一號目標飛行器失效無控飛行再入大氣層應對任務“天宮一號目標飛行器無控隕落預報及危害性分析”課題實施，分別開展服役期滿大型航天器無控隕落與受控再入解體模擬研究。

本文在過去20年開展Boltzmann方程碰撞積分物理分析與可計算建模氣體動理論統一算法（GKUA）［19-24］研究基礎上，采用轉動慣量描述氣體分子自旋運動，引入能量模式配分函數和非彈性碰撞松弛數，提出求解含能量輸運非平衡效應Boltzmann模型方程統一算法及其與結構熱力響應有限元算法耦合模擬技術，建立適于大型航天器服役期滿離軌再入跨流域氣動環境與結構熱力耦合響應變形/失效解體非線性力學行為一體化計算應用平臺。在分析驗證GKUA與DSMC、連續/近連續流區（滑移）N-S解算器、N-S/DSMC耦合算法、低密度風洞實驗模擬手段途徑正確性基礎上，將再入氣動融合結構變形失效解體算法分別應用于類天舟一號貨運飛船、天宮一號/二號空間實驗室再入解體數值預報，檢驗該算法的正確性。

2 求解含非平衡效應Boltzmann模型方程統一算法與大規模并行計算技術

基于轉動松弛特性Rykov模型研究思想［25］，在氣體分子速度分布函數演化更新求解中考慮內自由度影響，采用轉動慣量描述氣體分子自旋運動，利用分子總角動量守恒作為一個新的碰撞不變量。把分子內部能量作為分布函數自變量，引入能量模式配分函數將能量在各自由度平均分布。基于氣體分子速度分布函數f=f（r，V，t，e）（r、V分別是位置空間和速度空間的坐標，e為內能），在求解Boltzmann模型方程統一算法框架［19-22］下，基于權因子1和內能e對速度分布函數所依賴的速度空間進行無窮積分，引入能級約化速度分布函數f0與f1，去掉方程對內部能量連續依賴關系，確立含非平衡效應各流域統一Boltzmann模型方程，其無量綱形式為式（1）：

其中，上下標t、r分別表示彈性碰撞和非彈性碰撞，如νt和νr是彈性碰撞頻率和非彈性碰撞頻率；下標i表示空間維度方向，μ為粘性系數，n、T、p、q分別是氣體分子數密度、氣體流動溫度、壓力、熱流，ω0、ω1是與氣體特性相關的常數，δ對于分子間相互作用規律來說是一個常數，c是氣體分子熱運動速度，n是氣體分子數密度。反映粘性與熱傳導屬性的擴散時間與氣體分子平均碰撞時間可分別用于描述含內能激發的多原子分子彈性與非彈性碰撞松弛變化，其統一表達式如式（2）～（3），以表征氣體分子速度分布函數趨于當地平衡態分布碰撞松弛速率。

其中，B是與氣體相關的常數。氣體動力學各宏觀流動量，如密度、流動速度、溫度、壓力、應力張量、熱流矢量，均可通過式（4）所示氣體分子速度分布函數f0和f1對速度空間積分得到。

式中，u為流動的宏觀速度，其大小為u；τij為應力張量分量，δij表示i=j時其值為1、不相等時為0，i、j表示空間維度方向。

速度分布函數f0、f1對速度分量函數依賴關系屬指數型。需將速度空間離散降維，去掉分布函數對速度空間的連續依賴性。在離散數值求解過程中，通過構造自適應選取的離散速度坐標點來保證離散速度坐標點確定的速度分布函數值正定守恒。為此，經過速度空間離散化處理的速度分布函數方程在計算平面內矩陣形式如式（5）：

采用二階有限差分格式［26］，將時間分裂數值計算方法應用于方程（5）式中位置空間的三個對流運動方程和碰撞松弛源項方程，進行耦合迭代數值求解，構造在每個離散速度坐標點處直接求解分子速度分布函數演化更新的氣體動理論統一數值格式如式（6）：

式中，L為對應方向上的算符，n表示n時刻的值，下標S表示源項。對于（6）式，采用二階Runge-Kutta方法數值求解具有非線性源項的碰撞松弛方程；采用時間和位置空間均為二階精度基于原始變量（即相對于f0σ，δ，θ、f1σ，δ，θ）的NND-4（a）預測-校正兩步差分格式［26］數值離散對流運動方程。

針對顯式格式差分求解離散速度分布函數，計算時間步長受顯式格式穩定性條件制約［27］，為提高計算大型航天器再入解體過程非規則物形繞流計算效率，發展了Boltzmann模型方程隱式格式求解方法［28］。為方便起見，以一維Boltzmann模型方程為例構造隱式格式，簡單記為式（7）：

式中，A=?F/?U相應于Boltzmann模型方程對流輸運項的Jacobi矩陣，M=?S/?U相應于Boltzmann模型方程碰撞松弛源項線性化矩陣。為此，求解過程可分解為：①令（8a）右端第二項為零或取前一次迭代值，利用左邊的邊界條件由左向右掃描，依次求出；②將上一步求出的代入（8b）右端，利用右邊的邊界條件，從右向左進行掃描，求出；③以上左右各掃描依次完成一次迭代，將每次迭代求出的新值代入（8a）式右端并重復掃描，直到小于給定的某個值為止。一般只需進行兩次迭代就可滿足精度要求。由（8）式可看出前的系數為標量形式，計算過程不需塊矩陣求逆運算，計算實現較簡便。

根據數值格式（8），建立了考慮非平衡效應、基于分子速度分布函數特征化表征的基礎邊界條件的數學模型及數值處理方法［25，30］，一旦所有離散速度坐標點處的速度分布函數被數值求解，則任一時刻物理空間各點的宏觀流動量如氣體密度、流動速度、溫度、應力張量和熱流矢量等，便可通過速度分布函數乘以分子速度的某種函數形式再對離散速度空間數值積分而被演化更新。

適于航天再入氣動力熱問題的Boltzmann方程可計算建模氣體動理論統一算法，計算空間是由離散速度空間和位置空間組成的六維空間，并考慮內部能級分布，形成多相空間。從區域分解并行化原理出發，可將算法求解空間劃分為位置空間分解策略、速度空間分解策略。通過對各階段算法過程進行變量依賴關系分析看出，在各個離散速度坐標點 (Vxσ，Vyδ，Vzθ)處數值求解離散速度分布函數的差分格式中，位置空間Ωr各維方向都存在數據相關性，離散速度子空間ΩV各維方向則毫無數據相關。使用離散速度數值積分法對氣體分子離散速度分布函數進行矩積分確定飛行器繞流宏觀流動量的過程，在位置空間Ωr各維方向都不存在數據相關性，在離散速度空間ΩV各維方向以歸約形式體現了問題的相關性。因此，對于ΩV分解策略，計算宏觀流動量需要在ΩV內進行并行歸約計算，會產生數據通信。依據二叉樹方式并行歸約，可得單進程處理的數據量CV如式（9）：

式中，假設位置空間Ωr是NiNjNk網格陣列，處理器陣列為NpσNpδNpθ。通常可將 ΩV按三維方式等分分解：為了分析離散速度空間的網格陣列與處理器陣列的分布關系，將式（9）改寫為式（10）：

式中，V=NσNδNθNiNjNk，Nσ、Nδ、Nθ分別為3個計算方向離散速度點數。從式（10）看出，為獲得較小CV值，通常設置Npσ≤Nσ、Npδ≤Nδ、Npθ≤Nθ為宜。因此，基于離散速度空間ΩV并行分解策略能實現具有成千上萬CPU核甚至更大規模并行計算。將求解Boltzmann模型方程并行算法計算效率與賓西尼亞大學依靠高性能計算機開展求解BGK模型方程不同馬赫數一維激波結構演化問題的研究［31］進行比較，表1給出了本文基于256個CPU的GKUA并行算法用于計算三維繞流問題的并行加速比與文獻［20］BGK方程計算一維激波的并行加速比。可見，GKUA三維繞流并行計算加速比與文獻［31］一維BGK方程并行計算加速比相當，顯示出統一算法卓越的并行加速性能，是國際首創的求解各流域三維復雜繞流問題的氣體動理論高性能并行算法，且具有良好的并行可擴展性與并行效率。

表1 基于256CPU統一算法并行計算加速比與國際同行比較Table 1 Comparison of speedup ratio of parallel computation for GKUA on 256CPU w ith sim ilar international studies

圖2繪出本文GKUA程序經異構協同模式并行移植優化得到500～45000處理器核的并行計算加速比與并行效率，圖中橫坐標表示計算所用進程數。可看出，本文設計的適于眾核異構計算機的氣體動理論高性能并行算法加速比，仍屬擬線性加速比，并行效率80%以上，良好的并行加速性能足以保證在較高并行效率下，通過增加處理機數，即可大大增加計算規模，使航天領域依靠傳統計算條件難以解決從高稀薄自由分子流到連續流全飛行流域復雜高超聲速繞流問題計算求解變為現實成為可能。

3 建立結構動態熱力耦合響應變形行為有限元計算方法

圖2 統一算法在500～45000多核處理器進程并行計算加速比與并行效率Fig.2 Speedup ratio and parallel efficiency on 500～45000 processors for gas-kinetic parallel algorithm

直接求解Boltzmann模型方程統一算法基于氣體分子速度分布函數演化更新，使用（4）式數值積分得到所有宏觀流動物理量的做法，能準確捕捉物面不同位置的力、熱流、溫度、速度分布，這為進一步研究飛行器在極高速再入強氣動力熱環境下，其部件結構內部熱傳導致溫度梯度與彈塑性位移變形行為創造了條件。為研究航天器再入氣動力熱環境致結構變形/失效解體非線性力學響應行為，基于瞬態熱傳導方程與材料熱彈性動力學方程，需構造動態熱力耦合響應數學模型，建立再入強氣動力熱環境致材料結構熱力響應行為有限元算法。

設在空間區域Ω?R3中，不考慮阻尼影響，材料熱彈性動力學方程可表示為式（11）［32］：

式中，c表示材料比熱，kij為材料熱傳導張量，h為熱源項。利用變分原理推導［13，33-34］關于u與θ的動態熱力耦合控制方程的弱形式如式（13）：

該耦合弱形式屬于雙曲拋物耦合系統，其溫度場與位移場相互依賴，在有限元計算過程中兩方程需要聯立求解。

對于熱彈性動力學方程，使用Newmark隱式方法實施熱彈性動力學方程時間上的離散推進，計算速度與位移［13，34］；使用高精度Crank-Nicolson格式［13，35-36］耦合求解熱傳導方程。

4 高超聲速再入氣動環境致結構變形/失效非線性力學行為一體化模擬

為了準確可靠模擬航天器再入中高超聲速氣動力/熱環境致金屬（合金）桁架結構材料內部溫度分布與響應變形、熱損毀等非線性力學行為，需要將動態熱力耦合響應有限元計算方法（FEM）與求解Boltzmann模型方程氣體動理論統一算法（GKUA）相結合，發展一體化耦合模擬技術。為簡單起見，這里以圖3所示再入豎直平板含碳鋼結構繞流模型為例，把任一時刻GKUA得到的氣動熱與氣動力代入結構熱力耦合響應計算公式，作為有限元計算的邊界條件。對圖3所示鋼板厚度0.015 m、高度0.5 m，鋼板外部受再入近連續稀薄過渡流區高超聲速繞流沖擊，來流克努森數Kn∞=0.01、馬赫數Ma∞=8.366 6、γ=1.4。使用統一算法求解外流場，因高超聲速再入近連續過渡流，在豎直鋼板附近繞流場形成較厚的脫體激波層（如圖3（a）），鋼板表面受到強氣動熱與氣動力作用，在迎風區承受高溫與高壓，若不考慮鋼板振動，根據GKUA求解獲得鋼板表面溫度與壓力作為材料結構動態熱力耦合響應有限元算法界面條件，對鋼板內部溫度分布與變形情況進行計算分析。

圖3 高超聲速再入Kn∞=0.01、Ma∞=8.366 6豎直鋼板繞流場溫度等值線GKUA與內部溫度分布FEM計算驗證Fig.3 Tem perature contours and computation verification of GKUA and FEM along x axis between flow field and inner structure under hypersonic reentry when Kn∞=0.01，Ma∞=8.366 6

對鋼板截面內部結構進行網格剖分，固定鋼板下表面（y=0）。由于用于繞流物體外部流場計算的網格與用于熱力耦合響應計算的鋼板內部結構網格的不一致性，這里將GKUA計算得到的外部流場網格邊界的溫度與壓力進行線性插值施加到鋼板內部網格邊界中。假設鋼板處于初始溫度T0，外部流場處于穩定狀態，鋼板外表面溫度達到T=T0+θ，其中T為GKUA計算豎直鋼板外流場穩定時獲得的溫度值，由此得到鋼板物面邊界的溫度增量值θ=T-T0。

在外部流場計算中，為表征飛行器再入不同高度的跨流域繞流多尺度效應，通常采用無量綱模型化計算；而在物體內部結構場的熱力耦合有限元計算中，則需將其轉換成符合實際的有量綱溫度與壓力。如在流場中計算得到的溫度與壓力值分別為T*與p*，則實際溫度與壓力T與P分別為T=T*T∞、P=p*ρ∞V2∞。其中，T∞、ρ∞、V∞分別為參考溫度、來流密度與來流速度。作為算例，鋼板物性參數為E=2.06×1011Pa，ν=0.3，α=1.5×10-6/℃，ρ=7 850 kg/m3，cv=460 J/（kg·℃），k=50W/（m·℃）。初始參考溫度T0=500℃，取時間步長Δt=0.05 s，計算時間為5 s。圖3（b）繪出直接耦合法計算得到鋼板內部結構場溫度增量分布等值線云圖，圖3（c）繪出鋼板不同部位y=0.06、0.26、0.46外部流場與內部溫度場沿x方向變化曲線。可看出，鋼板外部使用GKUA計算得到的迎風面y=0.06、0.26、0.46處溫度與背風面溫度，均分別與鋼板內部熱力耦合有限元計算得到的鋼板表面溫度吻合一致，證實外流場求解Boltzmann模型方程統一算法與材料內部動態熱力響應耦合算法模型與程序設計實現的正確可靠性；物體內部靠近鋼板迎風面駐點（坐標原點）的溫度較附近其他點的溫度高，最高的溫度點則出現在鋼板的左上拐點，此處邊界曲率變化最大，流場熱流在此達到最大，由穩態熱傳導方程極大值原理可知，材料內部溫度在此也達到最大，故在鋼板迎風邊界由下到上溫度出現一個先降后升的變化，這種變化在材料內部也清楚地反映出來。

圖4繪出鋼板在高超聲速再入繞流環境受外部溫度與壓力雙重作用下t=1～5 s整體變形情況，可以看出，在受外部高超聲速強氣動力熱繞流迎風面與背風面較大溫度與壓力差雙重作用下，鋼板首先產生了較大幅度的彎曲（圖4（a）），這是迎風面壓力作用的結果，隨后氣動熱使得鋼板內部溫度逐漸升高，在鋼板內部形成了非均勻熱應力分布，高溫使得鋼板在迎風面產生較背風面更大的熱膨脹，所產生的熱應力方向與迎風面壓力方向相反，內外兩種力相互作用使得鋼板彎曲幅度變得相比于計算初始階段的彎曲幅度更小，在5 s時，鋼板內部溫度分布還未穩定，所以熱膨脹仍在持續，即5 s時鋼板的彎曲變形還將繼續，最終達到圖4（d）的穩態變形分布。

圖4 強氣動力熱環境致豎直平板不同時刻及穩態變形圖（變形放大尺度為50倍）Fig.4 Time history of plate deformation computed by FEM，displacement am plified by 50 times

5 航天器隕落再入過程跨流域氣動環境與結構熱力響應變形計算分析

5.1 求解Boltzmann模型方程統一算法對大型航天器跨流域氣動并行計算檢驗

為確認求解Boltzmann模型方程統一算法對大型航天器隕落再入過程跨流域氣動并行計算的精度和可靠性，對長度Lref=10.409 m、最大直徑Dref=3.452 m的類天宮飛行器兩艙體近連續過渡流區兩個繞流狀態H=65.2 km、Ma∞=12.5、Kn∞=5.1×10-5（case1）與H=62.1 km、Ma∞=12.79、Kn∞=3.37×10-5（case2），在流動介質氮氣（N2），迎角 α=2°、0o、2°、5o、10o、15o、20o、25o，取Pr=0.72、γ=1.4、Tw=300 K，分別開展GKUA與DSMC［37］、N-S/DSMC［17］、N-S方程計算和低密度風洞實驗測試結果間比較分析。圖5繪出GKUA采用4096CPU核并行計算得到的兩艙體隕落再入case1繞流場壓力、流線結構、俯仰力矩、壓心位置與DSMC、N-S/DSMC、滑移N-S（NS-Slip）、無滑移N-S（NS-NoSlip）、低密度風洞實驗（Exp.）結果比較曲線，可看出，對縱橫十數米大尺度復雜結構航天器，即使在較高飛行高度，但因極低的來流克努森數Kn∞=5.1×10-5，流場整體上呈現明晰的脫體激波、膨脹真空區與尾部流動分離的連續流區繞流面貌，因前體對接臺復雜結構，氣體繞流出現高壓、高溫區，在跨越前臺區進入背風圓柱面及尾部呈現較低密度、壓力分布；且前體對接臺、后體非規則處，稀薄氣體速度滑移效應表現嚴重，有迎角繞流出現迎風區與背風區物面流線匯聚，25°迎角繞流前后艙體上部均出現流線分離、非對稱渦結構，圖5（a）、（b）分別揭示了大型復雜結構航天器繞流特殊的連續流、近連續滑移、稀薄效應共存混合流動輸運現象及變化特點。圖5（c）、（d）分別繪出類天宮兩艙體隕落再入case1繞流狀態不同迎角俯仰力矩和壓心位置驗證比較情況，可看出，對case1（模擬高度H=65.2 km）的俯仰力矩，N-S方程及滑移N-S計算與其他方法計算結果均相差較大，說明該狀態下稀薄效應較為明顯，N-S方程已不再適用，而統一算法（GKUA）計算結果始終介于N-S/DSMC、DSMC和N-S方程所得結果之間，且更吻合于DSMC結果，誤差在合理范圍內，同時統一算法所得壓心位置較DSMC、N-S/DSMC與實驗值符合更好，表明統一算法結果更為準確可靠。

圖5 類天宮兩艙體H=65.2 km、Kn∞=5.1×10-5、Ma∞=12.5、α=25°流場及氣動力系數計算與實驗比較驗證Fig.5 Validation of computation and experiment for aerodynam ics on two-capsule vehicle of Tiangong type spacecraft w ith H=65.2 km，Kn∞=5.1×10-5，Ma∞=12.5 andα=25°

圖6繪出case2繞流狀態計算與實驗比較情況，由于飛行高度H=62.1 km更低，屬于來流Knudsen數Kn∞=3.37×10-5更小的連續流區，整體而言N-S方程計算所得結果與實驗值符合較好，其他幾種方法均有所偏離實驗值，但誤差均在5%以內。同時可看出統一算法結果與N-S方程結果吻合，表明統一算法對連續流區計算具備良好的適應性，計算結果準確可靠。幾種空氣動力學模擬手段對大型復雜結構航天器case1、case2繞流狀態模擬的分析比較表明，統一算法對更高稀薄效應H=65.2 km計算結果與DSMC、N-S/DSMC及實驗吻合更好，對更低飛行高度H=62.1 km計算與N-S方程及實驗符合較好，證實從H=65.2 km到H=62.1 km這樣的飛行高度段，大型航天器繞流流態從稀薄流到連續流發生明顯變化，而統一算法在跨流域氣動計算所具有的一致收斂性優勢。

圖6 類天宮兩艙體H=62.1 km、Kn∞=3.37×10-5、Ma∞=12.79氣動力系數計算與實驗比較驗證Fig.6 Validation of com putation and experiment for aerodynam ics on two-capsule vehicle of Tiangong type spacecraft w ith H=62.1 km，Kn∞=3.37×10-5 and Ma∞=12.79

5.2 類天宮飛行器隕落再入過程跨流域氣動環境與結構熱力響應計算分析

依靠跨流域空氣動力學模擬方法計算得到類天宮飛行器隕落再入飛行繞流物面壓力、熱流分布，使用第4節介紹的氣動環境與結構變形/失效解體非線性力學行為一體化模擬手段，就可實時計算得到結構溫度增量與位移變形，圖7（a）、（b）分別繪出類天宮飛行器隕落至H=120 km、110 km時的物面結構溫度等值線云圖，說明大型復雜結構航天器隕落至120 km高度時，因超高速氣流在物面附近強擾動，致太陽電池翼與對接臺肩部區域出現高溫區；隨著高度進一步降低，隕落至110 km，大氣密度增加，高超聲速氣流在飛行器前部形成脫體激波層，沖擊在靠近安裝支架的太陽電池翼面，形成更為強烈的高溫區，并與激波過渡帶對太陽電池翼面力熱沖擊形成雙重耦合作用，致太陽電池翼因結構的動態響應而變形彎曲，圖7（c）繪出結構動態熱力耦合響應有限元算法模擬該類天宮飛行器無控隕落至H=100 km、V∞=7.6 km/s強氣動力熱環境致結構整體穩態變形情況，計算表明強氣動力熱作用致太陽電池翼安裝支架變形達2 m，超過支架屈服強度，帆板聯接機構將失效、折斷損毀。

類天宮飛行器首次解體后兩艙體經空氣動力配平繞流隕落過程，將持續在外部流場氣動力/熱作用下，物面溫度不斷升高。圖8（a）、（b）分別繪出兩艙體隕落至H=90 km、α=0°與α=20°飛行繞流物面溫度分布，可看出，α=0°時，物面最高溫度在對接口達到1450 K，而α=20°時，物面最高溫度同時出現在對接口與前臺迎風面，達到1550 K。圖8（c）繪出類天宮飛行器兩艙體分別以α=0°與α=20°隕落H=120～90 km的過程中物面最高溫度隨飛行高度降低非線性增加的曲線，從圖中可以看出，當飛行器隕落至100 km左右，考慮與不考慮輻射效應計算得到的結構熱力響應溫度分布產生嚴重偏離，熱輻射開始在熱傳輸中產生顯著影響；當隕落飛行至90 km，兩艙體結構在0°迎角時駐點溫度達到1500℃左右，在20°迎角時兩艙體結構最高溫度超過1500℃，而且輻射效應引起的物面溫度計算偏差高達數百℃。這說明開展大型航天器隕落再入過程中因強氣動力熱環境和輻射傳熱耦合而導致結構熱力響應變形/失效解體的有限元計算研究是一個新的研究方向。

圖7 類天宮飛行器H=120～100 km、V∞=7.6 km/s強氣動力熱環境致結構溫度變化和穩態變形Fig.7 T Temperature distribution and steady state deformation of structure at H=120～100 km and V∞=7.6 km/s around Tiangong type spacecraft by strong aerodynam ic heating

6 服役期滿航天器離軌再入解體數值預報平臺研制與應用

針對非回收類航天器服役期滿離軌隕落回歸地球大氣層過程，遍歷外層空間到地面全飛行流域，單憑任何一種空氣動力學模擬手段難以勝任。作為跨流域空氣動力學模擬方法集成升華對接工程需求與應用，構建了Boltzmann方程可計算建模氣體動理論統一算法、DSMC、N-S/DSMC、（滑移）N-S解算器、低密度風洞實驗測試多種空氣動力學模擬手段［15，17，37］驗證結合，可靠模擬大型航天器從軌道衰降到再入跨流域空氣動力學一體化模擬平臺。以此為基礎，結合發展氣動熱環境與結構傳熱/復合材料熱解燒蝕計算方法［9］、航天器解體分離、多體繞流數值模擬方法［28］與數據庫、軟件工程手段，建立以運動彈道計算為主線，融合氣動力、氣動熱、金屬桁架結構熱力響應變形失效/軟化融熔、復合材料熱解/燒蝕及再入解體（部件級或碎片級）飛行航跡聯合計算分析機制。結合外測軌道信息，研制服役期滿航天器離軌隕落飛行軌道衰降與再入跨流域氣動/結構響應變形失效/燒蝕/解體數值預報模擬系統應用平臺［5，8，13，14，18，28，33］。

圖8 類天宮飛行器兩艙體隕落H=120～90 km強氣動力熱環境致結構溫度變化分布Fig.8 Tem perature distribution in structure during the reentry of H=120～100 km around Tiangong type two-capsule spacecraft by strong aerodynam ic heating

6.1 類天舟一號貨運飛船受控再入解體數值預報應用檢驗

以自制類天舟一號貨運飛船模型為試驗驗證對象，實施運行軌道自由分子流到再入解體近空間連續流各流域飛行繞流模擬預報，圖9繪出類天舟一號貨運飛船受控離軌再入341.2～100.9 km時飛行的繞流壓力分布與航跡預報結果。數值預報模擬表明，類天舟一號再入飛行高度115 km以上，氣動載荷不會使太陽電池翼結構損壞；到110 km時，氣動載荷作用，帆板結構有可能會損壞；100 km以下，氣動載荷作用急劇增大，帆板結構損壞，判斷發生首次解體高度在101.3～100.9 km太陽電池翼撕裂分離。隨后，類天舟一號兩艙體隕落飛行至90 km高度，頭部對接機構16 mm后退量；87～85 km實驗艙段后退量10mm燒穿毀壞。由此判斷類天舟一號發生二次解體高度90～87 km兩艙體實驗艙外殼熔融所致，此時中繼天線斷裂拋出、天舟一號與地面通信中斷，這與實際再入遙測信號消失點88.75 km吻合，偏差1.41%。其后，因天舟一號實驗艙外殼被燒穿，艙內容器或蓄電池模組爆炸可致部件級與碎片級多次解體。

6.2 天宮一號目標飛行器無控再入解體落區數值預報

針對天宮一號目標飛行器超期服役兩年半至2016年3月16日突發功能失效進入無控飛行軌道衰降，利用航天測站雷達成像、光學自適應圖像數據等，建立非合作外測估計天宮一號空間飛行姿態測量手段，實施二維成像三維重構空間飛行姿態及帆板姿態角識別，建立天宮一號空間飛行姿態估計方法，解決天宮一號飛行姿態聯測輸入氣動模擬的條件問題，結合外測軌道星歷數據，對沿軌道衰降飛行的天宮一號空氣動力三力（阻力、升力、側力）三矩（俯仰力矩、偏航力矩、滾轉力矩）開展大規模并行計算，并基于數值結果優化修正當地化橋函數關聯參數快速計算方法，融合彈（軌）道后實施高精度迭代計算。

根據2018年4月1日天宮一號外測姿態數據與2018年4月2日6時22分外測軌道星歷數據，使用系統各執行模塊組成的模擬平臺［18］數值預報：天宮一號于北京時間2018年4月2日7時54分41秒到達120 km再入大氣層，再入速度7492.13 m/s；首次解體高度：110～105 km（帆板鏈接框架斷裂），二次解體高度：資源艙小頭朝前飛行至100～95 km，中繼天線于97.5 km從資源艙飛出至95 km解體，主承力錐臺、軌控發動機83～56 km解體，未存留大的殘骸，多次解體殘骸碎片散布核心區縱向長約1200 km、橫向寬約100 km，圖10（a）繪出天宮一號隕落再入解體殘骸碎片落區數值預報結果與美國航空航天官方網站后續標注的天宮一號再入事后監測圖10（b）所示吻合一致，且落區數值預報結果和事后一天美國聯合太空作戰中心公布的圖10（c）所示的天宮一號墜落南太平洋海域極為相容，與美國太空司令部事后確認沒有大型殘骸落地的監測數據吻合很好，所建立的氣動融合結構/彈（軌）道再入解體數值預報算法模型的正確性與預報精度得到了檢驗。且分析表明，數值預報大型航天器繞地球一圈用時一刻鐘；停止姿態、外測軌道數據供給，6 h后地心慣性系位置矢量預報偏差2.3 km，這就開辟了大型航天器無控再入解體落區數值預報新領域發展方向。

圖9 類天舟一號貨運飛船離軌再入341.2～100.9 km解體飛行繞流壓力分布與航跡預報模擬結果Fig.9 Pressure distribution and track forecast simulation of Tianzhou-1-type cargo spacecraft during the reentry of 341.2～100.9 km

6.3 天宮二號空間實驗室受控再入解體落區數值預報

2019 年在前期研制天宮一號無控飛行軌降再入解體數值預報平臺基礎上，實施了天宮二號超期服役近一年后受控再入解體回歸地球大氣層過程，在2019年7月16日北京航天飛行控制中心提供的天宮二號第二次變軌制動發動機關機后的倒飛理論設計彈道參數基礎上，于7月18日18：00數值預報：天宮二號首次解體高度105～103 km（帆板鏈接框架斷裂）；二次解體高度95～92 km，中繼天線熔融解體高度95 km之后和地面通信中斷，與2019年7月19日晚天宮二號實際再入遙測信號消失點94.775 km吻合一致，偏差0.24%；多次解體殘骸碎片核心散布區域：南太平洋偏東南海域起點西經132.4°、南緯35.2°，終點西經126.7°、南緯37.5°，長約570 km，寬約100 km，較預定海域6000 km×1700 km，減小170余倍，證實數值預報模型能提前鎖定大型航天器受控再入解體殘骸碎片落區散布范圍，圖11（a）繪出天宮二號受控再入解體落區數值預報與軌道外推理論設計預定海域比較驗證情況。隨后為了進一步復核驗證該數值預報模型對大型航天器受控再入解體殘骸碎片落區散布范圍模擬結果正確性，在載人航天工程辦公室、北京航天飛行控制中心建議下，由北京中心軌道室提供天宮二號于2019年7月19日20：33第二次變軌制動發動機關機220～100 km實際再入飛行參數，使用數值預報模型開展復核復算，驗證了圖11（b）和（c）所示為天宮二號二次解體中繼天線毀壞通信中斷與天宮二號實際再入遙測信號消失點及數值預報復核復算落區與預定海域對照結果，證實數值預報彈道與天宮二號實測GPS彈道偏差為數米量級。

圖10 天宮一號無控飛行軌降再入解體落區數值預報結果與美國航空航天、太空司令部事后監測圖示吻合一致Fig.10 The present numerical forecast of uncontrolled reentry disintegration of Tiangong-1 was in good accordancew ith the USA monitoring plot of the falling area

7 結論

1）針對大型航天器服役期滿離軌再入墜毀飛行航跡數值預報問題，基于Boltzmann方程碰撞積分物理分析與可計算建模，建立了含非平衡效應各流域統一Boltzmann模型方程及隱式混合通量數值格式與穩定運行數千數萬CPU核大型航天器隕落再入跨流域氣動力熱繞流環境大規模可擴展并行算法，具有500～45000處理器核擬線性并行計算加速比與80%以上高性能并行計算效率。

圖11 天宮二號受控再入解體落區數值預報與實際再入遙測信號消失點、理論設計預定海域比較驗證Fig.11 Com parison and verification of the present numerical forecast，the actual reentry telemetry signal vanishing point，and the theoretical design of predeterm ined sea area for the controlled reentry of Tiangong-2

2）通過將求解Boltzmann模型方程統一算法對類天宮飛行器兩艙體跨流域氣動并行計算與DSMC、N-S/DSMC、N-S解算器、低密度風洞實驗測試結果間的比較分析，證實GKUA求解大型復雜結構航天器隕落再入解體過程跨流域高超聲速氣動問題精度與可靠性，揭示了大尺度非規則物形繞流多流區共存混合流動輸運現象與變化特點。

3）利用Newmark隱式方法設計熱彈性動力學方程時間上的離散推進，計算速度與位移；使用高精度Crank-Nicolson格式求解熱傳導方程，建立了適于航天器再入強氣動力熱環境致結構響應變形行為隱式有限元算法。通過開展求解Boltzmann模型方程統一算法對豎直平板、類天宮飛行器高超聲速流場計算與結構力-熱耦合響應變形毀壞非線性行為一體化模擬檢驗，證實所提出服役期滿航天器隕落再入跨流域氣動環境與結構響應變形/失效解體一體化算法可靠性，揭示了航天器隕落再入高超聲速繞流強氣動力熱沖擊致物面材料內部溫度分布與結構變形失效非定常演化機制。

4）通過研制服役期滿大型航天器離軌再入跨流域氣動融合結構響應變形失效/燒蝕/解體數值預報模擬系統應用平臺，對自制類天舟一號貨運飛船模型進行了數值預報，與實測數據對比表明所建立氣動融合結構/軌道衰降與再入解體落區數值預報模型算法正確性與高精度預報能力。對天宮二號空間實驗室受控再入解體過程氣動融合結構/彈道數值預報表明，數值預報模型能提前鎖定大型航天器受控再入解體殘骸碎片落區散布范圍。

5）天宮一號無控飛行與天宮二號受控再入解體落區數值預報比較分析表明，受控再入彈道傾角很大，如天宮二號較天宮一號無控再入的彈道傾角大13.5倍，致再入飛行距離、時間上天宮二號較天宮一號縮短近四倍，其殘骸/碎片核心落區散布范圍較無控再入解體落區范圍小若干倍，可根據殘骸碎片預定安全落區反演設計受控離軌控制策略，預期可大大縮小預定安全海域，如前期對外公布的天宮二號預定海域（西經160～90°、南緯30～45°）6000 km×1700 km，而實際數值預報天宮二號再入解體落區散布區域（西經132.4～126.7°、南緯35.2～37.5°）570 km×100 km，較預定區域小170多倍。數值預報模擬平臺可提前較長時間預報鎖定服役期滿大型航天器受控再入解體殘骸碎片落區散布范圍。

本文是階段性工作，搭建了服役期滿大型航天器離軌再入強氣動力熱環境致結構熱力響應變形失效解體數值預報模型算法及無控軌降與受控再入解體數值預報應用平臺框架，下一步擬研究材料變形熔融分子動力學可計算建模與熱力響應有限元算法、統一算法耦合雜交實時模擬，有望在提高服役期滿大型航天器再入解體數值預報精度方面取得更大突破。