航天器無控再入解體非規則碎片散布范圍分析研究

蔣新宇，黨雷寧，李志輝*，李四新，唐小偉

（1.中國空氣動力研究與發展中心超高速空氣動力研究所，綿陽621000）2.國家計算流體力學實驗室，北京100191）

1 引言

航天器被送入太空，在航天器服役期滿或壽命末期，部分大型航天器（如和平號空間站等）會采取主動離軌再入策略；部分航天器（如同步軌道衛星等）會進入墳場軌道，騰出寶貴的軌道資源，以免對正常軌道上的衛星構成威脅；而絕大部分失效航天器則是沿著螺旋型橢圓軌道衰降無控飛行，最終再入大氣層解體墜毀。失效航天器和因空間活動（如太空行走、航天器相互碰撞、燃料泄漏等）產生的小型零部件、殘骸等，都統稱為空間碎片。空間碎片會威脅在軌航天器的安全運行，需要對這類航天器做好監測防護或主動預報［1-3］，必要時進行主動軌道機動來躲避高風險碎片［4］。隨著空間碎片的增多，碎片清理也成為一個越來越受到重視的研究方向。Bérend等［5］給出了一種空間碎片清理策略的雙目標優化方法；Forshaw等［6］完成了一次在軌碎片清理試驗，分別對編網捕獲、魚叉捕獲、視覺導航交會、阻力帆技術進行了驗證，展示了未來可能的軌道碎片清理技術。受大氣阻力、空間環境和引力等因素影響，空間碎片的軌道高度都會逐漸衰降，特別是低軌道碎片受大氣阻力影響明顯，800 km以下的空間碎片一般會在十余年內再入大氣層［7］。目前每年都會發生100次左右的航天器再入事件，除去有特殊熱防護設計的返回式衛星/飛船等，大部空間物體會在再入環境的強氣動力/熱作用下熔融/燒蝕/解體，大部分質量熔融/燒蝕殆盡，但一般仍有10%～40%的殘余質量會到達地面，并可能對地面造成危害［8］。

個別較特殊的再入事件會人為主動處理，如美國軍方聲稱為避免有毒材料危害地球，采用標準-3導彈在約250 km高度將USA-193衛星擊毀。公開資料顯示大部分碎片會在24～48 h內墜落，余下的碎片則在40天內全部墜落。Pardini等［9］采用公開的觀測數據對此事件進行了模擬，研究了衛星被摧毀后的碎片云散布范圍和演化過程，并希望有空間監測能力的國家、機構對此類信息做到公開透明，以確保公共空間安全。

對于空間物體再入的研究主要分為2個方向，一個方向是再入環境下的材料毀壞，如Prevereaud等［10］對航天器再入過程的燒蝕破壞作了數值分析，并采用試驗方法研究了合金材料的氧化/熔融過程；另一個方向是復雜外形物體的彈道計算與落點散布分析，如Mehta等［11］建立了一種考慮參數不確定性的高維模擬方法，其模擬的落點較常規高斯分布和橢圓分布相去甚遠；Falsone等［12］提出了一種基于仿真的碎片落點散布分析方法，結果優于之前使用的協方差傳播方法，并認為以解體點為模擬起始點是一個有趣的研究方向。

目前主要的空間物體再入預測軟件包括美國的DAS［13］、ORSAT［14］和歐洲的SCARAB［15］等。胡銳鋒等［16-17］也開發了空間碎片方面的工程預測軟件DRAPS。空間物體再入預測軟件按分析方法可以分為面向物體法和面向航天器法，除歐洲的SCARAB外，上述軟件均采用了面向物體法［18］。面向物體法將航天器以及解體后產生的部件外形認為由簡單的基本幾何模型組成，且其幾何信息可用若干特征參數描述；面向航天器法則盡可能模擬真實的航天器外形，采用基于表面網格的飛行器和部件模型，理論上具有更高的預測精度，但建模較復雜并需要較長的計算時間。以上軟件主要是模擬了航天器部件的再入情況，分析可能到達地面的殘骸，分析對象都是已知且提前預設好的航天器結構部件，暫時還沒有考慮航天器在再入過程中因熔融、燒蝕、爆炸等因素突然產生的大量微小碎片。微小碎片因為數量大，散布范圍廣，可能造成更大的地面傷害概率，因此也需要研究。

