大型航天器再入隕落時太陽翼氣動力/熱模擬分析

梁 杰，李志輝，杜波強，方 明

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

0 引言

大型航天器如空間站或衛星在壽命即將結束時，地面控制中心利用剩余的推進劑進行多次剎車離軌控制，將其軌道高度逐漸減低，然后在無控狀態下再入大氣層。當再入到120～100 km左右高度時，在大氣的作用下開始氣動加熱。100 km高度左右，太陽翼電池帆板、大型天線等設備在氣動熱和氣動力的雙重作用下開始撕裂，并與航天器本體分離。隨著航天器表面溫度持續升高，航天器局部材料出現軟化、結構失效，艙體內外部設備陸續與航天器發生分離。解體后的部件在墜落過程中熔化或燃燒，質量逐漸減小，一些較大的部件則會殘留一部分，直到撞擊地面。所有的碎片沿航天器的軌跡方向散布在一個狹長的區域［1］。如果解體后的碎片散落在人口密集區，則會對人類的安全造成破壞［2］。這類危險目標的隕落跟蹤與落點預報均會引起世界各國的重視，如過去的美國天空實驗室、原蘇聯宇宙1402號核動力衛星、我國失控衛星以及俄羅斯“和平”號空間站的隕落，均受到世界各國的重視，提前進行了再入風險評估［3］。國外開展空間碎片再入預測及地面風險評估的研究已有十多年的歷史，并利用工程預測方法形成了相對成熟的軟件系統，如美國的碎片評估軟件(Debris Assessment Software，DAS)［4］和目標再入存活分析工具(Object Reentry Survival Analysis Tool，ORSAT)［5］，歐洲的航天器大氣再入與氣動加熱解體(Spacecraft Atmosphere Reentry and Aerothermal Breakup，SCARAB)［6］等。國內的清華大學利用工程預測手段初步建立了空間碎片再入分析方法和評估軟件［7］。我國的天宮一號目標飛行器由實驗艙和資源艙組成，長約7．8 m的太陽翼電池帆板位于資源艙的兩側，兩個艙體的總長約10．5 m，艙體最大直徑約3．4 m，設計在軌壽命2年。自2011年發射升空后，已經超期服役即將墜入大氣層損毀，需要預先對其隕落的風險進行詳細的評估，同時還要不斷完善這種大型航天器的再入解體和地面風險評估手段，為后續的貨運飛船、空間站、大型衛星的隕落風險預報做準備。

關于太陽翼在軌運行的熱分析，已有一些文獻開展過相關研究［8－9］，但對于大型航天器離軌隕落解體過程中太陽翼帆板的力熱分析國內尚屬空白。而在航天器再入隕落的初期，流動處于稀薄過渡流區域，求解N－S方程的連續流計算方法已經失效，需要采用基于分子碰撞理論的直接模擬蒙特卡洛(Direct Simulation Monte Carlo，DSMC)方法［10］。文獻［11］采用DSMC方法對鈍體返回艙外形在稀薄流域的高溫真實氣體效應進行了計算分析，文獻［12］則采用并行DSMC技術開展了過渡流區復雜噴流干擾流動的模擬，尚未針對帶太陽翼帆板的大型航天器存在復雜激波干擾環境下的氣動特性開展相關研究，因此本文在構建模擬大型復雜航天器高溫氣體熱化學非平衡特性的DSMC數值方法的基礎上，對壽命末期離軌隕落的低軌航天器再入稀薄流域120～90 km高度的氣動力、熱特性進行了模擬分析，重點分析了作用在水平和垂直放置的單側太陽翼帆板上的力、力矩和峰值壓力、熱流，為下一步的飛行器結構解體分析提供數據支撐。

1 DSMC數值模擬方法

DSMC方法是通過模擬仿真分子的運動、與邊界的碰撞以及分子之間的碰撞隨時間的變化，并存貯仿真分子的位置坐標、速度分量以及內能，最后通過統計網格內仿真分子的運動狀態實現對真實氣體流動問題的模擬。每個仿真分子代表1012～1014個真實氣體的分子，該方法的關鍵是在時間步長Δt內將分子的運動與碰撞解耦。在Δt時間內讓所有分子運動一定的距離并考慮在邊界的反射，然后計算此Δt內具有代表性的分子間的碰撞。

