一种复杂外形低轨航天器气动特性的快速预测方法

摘要

一种复杂外形低轨航天器气动特性的快速准确预测方法，首先，构建一个足够大的计算域，航天器位于该计算域内部，并在其边界产生一个试验粒子；然后，跟踪和模拟该试验粒子之后的运动轨迹和碰撞过程，直到其飞出计算域或撞到飞行器表面，若试验粒子和表面碰撞，计算其与物体表面的动量、能量交换并继续跟踪该试验粒子；最后，重复上述过程直至试验粒子数足够大，以保证计算结果收敛，统计出飞行器飞行环境下的气动特性规律，进而进行飞行器设计。

说明书

技术领域

本发明涉及一种复杂外形低轨航天器气动特性的快速预测方法，属于飞行器气动特性设计技术领域。

背景技术

为了得到精细的全波段地球重力场和厘米级大地水准面，在轨航天器，比如欧洲航天局设计的GOCE(Gravity field and steady state Ocean Circulation Explorer，重力场和稳态海洋环流探测者)，必须尽量接近地球表面，以最大化对地球重力场的敏感度。然而，航天器在地球地轨道飞行时会受到非保守力，如稀薄大气阻力、太阳光压和地球反照压等的极大影响，其中，大气阻力的影响最为明显。为了获得精确的地球重力场数据，必须抵消这些非保守力，尤其是稀薄大气阻力的影响，这就需要精确计算大气阻力等非保守力，进而根据计算结果设计出一套无阻力高度控制系统的卫星点火设备进行非重力补偿。GOCE等低轨卫星的飞行高度在240～450 km之间，流动区域属于稀薄程度最高的自由分子流区。

对于气体极端稀薄的自由分子流问题，气体分子与物体表面的碰撞占主导地位，由物体反射的分子只有在飞离物体表面很远之后才与其他分子碰撞，可以采用理论分析得出其分析解。这种理论分析方法近似认为Kn趋于无穷大，完全忽略分子间碰撞引起的气体速度分布函数的变化(此时称为无碰撞流动)，其基本方程是无碰撞项的Boltzmann方程。对于无碰撞流动物体定常扰流的问题，来流的速度分布函数是平衡态的分布，即Maxwell分布，来流分子对物体表面的动量传递和能量传递就可以通过对Maxwell分布求矩的方法计算出来，如果分子与表面的相互作用是已知的，那么反射分子从表面带走的动量和能量也可以计算出来，继而求出气动力和热。现在工程实际中广泛使用的面元积分方法就是把复杂外形分成很多个面元，利用上述理论分析给出的平板绕流结果对所有面元进行积分得到相应的气动力和热。然而，这种方法只适用于简单外形的流动问题，对于复杂外形由于忽略了各个实体部分之间的相互遮挡和多次反射效应，工程应用中存在一定误差。

近年来，直接模拟Monte Carlo(direct simulation Monte Carlo，DSMC)方法在稀薄气体动力学研究领域获得了广泛的应用，其算法也趋于成熟。DSMC方法能够模拟三维稀薄气体流动的复杂流场，而且能真实地模拟包含热化学非平衡反应等复杂的稀薄气体流动问题。然而，DSMC方法虽然计算精度高，却存在极端消耗计算机资源的问题，对计算机的运行速度和存储量要求都很高。另一方面，低轨航天器的无阻力高度控制系统需要长期的、精确的和实时的气动阻力输入，在当前的计算机硬件条件下，大型并行机上的DSMC程序也无法满足低轨航天器气动阻力计算的实时性要求。

综上分析可以发现，低轨航天器的正常工作要求无阻力高度控制系统长期对稀薄大气气动阻力的非重力补偿，而无阻力高度控制系统则需要精确实时的气动阻力输入。现有的低轨航天器气动特性计算方法，包括面元积分法和DSMC技术，前者因为忽略了复杂外形实体之间多次反射而无法满足精确性要求，后者则由于计算量太大而达不到实时性要求。

发明内容

本发明的技术解决问题是：克服现有技术的不足，根据现有复杂外形低轨航天器气动特性计算的需要，基于稀薄气体动力学中的自由分子流理论，结合真空技术中计算通过管道流率的随机模拟技术，提供了一种复杂外形低轨航天器气动特性的快速准确预测方法。

