由于重力场和稳态海洋环流精确测量的需要，超低地球轨道(200~500 km)航天器逐渐展现出广阔的应用前景。欧洲空间局(European Space Agency，ESA)分别于本世纪初相继发射了GRACE (Gravity Recovery and Climate Experiment)、CHAMP (CHAllenging Mission Payload)和GOCE (Gravity field and steady state Ocean Circulation Explorer)卫星^[1-3]。实现重力梯度精确测量的先决条件是为卫星提供一套无拖拽姿态控制系统(drag-free and attitude control, DFAC)，该系统采用一套离子发动机和磁力矩器来补偿大气阻力等非重力的影响，从而进行姿态控制^[4]，而这需要卫星大气阻力的实时准确预测。

此外，随着各国航天发射任务的增加，空间碎片逐年递增，计划中和正在进行的空间飞行任务不得不把与太空中较大块碎片的碰撞可能性纳入考虑范围。而且，较大碎片失控脱离轨道和意外进入地球大气层也对地面的生命和财产安全造成了巨大威胁^[5]。在超低轨轨道环境下，大气阻力是估计卫星轨道的最大误差来源^[6]，提高阻力计算精度，可以提高卫星和空间碎片轨道的预测精度，从而避免在轨卫星与空间碎片或计划发射卫星的碰撞。

Titov等^[7]把影响卫星和空间碎片大气阻力的所有因素划分为3类：卫星几何参数、卫星表面属性和大气参数。卫星阻力表达式为^[8]

$ D = {C_D} \cdot \left( {\frac{1}{2}\rho {{\left| {{\mathit{\boldsymbol{v}}_{{\rm{rel}}}}} \right|}^2}} \right) \cdot A. $

(1)

其中：v_rel为卫星相对于大气的速度，属于卫星几何参数；A为参考面积(通常为卫星沿-v_rel方向的投影面积)，也属于卫星几何参数；ρ为大气密度，属于大气参数；C_D为阻力系数，主要受卫星几何参数和卫星表面属性的影响，也受大气参数的影响。

根据阻力关系式(1)，给定卫星相对速度v_rel和投影面积A，卫星阻力的预测归结到阻力系数和大气密度的计算。阻力系数的计算属于自由分子流空气动力学的范畴，前人已经进行过大量的研究。Prieto等^[9]认为复杂外形航天器阻力系数的计算方法主要有4种：面元积分法^[10]、射线追踪面元积分法^[11]、直接模拟Monte Carlo (direct simulation Monte Carlo, DSMC)方法^[12]和试验粒子Monte Carlo (test particle Monte Carlo, TPMC)方法^[13-16]。其中，TPMC方法能准确模拟复杂外形导致的流动遮挡和多次反射现象，且计算速度较快，存储要求较低，是低轨卫星气动特性预测的理想方法。

影响卫星阻力预测的另外一个主要因素是大气参数，尤其是大气密度。大气参数建模的主要过程是基于大量卫星测量数据集，获得大气密度、温度和组分相对轨道位置、太阳辐射和地磁活动的变化^[17]。在航天工程领域，大气模型为卫星和空间碎片的气动阻力预测提供大气密度、温度和组分等大气参数^[18]。目前，广泛使用的大气模型有美国标准大气(US Standard Atmosphere, 以1976版即USSA-1976^[19]为代表)、Jacchia系列大气模型(以Jacchia-1977^[20]为代表)、质谱仪非相关散射模型(以NRLMSIS-00^[21]为代表)、卫星阻力总密度模型(以JB2008^[22]为代表)和阻力温度模型(以DTM-2013^[23]为代表)。