1．1計算網格及網格自適應

在DSMC方法中，流場中的網格是用來選取可能的碰撞分子對以及對宏觀流動參數取樣。“天宮一號”目標飛行器結構復雜，涉及到局部凸起物、薄壁和細支架等需要精細描述的復雜壁面，同時計算區域需要覆蓋飛行器本體、展開的太陽翼帆板以及激波干擾引起的復雜流動區域，因此網格的計算效率和高可靠性至關重要。在早期的研究中作者采用兩級直角網格結構對美國“獵戶座”飛船在稀薄流域的氣動力系數進行了計算，并與國際知名軟件DS3V和DAC的計算結果一致［13］。針對“天宮一號”復雜外形，本文采用計算效率高的流場直角網格和對飛行器幾何外形描述精確的表面三角形非結構網格(如圖1)，建立了混合網格結構的DSMC數值模擬方法。在描述物面幾何形狀的非結構網格建立以后，直接將其嵌入到直角網格的流場中，使DSMC計算對流場網格的依賴程度大大降低。同時通過判斷分子運動軌跡方程和物面三角形面元上任一點的位置方程，唯一確定出分子與物面的碰撞點坐標，解決了這種混合網格流場分子運動與物面碰撞的難題。對分子在物體三角形面元上碰撞、反射前后的流場參數進行統計取樣就可以獲得飛行器的整體氣動力特性以及表面力、熱載荷分布。

圖1 計算網格及網格自適應Fig．1 Computational grid and self-adaption

為了解決流場中因激波壓縮、激波干擾以及氣體膨脹后出現的密度梯度急劇變化的流動特征，計算中采用了網格自適應的策略［14］。即在背景網格的基礎上，根據流場中密度梯度的變化分別對碰撞網格和取樣網格進行細化(如圖1)。由于流場的梯度沿各個方向的變化是不同的，梯度變化大的方向網格細分的更密一些，因此沿三個坐標方向是各自獨立地進行自適應，碰撞分子則是在自適應后最小的亞網格內選取，保證了計算的空間精度。

1．2分子轉動和振動松弛模型

天宮飛行器再入隕落的速度達到7．6 km/s，需要考慮氣體分子的轉動和振動激發以及化學反應。連續流氣體分子的轉動松弛碰撞數和振動松弛碰撞數需要經過內自由度的修正才可用于DSMC的模擬［15］，修正如下:

式中:ζt、ζr和ζv分別是分子平動、轉動和振動自由度。Parker［16］給出了連續流轉動松弛碰撞數ZCr的表達式

式中:T*是氣體分子作用勢的特征溫度，(Zr)∞是實驗測定的極限值。

連續流的振動松弛碰撞數定義為

式中:τc為分子平均碰撞時間，τMW為振動松馳時間，τP為高溫修正項。

由Landau-Teller理論結果給出

Millikan和White［17］在高溫激波管內通過干涉儀觀測了氣體分子振動松馳過程，在104K的溫度范圍內，給出了如下擬合參數

式中:mr為分子折合質量，θv是振動特征溫度。Park［18］給出了高溫修正的 τP為

式中:n是數密度，σv是分子振動碰撞截面(=10－20m2)，ˉc是分子平均熱運動速度。

1．3化學反應模型及非彈性碰撞模擬策略

通常反應速率常數可以寫成下面的方程形式:

式中:Ea是反應中需要的活化能，k是Boltzmann常數，Λ和η是常數，由實驗定出。

化學反應幾率可以推導為下面的表達式［10］

式中:ε是對稱因子，對于不同類分子ε=1，對同類分子ε=2。ˉξ為碰撞分子的平均內自由度，ωAB為碰撞對的粘性溫度指數。當碰撞中的總碰撞能量Ec＞Ea時，反應截面σR與總碰撞截面σT的比值就是反應發生的概率，這種化學反應模型也稱為TCE(Total Collision Energy)模型。

