Journal of Space Science and Experiment >

2025 , Vol. 2 >Issue 2: 103 - 113

DOI: https://doi.org/10.19963/j.cnki.2097-4302.2025.02.012

Space Debris and Space Environment

Aerodynamic Drag Coefficients of Reentry Objects Based on Rarefied Gas Dynamics

  • Xu YANG , 1 ,
  • Yaping WANG , 2, * ,
  • Yue WANG 3 ,
  • Yao ZHANG 1 ,
  • Dawei LI 1
Expand
  • 1. National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China
  • 2. China Satellite Network System Research Institute Company Limited, Beijing 100070, China
  • 3. School of Astronautics, Beihang University, Beijing 102206, China

Online published: 2025-07-18

Fold

Abstract

The drag coefficient is a crucial parameter for calculating aerodynamic drag during space object re-entry prediction. In most cases, compact shape objects generally have a drag coefficient of 2.2, but there is a significant difference between this empirical value and the actual drag coefficient. This difference limits the accuracy of re-entry prediction to some extent. This paper calculates the aerodynamic drag coefficients of typical shaped objects with high precision, based on rarefied gas dynamics theory. The Schaaf-Chambre gas-surface interaction model is applied when the object is in the free molecular flow region. Upon entering the transition flow region with orbit attenuation, the Wilmoth bridging formula is used to calculate the drag coefficients, using the drag coefficients in the continuum and free-molecular limits. The study shows that the drag coefficients decrease as the orbital altitude decreases. Factors such as latitude, longitude, solar and geomagnetic activities indirectly affect the drag coefficient by influencing parameters such as the density and temperature of the incoming flow. These factors result in differences of less than 6 percent in the drag coefficient. The drag coefficients of several typical objects differ significantly from the empirical value of 2.2, confirming the necessity for this research.

Key words： rarefied gas; aerodynamics; reentry prediction; aerodynamic drag coefficient

Cite this article

Xu YANG , Yaping WANG , Yue WANG , Yao ZHANG , Dawei LI . Aerodynamic Drag Coefficients of Reentry Objects Based on Rarefied Gas Dynamics[J]. Journal of Space Science and Experiment, 2025 , 2(2) : 103 -113 . DOI: 10.19963/j.cnki.2097-4302.2025.02.012

0 引 言

随着人类航天活动的激增，在轨空间物体数量急剧增加。部分空间物体因受高层大气阻力的影响，轨道高度逐渐降低，最终无控再入大气层。大质量再入目标的烧蚀残留部件会对落点区域构成安全威胁。因此，必须提高空间目标的再入预报精度。
在再入预报末段，目标依次穿过地球高层稀薄大气的自由分子流区域和过渡流区域。随着轨道高度的降低，大气密度近似呈指数增长，气动阻力愈发显著，其强不确定性成为提高空间目标再入预报精度的主要障碍，而阻力系数是不确定性的主要来源之一。对于在自由分子流中的球体、椭球体等紧凑形状的目标，阻力系数通常取2.2，该值是库克在早期对轨道上气体-表面相互作用了解较少时的近似估计[1]，并被广泛应用于航天器的气动阻力计算[2-3]。但后续研究表明，该经验值与真实物理现象存在较大差距。随着对自由分子流区域气体与表面相互作用研究的深入，精确计算特定形状目标的阻力系数已成为可能。Moe等[1]根据Sentman模型计算得到了简单几何体以及S3-1卫星的阻力系数随轨道高度的变化曲线。Mehta等[4]在不同气体-表面相互作用模型下，分析了简单几何体阻力系数对来流参数及物体表面温度的敏感性，并研究了球体阻力系数在太阳极大期和极小期时随轨道高度的变化情况。Bullard[5]使用直接模拟蒙特卡罗（Direct Simulation Monte Carlo, DSMC）方法，计算得到了70°钝体在150~300 km轨道高度的阻力系数曲线，并进一步研究了基本几何形状的阻力系数随来流温度和分子数密度的变化。然而，目前的阻力系数计算主要集中于自由分子流区域，未考虑临再入前的过渡流区域，且部分研究仅分析了单一环境参数对阻力系数的影响，缺乏对大气参数、太阳和地磁活动、经纬度等环境参数综合影响的研究。
针对目前研究存在的不足，本文基于稀薄气体动力学理论，采用Schaaf-Chambre气体-表面相互作用模型，推导目标表面微元的压力和剪切应力系数，对目标表面受到来流影响的面积微元进行积分，最终得到目标的气动阻力系数。当目标随轨道衰减进入过渡流区域时，利用Wilmoth桥接公式进行计算。针对目标姿态未知的问题，将目标在不同姿态情况下阻力系数的平均值作为名义值。最后，基于上述计算结果，对比和分析不同外形目标在不同轨道高度的气动阻力系数变化规律，并讨论经纬度、太阳和地磁活动等环境参数的影响。本文研究采用了稀薄气体动力学理论，不仅可以克服目前广泛采用的经验型阻力系数的缺点，还对多种环境参数的影响进行综合分析，有望提高再入目标的气动阻力系数计算精度和再入预报精度。