20世纪80年代，Gaposchkin等^[24]对各种大气模型进行了对比分析，认为没有一种大气模型能够满足所有卫星型号或发射任务。时至今日，这种状况并没有得到有效改善，只不过有了更多的大气模型可以选择。国内外针对卫星大气阻力的研究大多只给出了无量纲阻力系数的结果^[17]，并且只是局限于某种大气模型。比如，尹凡等^[25]采用JB2008模型，联合考虑地球扁率和大气阻力摄动对CHAMP卫星的影响，模拟考察大气阻力引起的轨道高度衰减。刘卫等^[26]采用NRLMSISE-00大气模型，从小波分解的角度列出天宫一号卫星阻力系数与空间环境参量变化关系。汪宏波等^[27]在分析低轨卫星轨道预报误差发散规律的基础上，提出了一种大气阻力系数计算方法。

根据阻力关系式(1)，大气密度在确定阻力时与阻力系数同等重要，考虑到目前大气模型的诸多假设和半经验性质，对各种大气模型在卫星阻力预测时的表现进行对比研究显得更为迫切。因此，采用多种大气模型计算卫星大气阻力并进行对比分析，一方面可以为卫星的总体设计提供关键数据，对于卫星的外形和轨道设计尤其具有指导性意义；另一方面则能根据阻力预测结果对大气模型进行改进。

本文在自由分子流TPMC方法^[16]的基础上，通过构建对应的接口程序，嵌入5种国际主流的大气模型，包括USSA-1976^[19]、Jacchia-1977^[20]、NRLMSIS-00^[21]、JB2008^[22]和DTM-2013^[23]，发展了一套低地球轨道(low-earth orbit, LEO)航天器气动特性预测的通用的三维并行软件APSL (aerodynamic properties of satellites in LEO)，实现了任意复杂外形航天器气动力和力矩的快速准确预测。

1 APSL软件简介

APSL软件是一套低地球轨道航天器气动特性预测软件，可以根据航天器外形自动构建计算域，实现试验粒子的初始化和高效跟踪，包括了当下主流的大气模型和气体与表面相互作用(gas-surface interaction, GSI)模型，具备大规模并行处理的能力。此外，APSL软件还能沿着轨道以批处理的方式对航天器的气动特性进行预测，实现结果的可视化输出。该软件除气动计算模块之外，还包括前处理和后处理模块，如图 1所示。气动计算模块是APSL软件的核心模块，主要由自由分子流TPMC方法和大气模型构成。

1.1 TPMC方法

Davis^[14]提出TPMC方法并成功用于计算通过分子平均自由程远大于圆管特征尺度条件下的圆管分子流率。由于该方法忽略了气体分子之间的碰撞，因此只适用于自由分子流动。与DSMC方法相比，TPMC方法的特点是仿真分子是顺序而非同时产生的，一次只产生一个试验粒子，因此不会耗费太多的计算时间和存储空间，更适用于复杂边界导致的多次表面反射流动问题^[28]。TPMC方法不但具备准确模拟流动遮挡和多次反射效应的能力，还能直观地实现多种GSI模型，以便更准确地模拟气体分子与航天器表面之间的动量和能量交换^[29]。除了前述真空技术中的管道流率计算之外，TPMC方法已经成功应用于航天器表面出气导致的返回流污染问题^[30-31]和自由分子流区空气动力学问题^[15-16]。

TPMC方法的基本思路是构建一个包含航天器的计算域，通过跟踪和模拟计算域内部的试验粒子来模拟气体流动。每个试验粒子代表大量的真实气体分子，其初始速度为确定性质的来流速度与概率分布性质的热运动速度的矢量和^[29]。试验粒子逐个进入计算域，在运动过程中，要么在航天器表面发生碰撞并反射，要么直接飞出计算域。试验粒子在航天器发生反射时，按照给定的GSI模型与表面交换动量和能量。在试验粒子的采样数量大到足以保证真实物理过程的特征均得到准确模拟后，统计出气动力、热等空气动力学宏观量。

1.2 大气模型

