基于雷諾數的大氣阻力模型在飛行器再入預報中的應用*

劉勁宏 徐勁 杜建麗 潘建平

1(中國科學院紫金山天文臺 南京 210023)

2(重慶交通大學智慧城市學院 重慶 400074)

0 引言

大氣再入預報是指目標從當前軌道高度自然衰減到120 km 所需時間。大氣再入分為兩類，即受控再入和無控再入。受控再入以偏遠地區(qū)為目標，對地面威脅最小，而無控再入受多種因素影響，其預報的墜落時間、地點具有隨機性。部分目標結構堅固、耐高溫，在大氣再入過程中難以燒蝕和解體，對地面人員和財產安全構成潛在威脅。Pardini等[1]統(tǒng)計了2008－2017年完整航天器無控再入數量與質量的分布，共計448 個，總質量約900 t[1]。統(tǒng)計數據表明：過去10年內，平均每周至少有一個完整的飛行器（衛(wèi)星、火箭體和運載平臺）重返大氣層，平均質量約2 t，導致每年有50～100 t 的完整飛行器重返大氣層。隨著巨型衛(wèi)星星座的布設及碎片主動移除技術的論證與實現，未來再入事件的發(fā)生頻率將更高且質量更大，對人類空間活動及地面安全構成的威脅與日俱增。為此，亟需對在軌的大型航天器進行監(jiān)測、跟蹤和預警，以應對可能出現的緊急突發(fā)情況。

STK（Satellite Tool Kit）提供多種軌道預報模型，可選擇模型包括：TwoBody、J2 Perturbation、J4 Perturbation、HPOP（High Precision Numerical Orbit Propagator）等。本文選擇高精度數值軌道傳播器HPOP[2-4]和中國武漢大學提出的半解析軌道傳播器WHU-SST（Wuhan University Semi-Analytical Satellite Theory）[5]對低軌道目標，特別是臨近大氣層的目標進行再入預報。HPOP 采用COWELL 數值方法對力學模型進行積分，WHU-SST 是基于多尺度攝動理論的半解析軌道傳播方法，其核心思想是將軌道根數分為長周期變量和短周期變量，采用數值積分方法對目標的平均軌道根數進行大步長積分，利用解析算法重建短周期項，最后對計算得到的平均軌道根數和短周期項重組，獲得快速、精確的軌道。為便于后續(xù)分析討論，對以上兩種軌道預報模型設置相同的參數，采用Runge-Kutta-Fehlberg 7(8)積分，引力選擇JGM3 和70×70 階次，太陽輻射壓為球形，大氣密度模型選用NRLMSISE-00，大氣阻力系數CD=2.2，第三體引力考慮太陽和月球。

對于低軌目標的軌道預報，大氣阻力系數模型是影響再入大氣預報精度的關鍵。實際應用中，大氣阻力系數以三種方式取值，即常數值、擬合值和物理值[6,7]。一般軌道傳播器的大氣阻力系數采用常數值，通常設置為2.2，然而其不能反映目標阻力系數的真實情況，使得計算的大氣阻力存在很大偏差[7]。擬合大氣阻力系數通過軌道確定獲得，與選擇的大氣模型相關，同時也會吸收其他模型的誤差[8]。物理大氣阻力系數反映了航天器表面與自由分子流的能量和動量交換，主要受到覆蓋在航天器表面的氧原子吸附影響[9]，其數值精確，但是計算量大且耗時長。

事實上，已有較多研究者在實驗室通過模擬大氣環(huán)境獲得不同形狀、縱橫比的大氣動力經驗系數。針對縱橫比近似3∶1 的圓柱體飛行器，本文在簡化飛行器模型的基礎上引入基于雷諾數的大氣動力系數經驗模型，采用RK6 (7)對運動方程重新積分，探討該模型用于火箭體再入預報所需條件，為后續(xù)研究奠定基礎。該模型能夠同時考慮大氣阻力和升力的影響，適用雷諾數范圍為10～300[10]。

1 大氣動力模型

通常臨近大氣層的飛行器，其大氣阻力系數變化劇烈，可引起極大的預報誤差（見圖1）。在再入前的第10 天左右，飛行器運行高度開始陡然下降，表明大氣阻力系數取常數值或擬合值均不能真實模擬飛行器再入軌跡的改變。因此提出了基于雷諾數的大氣動力系數經驗模型。

