胡銳鋒, 龔自正, 吳子牛

（1.西安電子科技大學(xué) 機(jī)電工程學(xué)院，西安 710071；2．北京衛(wèi)星環(huán)境工程研究所 可靠性與環(huán)境工程技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100094；3．清華大學(xué) 航天航空學(xué)院，北京 100084）

0 引言

無控航天器再入地球大氣層時(shí)，在強(qiáng)烈的氣動(dòng)力和氣動(dòng)熱載荷作用下，會(huì)發(fā)生解體并產(chǎn)生大量碎片；其中一部分碎片在空中被完全燒毀，其余存活下來的碎片可能會(huì)對(duì)地面建筑、人群和生態(tài)系統(tǒng)構(gòu)成很大威脅。

人類有記錄的再入事件之一發(fā)生于1997年1月22日：美國(guó)俄克拉荷馬州圖爾薩市的一位女士在日出前陪朋友散步，她突然看到一道火球從北向南穿過天空，之后一塊碎片劃過這位女士的左肩并撞擊到地面。這塊碎片來自于1996年4月24日發(fā)射的德爾塔火箭的第二級(jí)，與此同時(shí)德爾塔火箭第二級(jí)上的一個(gè)更大的不銹鋼推進(jìn)劑儲(chǔ)箱落在德克薩斯州喬治城附近的一個(gè)農(nóng)場(chǎng)內(nèi)，此外還有一塊鈦制圓球落在了塞金鎮(zhèn)外。美國(guó)1997年發(fā)射的“銥星33”通信衛(wèi)星于2009年2月11日與俄羅斯 1993年發(fā)射的一顆衛(wèi)星（已報(bào)廢）在西伯利亞上空相撞，這是歷史上首次衛(wèi)星相撞事故，產(chǎn)生了大量碎片（約500～600 塊），部分碎片再入大氣層后降落在美國(guó)德克薩斯州等地。2011年9月，美國(guó)一顆名為“高層大氣研究衛(wèi)星”（UARS）的再入事件引起了全世界的關(guān)注和恐慌：NASA 相關(guān)專家預(yù)測(cè)UARS 墜落地球后，其碎片覆蓋的區(qū)域?qū)⑦_(dá)到800 km 左右，而在重返地球前大約2 h，NASA才能將墜落區(qū)鎖定在大約1萬km之內(nèi)；最終UARS于2011年9月24日墜落在南太平洋中，沒有對(duì)人類產(chǎn)生威脅。

提前預(yù)報(bào)航天器及碎片再入軌跡和殘存性是預(yù)防和減輕其危害的重要措施。近20年來，國(guó)內(nèi)外在無控航天器和空間碎片再入預(yù)測(cè)方法和地面危險(xiǎn)評(píng)估研究方面取得了很大進(jìn)展，并開發(fā)了若干工程預(yù)測(cè)軟件，如 NASA 的 DAS（Debris Assessment Software ）[1-2]和 ORSAT （ Object oRiented Surveillance Analysis Tool）[3]軟件，ESA的SCARAB（SpaceCraft Atmospheric Re-entry and Aerothermal Breakup）[4]軟件，以及國(guó)內(nèi)由清華大學(xué)和中國(guó)空間技術(shù)研究院共同開發(fā)的 DRAPS（Debris Reentry and Ablation Prediction System）[5-6]軟件等。

本文對(duì)航天器與空間碎片再入的工程預(yù)測(cè)方法研究現(xiàn)狀進(jìn)行跟蹤和總結(jié)。

1 預(yù)測(cè)模型

1.1 航天器模型

Lips 和Fritsche[7]將航天器再入工程預(yù)測(cè)方法分為兩類：基于集總參數(shù)模型的面向物體法和基于復(fù)雜表面網(wǎng)格模型的面向航天器法。

面向物體法中，航天器以及解體后產(chǎn)生的碎片外形是簡(jiǎn)單的基本幾何模型，并且是集總參數(shù)類型，即其幾何信息可用若干特征參數(shù)描述。工程軟件DAS 和ORSAT 提供了圓球、圓柱、平板、六面體4 種基本幾何模型[1-3]，DRAPS 則提供了15種基本幾何模型[5-6]。面向物體法可通過定義物體之間的包含邏輯關(guān)系來模擬航天器及部件之間的多層嵌套結(jié)構(gòu)。

在面向航天器法中，航天器最底層部件由基本幾何模型構(gòu)成，進(jìn)而組成子系統(tǒng)，最后集成為整個(gè)航天器?；編缀文Ｐ筒皇羌倕?shù)類型，而是劃分表面網(wǎng)格的三維模型。相比于面向物體法，面向航天器法不僅在幾何建模水平上更接近實(shí)際情況，并且還能夠準(zhǔn)確模擬航天器的質(zhì)量和慣量特性。

1.2 動(dòng)力學(xué)模型

航天器幾何建模方式?jīng)Q定所采用的動(dòng)力學(xué)模型。

面向物體法由于采用集總參數(shù)的幾何模型，所以采用僅確定質(zhì)心運(yùn)動(dòng)的三自由度彈道模型，而通過預(yù)先給定的運(yùn)動(dòng)模式考慮物體姿態(tài)的變化對(duì)所受氣動(dòng)力的影響。運(yùn)動(dòng)模式包括靜態(tài)模式和動(dòng)態(tài)模式，其中靜態(tài)模式中的物體再入姿態(tài)保持不變，而動(dòng)態(tài)模式中的物體可繞定軸翻滾或無序翻滾，但僅考慮翻滾時(shí)的平均氣動(dòng)力。

與面向物體法不同，面向航天器法獲得了沿物體表面的壓力分布從而得到三個(gè)方向的氣動(dòng)力和力矩，因此采用六自由度動(dòng)力學(xué)方程，即聯(lián)立求解質(zhì)心運(yùn)動(dòng)方程和姿態(tài)運(yùn)動(dòng)歐拉方程，計(jì)算中一般采用四元數(shù)代替歐拉角以避免求解時(shí)可能出現(xiàn)的奇異問題。

1.3 氣動(dòng)模型