美国标准大气1976(USSA-1976)^[19]作为USSA-1962的修订版，为满足人们对上层大气知识日益增加的需求而产生。在50 km海拔高度以上，该模型基于当时大量的火箭数据和大气中间层、热层底部的理论，以及超过一个完整太阳循环期间巨大的卫星观测数据。在50 km以下，该标准与ISO Standard (1973)一致。USSA-1976和早期的大气模型一般认为卫星的阻力系数近似为2.2，且独立于飞行高度，利用阻力数据，根据阻力和阻力系数的关系，给出密度的值。阻力数据来源于卫星所携带加速度计的测量或者轨道衰变的观测，因此，该模型给出的大气密度、温度和组分并不十分准确。

Jacchia^[20]推导出一系列多参数大气模型(Jacchia-1964、Jacchia-1970、Jacchia-1971和Jacchia-1977模型)，描述时间(世界时，即Greenwich时间)、空间位置和太阳活动、地磁干扰等物理现象对大气的影响。基于卫星阻力数据和大气组分、温度测量数据，这些模型可以给出大气密度、温度、压力和组分。同样，这一系列模型也假定阻力系数为2.2，且独立于飞行高度。通过引入更真实的阻力系数，Chao等^[32]通过研究发现，飞行高度大于250 km时，Jacchia-1971模型给出的大气密度偏高23%。

NRLMSISE-00模型^[21]是(美国)海军研究实验室(Naval Research Laboratory，NRL)在质谱仪非相干散射雷达(mass spectrometer incoherent scatter radar，MSIS)模型的基础上开发的模型。作为MSISE-90大气模型的升级版，新的NRLMSISE-00大气模型较好地结合了MSISE-90和Jacchia-1970两种模型的优点，覆盖了地球表面到外大气层底部的大气。相比美国标准大气1976和Jacchia-1977模型，该模型能够考虑时间、地理纬度、地理经度对大气参数的影响。

JB2008模型^[22]是JB2006的改进升级版，而后者的理论基础是大气动力学的Jacchia扩散方程。值得一提的是，JB2008模型的密度数据来源比较新，既包括2001—2005年期间CHAMP卫星的加速度密度数据，也包括2002—2005年GRACE卫星的加速度密度数据。该模型的输入参数除了前述的时间(世界时)、空间位置(包括海拔高度、地理纬度、地理经度)之外，还可以考虑不同年份对大气参数的影响。

DTM模型^[23]的全称是阻力温度模型，目前公开的最新版本是DTM-2013模型。作为描述地球大气温度、密度和组分的一个半经验模型，DTM-2013能够考虑时间(世界时)、海拔高度、地理纬度、地理经度、太阳活动和地磁干扰对大气参数的影响。开发该模型的初衷就是为了精确地预测卫星阻力和轨道，而且该模型的构建采用了通过CHAMP、GRACE和GOCE等卫星所测量加速度数据而推演出的密度数据。据构建该模型的研究者宣称^[23]，由于增加了GOCE卫星的测量数据，总体上DTM-2013是目前最精确的大气模型，在大多数高度下，其数据相对于DTM-2009有了大幅度的改善，且比JB2008精确。

2 方法校验 2.1 算例描述

GOCE卫星是世界上第一颗重力梯度测量卫星^[3]。GOCE卫星主要由卫星本体、太阳翼和稳定翼构成，本体呈八棱柱体形状，横截面是八边形，太阳翼面积较大，稳定翼面积较小，如图 2所示。图中的坐标系OXYZ原点位于卫星质心，X轴沿卫星轴向指向尾部，Y轴位于太阳翼平面，Z轴位于稳定翼平面。来流大气密度、温度和平均相对分子质量采用USSA-1976确定，具体的飞行和计算条件见表 1。

2.2 计算结果

在算例校验之前，先给出物理量的定义。在此定义速度比为卫星飞行速度v_∞与来流最可几分子热运动速度c_m的比值，即

$ s = {v_\infty }/{c_{\rm{m}}}, $

(2)

$ {c_{\rm{m}}} = \sqrt {2k{T_\infty }/m} . $

(3)

其中：k为Boltzmann常量，T_∞为大气来流温度，m为来流大气的平均相对分子质量。

显然，给定大气组分和温度，速度比与飞行速度成正比。阻力系数是表征阻力大小的无量纲系数，定义为阻力与来流动压的比值，即