圖1 火箭體40145 和41027 再入過程中的高度變化Fig.1 Altitude changes with time during the objects of 40145 and 41027 reentry atmosphere

雷諾數是用于表征流體流動狀態(tài)的無量綱數，其定義如下：

式中：ρ為流體密度，由NRLMSISE-00 大氣密度模型計算得到；v∞為自由流速度；D為等效球體直徑；μ為流體動力粘度系數。

為評估圓柱體受到的大氣阻力大小，Cao等[10]建立了基于雷諾數的阻力系數和升力系數表達式，討論了其隨流體入射角的變化。假定流體入射角為θ，則阻力系數CD,θ和升力系數CL,θ的計算公式如下[10]：

式（2）～（5）適用的雷諾數范圍為[10,300]。實驗表明：阻力系數平均相對誤差為1.5%，均方根誤差為1.5×10?3；升力系數平均相對誤差為5.77%，均方根誤差為1.9×10?3[10]。

假定火箭體的長度為L，直徑為D，則投影到入射流方向的面積可表示為

假定質量為m，彈道系數（Ballistic Coefficient，定義數學符號CB）可表示為

假定火箭體和殘存燃料質量共10 t，利用式（2）～（5）、式（6）和式（7）計算不同雷諾數下彈道系數隨入射角的變化情況，如圖2 所示，雷諾數變化范圍為[10,300]。從圖2 可以發(fā)現，盡管雷諾數不同，但是彈道系數隨入射角變化的峰值都在75°～80°內，四幅圖的變化趨勢完全一致，受火箭體尺寸影響，縱軸變化范圍有所不同。

圖2 不同縱橫比的彈道系數隨入射角度的變化Fig.2 Variation of ballistic coefficient with incident angle for different aspect ratio

2 火箭體長時間再入預報

火箭體的大氣再入過程具有如下特點：歷史再入數量多、質量大、頻率高且難燒蝕，結構保持完整。這里主要討論的對象為中國航天發(fā)射活動使用的長征系列運載火箭體。表1 給出了部分二級火箭體的尺寸和質量，以便估算與大氣阻力相關的物理參數，其中定義等效球體直徑為與火箭體體積相同的球體直徑*http://cn.cgwic.com/LaunchServices/LaunchVehicle/LM3 B.html。

表15 種長征系列火箭體的二級火箭尺寸和質量Table 1 Size and mass of two-stage rockets of five kinds of Long March series rocket carriers

假定火箭體在再入大氣層之前已處于長時間穩(wěn)定狀態(tài)，即姿態(tài)穩(wěn)定、質量穩(wěn)定，此時面質比及受到的流體入射角度應為常量。基于此，如何利用歷史TLE 數據確定合適的面質比和入射角，是采用基于雷諾數的大氣動力模型精確預報火箭體再入的關鍵。這里主要討論時間長度為一個月的再入預報。

在飛行坐標系下，大氣攝動加速度采用下式描述：

加速度的矢量方向分別為速度方向（阻力方向）、垂直于速度向上的方向（升力方向）以及與這兩個方向相互垂直的方向。其他攝動力包括EGM2008 20×20階引力場、太陽輻射壓、第三體（日、月）引力，大氣密度模型選用NRLMSISE-00 模型。采用RK6 (7)，即Runge-Kutta-Fehlberg 6 (7)對軌道運動方程以固定步長（10 s）向前積分，獲得再入時間的預報。

預報采用百分比誤差，其定義為[11]

式中，Tref為真實再入時間，Tpred為預報再入時間，Tinitial為預報的起始時刻。

收集2000年以來再入的火箭體數據，剔除因觀測值稀少及質量不佳的某些火箭體，經過數據質量清理、彈道系數估計和軌道定軌，獲得5 種火載體共計26 個預報結果，其中：CZ-4B R/B 12 個，CZ-3B R/B 6 個，CZ-2C R/B 5 個，CZ-2D R/B 2 個，CZ-3C R/B 1 個。圖3 為5 種火載體在入射角變化范圍為[1°,90°]時再入前30 天的預報結果?？梢钥闯觯S著入射角的增加，大氣阻力系數越大，預報的再入時間越短。圖3(f) 給出了26 個火箭體再入預報誤差分布，流體入射角為12°。

