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

（1.中國(guó)空氣動(dòng)力研究與發(fā)展中心超高速空氣動(dòng)力研究所，綿陽(yáng)621000）2.國(guó)家計(jì)算流體力學(xué)實(shí)驗(yàn)室，北京100191）

1 引言

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

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

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

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

本文旨在開展航天器再入解體后產(chǎn)生的大量非規(guī)則微小碎片地面散布范圍可計(jì)算建模，并研究其分布規(guī)律。根據(jù)碎片幾何形狀與尺度，將碎片分為有限類，利用空間分組策略實(shí)現(xiàn)對(duì)碎片運(yùn)動(dòng)隨機(jī)性的模擬，使用當(dāng)?shù)鼗焖俜椒ㄓ?jì)算不規(guī)則碎片的空氣動(dòng)力特性，采用六階顯式Adams積分方法求解彈道方程，最終獲得碎片的地面散布范圍。

2 非規(guī)則微小碎片散布范圍可計(jì)算建模

2.1 碎片的幾何分類與質(zhì)量分布模型

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

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

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

結(jié)合表1與式（1），可確定各類碎片的權(quán)重。例如假設(shè)碎片總質(zhì)量為200 kg，3個(gè)尺度取為0.1m，0.01 m和0.001 m，則各類碎片數(shù)量與質(zhì)量的關(guān)系如圖1所示。

2.2 碎片運(yùn)動(dòng)隨機(jī)性的模擬

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

圖1 各類碎片的數(shù)量與質(zhì)量關(guān)系Fig.1 Relationship between quantity and mass of various debris

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

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

3）對(duì)于具有N1個(gè)碎片的種類1，通過生成N1個(gè)隨機(jī)數(shù)R1，再按規(guī)則映射到某個(gè)方位角中去即可。同理進(jìn)行其它種類碎片的映射；

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

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

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

1）一次性生成Mφ個(gè)均勻隨機(jī)數(shù)，存放在數(shù)組Rφ中；

2）對(duì)Rφ中的元素進(jìn)行排序，記錄每個(gè)元素所對(duì)應(yīng)的序號(hào)，存放在數(shù)組Iφ中；

3）將飛行迎角均勻劃分成Mφ個(gè)組，Iφ即可根據(jù)序號(hào)映射到不同的飛行迎角；

4）同理，可對(duì)飛行測(cè)滑角進(jìn)行類似的映射；

5）分別從2個(gè)方向完成對(duì)Mα×Mβ碎片組的上述循環(huán)，即可確定各碎片組的飛行迎角和側(cè)滑角。

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

2.3 再入跨流域氣動(dòng)特性的當(dāng)?shù)鼗焖偎惴?/h3>

常用的數(shù)值方法如直接模擬蒙特卡羅（DSMC）方法、N-S方程解算器與求解玻爾茲曼（Boltzmann）模型方程、氣體動(dòng)理論統(tǒng)一算法（GKUA）等，由于計(jì)算量太大，不能滿足彈道計(jì)算中的氣動(dòng)特性快速獲取要求，本文采用了當(dāng)?shù)鼗焖偎惴ǎ?6］。

當(dāng)?shù)鼗焖偎惴ㄒ韵”怏w高超聲速繞流當(dāng)?shù)鼗嬖ㄓ?jì)算理論為基礎(chǔ)，將非規(guī)則解體物面劃分成若干塊小曲面，對(duì)于每個(gè)小曲面，選用一個(gè)小的平面三角形面元或四邊形面元來(lái)代替，利用面元逼近復(fù)雜解體物形。這樣，計(jì)算小曲面上的氣動(dòng)力就轉(zhuǎn)換成計(jì)算面元上的氣動(dòng)力，將這些面元上的氣動(dòng)力累加起來(lái)就可以得到整個(gè)非規(guī)則解體物的氣動(dòng)特性。當(dāng)面元分得極其細(xì)小時(shí)，由此引起的誤差就會(huì)變得很小。

當(dāng)?shù)孛嬖軞鈩?dòng)力系數(shù)依賴于來(lái)流和當(dāng)?shù)匦再|(zhì)，如當(dāng)?shù)赜?、?dāng)?shù)乇砻孀饔玫取?duì)于給定的入射角，面元上經(jīng)歸一化的壓力和摩擦力系數(shù)表達(dá)為式（2）：