$ {C_D} = D/\left( {{\rho _\infty }v_\infty ^2/2} \right). $

(4)

其中：D为阻力，ρ_∞为来流大气质量密度。

图 3为GOCE卫星在250 km高度阻力系数随着来流速度比的变化曲线，考虑2种GSI模型：完全漫反射(Maxwell适应系数σ=1.0)和部分漫反射部分镜面反射即Maxwell边界条件(σ=0.6，即60%的气体分子漫反射，40%的气体分子镜面反射)。显然，对于2种GSI模型而言，阻力系数都随着来流速度比(飞行速度)的增加而单调减小，并渐进地趋于稳定。此外，无论采用完全漫反射还是部分漫反射模型，本文计算结果都与Koppenwallner等^[33]发表的结果一致，验证了方法的可靠性和准确性。

3 GOCE卫星阻力影响因素

完成方法验证之后，研究飞行高度、纬度和经度等因素对GOCE卫星阻力的影响，重点考察阻力相对大气模型的敏感性。飞行时间为2009年5月1日0点，GSI采用完全漫反射模型，若无特殊说明，其他大气来流参数和计算条件如表 2所示。

3.1 飞行海拔高度

首选，考虑飞行高度对卫星大气阻力的影响，图 4为GOCE卫星大气阻力和阻力系数随着高度的变化曲线图，其中，飞行高度的变化范围是H∈[120,600] km。同时，为了对比分析，还给出了大气密度相对高度的变化曲线。

如图 4a所示，随着高度的增加，由于大气密度的急剧降低，阻力也随之呈指数型减小。通过对比图 4a和4c还可以看出，当飞行高度变化范围较大时，大气阻力主要受大气密度的影响，其他因素基本可以忽略不计。当然，大多数地球低轨道卫星的飞行高度变化范围不会超过100 km，所以其他因素的影响也得考虑。阻力系数相对于高度的变化则呈现不同的特征：随着高度的增加，阻力系数单调增大，这与大气密度相对于高度的变化明显不一致。原因是阻力随高度减小的速度比密度随高度减小的速度要慢。

在较低飞行高度(低于200 km)下，不同大气模型给出的大气密度、温度和组分差异较小，因此采用5种大气模型计算出的大气阻力和阻力系数近似一致。随着高度的升高，不同大气模型给出的大气密度、温度和组分差异增大，导致大气阻力和阻力系数的差异随之变大。这表明，随着飞行高度的升高，卫星大气阻力相对大气模型的敏感性增强。

就大气阻力和阻力系数而言，前3种大气模型即USSA-1976、Jacchia-1977和NRLMSISE-00表现相近；在合理的误差范围，后2种大气模型即JB2008和DTM-2013表现基本一致。对于大气阻力，采用前3种大气模型的计算结果较后2种的偏高；对于阻力系数，在较低高度(低于440 km)下，前3种大气模型的计算结果较高，在较高高度下(高于440 km)，后2种大气模型的计算结果较高。

3.2 轨道的地理纬度

卫星沿着轨道飞行时一般会跨越不同的纬度和经度，而轨道纬度或经度的变化会引起来流参数的变化，因此有必要分析轨道纬度和经度对卫星阻力的影响。考虑到前3种大气模型产生相似的结果，USSA-1976模型和Jacchia-1977程序包只考虑了高度对大气密度等参数的影响，故仅采用前3种模型中最新的NRLMSISE-00模型(该模型也是美国海军研究实验室推荐使用的大气模型)和后2种大气模型，即JB2008和DTM-2013进行计算。