1 气体动力学模型

区分地球高层稀薄大气各区域的主要依据是Knudsen数（Kn），定义为气体分子平均自由程与物体特征尺寸的比值。当Kn < 0.01时，大气被视为连续流体；而Kn > 0.01时，则为稀薄流体，并进一步划分为过渡流（0.01 < Kn < 10.00）和自由分子流（Kn > 10.00）[6]。
地球大气密度随轨道高度增加迅速减小，相应地，Kn随轨道高度而增大。如图1所示，Kn在轨道高度131 km时为10，这意味着该高度以上的大气可被视为自由分子流；在高度87 km为0.01，即该高度以下的大气被视为连续流体，87~131 km之间则为过渡流。
图 1 NRLMSISE-00大气模型下，Kn随海拔的变化曲线（特征长度1 m）

Fig.1 The variation of Kn with respect to altitude using the NRLMSISE-00 atmospheric model (characteristic dimension of 1 m)

近年来，随着大气环境数据的日益丰富，研究人员不断修正并提出新的大气模型，以更准确地模拟真实的大气环境，为低轨环境的相关研究提供更加有力的支撑。目前，使用比较广泛的大气模型有：空间研究委员会国际参考大气模型（Committee on Space Research International Reference Atmosphere，CIRA），美国标准大气模型（U.S. Standard Atmosphere，USSA），Jacchia系列大气模型，NRLMSISE-00大气模型（质谱仪非相关散射模型），JB2008大气模型（卫星阻力总密度模型）和阻力温度模型（Drag Temperature Model，DTM）[7-9]。本文采用NRLMSISE-00大气模型。该模型能够覆盖从地球表面到外大气层底部，并且考虑了地理位置、季节因素以及太阳和地磁活动对大气参数的影响，在超低轨环境的相关研究中应用广泛。利用该模型，可以得到大气温度、大气密度及各气体成分的分子数密度。

1.1 自由分子流中的气动力模型