本发明的技术解决方案是：等最后我补充。

本发明与现有技术相比的优点：

当前低轨航天器气动特性的计算方法，要么忽略了复杂实体之间的多次反射导致准确度降低，要么计算量太大无法满足实时性要求。基于此，本发明方法基于自由分子流理论，结合真空技术中计算通过管道流率的随机模拟技术，提供了一种复杂外形低轨航天器气动特性的快速准确预测方法。相对于现有的低轨航天器气动特性计算方法，本文的TPMC方法由于考虑了占主导地位的分子和表面之间的碰撞，不考虑次要的分子之间的碰撞，一次只产生和跟踪一个试验粒子，对计算存储要求低，计算速度快，且能够模拟复杂凹形几何体引起的多重反射效应和流动遮挡效应，同时满足准确性和实时性要求，是低轨航天器气动特性计算的理想方法。

附图说明：

图1：本发明方法流程图；

图2：本发明方法低轨卫星(GOCE)计算模型示意图；

图3：GOCE卫星阻力系数计算结果与文献结果对比图；

图4：本发明方法带太阳翼卫星简化模型示意图；

图5：带太阳翼简化卫星气动力计算结果((a)为法向力F_Z，b(为)俯仰力矩M_Y)与DSMC结果对比图；

图6：带太阳翼简化卫星流动遮挡和多次反射效应导致的附加气动力；其中(a)为卫星本体法向力F_ZSa，(b)为太阳翼侧向力F_YAr1，F_YAr2。具体实施方式

如图1所示，本发明提供了一种复杂外形低轨航天器气动特性的快速准确预测方法。首先，构建一个足够大的计算域，航天器位于该计算域内部，并在其边界产生一个试验粒子；然后，跟踪和模拟该试验粒子之后的运动轨迹和碰撞过程，直到其飞出计算域或撞到飞行器表面，若试验粒子和表面碰撞，计算其与物体表面的动量、能量交换并继续跟踪该试验粒子；最后，重复上述过程直至试验粒子数足够大，以保证计算结果收敛，统计出飞行器飞行环境下的气动特性规律，进而进行飞行器设计。

具体步骤为：

(1)构建一个足够大的圆柱体计算域，低轨航天器位于该计算域内部。设R_ref，L_ref分别为航天器的特征半径和长度，则计算域的半径R_c和长度L_c分别为

$>\left(\begin{array}{c}{R}_c=K·R_{ref}\\ L_c=K·L_{ref}\end{array}\right)$>

其中，K为足够大的正参数(例如，K＝4)以保证计算结果的稳定。

(2)根据自由来流条件在计算域边界产生试验粒子，对其初始位置r和速度c进行随机取样，具体实现过程如下：

引入来流速度比，其定义为来流气体宏观速度与最可几分子热运动速度的比值，即s＝c₀/c_m，其中，c0为来流气体宏观速度，为最可几分子热运动速度，k＝1.38×10^-23J·K^-1为Boltzmann常数，T为来流温度，m为来流分子质量，确定在计算域前面f、侧面s和后面b产生试验粒子的概率为

$>\left(\begin{array}{c}{P}_f=\frac{R_c·\chi \left(s\right)}{R_c·\chi \left(s\right)+2L_c+R_c·\chi (-s)}\\ P_s=\frac{2L_c}{R_c·\chi \left(s\right)+2L_c+R_c·\chi (-s)}\\ P_b=\frac{R_c·\chi (-s)}{R_c·\chi \left(s\right)+2L_c+R_c·\chi (-s)}\end{array}\right)$>

式中，函数χ(t)定义为

$>\chi \left(t\right)=\exp (-t^2)+\sqrt{\pi }t[1+erf\left(t\right)]$>

在(0，1)之间产生一个随机数R，根据随机数的大小结合上述概率确定产生试验粒子的计算域边界：若R≤P_f，则在计算域前面产生试验粒子；若P_f＜R＜P_f+P_s，则在侧面产生试验粒子，若R≥P_f+P_s，则在后面产生试验粒子。