图 5为250 km高度下GOCE卫星大气阻力、阻力系数和大气密度随轨道纬度的变化图，纬度W∈[-90°, 90°]。随着纬度的变化，卫星阻力、阻力系数和来流大气密度均呈现非单调的波动性变化。对于前3种大气模型而言，阻力和阻力系数均在北极地区取得最大值，而最小值位置则与大气模型有关，不同大气模型的最小值位置不同。整体而言，卫星位于北半球上空时受到的大气阻力和阻力系数都较位于南半球上空时的大。轨道纬度变化导致来流大气密度、温度和组分变化，从而影响卫星阻力和阻力系数。来流大气温度和组分影响阻力系数，而阻力系数和来流大气密度共同决定卫星阻力。当卫星跨越不同纬度时，来流大气温度和组分的变化导致阻力系数产生一定的波动，因此卫星阻力随纬度的变化曲线与大气密度随纬度的变化差异较小。总体而言，只有轨道纬度变化而其他参数恒定时，卫星阻力和阻力系数的变化都不太明显。就大气阻力和阻力系数而言，JB2008和DTM-2013这2种大气模型表现相近，而NRLMSISE-00模型计算出的阻力和阻力系数都明显较前两者偏高。

3.3 轨道的地理经度

研究轨道经度对卫星大气阻力和阻力系数的影响，同样仅利用NRLMSISE-00、JB2008和DTM-2013进行计算。与轨道纬度类似，轨道经度也是通过影响来流大气参数而间接影响阻力和阻力系数。

图 6是250 km高度下GOCE卫星大气阻力、阻力系数和大气密度随轨道经度的变化，经度的变化范围是J∈[-180°, 180°]。随着轨道经度的变化，3种大气模型计算出的卫星阻力、阻力系数和大气密度同样呈现非单调的波动变化特征，且变化规律相似，阻力系数的波动性最明显。整体而言，卫星位于较高经度地区上空时受到的大气阻力和阻力系数较处于较低经度地区上空时的大。3种大气模型的阻力和大气密度均在W150°附近取得最大值。对于NRLMSISE-00和DTM-2013模型，阻力系数在W120°附近达到最大；对于JB2008模型，阻力系数则大约在W/E180°处取得最大值。

需要注意的是，卫星阻力和阻力系数随着轨道经度的相对变化量明显比其随着轨道纬度的相对变化量更大，这表明沿着东西方向飞行的卫星受到阻力的波动较沿着南北方向飞行的大。同样，就大气阻力和阻力系数而言，JB2008和DTM-2013两种大气模型表现相近，而NRLMSISE-00模型计算出的阻力和阻力系数虽然规律与前两者类似，却明显较前两者偏高。

4 结论

本文在自由分子流TPMC方法的基础上，通过构建对应的接口程序，嵌入5种国际主流大气模型，发展了一套低地球轨道任意复杂外形航天器气动特性预测的通用的三维并行软件APSL。简要介绍TPMC方法和5种大气模型之后，以GOCE卫星为研究对象，计算该卫星的大气阻力，比较分析了大气模型、飞行高度、轨道纬度和经度对大气阻力的影响，重点考察大气阻力对大气模型的敏感性。结论如下：

1) 随着高度的增加，由于大气密度的急剧降低，阻力也随之呈指数减小，阻力系数却单调增大，卫星阻力对大气模型的敏感性增强；当飞行高度变化范围较大时，在卫星飞行速度保持不变的条件下，大气阻力主要受大气密度的影响，其他因素基本可以忽略不计。

2) 随着轨道地理纬度的变化，卫星阻力、阻力系数和来流大气密度均呈现非单调的波动性变化，在本文所考虑的计算条件下，阻力和阻力系数均在北极地区取得最大值，卫星位于北半球上空时受到的大气阻力和阻力系数都较位于南半球上空时的大。

3) 随着轨道经度的变化，卫星阻力、阻力系数和大气密度呈现非单调的波动变化特征，且变化规律相似，在本文所考虑的计算条件下，卫星位于较高经度地区上空时受到的大气阻力和阻力系数较处于较低经度地区上空时的大。

4) 就卫星阻力和阻力系数而言，JB2008和DTM-2013两种大气模型表现相近，而前3种大气模型即USSA-1976、Jacchia-1977和NRLMSISE-00计算出的阻力都明显较前两者偏高。