自由分子流中，分子数密度极小，气体分子的粒子特性开始显现，连续介质假设不再适用。自由分子流的主要假设是反射粒子对来流没有影响，对于具有一定攻角的平板，意味着平板前缘反射的分子流不会干扰平板后缘的入射分子流。该假设的另一个结果是，边界层非常微弱，不会影响气体分子在物体表面的反射[10]。
在自由分子流中建立气动力模型需要考虑三个相互关联的模型：高层大气模型、粒子流的力特性模型以及气体表面相互作用模型。现有的稀薄气体气动力模型的基础是气体的分子动力学理论，其中假定每个气体分子的运动由经典力学定律描述。气体动力学方程（玻尔兹曼方程）假设：自由流中的物体对来流方向无影响；物体表面反射的分子不会再次撞击物体表面；气体分子的平均自由程远大于物体的特征长度；气流上无外力作用[11]。除了入射流，确定作用在物体上的气动力还需要知道反射流的特性，而气体与表面相互作用的复杂机制尚未被充分理解。因此，自由分子流中气动力计算依赖于特定的气体-表面相互作用模型[12]。
气体-表面相互作用（Gas surface interaction，GSI）模型的关键在于对气体与表面之间由于碰撞和反射而产生的能量和动量的交换进行精确建模，模型中需要考虑许多因素，包括来流气体特性（密度、速度、组分、温度等），表面特性（材料、温度、清洁度、粗糙度等）以及来流与表面之间的入射角、入射速度等[12]。目前，使用比较广泛的GSI模型主要有：麦克斯韦（Maxwell）模型、尚贝格（Schamberg）模型、不完全调节漫反射（Diffuse Reflection with Incomplete Accommodation，DRIA）模型、Schaaf-Chambre模型、CLL（Cercignani-Lampis-Lord）模型）[12-15]。每个GSI模型的核心是量化气体和表面之间交换能量和（或）动量的调节系数[15]。本文采用Schaaf-Chambre模型，该模型引入法向调节系数$ {{\sigma }_N} $和切向调节系数${{\sigma }_T}$，用于描述法向和切向的动量传递到表面的程度，从而区分不完全法向和切向动量传递的影响。Schaaf-Chambre模型能够描述从完全镜面反射到完全漫反射之间的一系列复杂散射行为，能够准确反映不同表面条件和流体特性下气体分子散射的真实物理现象，相比一些简单模型更具优势。此外，目前基于Schaaf-Chambre模型已有简单凸几何形状的气动力系数封闭解析解，便于对比验证平板法的精度。
对于完全镜面反射，${{\sigma }_N} = {{\sigma }_T} = 0$；对于完全漫反射，则有${{\sigma }_N} = {{\sigma }_T} = 1$。
${{\sigma }_N}$、${{\sigma }_T}$的定义如下：
$ {{\sigma }_N} = \frac{{{p_i} - {p_r}}}{{{p_i} - {p_w}}} $
$ {{\sigma }_T} = \frac{{{{\text{τ} }_i} - {{\text{τ} }_r}}}{{{{\text{τ} }_i} - {{\text{τ} }_w}}} = \frac{{{{\text{τ} }_i} - {{\text{τ} }_r}}}{{{{\text{τ} }_i}}}，{{\text{τ} }_w} = 0 $
式中，$p$是法向动量通量；${\text{τ} }$是切向动量通量；下标$i$表示入射；下标$r$表示反射；下标$w$表示粒子与表面实现了完全热平衡并根据表面温度$({T_w})$的麦克斯韦分布重新发射。
各动量通量如图2所示，其中质量通量${{\mathrm{d}}}Q = \text{ρ}V\sin {\theta }$。通常认为这两个调节系数相等，在轨道高度200 km处且太阳活动强度较低时约为0.99，当轨道高度低于200 km时则取1[16]。
图 2 入射通量和反射通量

Fig.2 Incident and reflected fluxes

1.2 平板法

平板法是一种用于计算自由分子流中物体气动力的数值方法，适用于简单凸几何体的气动计算。该方法将物体理想地离散为多个面板，每个面板建立一个平板模型，并用所选的GSI模型计算迎风面板的压力系数和切向力系数。最后，将所有受来流影响的面板贡献值叠加，得到物体的整体动系数。
如图3所示，在一般表面上建立坐标系，${{V}_\infty }$是自由分子流的来流速度矢量，${V_\infty }$通常取轨道速度；${{\boldsymbol{n}}}$为表面法向量，指向表面外侧；${\boldsymbol{t}}$为表面切向量，其定义为$ {{\boldsymbol{t}}} = {{\boldsymbol{n}}} \times ( - {\text{ }}{{V}_\infty } \times {{\boldsymbol{n}}}) $；${\theta }$为入射角，由$\sin {\theta } = - {\text{ }}{{V}_\infty } \cdot {{\boldsymbol{n}}}$定义，攻角$ \alpha = \dfrac{{{\text π} }}{2} - {\theta } $。
图 3 表面的法向方向和切向方向

Fig.3 Normal and tangential directions of a surface