在DSMC模擬過程中為簡化碰撞算法，假設非彈性碰撞的化學反應過程和內能交換是相互獨立的事件。由于碰撞分子發生化學反應的幾率在所有非彈性碰撞事件中是最低的，所以首先判斷其是否發生化學反應。離解反應、交換反應等不同的化學反應過程在滿足反應幾率的條件下都有可能發生。在無法獲知反應分子遵循何種反應軌跡生成新的反應生成物時，可以假定每種反應過程都是獨立的。對所有可能發生化學反應的幾率求和得到總的反應幾率，如果Pr＞Rf，碰撞分子發生化學反應［19］。具體發生的是何種反應類型，則同樣根據幾率/Pr與隨機數Rf的比較來選取，然后按照相應的化學反應類型產生新的生成物(分子或原子)。

對于不發生化學反應的碰撞分子，即Rf＞Pr時，則根據分子振動抽樣概率，判斷是否發生振動自由度的激發;對于未發生振動激發的分子，再根據分子轉動抽樣概率，判斷是否發生轉動自由度的激發;對于發生內自由度激發的分子，按相應的能量交換模型進行能量的再分配，否則分子按彈性碰撞處理。

由于本文模擬的高度較高，僅考慮五組元空氣(O2，N2，O，N，NO)的離解反應和交換反應，燒蝕的影響暫不作考慮。

1．4 DSMC并行算法

DSMC并行算法采用區域分解的策略，根據計算的處理器數(或CPU核心數)的多少將計算區域劃分為等量的子區域，每個處理器在其分配的子區域內部獨立地計算模擬分子的碰撞和遷移，離開子區域的模擬分子把攜帶的信息傳遞給對應子區域的處理器。并行計算總的時間包括每個處理器計算碰撞、遷移的時間、處理器之間的通訊時間以及各處理器之間為同步而等待的時間。提高并行計算的效率主要是通過減少通訊和同步等待的時間來實現。

為了減少不同處理器計算時間的差別，通常采用負載平衡技術［20－22］，本文采用靜態隨機負載平衡方法用于解決不同處理器之間的計算時間同步問題［12］，該方法基于概率近似原理，將計算區域的全部網格平均分配給指定的所有處理器。由于采用相同的概率隨機選取流場網格，按照均分后的數量分配給每個處理器，因此當計算區域的網格數量較多時，每個處理器包含近似相等的物體邊界、高密度流動區域和稀薄氣體區域的網格數，這樣每個處理器的計算負載也非常接近。

2 計算結果分析

本文選取天宮一號目標飛行器的簡化模型作為研究對象，飛行器上的太陽翼帆板朝向考慮兩種情況，分別是板面平行于和垂直于飛行器中心軸線。計算的來流速度為7．6 km/s，高度范圍為120 km～90 km，每隔10 km一個狀態，攻角(Angle of Attack，AOA)為0°和20°，壁面溫度設置為500 K、壁面反射模型為Maxwell模型(80%的漫反射、20%鏡面反射)，僅計算對稱的半個流場，計算中考慮了分子的轉動、振動能量激發，90 km高度時還考慮了五組元空氣的離解反應和交換反應。

2．1太陽翼帆板平行于目標飛行器

圖2給出了太陽翼帆板平行于天宮一號目標飛行器情況下，0°攻角下不同再入高度時的壓力等值線分布云圖。從圖中可以看出120 km、110 km高度的大氣密度相對較低，盡管飛行速度達到了7．6 km/s，但飛行器的頭部并未形成明顯的脫體激波，而是逐漸增強的壓縮過程，氣體在整個太陽翼帆板的側面也產生明顯的壓縮。圖中也顯示出，飛行器頭部簡化的對接機構形狀對飛行器前方壓縮氣流的擾動影響，而隨著飛行高度的降低，這種影響不斷減弱，直至形成較強的頭部脫體激波。100 km和90 km的頭部激波作用在太陽翼帆板上產生了比駐點值高幾倍的峰值壓力和峰值熱流。

