跨流区超低轨航天器快速气动力计算方法

郭晨林，陈方，赵艳彬，尤超蓝，钱勇，李文龙，杨丽丽

（1.上海交通大学 航空航天学院，上海 200240；2.上海卫星工程研究所研发中心，上海 201109）

0 引言

超低轨飞行器由于其重访周期短、对地成像观测性能强的特点，在通信、对地探测等领域展现出了巨大的潜力［1］，但超低轨飞行器受到的气动摄动大，造成轨道维持与姿态控制的困难。同时，由于其运行高度在120～300 km 之间，属于过渡流区与自由分子流区，对于气动力的估算缺少经验理论与工具，造成了总体设计与动力选型上的困难。因此，当前急需有关超低轨卫星的快速气动力计算工具，以实现合理的动力选型，延长卫星工作寿命［2］与加速超低轨飞行器的初步设计迭代。

对于工作在自由分子流区中的飞行器的气动力特性求解，BIRD［3］提出的直接模拟蒙特卡洛法（Direct Simulation Monte Carlo，DSMC）求解分子之间的碰撞效应高精度求解其气动力特性，DAVIS［4］提出的测试粒子蒙特卡洛法（Test Particle Monte Carlo，TPMC）省略了分子间的碰撞效应，虽然降低求解精度，但是大幅减小了计算量，提升了计算效率。Maxwell 通过建立简单的面元与气体粒子的碰撞模型，求解自由分子流区中的面元受力状况，进而积分求解飞行器受力［5］。该方法成为自由分子流区中快速气动力计算的基本方法，后续许多学者根据分子与面元碰撞时的动量损失与漫反射效应对Maxwell 提出的模型进行优化，扩展了该模型的精度与适用范围［6-7］。

在过渡流区中，气体的运动规律更为复杂，难以抽象出简化的物理模型。当前求解过渡流区中的流动特性多依赖于DSMC 方法［8］。由于在该区域中大气密度较高，意味着DSMC 模拟的分子数量多，所需要的计算量大，因此，求解过渡流区气动特性的方法为：首先对连续流区的气动特性与自由分子流区中的气动特性进行求解，再使用桥函数加权来近似过渡流区的气动特性。常用的桥函数有sine-square 桥函数［9］与erf-log 桥函数［10］，通过假定多个待定系数函数拟合连续流区-过渡流区-自由分子流区的部分计算结果，实现对其他工况的结果进行估计。以上桥函数的拟合只对特定的气流来流条件与当前几何外形的飞行器适用，当外部条件发生变化时，需要新的DSMC 或试验结果进行标定。

目前，过渡流区的计算高度依赖于DSMC，TPMC 等高精度计算工具，快速气动力计算工具十分有限，常用的气动估算方法存在一定的局限性。同时，面元法应用的前景广泛，通过对面元法进行进一步改进，以上问题能被较好的解决。本文旨在解决超低轨飞行器的快速气动计算需求，提出一套适用于超低轨卫星从过渡流区到自由流区的快速气动估算方法。在面元法的基础上，对迎风单元筛选算法与桥函数进行了设计，达到了提高精度，增加适用范围的效果。

1 超低轨卫星快速气动力计算方法

本文提出的超低轨卫星快速气动力计算方法主要分为3 个部分，如图1 所示。包括输入数据、数据的预处理方式和气动力计算方法。其中，有2 个过程显著影响计算精度和计算效率：①数据预处理与迎风单元与背风单元的划分算法；②气动力的计算方法和桥函数模型。本文将对此进行较详细的阐述。

图1 超低轨卫星快速气动力计算方法流程Fig.1 Flow chart of the rapid aerodynamic calculation method for ultra-low orbit spacecrafts

1.1 数据预处理

在输入的模型STL 格式文件中，包含着卫星的几何尺寸信息，首先需要对其表面单元进行离散，得到其表面的三角形单元如图2 所示。在建模的过程中，常以飞行器的体坐标系进行建模，为方便后续的迎风背风单元算法及气动力计算，需要对离散单元在风轴系和体轴系下进行自由转换。通过以下旋转矩阵可以实现转换过程：

图2 迎风单元与背风单元筛选Fig.2 Screening diagram of windward and leeward units