根据式（1）和式（2），可以用入射动量通量和动量调节系数表示表面的反射动量通量：
$ {p_r} = (1 - {{\sigma }_N}){p_i} + {{\sigma }_N}{p_w} $
$ {{\text{τ}}_r} = (1 - {{\sigma }_T}){{\text{τ}}_i} $
而
$ {p_i} = \frac{{\text{ρ}{V_i}^2}}{{2{s^2}}}{\Gamma_1}(s\sin {\theta }) $
$ {{\text{τ} }_i} = \frac{{\text{ρ}{V_i}^2}}{{2s}}{\Gamma_2}(s\sin {\theta }) $
式中，${\text{ρ}}$为来流的密度；${V_i}$为来流的速度，即${V_i} = {V_\infty }$。
$ {\Gamma_1}(x) = \frac{1}{{\sqrt {{\text π}} }}\left[ {x\exp ( - {x^2}) + \frac{{\sqrt {{\text π}} }}{2}(1 + 2{x^2})(1 + {\mathrm{erf}}{\text{ }}x)} \right] $
$ {\Gamma_2}(x) = \frac{1}{{\sqrt {{\text π} } }}\left[ {\exp ( - {x^2}) + \sqrt {{\text π} } x(1 + {\mathrm{erf}}{\text{ }}x)} \right] $
$ {\text{erf }}x = \frac{2}{{\sqrt {{\text π} } }}\int_0^x {{{\mathrm{e}}^{ - {t^2}}}} {{\mathrm{d}}}t $
$ s = \frac{{{V_\infty }}}{{\sqrt {2\dfrac{{\mathrm{R}}}{{{M_{{\mathrm{mol}}}}}}{T_\infty }} }} $
式中，分子速率比$s$通常要求大于5；${{\mathrm{R}}}$是一般气体常数8.314 J/（mol·K）；${M_{{\mathrm{mol}}}}$是气体摩尔质量；${T_\infty }$是自由分子流的气体温度。参数$\text{ρ} $和${T_\infty }$可由NRLMSISE-00大气模型直接得到，而$ {M_{{\mathrm{mol}}}} $则需要进一步计算得到。
因此，可以得到暴露在自由流中的面板的压力和剪切系数的表达式：
$ {C_p} = \frac{1}{{{s^2}}}\left[ {\left( {2 - {\sigma _N}} \right){\Gamma_1}(s\sin {\theta}) + \frac{{{{\sigma }_N}}}{2}\sqrt {\frac{{{T_w}}}{{{T_\infty }}}} \sqrt {{\text π}} {\Gamma_2}(s\sin {\theta})} \right] $
$ {C_{\text{τ}}} = \frac{{{{\sigma }_T}\cos {\text{θ}}}}{s}{\Gamma_2}(s\sin {\theta}) $
式中，$ \theta $是入射角；${T_w}$为物体表面温度；通常取300 K。
最后，将所有受来流影响面板上压力和剪切系数按式（13）相加，得到自由流中物体的气动力系数$ {C_F} $的表达式：
$ {C_F} = \sum\limits_{i = 1}^n {\left\{ { - {C_{{p_i}}}{\boldsymbol{n}} - {C_{{{\text{τ}}_i}}}{\boldsymbol{t}}} \right\}} \frac{{{A_i}}}{{{A_{{\mathrm{ref}}}}}} $
式中，${A_i}$为面板的面积；${A_{{\mathrm{ref}}}}$为物体的参考面积。

1.3 过渡流区域中的气动力模型

在自由分子流中，需要考虑气体分子粒子特性以及分子之间或气体分子与表面的相互作用；而连续流则具有极短的弛豫时间，其宏观性质被认为是连续变化的。过渡流则处于两者之间，其气体分子的复杂运动开始受到分子间相互作用力的影响。在该区域内，气体分子-表面的相互作用和分子间的相互作用力同样重要。目前，过渡流区域内的复杂相互作用机制尚未被充分了解，因此该区域内的阻力系数计算通常采用桥接公式。
Wilmoth等[17]提出了简单的正弦平方形式的桥接公式，并用于预测星尘号探测器返回舱（Stardust Sample Return Capsule, SRC）在低密度过渡流中的气动力，计算结果能够较好地拟合高精确数据，公式的有效性和准确性已得到验证。Wilmoth桥接公式如下[6,17]：
$ {C_D} = C_D^{{\mathrm{CR}}} + \left( {C_D^{{\mathrm{FM}}} - C_D^{{\mathrm{CR}}}} \right){\sin ^2}\phi $
式中，$ C_D^{{\mathrm{CR}}} $和$ C_D^{{\mathrm{FM}}} $分别是考虑连续流-过渡流分界处、自由分子流-过渡流分界处的阻力系数，参数$\phi $由当地的Kn数给出：
$ \phi = {{\text π}}\left( {{a_1} + {a_2}\lg Kn} \right) $
$ {a_1} = - \frac{1}{2}\cdot \frac{{\lg K{n^{{\mathrm{CR}}}}}}{{\lg K{n^{{\mathrm{FM}}}} - \lg K{n^{{\mathrm{CR}}}}}} $
$ {a_2} = \frac{1}{2}\cdot \frac{1}{{\lg K{n^{{\mathrm{FM}}}} - \lg K{n^{{\mathrm{CR}}}}}} $
式中，$K{n^{{\mathrm{FM}}}} = 10$，$K{n^{{\mathrm{CR}}}} = 0.01$。
在本文研究中，自由分子流-过渡流分界处的阻力系数$ C_D^{{\mathrm{FM}}} $，可由第1.1~1.2节中介绍的方法计算得出。而连续流-过渡流分界处的阻力系数$ C_D^{{\mathrm{CR}}} $，则使用DSMC方法的软件包DS3V进行计算。Bird[18]于1994年提出了DSMC方法，并开发了原始的DSMC代码：DS2V和DS3V。DSMC方法直接模拟气体分子，每个模拟分子代表1012~1020个真实分子，分子之间可以相互碰撞。因此，该方法的计算所需时间与气体的密度成正比。此外，DSMC方法使用GSI模型处理气体-表面的相互作用，而物体则通常用表面网格表示[12]。在DS3V中，有两种GSI模型可供选择：完全漫反射模型和CLL模型。本文采用CLL模型，并在低海拔下的完全漫反射情况下，将法向能量调节系数和切向动量调节系数均设置为1。