本文旨在開展航天器再入解體后產生的大量非規則微小碎片地面散布范圍可計算建模，并研究其分布規律。根據碎片幾何形狀與尺度，將碎片分為有限類，利用空間分組策略實現對碎片運動隨機性的模擬，使用當地化快速方法計算不規則碎片的空氣動力特性，采用六階顯式Adams積分方法求解彈道方程，最終獲得碎片的地面散布范圍。

2 非規則微小碎片散布范圍可計算建模

2.1 碎片的幾何分類與質量分布模型

采用面向物體法的思想，將碎片抽象為具有簡單幾何形狀的有限類，對各類碎片進行彈道模擬，落點相近的類可以整合。根據大量模擬結果，本文將碎片按幾何形狀分為塊狀、片狀和桿狀等5類。同時注意到碎片面質比對彈道特性影響很大，所以在確定了幾何分類方案后，還要考慮碎片尺度的影響。本文將碎片按尺度分為3類，最終確定的碎片分類方案見表1所示。某類碎片的單個質量越小，則該類碎片的數量就越大，且質量和數量在取對數后基本呈線性關系［19］。為此，可設置n為某類碎片的數量，如式（1）所示：

式中，n為該類碎片的單個質量，C基于質量約束條件確定。k根據解體類型不同而不同，在本文中取為0.553。

表1 碎片有限分類方案*Table 1 Finite classifying scheme for debris

結合表1與式（1），可確定各類碎片的權重。例如假設碎片總質量為200 kg，3個尺度取為0.1m，0.01 m和0.001 m，則各類碎片數量與質量的關系如圖1所示。

2.2 碎片運動隨機性的模擬

服役期滿航天器再入多次解體生成大量碎片之前，航天器作為一個整體運動，具有確定的空間位置、速度、姿態角。在熔融、燒蝕、爆炸等因素下突然生成碎片的時間點，產生大量的不同形狀、不同尺度的碎片，在按形狀、尺度將碎片分為有限類以后，各類碎片的數目依然眾多。采用李志輝等［20］研究建立的箔條云隨機統計模擬方法來處理每一類碎片。同一類中的各碎片具有隨機的飛行姿態角，以相同的分離速度（絕對值）沿不同的空間方位角散開［21］，不同種類碎片的分布相互獨立。為實現對同類碎片的跟蹤模擬，可采用如下規則：

圖1 各類碎片的數量與質量關系Fig.1 Relationship between quantity and mass of various debris

1）將整個空間方位角均勻劃分為ng=Mφ×Mθ個組，每組具有確定的方位角 ( φ，θ)；

2）每組根據該組的方位角 ( φ，θ)和相對分離速度得到本組碎片生成后的初速度，初速度保證了各碎片組之間具有相互散開的趨勢；

3）對于具有N1個碎片的種類1，通過生成N1個隨機數R1，再按規則映射到某個方位角中去即可。同理進行其它種類碎片的映射；

4）若某個種類的碎片數目Ni?ng，由于在空間方位上是均勻分布，為減少隨機數的生成和映射，可以先將大部分數目（比如80%）的碎片平均分配到各空間方位角，再將余下部分按規則（3）進行映射。

為實時跟蹤有限的Mφ×Mθ組碎片運動飛行情況，需要將Mφ×Mθ組碎片一一映射分布到具有（αL，βM）飛行迎角與側滑角的各小組Mα×Mβ中，使得每組碎片具有若干碎片個體，且在每個d t時間內飛行姿態相對不變［21］。

依次處理每個組［21-22］：先生成一個隨機數，再把隨機數按規則映射到一個角度值，如果該角度值未使用過，就進行隨機篩選；如果已經使用過，則需要重新生成隨機數。針對大型航天器再入多次解體形成碎片問題，需要重復多次生成合適隨機數，采用如下規則進行碎片組飛行姿態角的隨機篩選確定：

1）一次性生成Mφ個均勻隨機數，存放在數組Rφ中；