錢學(xué)森[8]最早利用稀薄程度對(duì)流動(dòng)進(jìn)行分區(qū)。隨著克努森數(shù)Kn（Kn=λ/L，λ是分子平均自由程，L是特征長(zhǎng)度）的增加，流動(dòng)可粗略分為連續(xù)流、過渡流和自由分子流三個(gè)流區(qū)，而連續(xù)流區(qū)外的其他區(qū)域一般稱為稀薄流區(qū)。通常將Kn＜0.001 劃分為連續(xù)流區(qū)，0.001＜Kn＜10 為過渡流區(qū)，Kn＞10為自由分子流區(qū)。航天器的再入過程經(jīng)歷三個(gè)流區(qū)環(huán)境，即從自由分子流、過渡流到連續(xù)流。

在自由分子流區(qū)，流場(chǎng)基本方程是無碰撞的Boltzmann 方程，在Maxwell 平衡氣體分布假設(shè)下，可得到任意外形物體表面壓力、摩擦力和熱流分布的解析公式，其具體形式可參考文獻(xiàn)[9-11]。

在連續(xù)流區(qū)，航天器迎風(fēng)表面壓力分布可通過牛頓流模型、切楔法、切錐法、激波/膨脹波法等近似理論或工程方法計(jì)算[12]。其中牛頓流模型方法又包括經(jīng)典牛頓流模型、修正牛頓流模型等，適用于具有大撞擊角的鈍頭外形。切楔法、切錐法和激波/膨脹波法適用于撞擊角較小的楔面或錐面外形。此外，還有將牛頓流模型與切楔/切錐法組合使用的Dahlem-Buck 法，即在大撞擊角時(shí)采用牛頓流模型，在小撞擊角時(shí)采用切楔法或切錐法。對(duì)超聲速/高超聲速連續(xù)流區(qū)表面壓力分布計(jì)算還可采用激波層方法或CFD 方法（無黏或有黏），但這些方法計(jì)算量太大，并不適用于航天器再入燒蝕和解體的工程預(yù)測(cè)。

由于物體形狀的限制，面向物體法的氣動(dòng)模型給出的是整體阻力系數(shù)和平均熱流大小，而不是表面壓力和熱流分布。在連續(xù)流區(qū)，面向物體法一般利用典型外形（如圓球或圓柱）駐點(diǎn)熱流以及理論或?qū)嶒?yàn)得到的表面熱流分布形式，積分平均后得到物體表面平均熱流。高超聲速連續(xù)流區(qū)駐點(diǎn)熱流估算有許多成熟的近似理論或工程方法，如Lees 公式[13]、Fay-Riddell 公式[14]、Detra-Kemp-Riddell 公式[15]、Cohen 公式[16]、Sutton-Graves 公式[17]等。

在目前已有的面向航天器法應(yīng)用中，表面熱流分布通過與駐點(diǎn)參數(shù)以及撞擊角有關(guān)的近似經(jīng)驗(yàn)公式計(jì)算[4,18]。利用航天器表面網(wǎng)格，還可采用可信度更高的基于跟蹤表面流線的邊界層工程方法計(jì)算摩擦阻力和熱流分布，如Dejarnette 等[19-23]提出和發(fā)展的軸對(duì)稱比擬法。

在過渡流區(qū)，一方面，由于連續(xù)介質(zhì)假設(shè)逐漸失效從而無法采用上述連續(xù)流區(qū)氣動(dòng)估算方法；另一方面，分子間的碰撞不可忽略，因此也不能通過簡(jiǎn)化Boltzmann 方程得到氣動(dòng)力和氣動(dòng)熱的估算公式。工程中一般采用橋函數(shù)方法（Bridging Method），再利用連續(xù)流區(qū)和自由分子流區(qū)的結(jié)果插值得到過渡流區(qū)航天器的氣動(dòng)力和氣動(dòng)熱。常用的氣動(dòng)力橋函數(shù)有正弦平方橋函數(shù)[24]和誤差橋函數(shù)[25]。徐姍姝[26]研究了馬赫數(shù)對(duì)氣動(dòng)力橋函數(shù)適應(yīng)性的影響：在高馬赫數(shù)（Ma＞10）情況下，兩種橋函數(shù)都與DSMC 結(jié)果符合得較好；而在中低馬赫數(shù)（Ma＜5）情況下，正弦平方橋函數(shù)偏離DSMC 結(jié)果較大，而誤差橋函數(shù)與DSMC 結(jié)果仍符合得很好。

氣動(dòng)熱橋函數(shù)有Matting 橋函數(shù)[27]和Nomura橋函數(shù)[28]。Swaminathan 等[29]給出Matting 橋函數(shù)和Nomura 橋函數(shù)對(duì)某高超聲速導(dǎo)彈過渡區(qū)熱流分布預(yù)測(cè)結(jié)果與DSMC 計(jì)算結(jié)果的對(duì)比，結(jié)論是兩種橋函數(shù)均在導(dǎo)彈表面部分區(qū)域與DSMC 結(jié)果符合較好，但在某些區(qū)域與DSMC 結(jié)果存在明顯偏差（2 倍左右）。

1.4 燒蝕和解體模型

在氣動(dòng)加熱和輻射換熱作用下，再入過程中存在從航天器表面向內(nèi)部傳遞的凈熱流，該熱流隨著高度的降低而增大，在大約70～80 km 附近達(dá)到最大。航天器結(jié)構(gòu)在長(zhǎng)時(shí)強(qiáng)熱作用下材料溫度顯著升高，一方面溫度有可能達(dá)到材料熔點(diǎn)并使材料發(fā)生熔化，即發(fā)生燒蝕；另一方面在加熱作用下結(jié)構(gòu)強(qiáng)度降低，某些危險(xiǎn)截面可能由于局部應(yīng)力超過極限值而引起結(jié)構(gòu)破壞，即產(chǎn)生解體。燒蝕與解體均為在強(qiáng)氣動(dòng)熱效應(yīng)下航天器結(jié)構(gòu)的破壞現(xiàn)象。