2 典型形状目标的阻力系数

典型形状目标的阻力系数计算采用的部分参数如表1所示。对于不同的轨道高度，来流速度将根据相应的轨道速度进行选取，而大气温度、密度、分子自由程、分子数密度等大气参数则通过NRLMSISE-00大气模型计算得到。这些参数将用于后续的平板法和DS3V计算。平板法和DS3V将采用相同的网格尺寸，具体由后续的网格划分过程确定。DS3V的时间步长由软件根据计算需求动态调整，以确保计算过程的稳定性和结果的准确性。
表 1 典型形状目标阻力系数的计算参数

Table 1 Simulation conditions for drag coefficients of typical shaped objects

物理量 参数值
目标表面温度Tw / K 300
轨道纬度 40°N
轨道经度 120°E
时间 2022年9月6日00时
太阳活动F10.7 / SFU 125
地磁活动Ap / 2nT 19
参考面积Aref / m2 球：0.7854
平板：1.0000
正方体：1.0000
圆柱体：0.7854

2.1 球

针对再入目标为直径1 m的球体，本文采用平板法和DSMC方法进行网格划分，如图4所示。在本文算例中，设置最大棱长为0.1 m，由ANSYS软件自动进行网格划分。图5是DS3V计算球体在连续流-过渡流分界处阻力系数时，沿球体与来流方向平行的对称面上的来流密度分布情况。通过计算可知，球的阻力系数随轨道高度的变化如图6所示，其总体趋势是随轨道高度的降低而减小。具体地，当球处于自由分子流区域内，阻力系数随轨道高度降低而略微减小；而当球进入过渡流区域时，阻力系数呈现迅速减小的趋势。
图 4 球的网格划分

Fig.4 Mesh division of a sphere

图 5 球在连续流-过渡流分界处时，来流的密度分布

Fig.5 Density distribution of the incoming flow when the sphere is at the continuous-transitional flow boundary

图 6 球的阻力系数随轨道高度的变化

Fig.6 The variation of a sphere's drag coefficient with respect to the altitude

2.2 平板

当目标为1 m × 1 m × 1 mm的平板形物体时，平板法和DSMC方法采用的网格划分如图7所示。由于目标姿态未知，本文假设目标姿态运动均匀，计算目标在攻角分别为0°、15°、30°、45°、60°、75°和90°时的阻力系数，求其平均值作为平板形目标的名义气动阻力系数。在计算连续流-过渡流分界处不同攻角平板的阻力系数时，得到与来流平行的xOz面上来流的密度分布，如图8所示。
图 7 平板的网格划分

Fig.7 Mesh division of a panel

图 8 连续流-过渡流分界处，来流在平板攻角不同时的密度分布

Fig.8 Density distribution of incoming flow at the continuous-transition flow boundary with different angles of attack for a panel

如图9所示，平板的阻力系数随轨道高度降低呈现缓慢减小趋势，而当进入过渡流领域后，阻力系数随轨道高度降低而迅速减小。值得注意的是，轨道高度200 km降至190 km的过程中，阻力系数减小的速度略大于190 km降至180 km。这是因为当轨道高度从200 km降至190 km时，调节系数从0.99变为1，平板的表面上只发生漫反射，而漫反射传递给表面的法向动量小于镜面反射，从而对平板的阻力系数产生减小的作用。
图 9 平板的阻力系数随轨道高度的变化

Fig.9 The variation of a panel's drag coefficient with respect to the altitude

2.3 正方体