圖2 太陽翼水平放置時0°攻角不同高度壓力分布Fig．2 Pressure distribution of 0°AOA in different altitude with horizontal solar array

圖3給出了120 km和90 km在20°攻角下的壓力分布云圖。圖中顯示出，由于攻角的存在，除了在飛行器的頭部產生壓縮以外，整個飛行器艙體的迎風面也會產生較強的壓縮，使艙體和太陽翼帆板的迎風受力面增大，作用在太陽翼帆板上的力和力矩都有較大增加。

圖3 太陽翼水平放置時20°攻角不同高度壓力分布Fig．3 Pressure distribution of 20°AOA in different altitude with horizontal solar array

表1、表2分別給出了0°攻角和20°攻角不同高度情況下，作用在單側太陽翼帆板上的軸向力、偏航力矩和峰值壓力、峰值熱流。由于太陽翼帆板平行于軸向方向，0°攻角時只有帆板的側面為主要受力面，因此作用在帆板上面的軸向力和偏航力矩(力矩作用點取帆板與艙體的連接處)相對較小，但在100 km和90 km由于頭部激波的干擾作用，使干擾區域的峰值熱流非常高，100 km時的峰值熱流為150 kW/m2左右，而到90 km時則達到了1．2 MW/m2。20°攻角時，盡管作用在帆板上的峰值壓力和峰值熱流與0°攻角時相差不大，但作用的軸向力和偏航力矩明顯增大。因此，當太陽翼帆板平行放置時，激波干擾引起較強的氣動加熱可能最先破壞太陽翼帆板的結構，然后才是在氣動力的作用下將其撕毀并脫離飛行器艙體。

表1 作用在水平太陽翼上的氣動力和熱流峰值(0°攻角)Table 1 Aerodynamics and peak heat flux acting on horizontal solar array at 0°AOA

表2 作用在水平太陽翼上的氣動力和熱流峰值(20°攻角)Table 2 Aerodynamics and peak heat flux acting on horizontal solar array at 20°AOA

2．2太陽翼帆板垂直于目標飛行器

圖4給出了太陽翼帆板垂直于目標飛行器0°攻角、不同再入高度時的壓力等值線分布云圖。圖中顯示出的流場結構與圖2相似，只是此時太陽翼帆板有較大的迎風面，其前方的壓縮層更厚(120 km、110 km)，激波/激波、激波/邊界層的干擾更強、更復雜(100 km、90 km)，同時也與艙體產生了較強的相互作用，該區域的激波干擾使資源艙表面的壓力和熱流上升，尤其是在太陽翼帆板與艙體的連接處附近，更容易使太陽翼帆板遭到破壞，另外作用在帆板上的力和力矩明顯比水平放置時大得多，因此垂直放置的太陽翼帆板會比水平放置的太陽翼帆板更早地在氣動力、氣動熱的作用下被撕裂并脫離飛行器艙體。

圖5給出了90 km高度、0°攻角下的壓力等值線分布，可以更加詳細地分析復雜激波干擾下的流動特征。頭部脫體激波與太陽翼帆板的脫體激波相互作用并產生一道正激波對帆板表面形成非常強烈的壓縮，使壓力急劇上升。而在艙體表面附近區域，由于氣體的膨脹，壓力迅速下降，因此在頭部脫體激波、艙體和太陽翼帆板表面之間的空域形成較大的逆壓梯度，引起流動分離，同時產生的分離激波也與頭部脫體激波相互干擾。圖中還可以看出，從形成的正激波作用在帆板表面的位置開始，沿帆板左右各形成一股“噴流”，朝向艙體的“噴流”在較大的壓力梯度下加速至超聲速并直接撞擊在資源艙的表面，形成局部的脫體激波。由于在較大的區域存在著逆壓梯度，形成一個大的分離渦，而在帆板脫體激波影響的區域僅有微弱的流動分離，如圖6所示。