在面向物體法中，為了滿足對(duì)物體幾何形狀的限制要求，一般采用零維或一維熱傳導(dǎo)模型模擬燒蝕過程。零維模型是集總參數(shù)模型，假設(shè)物體內(nèi)部溫度均勻分布，不考慮熱傳導(dǎo)效應(yīng)，當(dāng)再入過程中的吸熱總量達(dá)到物體完全熔化所需的熱量后認(rèn)為物體完全燒蝕（質(zhì)量置零）。一維模型考慮沿物體厚度方向的熱傳導(dǎo)，將物體沿厚度方向劃分一定數(shù)量的層狀單元，通過建立每一層單元的熱傳導(dǎo)方程計(jì)算溫度分布并模擬徑向燒蝕過程。

面向物體法采用簡(jiǎn)單解體準(zhǔn)則，包括高度準(zhǔn)則、溫度準(zhǔn)則等。高度準(zhǔn)則是基于對(duì)歷史上再入事件的統(tǒng)計(jì)分析直接指定解體高度，即認(rèn)為當(dāng)航天器降到一定高度時(shí)則發(fā)生解體。IADC 空間碎片減緩指南[30]指出，航天器通常在75 km 高度處解體，實(shí)際觀測(cè)發(fā)現(xiàn)航天器解體的典型高度為(80±10) km，ORASAT 軟件中采用的解體高度為78 km。溫度準(zhǔn)則根據(jù)物體表面溫度進(jìn)行解體判斷，即認(rèn)為當(dāng)溫度升高到某一臨界溫度時(shí)航天器會(huì)發(fā)生解體。臨界溫度與物體的材料特性（如熔點(diǎn)）密切相關(guān)，不同部件的臨界溫度不同。對(duì)薄壁物體或物體壁面較厚并且再入角較大的情況，Baker 等[31]基于輻射平衡假設(shè)或考慮物體升高到臨界溫度所吸收的熱量給出再入解體的失效溫度判據(jù)表達(dá)式，其本質(zhì)是臨界溫度為材料熔化溫度的溫度準(zhǔn)則?，F(xiàn)有的面向物體法工程軟件尚未能考慮基于結(jié)構(gòu)破壞的解體準(zhǔn)則。

面向航天器法假定結(jié)構(gòu)為均勻厚度的殼體，可以采用一維或三維熱傳導(dǎo)模型模擬燒蝕過程。一維燒蝕模型只考慮沿厚度方向的熱傳導(dǎo)，三維熱傳導(dǎo)模型則考慮沿周向網(wǎng)格單元與當(dāng)前網(wǎng)格單元存在溫差時(shí)所引起的熱傳導(dǎo)效應(yīng)。與面向物體法類似，當(dāng)某個(gè)單元的吸熱總量達(dá)到完全熔化所需的熱量后認(rèn)為完全燒蝕，質(zhì)量置零并改變局部幾何外形?，F(xiàn)有面向航天器法工程軟件能夠通過采用監(jiān)視危險(xiǎn)截面局部應(yīng)力[4]或自定義物理量[18]的方式建立解體準(zhǔn)則，一旦結(jié)構(gòu)發(fā)生解體成兩個(gè)或多個(gè)碎片，將分別計(jì)算碎片的再入過程。

2 國(guó)外軟件

2.1 DAS

面向物體法工程軟件DAS[1-2]由美國(guó)Lockheed公司于1998年開發(fā)，用于滿足NASA 對(duì)航天器再入風(fēng)險(xiǎn)快速分析的需求。

DAS 中提供的航天器模型包括圓球、圓柱、立方體和平板，但只能是實(shí)心而不能是空心，需定義的參數(shù)包括物體形狀、尺寸、質(zhì)量和材料特性。DAS 的再入動(dòng)力學(xué)模型采用三自由度彈道模型，燒蝕模型采用集總參數(shù)的零維熱傳導(dǎo)模型，因此不能模擬部分燒蝕效應(yīng)以及解體。

DAS 的主要計(jì)算結(jié)果是再入物體完全燒蝕的消亡高度或?qū)Φ孛娴臍娣e，以表格形式顯示。DAS 可用于航天器再入風(fēng)險(xiǎn)的初步分析，如果對(duì)其預(yù)測(cè)結(jié)果存在疑問，應(yīng)采用更為精確的方法和工程軟件作進(jìn)一步計(jì)算和分析。

2.2 ORSAT

ORSAT 由NASA 約翰遜空間中心開發(fā)，是目前NASA 用于空間碎片再入分析的主要工具，也是面向物體法的工程軟件。

與DAS 類似，ORSAT 的航天器模型也包括圓球、圓柱、立方體和平板，但可以模擬空心物體并提供總共10 種運(yùn)動(dòng)模式。ORSAT 的動(dòng)力學(xué)模型采用三自由度彈道模型，燒蝕模型提供了零維熱傳導(dǎo)模型或考慮徑向溫度分布的一維熱傳導(dǎo)模型。ORSAT 采用高度準(zhǔn)則判斷是否解體，解體高度給定為78 km。

在連續(xù)流區(qū)，ORSAT 采用牛頓流模型計(jì)算物體的阻力系數(shù)，采用Detra-Kemp-Riddell 公式或Fay-Riddell 公式計(jì)算駐點(diǎn)熱流從而得到平均熱流。在過渡流區(qū)，ORSAT 采用正弦平方橋函數(shù)計(jì)算阻力系數(shù)，采用對(duì)數(shù)或Matting 橋函數(shù)計(jì)算氣動(dòng)熱。此外，ORSAT 中還考慮了物體表面氧化熱流效應(yīng)。

ORSAT 已被NASA 用于許多航天器再入事件的預(yù)測(cè)分析以及標(biāo)準(zhǔn)算例驗(yàn)證中，比如用于與ESA的SCARAB 對(duì)比的空心圓球、德爾塔火箭第二級(jí)、鋇燃料棒、SPARTAN 航天器以及UARS 航天器 等[3,32-34]。

2.3 SCARAB

ESA 的SCARAB 是目前唯一的面向航天器法工程軟件，由德國(guó)宇航院 HTG（Hypersonic Technology G?ttingen）中心于1995年開始開發(fā)。