正方体的外形尺寸和网格如图10所示。图11展示了正方体再入目标的3种可能姿态：面指向、棱指向和角指向。通过计算目标在这些姿态情况的阻力系数，求其加权平均值作为正方体目标的名义气动阻力系数。在连续流-过渡流交界处，正方体处于3种不同姿态情况下表面的压力分布如图12所示。
图 10 正方体的网格划分

Fig.10 Mesh division of a cube

图 11 正方体目标的三种姿态情况

Fig.11 Three attitudes of a cube

图 12 连续流-过渡流分界处，正方体不同姿态时表面的压力分布

Fig.12 Pressure distributions on the cubic surface with different attitudes at the continuous-transition flow boundary

图13显示，正方体的阻力系数在自由分子流区域内随轨道高度降低略微减小；而在过渡流领域内，阻力系数则从略大于3快速降低至略小于2。
图 13 正方体的阻力系数随轨道高度的变化

Fig.13 The variation of a cube's drag coefficient with respect to the altitude

2.4 圆柱体

假设目标为长细比$l{/}d = 1$的圆柱形物体，平板法和DSMC方法采用的网格划分如图14所示。根据来流方向建立坐标系，使得来流方向在xOz面内，此时仅存在攻角，侧滑角为0。由于目标姿态未知，计算圆柱在攻角分别为0°、15°、30°、45°、60°、75°和90°时的阻力系数，取其平均值作为目标的名义气动阻力系数。图15展示了连续流-过渡流交界处，与来流平行的xOz面上，来流在不同攻角下圆柱体的密度分布情况。
图 14 圆柱的网格划分

Fig.14 Mesh division of a cylinder

图 15 连续流-过渡流分界处，来流在圆柱体攻角不同时的密度分布

Fig.15 Density distribution of incoming flow at the continuous-transition flow boundary with different angles of attack for a cylinder

图16展示了圆柱阻力系数随轨道高度的变化情况。在过渡流区域内，圆柱体的阻力系数随轨道高度的降低而迅速减小；而在自由分子流区域内，其阻力系数呈现轻微的减小趋势。
图 16 圆柱的阻力系数随轨道高度的变化

Fig.16 The variation of a cylinder's drag coefficient with respect to the altitude

2.5 小结

如图17所示，随着轨道高度的降低，各典型形状目标的阻力系数呈现减小趋势。其中，平板的阻力系数最小，球体次之，而正方体和圆柱体的阻力系数明显高于球体和平板。与经验值2.2相比，正方体和圆柱体的阻力系数总体偏大，平板则偏小。此外，球体在过渡流区域内随着轨道高度降低，与经验值的偏差逐渐增大。这证实了本文基于稀薄气体动力学理论研究阻力系数的必要性。
图 17 典型形状目标的阻力系数

Fig.17 Drag coefficients of typical shaped objects at different altitudes

3 阻力系数的其他影响因素

目标在再入大气过程中通常会跨越不同的经纬度区域，并且不同时期的太阳和地磁活动强度也不同。这些因素的变化会导致来流参数发生变化。因此，本节以第2.1节的球形目标为例，分析轨道经纬度、太阳和地磁活动强度对再入目标阻力系数的影响。

3.1 经纬度

3.1.1 经度

图18是纬度为40°N，而经度分别为0°E、60°E、120°E和180°E时，球的阻力系数随轨道高度变化的曲线。经度不同会导致来流的温度、密度以及各组分的占比不同，从而间接影响目标的阻力系数。图中可以看出，在自由分子流区域内，目标处于不同经度时阻力系数的差异小于0.1%，但120°E的过渡流区域内，轨道高度110 km处的偏差接近4.0%。
图 18 在纬度40°N、不同经度时，球的阻力系数随轨道高度的变化

Fig.18 The variation of a sphere's drag coefficient with respect to the altitude at latitude 40°N with different longitudes

3.1.2 纬度

当经度为120°E，纬度分别为0°、30°N、60°N和90°N时，球的阻力系数随轨道高度变化的情况如图19所示。与经度类似，纬度的不同也会导致来流的参数的变化，进而间接影响目标的阻力系数。从图19可以看出，不同纬度的球形目标在自由分子流区域内的阻力系数差异较小，但在过渡流区域内的最大偏差超过2.5%，且偏差的最大值通常出现在轨道高度100~110 km处。
图 19 在经度120°E、不同纬度时，球的阻力系数随轨道高度的变化