圖3 入射角對火箭體再入預報的影響和入射角為12°時火箭體的預報誤差百分比（紅色實橫線是預測誤差為0 的基線，即第30 天結果，紅色點線為預報誤差±10%的范圍線，品紅點線為預報誤差±20%的范圍線）Fig.3 Influence of the incident angle on the reentry prediction of the rocket body and the prediction error percentage of the rocket body when the incident angle is 12° (Red solid line is baseline with 0 prediction error,the red dotted line deviates ±10% from the baseline,and the magenta dotted line deviates ±20%）

已知實際再入時間在第30 天。因此，CZ-4B 火箭體對應的入射角應在10°～20°內，CZ-3B 火箭體對應的入射角在10°～30°內，CZ-2C 火箭體對應的入射角在12°～24°內，CZ-2D 火箭體對應的入射角在10°～23°內，CZ-3C 火箭體的對應的入射角在2°～3°內。通過分析，先驗入射角范圍為12°～14°，可以使得84.62% 的火箭體再入時間預報誤差小于20%，53.84%的火箭體預報誤差小于10%。

為了與高精度數值軌道傳播器HPOP 和半解析軌道傳播器WHU-SST 的預報結果進行比較，選擇8 個火箭體作為目標。表2 列出了8 個目標的軌道信息和預報誤差。WHU-SST 的預報誤差整體比HPOP預報誤差低30% 以上，而新模型的預報誤差最高為16.1%，遠小于WHU-SST 的最高預報誤差83.61%。因此，利用基于雷諾數的大氣阻力模型不僅能夠有效降低再入預報誤差，還可以提高預報成功率。

需要說明的是，表2 中目標40550 屬于地球同步轉移軌道目標，而本文使用的基于雷諾數的大氣模型僅適用于200 km 左右高度的火箭體目標，因此沒有進行再入預報計算。

3 在隕落位置預報中的應用

將基于雷諾數的大氣阻力系數模型應用于隕落預報，即高度范圍為10～120 km。目標22659、24770、25240、26607、38086 和41107 已墜落于地球表面，圖4 中黑色三角形為殘骸找到的真實位置，這里將其作為隕落位置，紅色點為隕落預報位置。隕落預報考慮了定軌誤差和大氣密度模型誤差，因此呈現散點分布。本文基于TLE 數據進行隕落預報，僅能實現實際隕落位置位于預報散點中的某一點。圖4 所列6 個目標（形狀僅包括球體和圓柱體）的隕落位置預報受各種誤差影響，呈現散點或者近似多條曲線分布，真實隕落位置位于其中的預報軌跡中。

圖4 基于雷諾數的大氣阻力系數模型的隕落預報位置分布（黑色三角形為殘骸找到的真實位置，這里將其作為隕落位置；紅色點為隕落預報位置）Fig.4 Distributions of decayed locations predicted by new atmospheric drag coefficient model based on Reynolds number (Black triangle is for the real falling location and the red dot is for the falling location of prediction)

4 結論

針對現有軌道傳播器對低軌目標大氣阻力系數描述不準確的問題，在簡化航天器模型的基礎上，引入基于雷諾數的大氣動力系數經驗模型，確定了適用于低軌火箭體的最佳流體入射角，重新對運動微分方程進行軌道數值積分，所獲結果與現有軌道傳播器的預報結果進行對比，結論如下。

（1）與已有的軌道傳播器相比，基于雷諾數的大氣阻力系數模型在中國火箭體的再入預報方面具有明顯優(yōu)勢，預報誤差顯著降低，預報成功率上升。例如目標41027 的預報誤差從75%～96%下降至7.6%，目標38253 和40120 預報成功。通過6 個目標的隕落位置預報，能夠確保真實隕落位置位于預報的統(tǒng)計隕落位置中，間接證明了基于雷諾數的大氣阻力系數模型的有效性。

（2）基于雷諾數的大氣動力系數模型雖然在火箭體的預報中取得了較好的預報結果，但是受限于雷諾數的應用范圍，該模型不適用于地球同步轉移軌道目標的再入情況。由于收集到的火箭體數量有限，未對再入預報進行外推核驗。未來會利用新的再入火箭體，驗證基于雷諾數的大氣動力系數模型的可靠性。

表28 個火箭體軌道信息和預報誤差Table 2 Orbit information and prediction errors of 8 rocket bodies