如前所述，SCARAB 中航天器幾何模型采用劃分表面網(wǎng)格的三維模型，其對(duì)BeppoSAX 衛(wèi)星的建模見圖1[35]。由于氣動(dòng)模型能夠提供航天器表面的壓力分布，因此SCARAB 軟件采用六自由度動(dòng)力學(xué)模型，直接求解三個(gè)方向的平動(dòng)和轉(zhuǎn)動(dòng)。該軟件具有圖形化操作界面，易于操作并且能夠直觀顯示計(jì)算結(jié)果。采用一維或二維熱傳導(dǎo)模型得到結(jié)構(gòu)內(nèi)部的溫度分布，并通過監(jiān)視預(yù)先定義截面的局部 等效應(yīng)力來判斷結(jié)構(gòu)在危險(xiǎn)點(diǎn) 是否發(fā)生解體。

圖1 SCARAB 對(duì)BeppoSAX 衛(wèi)星的幾何建模Fig.1 BeppoSAX satellite model by SCARAB

SCARAB 已被用于許多大型復(fù)雜航天器的再入預(yù)測(cè)問題，包括ATV、德國(guó)X 射線衛(wèi)星ROSAT、阿里安娜5 火箭、意大利BeppoSAX 衛(wèi)星以及Mir空間站等[35-37]。

2.4 其他軟件

除了上述介紹的軟件外，目前公開報(bào)道的航天器再入預(yù)測(cè)軟件還有ESA 的DRAMA（Debris Risk Assessment and Mitigation Analysis）[38]、美國(guó)宇航公司（The Aerospace Corporation）的 AHaB（Atmospheric Heating and Breakup tool）[39]以及韓國(guó)的SAPAR（Survivability Analysis Program for Atmospheric Reentry）[40]等。

DRAMA 是歐洲空間碎片減緩標(biāo)準(zhǔn)（European Space Debris Mitigation Standard，ESDMS）的支持工具，其再入預(yù)測(cè)分析模塊名為 SESAM（Spacecraft Entry Survival Analysis Module）。SESAM 采用與ORSAT 類似的模型和方法，但只采用零維熱傳導(dǎo)模型。AHaB 目前沒有可獲取的詳細(xì)介紹資料。SAPAR 由韓國(guó)首爾國(guó)立大學(xué)開發(fā)，基于ORSAT 的模型和方法發(fā)展而來，主要區(qū)別是考慮了空心圓柱物體底部的輻射熱損失效應(yīng)，從而與德爾塔火箭第二級(jí)圓柱儲(chǔ)箱再入觀測(cè)結(jié)果符合得更好。

總的來說，目前基于面向物體法的航天器再入預(yù)測(cè)工程軟件占絕大多數(shù)，而基于面向航天器法的僅有SCARAB。其原因一方面在于面向物體法軟件易于實(shí)現(xiàn)且計(jì)算量小，而面向航天器法軟件難于實(shí)現(xiàn)且計(jì)算量大；另一方面建模更加精細(xì)的面向航天器法并不一定比面向物體法更準(zhǔn)確，因?yàn)榍罢咚捎玫念A(yù)測(cè)模型仍然經(jīng)過了很大簡(jiǎn)化，與實(shí)際情況相比存在差別。因此，計(jì)算量較小的面向物體法軟件在工程應(yīng)用中更受歡迎，而且通過不確定性分析能夠在一定程度上給出模型近似所帶來的影響，對(duì)工程實(shí)際也有重要指導(dǎo)意義。

3 國(guó)內(nèi)軟件

在國(guó)家國(guó)防科技工業(yè)局空間碎片專項(xiàng)資助下，清華大學(xué)與中國(guó)空間技術(shù)研究院合作開發(fā)了目前國(guó)內(nèi)唯一一款面向物體法工程軟件DRAPS。

3.1 DRAPS 簡(jiǎn)介

DRAPS 提供15 種航天器模型、51 種運(yùn)動(dòng)模式及其相應(yīng)氣動(dòng)模型，從而能夠?qū)Ω鄰?fù)雜外形的航天器再入進(jìn)行模擬，其部分基本幾何模型如圖2所示。

圖2 DRAPS 中提供的基本幾何模型Fig.2 The basic geometry models in DRAPS

DRAPS 的動(dòng)力學(xué)模型采用三自由度彈道模型，燒蝕模型也采用零維熱傳導(dǎo)模型或考慮徑向溫度分布的一維熱傳導(dǎo)模型。其所采用的氣動(dòng)模型與ORSAT 類似，并通過CFD 和DSMC 模擬對(duì)部分氣動(dòng)模型進(jìn)行了驗(yàn)證[5]。在解體模型方面，DRAPS除了高度準(zhǔn)則外還提供溫度準(zhǔn)則、燒蝕準(zhǔn)則（燒蝕到一定程度發(fā)生解體）和綜合準(zhǔn)則（達(dá)到上述任何一個(gè)解體條件就認(rèn)為解體）。此外，DRAPS 中建立了多層嵌套碎片結(jié)構(gòu)模型，能夠模擬航天器部件間的包容關(guān)系。

DRAPS 的地面風(fēng)險(xiǎn)評(píng)估模型與ORSAT 相同，人口密度數(shù)據(jù)采用Columbia大學(xué)的2005年度地球人口密度分布數(shù)據(jù)，大氣模型采用國(guó)際標(biāo)準(zhǔn)大氣模型。

DRAPS 軟件采用C++語言編寫，GUI 界面通過Borland C++ 6.0 建立，其界面包括衛(wèi)星碎片定義、初始狀態(tài)、解體模型設(shè)置、軌道衰減仿真、再入結(jié)果與地面風(fēng)險(xiǎn)評(píng)估等。

3.2 典型仿真結(jié)果

目前DRAPS 已被用于多個(gè)簡(jiǎn)單或復(fù)雜外形航天器再入算例驗(yàn)證與測(cè)試[6]，包括標(biāo)準(zhǔn)空心圓球、德爾塔火箭第二級(jí)、SPARTAN 航天器、Intelsat VII-A 衛(wèi)星、某航天器部分艙段等，通過對(duì)比DRAPS 與ORSAT 計(jì)算結(jié)果發(fā)現(xiàn)二者精度相近。