Fig.19 The variation of a sphere's drag coefficient with respect to the altitude at longitude 120°E with different latitudes

3.1.3 经纬度

图20展示了轨道高度150 km处，球的阻力系数随经纬度变化的分布情况。从图20可以看出，球的目标阻力系数在低纬度的不规则带状区域内较大，在（150°W，0°）附近区域达到最大值，而两极区域的阻力系数偏小，尤其是在北极的高纬度区域更为明显。这是因为轨道高度150 km处，北极附近区域的大气温度显著升高。由式（11）和式（12），可得压力系数会随着大气温度的升高而减小，而剪切系数则变化不大，因此北极区域的阻力系数偏小。图21则显示了以（0°，0°）处阻力系数为参考的阻力系数偏差随经纬度变化的分布情况。结果表明，最大值与最小值之间的偏差略大于0.2%。轨道高度150 km处属于自由分子流区域，图21的计算结果与图18和图19的结果一致。
图 20 球在轨道高度150 km处的阻力系数

Fig.20 Drag coefficient of the sphere at 150 km altitude

图 21 轨道高度150 km处，球在不同经纬度的阻力系数与（0°，0°）处的偏差

Fig.21 Deviation of a sphere's drag coefficient at 150 km altitude with respect to the value at (0°, 0°)

3.2 太阳和地磁活动

太阳活动强度呈现周期性变化，其平均周期约为11年。大气的温度和密度通常会受到太阳活动的影响也会呈现周期性变化。此外，地磁活动的强度同样受太阳活动的影响。太阳活动和地磁活动的强度通常分别使用F10.7和Ap两个指数来描述，具体取值如表2所示。在200 km高度上，高太阳活动期间的大气密度大约是低太阳活动期间的3.5倍。
表 2 太阳和地磁活动强度参数设置

Table 2 Parameters for low, mean and high solar and geomagnetic activities

太阳和地磁活动强度 F10.7/SFU Ap/2nT
低（low） 70 1
一般（mean） 125 19
高（high） 250 300
图22展示了不同太阳和地磁活动强度下，球体阻力系数随轨道高度的变化及其偏差曲线。结果显示，在太阳和地磁活动活跃期，阻力系数相较于一般情况会偏大，而在低太阳和地磁活动期间则偏小。太阳和地磁活动对自由分子流区域的阻力系数影响较小，其偏差小于0.5%；而在过渡流区域，影响较大，最大偏差可达6.0%。太阳活动较为活跃的时期，同一轨道高度上的大气密度增大，温度升高。由式（10）可知，分子速率比$s$减小，且$ s \propto T_\infty ^{ - 0.5} $，将其代入式（11）和式（12），可得压力系数和剪切系数与$s$负相关关系。因此，$s$减小，导致阻力系数增大。
图 22 不同太阳和地磁活动强度下，球体阻力系数随轨道高度的变化

Fig.22 The variation of a sphere's drag coefficient with respect to the altitude in cases of low, mean, and high solar and geomagnetic activities

4 结 语

针对目前再入预报中广泛使用的大气阻力系数经验值与实际气动阻力系数差异较大的问题，本文基于稀薄气体动力学理论，研究了球、平板、正方体和圆柱体等典型形状的再入目标。利用自由分子流气动计算的平板法和应用于过渡流区域的Wilmoth桥接公式，计算了不同姿态下的阻力系数并取平均值，以此作为再入目标的名义阻力系数。同时，比较分析了典型形状的阻力系数随轨道高度的变化情况，研究了经纬度、太阳和地磁活动对目标阻力系数的影响。以下是研究的主要结论：
（1）随着轨道高度的降低，各典型形状目标的阻力系数均呈现单调递减的趋势，且在过渡流区域内的减小速度高于在自由分子流区域。
（2）与阻力系数的经验值相比，除了球形目标在自由分子流区域内的阻力系数外，几种典型形状的阻力系数与经验值2.2存在明显差异。例如，在轨道高度87~200 km内，平板的阻力系数仅为经验值2.2的43%~60%。如果采用经验值2.2，正方体和圆柱体在自由分子流区域内的阻力系数的误差约为25%。此外，球形目标的阻力系数随轨道高度降低，从不足5%逐渐增大至76%，这验证了基于气动力模型的阻力系数研究的必要性。
（3）目标所在经纬度的变化会导致来流的温度、密度及各成分占比变化，从而导致再入目标的阻力系数出现差异。在自由分子流区域，这种偏差小于0.3%，而在过渡流区域，最大偏差接近4.0%。
（4）太阳和地磁活动也会对大气参数产生影响，导致目标阻力系数在太阳和地磁活动活跃期偏大，而在安静期偏小。类似地，在过渡流区域的偏差可达6.0%，而在自由分子流区域的偏差不足0.5%。
基于本文的研究结果，在空间目标再入预报中，利用目标的几何外形参数甚至姿态运动信息，可以更加准确地计算目标气动阻力，从而预期提高轨道寿命估计和再入预报的精度。