试验粒子初始位置取样的基本思想是试验粒子在计算域边界均匀随机分布。如果试验粒子在计算域前面产生，则其初始位置为

$>\left(\begin{array}{c}x=-L_c/2\\ y=R_c\sqrt{R_1}cos\theta ,\theta =2\pi ·R_2\\ z=R_c\sqrt{R_1}sin\theta \end{array}\right)$>

同理，计算域侧面和后面产生试验粒子的初始位置为

θ＝2π·R₂和

试验粒子的初始速度c为来流宏观速度和分子热运动速度的矢量和，即

c＝c₀+c_m·ξ

式中，无量纲热运动ξ在柱坐标中前两个分量的随机取样为

θ_r＝2πR₄

第三个分量ξ_z可通过多项式拟合给出，即

ξ_z＝a₀+a₁η+a₂η²+a₃η³+a₄η⁴+a₅η⁵+a₆η⁶+a₇η⁷

其中多项式系数a₀～a₇可由相关文献给出。

(3)根据步骤(2)给出的试验粒子初始位置r和速度c，计算该试验粒子的运动轨迹；

(4)判断运动轨迹是否与航天器表面发生碰撞：如果没有发生碰撞，则认为该试验粒子直接飞出计算域，转到步骤(2)；如果发生碰撞，计算碰撞发生的位置并继续执行步骤(5)；

(5)按照指定的气体与表面相互作用模型计算试验粒子与航天器表面碰撞后的速度，计算单次碰撞的动量和能量交换，之后按照新的位置和速度继续跟踪该试验粒子(即重新计算其运动轨迹后执行步骤4)；

当气体与航天器表面相互作用为漫反射时，反射后的分子的速度分布是以物面温度T_w为参考值的Maxwell平衡分布；气体与表面相互作用为镜面反射时，反射后的分子在法向的速度分量反向，在切向的速度分量保持不变。

(4.1)计算试验粒子碰撞前的动量P_i和总能E_i，具体公式如下

$>P_i={\mathrm{mc}}_i,E_i=\frac{1}{2}m{\left|c_i\right|}^2+\frac{5-3\gamma }{\gamma -1}\frac{{\mathrm{kT}}_{\infty }}{2}$>

式中，c_i为碰撞前试验粒子的速度矢量，γ为来流气体比热比，T_∞为来流温度。

(4.2)计算试验粒子碰撞后的动量P_r和总能E_r，具体公式如下

$>P_r={\mathrm{mc}}_r,E_r=\frac{1}{2}m{\left|c_r\right|}^2+\frac{5-3\gamma }{\gamma -1}\frac{{\mathrm{kT}}_w}{2}$>

式中，c_r为碰撞前试验粒子的速度矢量，T_w为航天器表面温度。

(4.3)计算单个试验粒子碰撞过程的动量交换ΔP和能量交换ΔE

$>\left(\begin{array}{c}\Delta P={\mathrm{mc}}_i-{\mathrm{mc}}_r\\ \Delta E=E_i-E_r\end{array}\right)$>

(6)跟踪足够多的试验粒子后，统计空气动力学力学宏观量，从而得到低轨航天器飞行环境下的气动特性规律。

一般跟踪的试验粒子的数量大于等于10⁷时，开始统计空气动力学宏观量。根据统计误差和计算时间情况考虑，可以在10⁷～10⁸内取值。

设N为所有试验粒子与该表面碰撞的总次数，则作用于航天器表面的气动力F、力矩M和热交换量Q为

$>F=A\underset{j=1}{\overset{N}{\Sigma }}\Delta P_j,\mathrm{MA}\underset{j=1}{\overset{N}{\Sigma }}(r_j\times {\mathrm{\Delta P}}_j),Q=A\underset{j=1}{\overset{N}{\Sigma }}{\mathrm{\Delta E}}_j$>

式中，r_j，ΔP_j，ΔE_j分别为第j次碰撞时力矩参考点到碰撞位置的矢量、动量交换量和能量交换量，气动标准化参数A的表达式为

$>A=\frac{n_{\infty }(Q_f+Q_s+Q_b)}{N_{tp}}$>

式中，n_∞为来流气体分子数密度，N_tp为试验粒子总数，Q_f、Q_s和Q_b分别为单位时间通过计算域前面、侧面和后面单位面积的气体分子数，其表达式为