式中：Cbo为体坐标系ObXbYbZb至风轴坐标系OwXwYwZw下的旋转矩阵；α、β分别为当前飞行器的迎角、侧滑角。

进行坐标转换后，需对表面的单元进行分类，分为迎风单元与背风单元。根据单元属性的不同，在后续的气动力计算中使用不同的计算方式。本文设计了可用于非凸表面的迎风单元筛选算法。通过验证发现，随着网格密度达到一定要求时，其具有较高的计算精度，满足对复杂外形的超低轨卫星气动力估算的需求。迎风单元和背风单元筛选算法的输入包含如下数据：

1）节点位置Pi=[xi，yi，zi]，为编号为i节点其对应在风轴系下的空间坐标位置；

2）单元列表Un=[u，v，w]，n为单元编号，u，v，w为三角形单元的3 个顶点节点编号。

算法的输出包含以下几个部分：

1）迎风单元列表Wn=[u'，v'，w']，u'，v'，w'为三角形单元的3 个顶点节点编号；

2）背风单元列表Ln=[u'，v'，w']，u'，v'，w'为三角形单元的3 个顶点节点的编号；

3）迎风节点列表PWi=[xi，yi，zi]，xi，yi，zi为节点在风轴系下的坐标位置。

具体算法见表1。

表1 迎风单元筛选算法流程Tab.1 Flow chart of the screening algorithm for windward units

从迎风方向依次选取节点，判断其是否在迎风单元列表的投影内，若不在任意投影内，则将该节点加入迎风节点列表中，并将与其连接的单元加入迎风单元列表。最后将所有总单元列表Un减去筛选出的迎风单元列表，则可得到背风单元列表。具体筛选过程如图2 所示，在对表面进行离散后的卫星网格中，黑色单元为筛选出的迎风单元，灰色的是待筛选的单元列表。红色节点与蓝色节点为两个待筛选节点。在图示的沿着来流方向的投影平面上，红色节点在黑色单元投影外，所以其为新的迎风单元；而蓝色节点在当黑色单元的投影内，其为背风节点。

为验证该算法的准确性，以及提出在后续计算中划分网格的精细程度，本文设计了如下算例，如图3（a）所示，该图所示的梨状外形，为一个非凸的几何结构，其在轴向投影方向的最大投影半径R=0.5 以此作为模型的特征尺度，因此可以得到其准确的投影面积Sproj=πR2=0.25π。本文依次增加特征尺度与三角形网格大小R/δx之间的比值，观察通过迎风单元与背风单元算法计算出的投影面积与实际投影面积的差距，验证算法精度及求解合适的特征尺度与网格大小的比值关系R/δx。计算的结果如图3（b）所示，发现在比值为10 时，算法计算出的投影面积与实际的投影面积十分接近，进一步缩小网格尺度，增加特征尺度与网格大小的比值，精度不再改善。因此，本文认为该算法的准确性良好，在对后续的计算模型进行表面单元的离散过程中，将选取特征尺度与网格大小的比值R/δx=10。

图3 迎风单元与背风单元筛选算例模型Fig.3 Example model for screening windward and leeward units

1.2 气动力计算

在气动力的计算过程中，需要从大气模型中根据低轨卫星的飞行高度和经纬度获得当前的流动条件，获得准确的流动条件对算法的精度与可靠性也有较大的影响。因此，选取合适的大气模型也至关重要。当前常用的大气模型有以下几种，US Standard 1976、MSISe-90、MSISe-00 和Jacchia-71等［11-13］。本文选取NRLMSISE-00 作为计算使用的大气模型，主要基于其拥有以下良好性质［14］：①使用广泛，有丰富的程序接口可调用；②能够实现努森数（Knudsen Number，Kn）与海拔高度的转换；③其数据丰富，满足本文的计算需求。

本文选取的模型计算条件见表2。

表2 大气模型参数Tab.2 Parameters of the atmospheric model

为方便后续计算，本文固定部分大气模型参数（见表2），只保留高度ℎ 作为变量。以上数据取太阳丰度为平均值时的大气系数，能反映当前的大气较为普遍的状态。选取以上参数后，本文假设大气模型为函数Fa，其输入变量为高度ℎ，则其可输出Kn、T、m、ρ分别为当前高度下的努森数、大气温度、大气的摩尔质量与大气密度，即：

根据不同高度对应的流区的不同，本文将计算分为3 个部分，见表3。

表3 对应不同流区适用的不同气动力计算方式Tab.3 Aerodynamic calculation methods for different flow regions

1.2.1 自由分子流模型

在自由分子流中，分子之间的碰撞概率很小，经碰撞后反射分子对流动的影响可以忽略不计［6］。因此，主要关注的过程是气体分子与卫星表面发生碰撞（见图4），分子与卫星表面发生能量与动量的交换情况。因此，自由分子流动时卫星的受力状况是由气体-表面相互作用性质决定的。在假设气体仅与卫星表面发生动量交换的前提下，作用在卫星表面上的力等于气体动量的变化率。对于一个外凸形状的卫星，其单位面积受力大小与气体的入射动量与反射动量的大小有关。以一个卫星的表面单元为例，其受到分子的作用力可表示为

图4 气体分子与单元表面碰撞模型Fig.4 Model for the collision of gas molecule and unit surface（a）Maxwell 模型（b）牛顿理论

式中：f为分子的作用力；A为卫星表面的单位面积；p为法向的动量通量；τ为切向的动量通量。下标a、b 分别对应着入射通量和反射通量，法向和切向入射动量流pa、τa取决于入射速度（V）和质量流（dQ）；其反射动量通量为-pbn+τbt；n为在物面表面的法向分量；t为在物面表面的切向分量。

在自由分子流中，可以认为

式中：dQ=ρVcosδ，δ为粒子入射方向与平面之间的夹角，ρ为分子流密度；V为分子流速度大小。

在流动中气体分子的速度服从Maxwell 速度分布函数［15］。则入射动量通量可以表示为

式中：m为单个气体分子的质量；F(u)为麦克斯韦分布函数。

式中：s为分子速度比；V为在速度分布函数下最可能的速度大小；Rc为气体常数；T∞为当前高度下的来流温度。

目前，如何求解确定反射通量是个难点。有许多采用简化的气体-表面作用模型来对反射动量通量进行建模，如Maxwell［15］、Schamberg、Schaaf&Chambre［6］等模型。本 文选取Sentman 模型［16］进行建模，假设一个粒子在撞击后会依据一定的概率分布向不同方向进行漫反射，定义分子与平面发生碰撞之后，平面的法向力系数与切向力系数分别为

式中：Tk，b、Tk，a为

式中：αacc为能量调节系数，根据飞行器的高度确定［16］，随着飞行高度的增加，能量调节系数的取值不断减小；Tw为飞行器的表面壁温。

1.2.2 连续流区模型

牛顿理论（牛顿正弦平方律）阐述了一种推导平面与流体相互作用的理论模型，其假设气体分子为质点且互相孤立，气体与平面撞击后，其动量在物面的法向上发生完全的动量交换，分子沿着物面的切向动量完全保留。其理论在低速流动条件下偏差较大，但物体在高速连续流区运动中，其对气动力的估算效果良好。物面受到的法向力系数CP为

牛顿理论提出的模型仅仅与气流方向与物面的夹角有关，在马赫数不够大时，会出现明显的偏差。为了对该偏差进行修正，Lester Lee 提出了新的计算方式：

其中：

式中：Ma∞为自由来流的马赫数。

在牛顿理论中，背风单元不受气动力的作用，但在连续流区中，忽略背风单元的气动力会导致计算出现明显的误差。因此，本文使用以下普朗特-迈耶膨胀波理论估算背风单元的气动系数［17］

式中：γ为绝热气体指数，在量热完全气体的假设下γ为1.4。

在实际估算中，以述方法能够较好地估计连续流区中的气动力系数值，但由于缺少黏性力的估算，依然会有一定的偏差，为做进一步修正，本文使用可压缩流的摩擦系数计算公式对摩擦力系数进行估算。

1.2.3 过渡流区模型

在过渡流区的气动力计算中，本文使用桥函数的方式进行计算，即在计算出连续流区与自由分子流区的气动特性参数后，根据在过渡流区中的高度h由大气模型推导出当前高度下的努森数Kn对连续流区与自由分子流区的气动特性参数进行权重分配，得到当前位置的气动特性参数。

式中：Cfm为在自由分子流中求得的气动力参数；Ccont为在连续流区中求得的气动力参数；Pb为桥函数。当前常用的桥函数有sine-squared 桥函数与log-erf 桥函数。

在设定合理的参数条件下，2 种函数均表现出良好的拟合性能。以上2 种函数都需要根据求解物体的几何外形的变化与来流条件的改变对相关参数进行修正。在实际使用过程中，需要使用部分实验数据或高精度的求解结果对参数进行调整，这增加了工作量。

为了对上述问题进行改善，本文提出了以logistic 函数为基础的桥函数，并在函数中加入了形状修正因子，拓展了桥函数的适用范围。本文设计的桥函数如下：

式中：q为形状因子，其与气体来流方向与物面间的夹角有关；桥函数为与当地Kn数和形状有关的函数。

2 模型验证

为验证提出的快速算法的有效性及正确性，本文结合DSMC 的计算结果，对以上方法进行验证。吴子牛等［10］曾对圆柱体与钝头双楔体进行DSMC与桥函数的研究，其研究对象如图5 所示，2 个模型展向上的长度均为1 000 mm。

图5 几何外形Fig.5 Geometric shape

计算条件见表4，给定模型壁温分别500 K、300 K；自由来流温度条件给定为300 K。选取圆柱模型的速度工况为马赫数4 与马赫数16，钝头双楔体模型速度工况为马赫数4 与马赫数8；选取计算努森数Kn为10-3～102，在这个范围内大气特性变化显著。其余输入条件由Nrlmsise-00 大气模型给出。

表4 模型验证的计算条件Tab.4 Calculation conditions for model verification

为了进行对比，本文对DSMC 结果及使用本文提出的logistic-log 函数与sine-squared 函数进行比较。图6 表示了在本节所研究两种不同模型中，本文所提出的logistic-log 桥函数与sine-squared 桥函数在2 种不同模型中不同条件下的其努森数（Kn）与阻力系数（CD）之间的关系如图6 所示。在图6（a）中，红色线代表的是Ma=16 结果，蓝色线代表的是Ma=4 结果。其中，粗虚线为logistic-log 桥函数的计算结果，而细虚线为sine-squared 桥函数的计算结果，三角符号代表着使用自由分子流算法（Free Molecular Fuction，FMF）的计算结果，其余图像中的散点代表的是DSMC 计算数据。由图6 可知，对于圆柱的计算工况，在Ma=16 的高速工况下，logisticlog 桥函数与sine-squared 桥函数与DSMC 结果吻合都较好。在Ma=4 的工况下，sine-squared 桥函数的计算结果在Kn数较高时与实验数据产生了较为明显的差异。在图6（b）中，红色线代表的是Ma=8结果，蓝色线代表的是Ma=4 结果。其中，粗虚线为logistic-log 桥函数的计算结果，而细虚线为sinesquared 桥函数的计算结果，三角符号代表着使用自由分子流算法（Free Molecular Fuction，FMF）的计算结果，其余图像中的散点代表的是DSMC 计算数据。同样可以发现，对于钝头双楔体的计算工况，在速度较高的情况下（Ma=8），2 种桥函数与实验结果基本吻合。当Ma=4 时，使用sine-squared 桥函数在Kn=12 时出现了最大值，与DSMC 数据相差较大，也违背了随着Kn数的增加，气动阻力系数也随之增加的基本规律。以上结果说明了在固定系数的情况下，sine-squared 桥函数无法满足在不同外形，多种工况下均保持良好的估算精度；而logistic-log 桥函数可以实现在不改变系数的情况下满足对多种工况的气动系数估算。

图6 不同桥函数对应的气动特性估算结果Fig.6 Results of the aerodynamic characteristics corresponding to different bridge functions

为验证本文提出的形状因子对其结果的影响，本 文对q=1+0.632 6cosδ+6.124 9与q=1 进行对比，计算结果如图7 所示。

在图7（a）中，红色线图代表的是Ma=16 结果，蓝色线代表的是Ma=4 结果。粗虚线为形状因子q=1+0.632 6cosδ+6.124 9 的计算结果，而细虚线为形状因子q=1 的计算结果，图像中的散点代表的是DSMC 计算数据。由图7 可知，对于圆柱的计算工况，在Ma=16 的情况下，2 种形状因子与实验结果的吻合情况较好。在Ma=4 的计算条件下，不使用q=1+0.632 6cosδ+6.124 9 的形状因子则会出现明显的误差。对于钝头双楔体的计算工况，发现形状因子在不同工况下的影响不大。因为在钝头双楔体中，表面外形的倾角相对固定，形状因子对其影响则不明显；而在圆柱工况中，圆柱各个位置的倾角均有变化，则形状因子在速度较低的情况下起到了明显的修正作用。在速度较低时，形状因子起到的修正效果明显，在高速情况下形状因子的修正作用不显著。为了进一步量化计算结果与DSMC 数据的吻合程度，本文引入了均方误差对计算结果进行了比较，其计算公式如下：

图7 形状因子对计算结果的影响Fig.7 Effects of the shape factor on the calculation results

式中：ε为误差值为DSMC 的采样数据为使用不同桥函数得到的数值结果数据；k为DSMC 的数据维度。

通过对均方误差进行比较，计算出不同方法的误差，见表5。可见本文提出的计算方式误差最小。

表5 不同桥函数计算结果均方误差（ε）Tab.5 Mean square errors of the calculation results of different bridge functions

综上所述，本文提出logistic-log 桥函数与当前常用的桥函数相比，具有适用范围广，估算精度高的特点。本文提出的形状因子q能够根据计算外形对计算结果进行修正。

3 典型超低轨飞行器气动力计算

欧洲空间局（European Space Agency，ESA）的地球重力与海洋环流探测卫星GOCE 为典型的超低轨飞行器［18］，有关学者也使用了DSMC 或TPMC方法对其气动特性进行了研究［19-20］，本文第2 章提出的方法对其在工作高度范围内的阻力系数进行测算。

在不同侧滑角下，不同攻角下随着工作高度的变化GOCE 的阻力系数与阻力大小的变化趋势如图8 所示。随着工作高度的增加，阻力系数的变化范围越小，GOCE 在俯仰与偏航时所受的气动力变化相对较低。但是随着运行高度降低，大气密度显著增加，阻力呈几何量级增大。在平飞状态下，阻力的最大值为42.526 N，最小值为0.005 41 N。因此，超低轨卫星在进行工作时，需要及时进行姿轨控制维持其高度。

图8 不同高度下GOCE 的阻力与阻力系数大小与其姿态角度关系Fig.8 Relationships between the drag or drag coefficient of GOCE and its attitude angle at different altitudes

续图8 不同高度下GOCE 的阻力与阻力系数大小与其姿态角度关系ContinuedFig.8 Relationships between the drag or drag coefficient of GOCE and its attitude angle at different altitudes

GOCE 形态细长，姿态变化导致的气动力变化明显，在运行高度在100 km 时，其阻力系数变化范围在3.2～9.8；在高度为260 km 时其阻力系数的变化范围在4.1～9.5。其太阳能帆板占主要浸湿面积，在偏航时，迎风面积变化明显，导致了其在侧滑时的阻力系数与阻力变化明显高于俯仰时的情况。

4 结束语

1）本文提出了一种超低轨飞行器的快速气动力计算工具方法。使用面元法进行求解，通过将飞行器进行便面单元离散为迎风单元与背风单元，分别求解不同种类单元的气动力系数，通过对离散单元的数值结果进行积分得到了飞行器的气动力特性。通过重新设计桥函数，降低计算方法的数据依赖程度，提高了应用范围。

2）使用2 种桥函数计算，并对比钝头双楔体与圆柱体的DSMC 结果，说明了算法的可靠性和准确性，解释了本文提出的在计算不同工况气动力时，形状因子对计算结果有修正作用。

3）分析阐述了典型超低轨飞行器GOCE 在100～260 km 工作范围内随着俯仰角与偏航角变化气动特性的变化趋势。发现在发生姿态变化时的阻力系数变化显著，其变化范围在3.2～9.8，随着高度的降低阻力系数变化的范围越大。因此，飞行器在低轨工作时，需要为飞行器设计舵面或陀螺仪对其姿态进行控制。在过渡流区内，随着高度增加，气动阻力减小，阻力最大值为42.526 N，最小值为0.005 41 N。由此可见，超低轨飞行器需要在宽域过渡流区内运行，需要额外的动力装置。