式中，CP和Cf分別是過渡流區(qū)域的壓力和摩擦力系數(shù)，下標(biāo)FM和Cont分別表示自由分子流和連續(xù)流；FP，b和Ff，b分別是壓力和摩擦力橋函數(shù)，下標(biāo)b表示橋函數(shù)，依賴于獨(dú)立參數(shù)當(dāng)?shù)嘏瓟?shù)Kn、壁溫總溫比TW/T0和入射角θ。橋函數(shù)的意義在于使氣動(dòng)系數(shù)在自由分子流區(qū)和連續(xù)流區(qū)分別趨近于理論值，而在中間過渡流區(qū)實(shí)現(xiàn)光滑過渡。

分別使用氣體動(dòng)理論統(tǒng)一算法（GKUA）［23-24］、DSMC方法［25，27］或N-S方程解算器對(duì)殘骸碎片簡(jiǎn)化外形典型繞流狀態(tài)進(jìn)行精細(xì)計(jì)算，獲取典型狀態(tài)的精確氣動(dòng)數(shù)據(jù)。調(diào)試修正非規(guī)則碎片的橋函數(shù)關(guān)聯(lián)參數(shù)，使得當(dāng)?shù)鼗焖偎惴ǖ慕Y(jié)果在典型狀態(tài)與數(shù)值算法吻合后，則其它任意狀態(tài)的氣動(dòng)數(shù)據(jù)可由式（2）直接計(jì)算得出。

2.4 碎片運(yùn)動(dòng)彈道計(jì)算

2.2 節(jié)的隨機(jī)分組方法已經(jīng)考慮了飛行姿態(tài)的影響，則碎片運(yùn)動(dòng)的彈道計(jì)算中只需計(jì)算三自由度彈道。矢量形式的彈道方程如式（3）所示：

式中，n為碎片的瞬時(shí)質(zhì)量；V為碎片的瞬時(shí)速度矢量；R為碎片空氣動(dòng)力；G為地球引力，G=n▽U；U為地球的引力勢(shì)；Fe為離心慣性力，F(xiàn)e=n rΩ×( Ω×r)；Ω為地球的自轉(zhuǎn)角速度矢量，r是地心慣性坐標(biāo)系中質(zhì)心的矢徑；Fco為科氏慣性力，F(xiàn)co=2 n( Ω×V)。

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

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

式中，h為時(shí)間步長(zhǎng)，un+1為待求的下一時(shí)間步函數(shù)值，un為當(dāng)前時(shí)刻函數(shù)值，fn，fn-1，…，fn-5為當(dāng)前時(shí)刻及當(dāng)前時(shí)刻之前1～5時(shí)間步的右函數(shù)值。在計(jì)算啟動(dòng)時(shí)刻，由于時(shí)間步結(jié)果不足，無(wú)法直接采用此數(shù)值格式，需要采用龍格-庫(kù)塔等其它數(shù)值方法來(lái)計(jì)算初始的幾個(gè)時(shí)間步。

3 碎片散布范圍計(jì)算分析

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

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

圖5繪出特征尺度為0.01 m的立方體狀碎片飛行軌跡與碎片云質(zhì)心彈道的關(guān)系，其中藍(lán)色點(diǎn)劃線為碎片云質(zhì)心運(yùn)動(dòng)彈道，彩色圓圈為該組碎片不同時(shí)刻的空間位置，其顏色代表飛行馬赫數(shù)。由圖中可看出碎片飛行軌跡與質(zhì)心運(yùn)動(dòng)彈道的一致性很好。

圖3 解體碎片在地面的質(zhì)量分布（φ=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立方體狀碎片飛行軌跡與質(zhì)心彈道Fig.5 Flight path of cube-like debris w ith the characteristic scale of 0.01 m

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

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

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

4 結(jié)論

1）對(duì)碎片按幾何形狀和空間方位角的分組方案，可以實(shí)現(xiàn)對(duì)碎片外形非規(guī)則性和運(yùn)動(dòng)隨機(jī)性的統(tǒng)計(jì)模擬，進(jìn)而實(shí)現(xiàn)對(duì)航天器再入解體生成的大量非規(guī)則碎片的整體分析。

2）在一次解體過程生成的碎片，單個(gè)碎片質(zhì)量越大的碎片組，其碎片數(shù)量越少，航程越遠(yuǎn)。碎片在地面的散布整體呈細(xì)長(zhǎng)條狀，縱向可達(dá)數(shù)千公里，橫向僅約百余公里。

3）碎片飛行軌跡與三自由度質(zhì)心運(yùn)動(dòng)彈道一致。起始點(diǎn)彈道傾角對(duì)落區(qū)散布范圍影響很大，在實(shí)際應(yīng)用中應(yīng)確保起始點(diǎn)彈道傾角的高精度獲取。

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