圖4 太陽翼垂直放置時0°攻角不同高度壓力分布Fig．4 Pressure distribution of 0°AOA in different altitude with perpendicular solar array

圖5 90 km、0°攻角的壓力等值線分布Fig．5 Pressure contours distribution of 90 km at 0°AOA

圖7、圖8分別給出了90 km高度、0°和20°攻角下的太陽翼帆板表面的壓力和熱流分布的對比，0°攻角作用在帆板表面的高壓力和高熱流區域近似為條狀，而20°攻角時則近似為弧形，并且0°攻角時作用在帆板上的壓力和熱流強度高于20°攻角情況。從表3和表4給出的作用在帆板上的力和力矩表明，兩種攻角下的力和力矩相差并不明顯。在相同高度，作用在帆板的力和力矩要比水平放置時大一個量級。盡管100 km高度帆板表面的峰值熱流比帆板水平放置時有較大降低，但由于100 km高度時太陽翼帆板所受的力和力矩是帆板水平放置時在90 km高度時的近3倍，超過了600 N和2700 N·m，到90 km時更是達到了3950 N和17000 N·m。因此太陽翼帆板垂直放置時，在100 km高度就有可能在氣動力、氣動熱的雙重作用下被撕裂并脫離飛行器本體。

圖6 90 km、0°攻角的流線分布Fig．6 Stream lines distribution of 90 km at 0°AOA

圖7 太陽翼帆板表面壓力分布Fig．7 Pressure distribution of solar panels surface

圖8 太陽翼帆板表面熱流分布Fig．8 Heat flux distribution of solar panels surface

表3 作用在垂直太陽翼上的氣動力和熱流峰值(0°攻角)Table 3 Aerodynamics and peak heat flux acting on perpendicular solar array at 0°AOA

表4 作用在垂直太陽翼上的氣動力和熱流峰值(20°攻角)Table 4 Aerodynamics and peak heat flux acting on perpendicular solar array at 20°AOA

3 結論

通過上面的計算分析，得出以下幾點初步結論:

1)100 km以上高度，來流氣體對目標飛行器及太陽翼帆板是漸進的壓縮過程，作用在太陽翼帆板上的力和力矩、峰值壓力和熱流相對較低。

2)100 km及以下高度，由于飛行器頭部脫體激波與太陽翼帆板產生比較強的干擾，帆板所受的氣動力、力矩以及表面上的峰值壓力、熱流隨高度降低明顯增大，且增大的速率加快。

3)太陽翼帆板水平放置時，激波干擾引起的側緣峰值熱流達到了兆瓦量級，較高的氣動加熱可能最先引起帆板結構的破壞，然后在氣動力的共同作用下將其撕毀并脫離飛行器本體。

4)太陽翼帆板垂直放置時，飛行器頭部脫體激波與帆板脫體激波產生更強烈、更復雜的激波/激波、激波/邊界層的干擾。在相同高度，作用在帆板上的氣動力和力矩要比水平放置時大一個量級，因此太陽翼帆板垂直放置時要比水平放置時更快地在氣動力、氣動熱的雙重作用下被撕裂并脫離目標飛行器。

［1］ 朱魯青．美國“高層大氣研究衛星”隕落事件分析［J］．國際太空，2011，11:27－33．［Zhu Lu-qing．Analysis of UARS reentry event［J］．Space International，2011，11:27－33．］

［2］ 都亨，張文祥，龐寶軍，等．空間碎片［M］．北京:中國宇航出版社，2007:1－128．

［3］ 吳連大，王昌彬．危險目標再入跟蹤預報研究［J］．紫金山天文臺臺刊，2000，19(2):100－105．［Wu Lian-da，Wang Chang-bin．Research on tracking and prediction of risk object reentry［J］．Publications of Purple Mountain Observatory，2000，19(2):100－105．］

［4］Opiela JN，Hillary E，Whitlock D O，et al．DASuser's guide［R］．Houston:Orbital Debris Program Office in NASA Johnson Space Center，TX，2007