References

Publishing order | Descend order by publishing year | Descend order by cited within
1
MOE K, MOE M M. Gas–surface interactions and satellite drag coefficients[J]. Planetary and Space Science, 2005, 53, 793- 801.

DOI

2
LIU J H, XU J, DU J L, et al. Application of atmospheric drag model based on Reynolds number in reentry prediction of rocket bodies[J]. Chinese Journal of Space Science, 2022, 42 (2): 277- 283.

DOI

3
NWANKWO V U J, DENIG W, CHAKRABARTI S K, et al. Atmospheric drag effects on modelled low Earth orbit （LEO） satellites during the July 2000 Bastille Day event in contrast to an interval of geomagnetically quiet conditions[J]. Annales Geophysicae, 2021, 39 (3): 397- 412.

DOI

4
MEHTA P M, WALKER A, MCLAUGHLIN C A, et al. Comparing physical drag coefficients computed using different gas-surface interaction models[J]. Journal of Spacecraft and Rockets, 2014, 51 (3): 873- 883.

DOI

5
BULLARD J. Satellite Drag Analysis Using Direct Simulation Monte Carlo (DSMC)[D]. Herfordshire:University of Herfordshire,2018.

6
CARNA S F R, BEVILACQUA R. High fidelity model for the atmospheric re-entry of CubeSats equipped with the drag de-orbit device[J]. Acta Astronautica, 2018, 156, 134- 156.

7
JIN X H, HUANG F, CHENG X L, et al. Atmospheric drag on satellites flying in lower low-earth orbit[J]. Journal of Tsinghua University: Science and Technology, 2020, 60 (3): 219- 226.

8
CUI W H,QU W,LI H,et al. DTM2013 model parameter inversion and correlation analysis between its accuracy[C]. China Satellite Navigation Conference(CSNC) 2020 Proceedings. Singapore:Springer,2020:36-46.

9
HUANG J, CHANG L, DONG B Y, et al. Development on space environment and its dynamic and thermal problems of Ultra-LEO satellites[J]. Chinese Journal of Space Science, 2023, 43 (4): 711- 723.

DOI

10
HART K A, SIMONIS K R, STEINFELDT B A, et al. Analytic free-molecular aerodynamics for rapid propagation of resident space objects[J]. Journal of Spacecraft and Rockets, 2018, 55 (1): 27- 36.

DOI

11
MASLOVA A I, PIROZHENKO A V. Modeling of the aerodynamic moment acting upon a satellite[J]. Cosmic Research, 2010, 48 (4): 362- 370.

DOI

12
PRIETO D M, GRAZIANO B P, ROBERTS P C E. Spacecraft drag modelling[J]. Progress in Aerospace Sciences, 2014, 64, 56- 65.

DOI

13
LIVADIOTTI S, CRISP N H, ROBERTS P C E, et al. A review of gas-surface interaction models for orbital aerodynamics applications[J]. Progress in Aerospace Sciences, 2020, 119, 100675.

DOI

14
PADILLA J F, BOYD I D. Assessment of gas-surface interaction models for computation of rarefied hypersonic flow[J]. Journal of Thermophysics and Heat Transfer, 2009, 23 (1): 96- 105.

DOI

15
WALKER A, MEHTA P, KOLLER J. Drag coefficient model using the Cercignani-Lampis-Lord gas-surface interaction model[J]. Journal of Spacecraft and Rockets, 2014, 51 (5): 1544- 1563.

DOI

16
PILINSKI M D, ARGROW B M, PALO S E. Semiempirical model for satellite energy-accommodation coefficients[J]. Journal of Spacecraft and Rockets, 2015, 47 (6): 951- 956.

17
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.

DOI

18
BIRD G A. Molecular gas Dynamics and the Direct Simulation of Gas Flows[M]. Oxford:Clarendon Press,1994.

Options
Outlines

/

〈 〉