下面介紹采用DRAPS 對(duì)某典型航天器部分艙段再入仿真分析結(jié)果。該航天器內(nèi)部分系統(tǒng)組成包括推進(jìn)分系統(tǒng)、電源分系統(tǒng)、環(huán)控生保分系統(tǒng)、測(cè)控通信分系統(tǒng)、熱控分系統(tǒng)和數(shù)管分系統(tǒng)，主要部件包括燃燒室、擴(kuò)張噴嘴、推進(jìn)劑儲(chǔ)箱、氣瓶、管路、閥門、輻射器、天線、配電器、數(shù)傳設(shè)備、太陽電池板等共38 種。圖3和圖4給出了DRAPS對(duì)該航天器仿真得到的各部件再入軌跡、消亡高度以及落點(diǎn)位置。

圖3 DRAPS 對(duì)某航天器部分艙段再入仿真結(jié)果Fig.3 Reentry prediction results of certain spacecraft by DRAPS

圖4 DRAPS 對(duì)存活碎片落點(diǎn)仿真結(jié)果Fig.4 Landing locations of survived debris by DRAPS

3.3 下一步發(fā)展方向

盡管DRAPS 已具備航天器與空間碎片再入預(yù)測(cè)能力并與ORSAT 的預(yù)測(cè)精度相近，但還需在以下方面進(jìn)一步改進(jìn)與完善：

1）考慮更復(fù)雜的航天器外形。在航天器解體前采用面向航天器法模擬其再入過程，而解體后由于產(chǎn)生大量碎片，采用計(jì)算速度更快的面向物體法顯然更合適。

2）采用更精確的氣動(dòng)模型。所采用氣動(dòng)力和氣動(dòng)熱模型十分古老并經(jīng)過很大簡(jiǎn)化，對(duì)其精度和適用性的驗(yàn)證還很不夠；而隨著計(jì)算手段的快速發(fā)展，利用精細(xì)流場(chǎng)仿真結(jié)果對(duì)氣動(dòng)模型進(jìn)行校驗(yàn)和修正已成為可能。

3）設(shè)計(jì)更友好的用戶界面?？刹捎媚壳傲餍械能浖缑骈_發(fā)工具（如QT）對(duì)用戶界面重新設(shè)計(jì)和美化，增加三維視景顯示功能，虛擬再現(xiàn)航天器與空間碎片再入過程。

此外，還可采用更新的大氣模型、地球人口密度分布模型等，以及建立航天器與空間碎片再入數(shù)據(jù)庫(kù)等。

4 其他問題分析

4.1 碎片存活性

再入分析需要得到的一個(gè)重要結(jié)論是碎片能否存活至地球表面，這與再入物體的材料、幾何形狀、尺寸、質(zhì)量等因素緊密相關(guān)。采用解析手段研究簡(jiǎn)單外形再入物體能夠得到關(guān)于碎片存活規(guī)律的定性認(rèn)識(shí)。

Koppenwallner 等[41]采用Allen-Eggers 彈道方程[42]以及Lees 駐點(diǎn)熱流公式[13]研究實(shí)心圓球、圓柱和圓盤再入存活性問題。通過研究發(fā)現(xiàn)，碎片存活至地面源于兩種機(jī)理：第一種是輻射，即當(dāng)碎片表面熱輻射能力足夠強(qiáng)時(shí)，其溫度始終低于熔點(diǎn)從而可以存活至地面；另一種是熱容效應(yīng)，即當(dāng)物體容納熱量能力足夠強(qiáng)時(shí)，高速氣流的氣動(dòng)加熱并不能使其完全燒蝕或熔化。Koppenwallner 等[41]給出不同材料以及不同外形碎片存活性隨尺寸的變化，如圖5所示。從圖中可以看到，碎片在一定尺寸范圍內(nèi)才能存活。當(dāng)d＜dmin，且碎片的面質(zhì)比（A/m∝1/d）較小，再入速度較低，并且具有較強(qiáng)的輻射能力，從而能夠存活至地面。當(dāng)d＞dmax時(shí)，碎片的熱容能力很強(qiáng)，從而能夠存活至地面。圖5上顯示了實(shí)心圓球再入存活性隨其直徑及材料的變化關(guān)系。從中可以看出，鈦質(zhì)圓球具有最高的存活可能性，然后依次是不銹鋼、鎳鉻合金和銅的，而具有低熔點(diǎn)和較差輻射特性的鋁質(zhì)圓球具有最低的再入存活性。圖5(b)給出的是不同形狀（圓球、圓柱、圓盤）鈦金屬碎片的再入存活性與其直徑的關(guān)系。由于具有較小的面質(zhì)比，圓盤的再入存活性最高，而在尺寸較大時(shí)圓球更容易存活至地面。Fritsche 等[43]在此基礎(chǔ)上進(jìn)一步研究了空心碎片（圓球、圓柱和立方體）再入的存活性問題，并得到與碎片尺寸和壁厚有關(guān)的存活和消亡區(qū)域（圖6）。

圖5 碎片再入存活性與幾何尺寸、材料的關(guān)系 Fig.5 Debris survivability versus size and material

圖6 空心碎片存活與消亡區(qū)域示意圖Fig.6 Schematic diagrams for survival and demise regimes of hollow debris

通過上述研究可知，只有尺寸很小或很大的碎片能夠存活至地面。小尺寸的存活碎片直徑一般在1 mm 甚至更小，其對(duì)地面的影響可以忽略不計(jì)。工程中更關(guān)心的是大尺寸的再入存活碎片及其對(duì)地面的威脅程度。目前研究碎片存活性所考慮的因素只有碎片幾何尺寸和壁厚，進(jìn)一步可研究解體高度、再入角、方位角、運(yùn)動(dòng)模式對(duì)不同外形和材料碎片的存活性的影響，并形成數(shù)據(jù)庫(kù)。

4.2 預(yù)測(cè)不確定性

盡管上述航天器再入預(yù)測(cè)方法均為確定性方 法，即給定相關(guān)參數(shù)和初始條件便能得到唯一確定的解，但是由于所采用的數(shù)學(xué)模型經(jīng)過很多簡(jiǎn)化并包含許多不確定性因素，因此最終給出具有一定概率分布形式的預(yù)測(cè)結(jié)果應(yīng)更為合理。