2）對Rφ中的元素進行排序，記錄每個元素所對應的序號，存放在數組Iφ中；

3）將飛行迎角均勻劃分成Mφ個組，Iφ即可根據序號映射到不同的飛行迎角；

4）同理，可對飛行測滑角進行類似的映射；

5）分別從2個方向完成對Mα×Mβ碎片組的上述循環，即可確定各碎片組的飛行迎角和側滑角。

使用上述隨機動力學篩選方法，可將Mφ×Mθ組碎片一一映射到具有（αi，βj）飛行迎角與側滑角的各小組Mα×Mβ中，使得同一類碎片能被模擬成有限的碎片組，各碎片組在每個d t時間內具有確定的碎片個體數Nij，以確定的飛行速度Vij、飛行迎角αi和側滑角βj飄落飛行，于是可研究發展相關空氣動力學理論、方法，計算確定此d t時間內碎片個體所受的空氣動力學特性［23-25］。

2.3 再入跨流域氣動特性的當地化快速算法

常用的數值方法如直接模擬蒙特卡羅（DSMC）方法、N-S方程解算器與求解玻爾茲曼（Boltzmann）模型方程、氣體動理論統一算法（GKUA）等，由于計算量太大，不能滿足彈道計算中的氣動特性快速獲取要求，本文采用了當地化快速算法［26］。

當地化快速算法以稀薄氣體高超聲速繞流當地化面元法計算理論為基礎，將非規則解體物面劃分成若干塊小曲面，對于每個小曲面，選用一個小的平面三角形面元或四邊形面元來代替，利用面元逼近復雜解體物形。這樣，計算小曲面上的氣動力就轉換成計算面元上的氣動力，將這些面元上的氣動力累加起來就可以得到整個非規則解體物的氣動特性。當面元分得極其細小時，由此引起的誤差就會變得很小。

當地面元所受氣動力系數依賴于來流和當地性質，如當地迎角、當地表面作用等。對于給定的入射角，面元上經歸一化的壓力和摩擦力系數表達為式（2）：

式中，CP和Cf分別是過渡流區域的壓力和摩擦力系數，下標FM和Cont分別表示自由分子流和連續流；FP，b和Ff，b分別是壓力和摩擦力橋函數，下標b表示橋函數，依賴于獨立參數當地努森數Kn、壁溫總溫比TW/T0和入射角θ。橋函數的意義在于使氣動系數在自由分子流區和連續流區分別趨近于理論值，而在中間過渡流區實現光滑過渡。

分別使用氣體動理論統一算法（GKUA）［23-24］、DSMC方法［25，27］或N-S方程解算器對殘骸碎片簡化外形典型繞流狀態進行精細計算，獲取典型狀態的精確氣動數據。調試修正非規則碎片的橋函數關聯參數，使得當地化快速算法的結果在典型狀態與數值算法吻合后，則其它任意狀態的氣動數據可由式（2）直接計算得出。

2.4 碎片運動彈道計算

2.2 節的隨機分組方法已經考慮了飛行姿態的影響，則碎片運動的彈道計算中只需計算三自由度彈道。矢量形式的彈道方程如式（3）所示：

式中，n為碎片的瞬時質量；V為碎片的瞬時速度矢量；R為碎片空氣動力；G為地球引力，G=n▽U；U為地球的引力勢；Fe為離心慣性力，Fe=n rΩ×( Ω×r)；Ω為地球的自轉角速度矢量，r是地心慣性坐標系中質心的矢徑；Fco為科氏慣性力，Fco=2 n( Ω×V)。

為求解彈道方程（3），采用六階顯式Adams-Bashforth格式，在J2000.0地心慣性坐標系下進行數值積分。對于微分方程式（4）

其六階顯式Adams-Bashforth格式為式（5）：

式中，h為時間步長，un+1為待求的下一時間步函數值，un為當前時刻函數值，fn，fn-1，…，fn-5為當前時刻及當前時刻之前1～5時間步的右函數值。在計算啟動時刻，由于時間步結果不足，無法直接采用此數值格式，需要采用龍格-庫塔等其它數值方法來計算初始的幾個時間步。