［5］Bouslog SA，Ross B P，Madden C B．Space debris reentry risk analysis［R］．AIAA 1994－591，1994．

［6］Fritsche B，Klinkrad H，Kashkovsky A，et al．Spacecraft disintegration during uncontrolled atmospheric re-entry［J］．Acta Astronaut，2000，47(2－9):513－522．

［7］Wu Z N，Hu R F，Qu X，et al．Space debris reentry analysis method and tools［J］．Chin J．Aeronaut，2011，24:387－395．

［8］Li JL，Yan SZ．Thermally induced vibration of composite solar array with honeycomb panels in low earth orbit［J］．Applied Thermal Engineering，2014，71:419－432．

［9］Li JL，Yan SZ，Cai R Y．Thermal analysis of composite solar array subjected to space heat flux［J］．Aerospace Science and Technology，2013，27(1):84－94．

［10］Bird G A．Molecular gas dynamics and the direct simulation of gas flows［M］．Oxford:Clarendon Press，1994．

［11］ 黃飛，趙波，程曉麗，等．返回器高空稀薄氣動特性的真實氣體效應研究［J］．宇航學報，2014，35(3):283－290．［Huang Fei，Zhao Bo，Cheng Xiao-li，et al．Real gas effects of reentry vehicle aerodynamics under hypersonic［J］．Journal of Astronautics，2014，35(3):283－290．］

［12］ 梁杰，閻超，楊彥廣，等．過渡區側向噴流干擾的并行DSMC數值模擬研究［J］．宇航學報，2011，32(5):1012－1018．［Liang Jie，Yan Chao，Yang Yan-guang，et al．Parallel DSMC simulation of lateral jet interaction in rarefied transitional region［J］．Journal of Astronautics，2011，32(5):1012－1018．］

［13］ 梁杰，閻超，杜波強．基于兩級直角網格結構的三維DSMC算法研究［J］．空氣動力學學報，2010，28(4):466－471．［Liang Jie，Yan Chao，Du Bo-qiang．An algorithm study of three-dimensional DSMC simulation based on two-level Cartesian coordinates grid structure［J］．Acta Aerodynamica Sinica，2010，28(4):466－471．］

［14］ 梁杰，閻超，李志輝，等．稀薄過渡流區橫向噴流干擾效應數值模擬研究［J］．空氣動力學學報，2013，31(1):27－33．［Liang Jie，Yan Chao，Li Zhi-hui，et al．Numerical investigation of lateral jet interaction effects in rarefied transition flow regime［J］．Acta Aerodynamica Sinica，2013，31(1):27－33．］

［15］Lumpkin F E，Hass B L，Boyd I D．Resolution of differences between collision number definition in particle and continuum simulations［J］．Physics of Fluids，1991，A3:2282－2284．

［16］Parker J G．Rotational and vibrational relaxation in diatomic gases［J］．Physics of Fluids，1959，2(4):449－462．

［17］Millikan R C，White D R．Systematics of vibrational relaxation［J］．The Journal of Chemical Physics，1963，39(12):3209－3213．

［18］Park C．Problems of rate chemistry in the flight regimes of aeroassisted orbital transfer vehicles［C］．Progress in Astronautics and Aeronautics，AIAA，New York，USA，Jule9－14，1985:511－537．

［19］Gimelshein N E，Gimelshein S F，Levin D A，et al．Reconsideration of DSMC models for internal energy transfer and chemical reactions［C］．23rd International Symposium of Rarefied Gas Dynamics，Whistler，Canada，July 20－25，2002:349－357．

［20］Ivanov M，Markelov G，Taylor S，et al．Parallel DSMCstrategies for 3D computation［C］．Proc．Parallel CFD’96，Capri，Italy，July 13－16，1996:485－492．

［21］LeBeau G J．A parallel implementation of the direct simulation Monte Carlo method［J］．Comput．Methods Appl．Mech．Engrg．1999，174:319－337．

［22］Wu J S，Tseng K C．Parallel DSMC method using dynamic domain decomposition［J］．Int．J．Numer．Meth．Engng．2005，63:37–76．