對(duì)于航天器在軌運(yùn)動(dòng)的長(zhǎng)期預(yù)測(cè)，不確定性來源于軌跡追蹤、航天器姿態(tài)、大氣密度、太陽活動(dòng)以及地磁場(chǎng)等。對(duì)于航天器從進(jìn)入大氣層至地球表面的再入預(yù)測(cè)，不確定性來源于再入初始條件（高度、速度、再入角等）、大氣模型（大氣參數(shù)、風(fēng)場(chǎng)等）、氣動(dòng)模型（阻力系數(shù)、升力系數(shù)、氣動(dòng)熱等）、解體高度等[5]，上述誤差源的綜合影響可能導(dǎo)致很大的再入碎片落點(diǎn)散布預(yù)測(cè)范圍。

蒙特卡羅模擬法是考慮各個(gè)誤差源對(duì)碎片落點(diǎn)影響的常用方法，如Rao 和Woeste[44]采用該法研究了衛(wèi)星碎片的落點(diǎn)散布問題，并且該法已在返回艙再入落點(diǎn)散布估計(jì)中得到大量應(yīng)用[45-48]。DRAPS 中提供對(duì)單個(gè)碎片再入落點(diǎn)散布范圍進(jìn)行計(jì)算的蒙特卡羅模擬功能[5]，并能給出落點(diǎn)散布橢圓形區(qū)域的長(zhǎng)軸和短軸的3-σ值。

蒙特卡羅法的缺點(diǎn)是計(jì)算量太大，因此近年來人們致力于研究能夠減小計(jì)算量的不確定性量化方法。Purcell 和 Barbery[49]采用溫度協(xié)方差法（Temperature Covariance Method）研究了大氣參數(shù)不確定性對(duì)再入落點(diǎn)散布的影響，獲得與蒙特卡羅法相當(dāng)?shù)慕Y(jié)果但計(jì)算量減小很多。Frank 等[50]采用一種概率模型方法研究了航天器再入解體不確定性對(duì)落點(diǎn)散布的影響。Reyhanoglu 和Alvarado[51]采用協(xié)方差傳播法（Covariance Propagation Method）來估計(jì)存在解體時(shí)碎片再入落點(diǎn)散布范圍。此外，近年來得到大量研究和應(yīng)用的不確定性量化方法還包括多項(xiàng)式混沌法（Polynomial Chaos）[52-55]、響應(yīng)面法（Response Surface Method）[56-57]、隨機(jī)劉維爾方程（Stochastic Liouville Equation）[58]等。進(jìn)一步研究可將上述方法應(yīng)用到碎片再入落點(diǎn)散布的不確定性分析中來。

4.3 氣動(dòng)模型精度

無論面向物體法還是面向航天器法，目前所采用的氣動(dòng)力和氣動(dòng)熱模型均基于大量假設(shè)和簡(jiǎn)化，并且面向物體法中對(duì)碎片運(yùn)動(dòng)模式的人為假定為氣動(dòng)力和氣動(dòng)熱估算結(jié)果的可靠性帶來很多質(zhì)疑。因此，為獲得氣動(dòng)估算的不確定性范圍，對(duì)各種情況下、各種外形物體的氣動(dòng)力和氣動(dòng)熱模型進(jìn)行驗(yàn)證與校正十分必要。

本文作者利用CFD 和DSMC 方法對(duì)圓柱在自由分子流、過渡流和連續(xù)流中的氣動(dòng)力與氣動(dòng)熱模型進(jìn)行了初步驗(yàn)證研究[5]，圖7是圓柱再入至不同高度處的流場(chǎng)數(shù)值計(jì)算結(jié)果并用于氣動(dòng)模型驗(yàn)證，其中H=30 km 流場(chǎng)通過CFD 方法計(jì)算得到，而H=110 km 和150 km 流場(chǎng)通過DSMC 方法計(jì)算得到。鈍錐氣動(dòng)熱估算結(jié)果見圖8[59]。本文還研究了基于表面網(wǎng)格的高超聲速氣動(dòng)力/熱快速估算方法并開發(fā)了相應(yīng)軟件，對(duì)于典型外形，氣動(dòng)力和表面氣動(dòng)熱的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果以及CFD 結(jié)果吻合很好。進(jìn)一步還需對(duì)各外形物體在不同馬赫數(shù)、雷諾數(shù)、克努森數(shù)和姿態(tài)角條件下的氣動(dòng)模型進(jìn)行更全面的驗(yàn)證研究。

圖7 圓柱再入流場(chǎng)計(jì)算結(jié)果：(a)～(c) H=30、110、150 km 的流場(chǎng)壓力云圖；(d)～(f) H=30、110、150 km 的流場(chǎng) 溫度云圖Fig.7 Numerical simulations of flow field for reentry cylinder∶(a)-(c) Pressure contours at H=30 km, 110 km and 150 km;(d)-(f) Temperature contours at H=30 km, 110 km and 150 km

圖8 鈍錐氣動(dòng)熱估算結(jié)果Fig.8 Aerodynamic heating estimation of a blunt cone

此外，在過去的半個(gè)世紀(jì)里，研究者們利用風(fēng)洞實(shí)驗(yàn)、激波管實(shí)驗(yàn)或流場(chǎng)高精度數(shù)值模擬手段，對(duì)圓球、圓柱等簡(jiǎn)單外形物體在高速連續(xù)流或稀薄流中的氣動(dòng)問題開展了大量研究，并得到許多估算阻力系數(shù)的近似理論或經(jīng)驗(yàn)關(guān)聯(lián)式[60-66]?？沙浞掷靡延醒芯繑?shù)據(jù)來驗(yàn)證現(xiàn)有碎片氣動(dòng)模型。

5 結(jié)束語