3 碎片散布范圍計算分析

根據第2節航天器再入多次解體后產生的大量碎片存活墜落地面散布范圍可計算模型方法，對航天器解體后的碎片進行氣動融合彈道的全程模擬。本文以數值預報某大型航天器再入為例，其解體點位置為北緯42°，東經0°，海拔高度100 km；解體點速度為7500 m/s，朝向正東；彈道傾角為0.1°；碎片總質量為200 kg。對該航天器無控再入解體碎片散布范圍計算分析，圖2、圖3分別繪出碎片在地面散布的數量與質量分布。可以看出，在該給定初始狀態下，解體碎片的縱向散布范圍可達5 000多公里，而且碎片質量越大，其數量也就越少，航程也就越遠。這是因為在再入環境主要受重力和空氣動力影響，空氣動力與尺度平方（橫截面積）正相關，而重力與尺度立方（體積）正相關。在尺度減小時，重力下降得更快，相對而言，空氣動力對碎片運動的影響就逐漸增大。圖4繪出特征尺度為0.1 m的片狀碎片在地面的分布范圍。計算表明，由于考慮到了飛行姿態的隨機性，碎片在橫向也會散開，但橫向航程較縱向航程小很多，寬度約為100 km。

圖2 碎片在地面散布的數量分布（φ=0.1°）Fig.2 Quantity distribution of debris on the ground（φ=0.1°）

圖5繪出特征尺度為0.01 m的立方體狀碎片飛行軌跡與碎片云質心彈道的關系，其中藍色點劃線為碎片云質心運動彈道，彩色圓圈為該組碎片不同時刻的空間位置，其顏色代表飛行馬赫數。由圖中可看出碎片飛行軌跡與質心運動彈道的一致性很好。

圖3 解體碎片在地面的質量分布（φ=0.1°）Fig.3 M ass distribution of debris on the ground（φ=0.1°）

圖4 特征尺度為0.1m的片狀碎片在地面的分布Fig.4 Distribution of sheet-like debris（characteristic scale is 0.1 m）on the ground

圖5 特征尺度0.01 m立方體狀碎片飛行軌跡與質心彈道Fig.5 Flight path of cube-like debris w ith the characteristic scale of 0.01 m

為了研究航天器再入解體彈道傾角對落區散布影響，采用前述擬定計算狀態，僅將彈道傾角變更為1.1°，所得落區散布見圖6、圖7所示。對比分析圖2、圖3可以發現，碎片的縱向散布范圍大大減小，約為2500 km。該結果表明彈道傾角對落區散布范圍影響很大，彈道傾角越大，再入飛行時間越短，落區散布范圍越小，在后續服役期滿大型航天器受控再入解體數值預報與離軌控制策略的制定與分析中，應特別注意這方面的研究。

圖6 碎片在地面散布的數量分布（φ=1.1°）Fig.6 Quantity distribution of debris on the ground（φ=1.1°）

圖7 碎片在地面散布的質量分布（φ=1.1°）Fig.7 M ass distribution of debris on the ground（φ=1.1°）

4 結論

1）對碎片按幾何形狀和空間方位角的分組方案，可以實現對碎片外形非規則性和運動隨機性的統計模擬，進而實現對航天器再入解體生成的大量非規則碎片的整體分析。

2）在一次解體過程生成的碎片，單個碎片質量越大的碎片組，其碎片數量越少，航程越遠。碎片在地面的散布整體呈細長條狀，縱向可達數千公里，橫向僅約百余公里。

3）碎片飛行軌跡與三自由度質心運動彈道一致。起始點彈道傾角對落區散布范圍影響很大，在實際應用中應確保起始點彈道傾角的高精度獲取。

本文工作屬初步研究，尚未考慮碎片生成后的進一步熔融/燒蝕效應，如果考慮，碎片質量將在飛行過程中逐漸減小，并有相當一部分碎片根本無法到達地面。本文所得碎片地面散布結果較實際數值預報散布范圍偏大，有待進一步開展碎片生成后的氣動熱軟化熔融/熱解燒蝕模擬、多次解體綜合分析等工作。