$>\left(\begin{array}{c}{Q}_f=\frac{n_{\infty }}{2\sqrt{\pi }}{c}_m·\chi \left(s\right)\\ Q_s=\frac{n_{\infty }}{2\sqrt{\pi }}{c}_m·\chi (-s)\\ Q_b=\frac{n_{\infty }}{2\sqrt{\pi }}{c}_m\end{array}\right).$>

实施例

低轨航天器气动特性预测的具体求解实例如下：

首先考虑欧洲空间局(European Space Agency，ESA)的重力场和稳态海洋环流探测卫星(Gravity field and steady state Ocean Circulation Explorer，GOCE)，图2为其模型示意图。为了最大化对地球重力场的敏感度，该卫星的飞行高度低至250-300 km，表1为计算条件，气体分子反射模型为漫反射。

表1 GOCE卫星验证算例计算条件

图3给出了GOCE卫星阻力系数随速度比(卫星飞行速度大小与来流最可几分子热运动速度的比值)，并与德国HTG(哥廷根高超声速研究所)发表的结果进行对比，实线代表程序计算结果，散点代表文献结果。速度比的变化范围为[2，20]，覆盖了GOCE卫星的飞行速度范围，同时给出完全漫反射(σ＝1)和部分漫反射(σ＝0.6)这2种情况的结果。根据图3，随着速度比的增加(可以理解为飞行速度大小的增加)，阻力系数呈指数型减小；速度比较小时，完全漫反射对应的阻力系数大于部分漫反射的，速度比较大时则刚好相反。需要注意的是，完全漫反射和部分漫反射情况下，计算结果和文献结果都符合得很好，验证了本发明方法的正确性。

完成方法验证之后，考虑运行在地球低轨道的带太阳翼简化卫星(图4)的气动特性，以研究本发明方法的计算效率和优势，表2为计算条件。

表2 带太阳翼简化卫星计算条件

图5给出了法向力、俯仰力矩随来流攻角(α∈[-10°，10°])的变化曲线，并和不考虑分子之间碰撞DSMC的计算结果进行对比。DSMC的计算网格数为4.5万，仿真粒子数约1300万，50个CPU并行计算，时间步长10^-6s，16万步开始采样，共采样60万步，共耗时约40小时；而本发明方法的试验粒子数1亿，1个CPU计算，耗时20s。显然，本发明方法的结果和DSMC结果符合较好，且前者的统计误差明显小于后者，通过两者计算消耗时间的对比，容易看出本发明方法的计算速度比DSMC快几个量级。因此，在低轨卫星气动计算领域，本发明方法具有计算效率高，存储要求低等DSMC方法不具有的优势。

根据几何关系，卫星本体的某些表面与2个太阳翼的某些表面之间存在遮挡和多次反射效应(下文统称为干扰效应)。由于干扰效应引起的气动力变化比较微弱，容易被统计误差掩盖，图6给出了相同条件下20次计算卫星本体法向力、太阳翼侧向力的结果，并进行平均。显然，如果不考虑卫星本体与太阳翼之间的干扰效应，由于几何对称性，卫星本体的法向力F_ZSa≡0，2个太阳翼的侧向力F_YAr1≡0，F_YAr2≡0。但是，卫星本体展向(Y轴方向)2个侧面与太阳翼迎风面之间的干扰效应(主要是多次反射效应)导致卫星卫星本体受到一个负向的法向力。同样，卫星本体2个侧面和太阳翼迎风面、内侧面的干扰效应在2个太阳翼上产生了额外的侧向力，而且，经过平均之后，2个太阳翼受到的侧向力几乎大小相等，方向相反，故整个卫星受到的侧向力为0。

综合以上GOCE卫星和带太阳翼简化卫星气动特性的分析可以得到以下结论：本文的试验粒子Monte Carlo方法对计算存储要求低，计算速度快，且能够模拟复杂几何外形引起的流动遮挡效应和多次反射效应，同时满足低轨卫星气动计算准确性和实时性要求，是复杂外形低轨航天器气动特性计算的理想方法。

本发明未公开技术属本领域技术人员公知常识。