隨著人類航天活動(dòng)的不斷增多，航天器與碎片再入問題在國(guó)內(nèi)外航天界受到越來越多的關(guān)注和重視。但由于其復(fù)雜性和多學(xué)科交叉性，相關(guān)研究和應(yīng)用水平仍十分有限，一些基本問題尚未得到很好解決。在此背景下，本文對(duì)無控航天器與空間碎片再入問題工程預(yù)測(cè)方法與軟件開發(fā)的研究現(xiàn)狀進(jìn)行了跟蹤與總結(jié)，并討論了若干關(guān)鍵問題和研究設(shè)想，希望能對(duì)此領(lǐng)域研究與工程應(yīng)用起到拋磚引玉的作用。

（References）

[1]Reynolds R C, Soto A.Debris assessment software∶operators manual[G].NASA Johnson Space Center, 2001

[2]Opiela J N, Hillary E, Whitlock D O, et al.DAS user’s guide[G].Version 2.0.Orbital Debris Program Office, NASA Johnson Space Center.Houston, TX, 2007

[3]Bouslog S A, Ross B P, Madden C B.Space debris reentry risk analysis, AIAA 1994-591[R]

[4]Fritsche B, Klinkrad H, Kashkovsky A, et al.Spacecraft disintegration during uncontrolled atmospheric re-entry[J].Acta Astronaut, 2000, 47(2-9)∶513-522

[5]Wu Z N, Hu R F, Qu X, et al.Space debris reentry analysis methods and tools[J].Chin J Aeronaut, 2011, 24∶387-395

[6]胡銳鋒, 吳子牛, 曲溪, 等.空間碎片再入燒蝕預(yù)測(cè)與地面安全評(píng)估軟件系統(tǒng)[J].航空學(xué)報(bào), 2011, 32(3)∶390-399 Hu Ruifeng, Wu Ziniu, Qu Xi, et al.Debris reentry and ablation prediction and ground risk assessment software system[J].Acta Aeronaut ET Astronaut Sinica, 2011, 32(3)∶390-399

[7]Lips T, Fritsche B.A comparison of commonly used re-entry analysis tools[J].Acta Astronaut, 2005, 57(2-8)∶312-323

[8]Tsien S H.Similarity laws of hypersonic flows[J].J Math Phys, 1946, 25∶105-106

[9]Bird G A.Molecular gas dynamics and the direct simulation of gas flows[M].New York∶Oxford University Press, 1994∶162-172

[10]Cercignani C.Rarefied gas dynamics[M].Cambridge∶Cambridge University Press, 2000∶1-30

[11]沈青.稀薄氣體動(dòng)力學(xué)[M].北京∶國(guó)防工業(yè)出版社, 2003∶141-156

[12]Anderson J D.Hypersonic and high temperature gas dynamics[M].New York∶McGraw-Hill Book Company, 1989∶45-75

[13]Lees L.Laminar heat transfer over blunted-nosed bodies at hypersonic flight speeds[J].Jet Propul, 1956, 26(4)∶259-269

[14]Fay J A, Riddell F R.Theory of stagnation point heat transfer in dissociated air[J].J Aeronaut Sci, 1958, 25(2)∶73-85

[15]Detra R W, Kemp N H, Riddell F R.Addendum to ‘Heat transfer to satellite vehicles re-entering the atmosphere’[J].Jet Propul, 1957, 27(12)∶1256-1257

[16]Cohen N B.Boundary-layer similar solutions and correlation equations for laminar heat-transfer distribution in equilibrium air at velocities up to 41,100 feet per second, NASA TR R-118[R], 1961

[17]Sutton K, Graves R A.A general stagnation-point convective- heating equation for arbitrary gas mixtures, NASA TR R-376[R], 1971

[18]Tewari A.Entry trajectory model with thermomechanical breakup[J].Journal of Spacecraft and Rockets, 2009, 46(2)∶299-306

[19]Dejarnette F R, Davis R M.A simplified method for calculating laminar heat transfer over bodies at an angle of attack, NASA TN D-4720[R], 1968

[20]Dejarnette F R, Hamilton H H.Inviscid surface streamlines and heat transfer on shuttle-type configurations[J].Journal of Spacecraft and Rockets, 1973, 10(5)∶314-321

[21]Hamilton H H, Dejarnette F R, Weilmuenster K J.Application of axisymmetric analog for calculating heating in three-dimensional flows[J].Journal of Spacecraft and Rockets, 1987, 24(4)∶296-302

[22]Hamilton H H, Greene F A, Dejarnette F R.Approximate method for calculating heating rates on three- dimensional vehicles[J].Journal of Spacecraft and Rockets, 1994, 31(3)∶345-354

[23]Hamilton H H, Weilmuenster K J, Dejarnette F R.Approximate method for computing laminar and turbulent convective heating on hypersonic vehicles using unstructured grids, AIAA 2009-4310[R]

[24]Wilmoth R G, Mitcheltree R A, Moss J N.Low-density aerodynamics of the stardust sample return capsule[J].Journal of Spacecraft and Rockets, 1999, 36(3)∶436-441

[25]Ivanov M S, Markelov G N, Gimelshein S F, et al.High-altitude capsule aerodynamics with real gas effects[J].Journal of Spacecraft and Rockets, 1998, 35(1)∶16-22

[26]徐姍姝.過渡區(qū)飛行器流場(chǎng)的數(shù)值模擬研究[D].北京∶清華大學(xué), 2008

[27]Matting F W.Approximate bridging relations in the transitional regime between continuum and free- molecule flows[J].Journal of Spacecraft and Rockets, 1971, 8(1)∶35-40

[28]Nomura S.Correlation of hypersonic stagnation point heat transfer at low Reynolds number[J].AIAA Journal, 1983, 21(11)∶1598-1600

[29]Swaminathan P K, Taylor J C, Rault D F G, et al.Transition regime aerodynamic heating of missiles[J].Journal of Spacecraft and Rockets, 1996, 33(5)∶607-613

[30]Standard N S.Guidelines and assessment procedures for limiting orbital debris, NASA NSS 1995-1740[R]

[31]Baker R L, Weaver M.Orbital spacecraft reentry breakup[C]//50thInternational Astronautical Congress.Amsterdam, The Netherlands, 1999∶1-14

[32]Rochelle W M C, Kinsey R, Reid E, et al.Spacecraft orbital debris reentry aerothermal analysis[C]// Proceedings of the 8thAnnual Thermal and Fluids Analysis Workshop on Spacecraft Analysis and Design.Houston, TX, USA, 1997∶1-14

[33]Rochelle W C, Kirk B S, Ting B C, et al.Modeling of space debris reentry survivability and comparison of analytical methods[C]//50thInternational Astronautical Congress.Amsterdam, The Netherlands, 1999∶351-364

[34]Rochelle W C, Marichalar J J, Johnson N L.Analysis of reentry survivability of UARS spacecraft[J].Adv Space Res, 2004, 34(5)∶1049-1054

[35]Portelli C, Salotti L, Anselmo L, et al.BeppoSAX equatorial uncontrolled re-entry[J].Adv Space Res, 2004, 34(5)∶1029-1037

[36]Fritsche B, Koppenwallner G.Computation of destructive satellite re-entries[C]//Proceedings of the Third European Conference on Space Debris.Darmstadt, Germany, 2001∶527-532

[37]Fritsche B.Note on the application of SCARAB to the MIR re-entry[C]//Proceedings of the International Workshop MIR Deorbit.Darmstadt, Germany, 2001∶99-102

[38]Klinkrad H, Fritsche B, Lips T.A standardized method for re-entry risk evaluation[C]//55thInternational Astronautical Congress.Vancouver, Canada, 2004∶141-155

[39]The Aerospace Corporation.Atmospheric heating and breakup (AHaB) tool[G]

[40]Sim H S, Kim K H.Reentry survival analysis of tumbling metallic hollow cylinder[J].Adv Space Res, 2011, 48∶914-922

[41]Koppenwallner G, Fritsche B, Lips T.Survivability and ground risk potential of screws and bolts of disintegrating spacecraft during uncontrolled re-entry[C]// Proceedings of the Third European Conference on Space Debris.Darmstadt, Germany, 2001∶533-539

[42]Allen H J, Eggers A J.A study of the motion and aerodynamic heating of missiles entering the Earth’s atmosphere at high supersonic speeds, NACA R-1381[R], 1957

[43]Fritsche B, Lips T, Koppenwaller G.Analytical and numerical re-entry analysis of simple-shaped objects[J].Acta Astronaut, 2007, 60(8/9)∶737-751

[44]Rao P P, Woeste M A.Monte Carlo analysis of satellite debris footprint dispersion, AIAA 1979-1628[R]

[45]Spencer D A, Braun R D.Mars Pathfinder atmospheric entry∶trajectory design and dispersion analysis[J].Journal of Spacecraft and Rockets, 1996, 33(5)∶670-676

[46]Desai N P, Braun R D, Powell R W, et al.Six- degree-of-freedom entry dispersion analysis for the METEOR recovery capsule[J].Journal of Spacecraft and Rockets, 1997, 34(3)∶334-340

[47]Queen E M, Cheatwood F M, Powell R W, et al.Mars polar lander aerothermodynamic and entry dispersion analysis[J].Journal of Spacecraft and Rockets, 1999, 36(3)∶421-428

[48]Desai N P, Cheatwood F M.Entry dispersion analysis for the Genesis Sample Return Capsule[J].Journal of Spacecraft and Rockets, 2001, 38(3)∶345-350

[49]Purcell E W, Barbery T B.Dispersions resulting from atmospheric variations[J].Journal of Spacecraft and Rockets, 1968, 5(8)∶1005-1007

[50]Frank M V, Weaver M A, Baker R L.A probabilistic paradigm for spacecraft random reentry disassembly[J].Reliab Eng Syst Saf, 2005, 90∶148-161

[51]Reyhanoglu M, Alvarado J.Estimation of debris dispersion due to a space vehicle breakup during reentry[J].Acta Astronaut, 2013, 86∶211-218

[52]Wiener N.The homogeneous chaos[J].Am J Math, 1938, 60(4)∶897-936

[53]Eldred M S.Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design, AIAA 2009-2274[R]

[54]Sahai T, Pasini J M.Uncertainty quantification in hybrid dynamic systems[J].J Comput Phys, 2013, 237∶411-427

[55]Luchtenburg D M, Brunton S L, Rowley C W.Long- time uncertainty propagation using generalized polynomial chaos and flow map composition[J].J Comput Phys, 2014, 274∶783-802

[56]Myers R H, Montgomery D C, Anderson-Cook C M.Response surface methodology∶process and product optimization using designed experiments[M].Hoboken, NJ∶John Wiley &Sons, 2009∶1-680

[57]Riley M E, Grandhi R V.Quantification of modeling- induced uncertainties in simulation-based analyses[J].AIAA Journal, 2014, 52(1)∶195-202

[58]Halder A, Bhattacharya R.Dispersion analysis in hypersonic flight during planetary entry using stochastic Liouville equation[J].J Guide Control Dyn, 2011, 34(2)∶459-474

[59]胡銳鋒.高速飛行氣動(dòng)環(huán)境、氣動(dòng)特性快速預(yù)測(cè)與應(yīng)用[D].北京∶清華大學(xué), 2013

[60]Kinslow M, Potter J L.Drag of spheres in rarefied hypervelocity flow[J].AIAA Journal, 1963, 1(11)∶2467-2473

[61]Bailey A B, Hiatt J.Sphere drag coefficients for a broad range of Mach and Reynolds numbers[J].AIAA Journal, 1972, 10(11)∶1436-1440

[62]Henderson C B.Drag coefficients of spheres in continuum and rarefied flows[J].AIAA Journal, 1976, 14(6)∶707-708

[63]Loth E.Compressibility and rarefaction effects on drag of a spherical particle[J].AIAA Journal, 2008, 46(9)∶2219-2228

[64]Klett R D.Drag coefficients and heating ratios for right circular cylinders in free-molecular and continuum flow from Mach 10 to 30, SC-RR-1964-2141[R].Sandia Corp, Albuquerque, NM

[65]Randall D E.Method for estimating the drag coefficient of a tumbling circular cylinder, SC-TM-1964-0528[R].Sandia Corp, Albuquerque, NM

[66]Lutz S A.Heating correlations for bluff cylinder hypersonic rarefied flows, AIAA 2003-4060[R]