天宫飞行器低轨控空气动力特性一体化建模与计算研究

李志辉，吴俊林，彭傲平，唐歌实（.中国空气动力研究与发展中心超高速所，绵阳62000；2.国家计算流体力学实验室，北京009；.北京航天飞行控制中心航天飞行动力学技术重点实验室，北京00094）

·工程技术·

天宫飞行器低轨控空气动力特性一体化建模与计算研究

李志辉1，2，吴俊林1，彭傲平1，2，唐歌实3
（1.中国空气动力研究与发展中心超高速所，绵阳621000；2.国家计算流体力学实验室，北京100191；3.北京航天飞行控制中心航天飞行动力学技术重点实验室，北京100094）

对非规则板舱组合体天宫飞行器300～200 km低轨道飞行过程空气动力特性一体化计算建模，提出考虑复杂构型物面遮盖效应面元解析法与经修正的Boettcher/Legge非对称桥函数，发展基于三角形面元逼近复杂外形通用处理方法，建立适于天宫飞行器复杂物形处理与面元气动力系数计算规则；将DSMC方法与求解Boltzmann模型方程气体运动论统一算法应用于天宫飞行器简化外形，进行气动力当地化关联参数计算修正，建立针对大型复杂结构天宫飞行器低轨道飞行控制过程空气动力特性一体化快速算法与程序软件。对大尺度圆柱体外形与天宫飞行器300～200 km不同高度变轨飞行过程不同迎角与侧滑角及帆板平面与本体主轴不同夹角复杂构型气动力特性计算分析验证，表明天宫飞行器在200 km以上低轨道飞行控制过程中所受空气动力系数随飞行高度发生显著变化（8%～50%），证实长期在轨运行的大型航天器若采用统一固定的气动力系数，误差累积巨大，需要采取防护措施，低轨道飞控大气阻力仍是制约航天器定轨预报精度最关键因素。

天宫飞行器；低轨控空气动力特性；当地化桥函数；关联参数确定；工程计算；DSMC方法；统一算法

1　引言

航天器从离地面数百公里外层空间往地球低轨道再入飞行过程，是一个自数百到数十量级Knudsen（克努森）数高稀薄自由分子流动状态多物理场复杂构型极高超声速流动问题［1］。低地球轨道环境，是对地观测卫星、气象卫星、通讯遥感卫星、载人飞船、空间站等大型复杂结构航天器运行区域。距离地面200～1000 km的低地球轨道环境大气压力很低、气体非常稀薄，分子间断粒子效应与稀疏气体分子分布变化通过对航天器产生力/热冲击作用，直接影响航天器运行轨道、姿态和寿命。为了研究和控制航天器低地球轨道环境飞行过程，在过去计算机条件不够先进情况下，传统做法是近似处理，不同流动区域发展相应的计算方法。对于Knudsen数很大的稀薄气体自由分子流区，气体流动由无碰撞或近于无碰撞的自由分子流理论控制。为了进行航天器从340～200 km低轨道飞行试验，通常认为飞行器在200 km以上，气体分子平均自由程达两百多米以上，其Knudsen数数十以上，飞行绕流已属于完全自由分子流，阻力、升力等空气动力不再随飞行高度发生变化，而采用固定的气动力系数参与飞行力学轨道计算［2］。于是大型航天器开展低轨道飞行试验采用固定阻力系数2.2，发现飞行偏差有时越来越大，以致使用其上RCS姿控系统强制控回标称轨道，造成燃料消耗，甚至有时难以准确有效控制［3，4］。为此，提出这样的问题：我国天宫飞行器低轨道飞行所受大气阻尼系数是否会变化？在轨天宫飞行器本体直径数米、长度达十余米、其携带的太阳电池翼帆板展开宽度数十米，对如此复杂巨大不规则板舱组合体飞行器实施低轨道飞行控制、轨道参数计算时，能否采用统一、固定的空气动力系数？如何开展天宫飞行器低轨道飞行过程空气动力特性一体化建模与计算研究［4］，剖析天宫飞行器气动特性随飞行高度变化规律？为此，本文在研究建立航天空气动力学统一计算理论与各流动区域多物理场复杂流动机理模拟方法研究基础上，尝试开展上述问题计算研究。

连接微观分子动力学与宏观流体力学的介观Boltzmann（玻尔兹曼）速度分布函数方程［5］本身可描述各个流域气体分子输运现象，该方程自1872年提出以来，一直是国际学术界追踪研究的问题，它作为一个高度复杂高维积分、微分、多相空间多尺度非线性刚性问题［5，6］，精确求解描述各流域气体流动特征的Boltzmann方程至今未成现实。为此，近五十年来世界航天发展了基于微观分子动力学将分子运动与碰撞解耦进行随机统计模拟的DSMC（直接模拟Monte Carlo）方法，该方法自1963年Bird［7］将其发展用于稀薄气体流动模拟至今，已在稀薄气体动力学领域获得了广泛的应用与检验，国内外学者将DSMC方法用于模拟包括内能、化学反应、电离、辐射等多种作用物理模型的稀薄流区气体流动问题研究，已取得较好的应用发展［8-17］，并利用巨型计算机研究并行算法，引入中枢网格系统开展并行计算［18，19］。在DSMC方法对较高Knudsen数稀薄气体流动仿真取得巨大成功同时，该方法因受网格划分、时间步长和模拟分子数等方法本身模拟准则所限，难以对低Knudsen数复杂大尺度飞行器气动问题数值仿真［14，20，21］。

为了探索跨流域气体流动问题一体化模拟方法，通过跟踪国际上关于Boltzmann方程研究现状与发展趋势［22-27］，结合从事DSMC方法与计算流体力学有限差分法研究基础，过去十余年本文第一作者等从研究Boltzmann方程碰撞松弛间隔理论出发，确立描述各流域微观分子输运现象统一的Boltzmann模型方程，提出并应用离散速度坐标法对分子速度分布函数数值离散，发展可用于速度空间宏观流动取矩的离散速度数值积分法。将计算流体力学有限差分方法推广拓展到基于时间、位置空间与速度空间的Boltzmann模型方程数值求解，先后建立起模拟各流域一维、二维、三维绕流问题气体运动论统一算法（GKUA）理论与系列计算技术［28-33］，并开展了飞船返回舱、近空间飞行器、箔条云等跨流域气动问题研究［33-37］，解决了系列工程研制需求问题。由于上述求解Boltzmann模型方程统一算法将气体分子速度分布函数基于位置空间与速度空间的信息均保存起来数值求解，需要在位置空间与速度空间组成的六维多相空间进行计算，由此占用大量计算机内存资源。如果进行三维高超声速气体流动问题数值计算，如天宫飞行器跨流域绕流问题研究，必须开展大规模并行计算。

另一方面，航天器轨道计算可靠性直接依赖于轨道动力学方程空气动力项求解的准确性，需要与飞行力学数值推进方法耦合，计算确定航天器任意时刻空间位置与速度，需要发展跨流域空气动力特性快速计算方法。为了能在有限时间，建立天宫飞行器低轨道飞行过程空气动力一体化计算平台，拟以稀薄气体高超声速绕流当地化桥函数理论为基础，开展天宫飞行器低轨道飞行控制过程300～200 km空气动力特性当地化桥函数工程计算建模，研究使用DSMC方法、求解Boltzmann模型方程的气体运动论统一算法对天宫飞行器简化外形典型绕流状态计算，修正适于大尺度复杂结构航天器当地化关联参数，建立适于天宫飞行器低地球轨道环境不同高度、马赫数与不同迎角、侧滑角空气动力特性快速计算分析应用研究框架。

2　适于天宫飞行器低轨控不同高度空气动力特性当地化快速算法

稀薄气体流动领域气动力系数计算所用的当地化工程计算方法，是基于半经验桥函数理论的当地桥公式法［38，39］。自由分子流和牛顿连续流两个极限流态都可以用当地气动力的解析来表述，当地表面面元上的气动力系数只依赖于来流和当地性质，如当地迎角、表面作用等。当地连续流和自由分子流系数之间的间隙可以用半经验的桥公式光滑搭接起来。这是一种加权函数，依赖于独立的关联或相似参数（如克努森数），通常由实验结果与数值计算分析确定。对连续流区域的压力系数计算，采用修正的牛顿无粘流理论以及背风真空效应法进行估算；对自由分子流区，采用基于不同模型材料修正的Nocilla壁面反射模型进行压力系数和摩擦力系数计算；对过渡流区域，基于修正Boettcher/Legge非对称桥函数理论［40］，使用可分段描述的非对称压力与摩擦力系数关联桥函数。通常针对轴对称钝体飞行器设计，能预测估算钝柱形高超声速飞行器的稀薄气动力特性，但本文研究大型复杂结构航天器如天宫飞行器，其空气动力特性计算不能沿用常规旋成体当地化计算方法［41］。可研究使用DSMC方法、求解Boltzmann模型方程从稀薄流到连续流统一算法［28-37］，基于准确计算进行当地化桥关联系数的选择调试与检验修正，确定各关联参数的最佳组合［35，42］，发展针对纵横数十米大尺度复杂结构航天器，跨越高稀薄自由分子流、过渡流到大气层内低高度各流区、高低不同马赫数、不同迎角与侧滑角气动力系数快速计算方法。

2.1高超声速空气动力特性当地化计算方法

当地化桥公式法首先考虑的是在两个边界流域连续流和自由分子流区域，气动力特性可用两个不同的数学函数称为连续流函数FCont和自由分子流函数FFM来描述，而在过渡流区域需要给出一个新的函数Fb，亦即桥函数［4，35，38，39］，使其在连续流或自由分子流分别趋近FCont和FFM。对给定入射角θ，当地面元上经归一化压力和摩擦力系数表达为式（1）：

式中，Fb，p和Fb，τ分别是压力和摩擦力桥函数，依赖于独立参数Kn、TW/T0和θ。

对连续流区压力系数计算，可采用诸如修正牛顿无粘流理论及背风真空效应法计算；对自由分子流区，可使用基于不同模型材料修正的Nocilla壁面反射模型进行压力与摩擦力系数计算；对过渡流区，研究修正的Boettcher/Legge非对称桥函数理论［39］，发展可分段描述的非对称压力与摩擦力系数关联桥函数，其中，非对称压力桥函数［38，39，42］可表示为式（2）：

并且Knm，τ、ΔKnτ1、ΔKnτ2为可调参数［40］。

当地化方法关于飞行器物形处理［38，43，44］有两种方法，一种是用解析式分别描述飞行器的各部分，这种方法适用于外形简单且易于解析表达的轴对称旋成体飞行器［38，41，44］，但大多数飞行器外形复杂且很不规则，难于进行解析表示。于是当地化计算方法所依赖的物形处理常常使用另一种较为通用而近似的面元处理方法［4，39，43］，该方法将飞行器物面划分为若干个面元，利用四边形或三角形来逼近复杂外形，计算出每个面元上中心点坐标、中心点处的法向及面元面积。

飞行器外形曲面分成若干块小的曲面，对于每一个小曲面，选用一个小的平面四边形来代替。这样，计算小曲面上的气动力就转换成计算面元上的气动力，将这些面元上的气动力加起来就可以得到整个飞行器的气动力。当面元分得很细小时，由此引起的误差就更小［35，44］。

2.2天宫飞行器计算模型构建及表面非结构网格物形表征技术

由于天宫飞行器与其左右宽度达数十米的太阳电池翼帆板及连接装置构成一个相当复杂而极不规则巨大板舱组合体。为了解决天宫飞行器及附件复杂物形表征的困难，研究引入分区网格与非结构网格生成技术，以此获得能准确表征复杂物面边界的非结构网格，以便进行复杂物形表征处理。模型表面采用非结构贴体网格，通过网格处理得到计算需要的表面网格，具体过程见文献［4］，其中基于三角非结构网格表征天宫飞行器物形，图1绘出天宫飞行器计算模型表面三角非结构网格布局。

图1　天宫飞行器计算模型表面非结构网格布局Fig.1 Com putational m odel and sur face grid layout of the TG spacecraft

2.3天宫飞行器跨流域气动特性快速算法

根据稀薄过渡流流动特征，当地化快速算法将稀薄过渡流区压力桥函数（2）式与摩擦力桥函数（3）式用非对称误差函数形式关联，两式中Knm，p、ΔKnp1、ΔKnp2及Knm，τ、ΔKnτ1、ΔKnτ2应根据风洞试验数据或数值计算结果进行调试确定，由于天宫飞行器的大尺度复杂结构与在轨飞行高稀薄自由分子流区绕流特点，风洞实验条件并不现实，可依靠高性能计算机开展DSMC方法［4，13，17］或求解Boltzmann模型方程统一算法［30-34］典型飞行状态数值计算，调试选择与计算修正上述六个关联参数，使当地化工程计算结果最大程度地满足实际计算精度。由于本文旨在有限时间搭建天宫飞行器低轨道飞行过程300～200 km空气动力特性初步计算平台，作为算法验证，桥函数中当地化可调关联参数的确定采用250 km高度与天宫飞行器具有同样尺寸的圆柱体绕流数值计算结果进行修正确定。首先利用数值计算得到该圆柱体绕流典型状态所受到的气动阻力系数为参考对象；改变上述六个关联参数，使用当地化快速算法进行调试计算修正；最终获得一套针对天宫飞行器计算适用的六个关联参数的最佳组合，使得经关联参数修正后的气动特性当地化计算方法能满足本文关于天宫飞行器在不同高度、不同马赫数、不同迎角与侧滑角的气动力系数计算要求。图2给出了250 km高度0°迎角上述大尺度圆柱体绕流状态压力等值线云图分布，可看出圆柱体前端表面压力值远大于后表面，气流完全附着物面强扰动，流场呈现出高稀薄自由分子流动特征。

图2　 H=250 km，α=0°大尺度圆柱体绕流压力等值线分布Fig.2 Pressure distribution contours around largescale cylinder body w ith H=250 km，α=0°

数值计算得到0°迎角时圆柱体模型阻力系数为3.169025，以此为基础得到当地化快速算法中桥函数的可调关联参数如式（5）～（6）：

以上述六个当地化关联参数为基础，发展可用于确定大尺度飞行器在不同高度、马赫数、迎角与侧滑角的气动力系数快速算法。由此计算得到上述圆柱体绕流阻力系数3.169013，与数值模拟结果相当一致，彼此偏差仅0.000379%。表1列出了该状态数值计算与当地化快速算法分别得到4°迎角上述圆柱体模型气动力系数比较情况，可看出二者吻合很好，偏差范围0.063%～2.85%，验证了本文提出发展低轨道飞行过程大尺度飞行器300 km到200 km空气动力特性快速计算平台可行性。

表1　当地化快速算法计算H=250 km，α=4°圆柱体绕流气动力系数验证比较Table 1 Validation and comparison of aerodynam ic coefficients around the large-scale cylinder w ith H=250 km，α=4°

3　天宫飞行器简化外形200 km至300 km空气动力特性计算分析

3.1算法验证

为了进一步验证本文发展的当地化快速算法模型正确性以及200～300 km高度飞行的气动力系数变化规律，分别采用数值模拟方法和当地化快速算法计算了200 km和300 km上述大尺度圆柱体0°迎角气动力系数，表2给出了两类方法计算0°迎角200 km和300 km阻力系数Cd比较情况。可看出，两种计算方法得到的气动力系数完全吻合，最大偏差不超过0.48%，证实本文发展的当地化快速算法正确可靠；在300 km飞行高度阻力系数比200 km阻力系数均高出至少7.43%，说明近地轨道飞行过程200～300 km不同高度飞行器阻力系数是随飞行高度增加而增大，这反映出传统工程处理手段认为200～300 km气体分子平均自由程λ∞达235.2～2595.37m，远大于飞行器特征尺寸，来流气体处于自由分子流态，飞行器绕流使用固定不变阻力系数统一处理方式对类天宫飞行器两舱结构大尺度圆柱体绕流会产生至少7.43%偏差。

3.2天宫飞行器200～300 km气动力特性计算分析

为使本文当地化快速算法计算300～200 km期间天宫飞行器气动特性结果能尽可能真实可靠，本文对任一电池翼帆板平面与天宫飞行器本体主轴夹角0°、22.5°、45°、67.5°、90°确定的计算模型，进行300～200 km不同高度300 km、280 km、250 km、220 km、200 km，不同迎角α= -13.5°，0°，13.5°，22°，38°与侧滑角β=5.5°，2.6°，0°，-2.6°，-4.1°的天宫飞行器气动特性计算，提供了大量可供轨道控制弹道飞行力学所用空气动力系数计算数据。下面仅列两组典型状态。

表2　 200 km和300 km飞行高度0°迎角圆柱体绕流阻力系数验证比较Table 2 Com parison of drag coefficients around the large-scale cylinder w ith H=200、300 km，α=0°

3.2.1帆板平面与天宫飞行器主轴夹角为0°气动特性计算

拟定天宫飞行器250 km、高度13.5°迎角、0°侧滑角外形，根据数值计算调试确定的当地化桥函数关联参数为Knm，p=0.9、ΔKnp1=3.5、ΔKnp2=3.5，Knm，τ=0.85、ΔKnτ1=3.5、ΔKnτ2=3.0。图3绘出天宫飞行器简化外形展向水平面与纵向竖直平面内压力等值线分布，可看出绕流流场附着物面流动，仅在物面附近出现强扰动，逐渐扩展到远前方，表现出高稀薄自由分子流动特征，等值线分布范围相当广，达天宫飞行器本体尺寸2～3倍。

图3　天宫飞行器简化外形帆板与主轴夹角0°高度250 km迎角13.5°、侧滑角0°轴对称面压力等值线分布Fig.3 Pressure distribution contours around the sim p lified TG-spacecraft w ithθ=0°，H=250 km，α= 13.5°，β=0°

表3列出了帆板平面与本体主轴夹角θ=0°时，考虑帆板与天宫飞行器两舱本体三脚架连接（见图1）天宫一号目标飞行器计算模型，采用本文算法计算得到α=-13.5°、β=5.2°条件下，不同飞行高度气动力系数，包括轴向力Ca、法向力Cy、侧向力Cz、阻力Cd、升力Cl、侧力Czl，其中“最大差异”指300 km的气动力系数与200 km相应气动力系数值之间的相对偏差。可看出，轴向力系数与法向力系数最大差异在5%～10%，而侧向力系数最大差异在7%～13%。但是，由轴向力、法向力和侧向力系数经坐标变换得到的风轴系下阻力、升力和侧力系数的最大差异分别达到了10%、30%和80%。由此可见，在200～300 km的大型复杂结构航天器飞行过程，若采用固定统一的气动力系数进行飞行控制会产生较大误差。

表3　天宫飞行器帆板与本体0°夹角、-13.5°迎角、5.2°侧滑角200～300 km各高度气动力系数Table 3 Aerodynam ic coefficients around the TG-spacecraft from 200～300 km withθ=0°，α=-13.5°，β=5.2°

3.2.2帆板平面与天宫飞行器主轴夹角为45°气动特性计算

对θ=45°，天宫飞行器与帆板形成计算模型见图4（a），根据数值计算确定Knm，p=0.9、ΔKnp1=3.2、ΔKnp2=3.5，Knm，τ=0.65，ΔKnτ1=3.1，ΔKnτ2=3.0。图4（b）绘出天宫飞行器简化外形绕流展向水平面与纵向竖直平面内压力等值线分布云图，同样看出绕流流场附着物面流动，仅在物面附近出现强扰动，并逐渐扩展到远前方，表现出高稀薄自由分子流动特征。

表4介绍了本文方法计算天宫飞行器θ= 45°、α=0°、β=0°条件下，200～300 km间隔20 km各飞行高度气动力系数。可看出，300 km与200 km的气动力系数差别很明显，随飞行迎角状态不同，这种差异也不同。对于某些飞行状态，出现侧力系数正负符号变化，法向力系数“最大差异”达25%以上。

图4　天宫飞行器简化外形帆板与主轴夹角θ= 45°，H=250 km，α=0°，β=0°模型与对称面压力分布Fig.4 Pressu re distribution contours aroundthe simp lified TG-spacecraft withθ= 45°，H=250 km，α=0°，β=0°

表4　θ=45°，α=0°，β=0°飞行过程200～300 km各高度气动力系数Table 4 Aerodynam ic coefficients around the TG-spacecraft from 200～300 km w ithθ=45°，α=0°，β =0°

4　结论

针对在轨服役天宫飞行器板舱组合体大型复杂构型特点，研究修正的Boettcher/Legge非对称桥函数，发展了可靠模拟不同高度、不同迎角与侧滑角及不同复杂外形高超声速空气动力特性快速计算方法；利用三角形逼近复杂物形，研制适于天宫飞行器复杂物面表征的非结构网格生成技术；研究通用面元处理法，发展适于天宫飞行器复杂物形处理与面元气动力系数计算规则；建立了适于天宫飞行器近地轨控200～300 km空气动力特性一体化快速算法应用研究平台，揭示大尺度天宫飞行器空气动力变化规律，并得到如下结论：

1）天宫飞行器在300 km到200 km变轨飞行过程中，气动力系数随飞行高度变化幅度较大：体轴系中轴向力系数与风轴系中阻力系数变化幅度5%～9%；法向力、侧向力等其它气动力系数变化幅度更大，有的甚至超过50%，某些状态（如帆板平面与本体主轴夹角67.5°，以0°迎角飞行）法向力系数“最大差异”达48%，而一些飞行状态升力和侧力系数出现了正负符号与方向差别。对长期在轨飞行大型航天器，这种误差累积后果难以控制，须采取防护措施；

2）对天宫飞行器这类由太阳电池翼与主体构成大型复杂板舱结构航天器低轨道飞行300～200 km气动力系数随飞行高度、姿态与太阳电池翼绕连接轴转动外形变化而改变，变化幅度与飞行状态相关；不同迎角、侧滑角气动力系数存在较大差异；同时，太阳能电池翼围绕主体连接轴旋转对天宫整器气动力特性也产生较大影响；

3）天宫飞行器近地轨道飞行过程，气动力系数不仅发生变化，且变化幅度与方向不同。如轴向力和阻力系数随飞行高度300～200 km变化的“最大差异”在5%～10%，该类气动力系数本身基数大，随高度变化量绝对值较大（一般在0.18左右）；与之相应法向力系数“最大差异”有时达48%，但其随高度变化量绝对值为0.06；而某些状态，如太阳电池翼平面与本体主轴夹角45°，以0°迎角、2.6°侧滑角所受侧力系数200 km为-0.0001，而300 km变为0.0032，出现正负侧力方向的改变，它们之间“最大差异”就会出现异常严重，而该侧力系数随高度变化量绝对值仅0.0033。这种由太阳电池翼旋转所致气动力方向上的变化会直接危及航天器精细轨控。

（致谢：本文写作过程中和李革非研究员、李勰工程师及课题组蒋新宇助理工程师等进行了交流讨论，部分计算在总参三部北方计算中心完成，特此感谢！）

［1］PilinskiM D，Argrow BM，Palo SE.Drag coefficientsofsatellites with concave geometries：comparingmodels and observations［J］.Journal of Spacecraft and Rockets，2011，48（2）：312-325.

［2］Moe K，Bowman B R.The effects of surface composition and treatment on drag coefficients of spherical satellites［C］// AAS/AIAA Astrodynamics Specialist Conference，Americal Astronautical Soc.，Paper 05-258，Springfield，VA，2005.

［3］Moe K，Moe M M，Wallace SD.Improved satellite drag coefficient calculations from orbitalmeasurements of energy accommodation［J］.Journal of spacecraft and rockets，1998，35（3）：266-272.

［4］李志辉，吴俊林，彭傲平.天宫飞行器低轨道飞行控制过程空气动力建模与计算研究［R］.中国国防科技报告，绵阳：中国空气动力研究与发展中心超高速研究所，2013. Li Z H，Wu JL，Peng A P.Aerodynamicsmodeling and calculation during low-orbit flying control of TG-1 target spacecraft［R］.Technical Report，Hypersonic Aerodynamics Institute，China Aerodynamic Research and Development Centre，2013.（in Chinese）

［5］Cercignani C.The Boltzmann Equation and its Applications［M］.Springer Verlag，Berlin，1988：55-58.

［6］应纯同.气体输运理论及应用［M］.北京：清华大学出版社，1990：123-136. Ying C T.Gas transport theory and application［M］.Tsinghua Univ.Press，Beijing，1990：123-136.（in Chinese）

［7］Bird G A.Approach to translational equilibrium in a rigid sphere gas［J］.Physics of Fluids（1958-1988），1963，6（10）：1518-1519.

［8］Pham-Van-Diep G，Erwin D，Muntz E P.Nonequilibrium molecular motion in a hypersonic shock wave［J］.Science，1989，245（4918）：624-626.

［9］Carlson A B，Hassan H A.Direct simulation of reentry flows with ionization［R］.AIAA 90-144，1990.

［10］Haas B L，Boyd ID.Models for direct Monte Carlo simulation of coupled vibration-dissociation［J］.Physics of Fluids A：Fluid Dynamics，1993，5（2）：478-489.

［11］Bird G A.Molecular Gas Dynamics and Direct Simulation of Gas Flows［M］.Oxford University Press，1994：78-189

［12］樊菁，沈青.过渡领域高超声速圆柱绕流的直接模拟［J］.空气动力学学报，1995，13（4）：405-413. Fan J，Shen Q.Direct simulation of hypersonic transitional flow around circular cylinder，Acta Aerodynamica Sinica，1995，13（4）：405-413.（in Chinese）

［13］李志辉，吴振宇.阿波罗指令舱稀薄气体动力学特征的蒙特卡罗数值模拟［J］.空气动力学学报，1996，14（2）：230-233. Li Z H，Wu Z Y.DSMC simulation of hypersonic rarefied flow past Apollo-CM，Acta Aerodynamica Sinica，1996，14（2）：230-233.（in Chinese）

［14］Ivanov M S，Gimelshein S F.Computational hypersonic rarefied flows［J］.Annual Review of Fluid Mechanics，1998，30（1）：469-505.

［15］吴其芬，陈伟芳.高温稀薄气体热化学非平衡流动的DSMC方法［M］.长沙：国防科技大学出版社，1999：165-176. Wu Q F.Chen W.F.DSMC method for high-temperature thermochemical nonequilibrium flow of rarefied gas［M］.National University of Defense Technology Press，ChangSha，1999：165-176.（in Chinese）

［16］Sun Q，Fan J，Boyd ID.Improved sampling techniques for the direct simulation Monte Carlomethod［J］.Computers& Fluids，2009，38（2）：475-479.

［17］梁杰，阎超，李志辉，等.稀薄过渡流区横向喷流干扰效应数值模拟研究［J］.空气动力学学报，2013，31（1）：27-33. Liang J，Yan C.，Li Z H，et al.Numerical investigation of lateral jet interaction effects in rarefied transition flow regime，Acta Aerodynamica Sinica，2013，31（1）：27-33.（in Chinese）

［18］Wong B C，Long L N.Direct simulation Monte Carlo（DSMC）on the connection machine［R］.AIAA 92-0564，1992.

［19］Scanlon T J，RoohiE，White C，etal.An open source，parallel DSMC code for rarefied gas flows in arbitrary geometries［J］.Computers&Fluids，2010，39（10）：2078-2089.

［20］李志辉，杨彦广，梁杰，等.飞船返回舱高空稀薄气动特性研究［C］.全国载人航天工程气动工作总结暨关键问题研讨会，6月，北京，2004. Li Z H.Yang Y G.，Liang J，etal.Study of rarefied aerodynamics for spacecraft capsule in high-altitude，Symposium of front key problems for national manned spaceflight，Beijing，June，2004.（in Chinese）

［21］Yoon S，Gnoffo P A，White JA，et al.Computational challenges in hypersonic flow simulations［R］.AIAA 2007-4265，2007.

［22］Bhatnagar P L，Gross E P，Krook M.A model for collision processes in gases.I.Small amplitude processes in charged and neutral one-component systems［J］.Physical review，1954，94（3）：511-525.

［23］Holway Jr L H.New statistical models for kinetic theory：methods of construction［J］.Physics of Fluids（1958-1988），1966，9（9）：1658-1673.

［24］Shakhov E M.Generalization of the Krook kinetic relaxation equation［J］.Fluid Dynamics，1968，3（5）：95-96.

［25］Yang JY，Huang JC.Rarefied flow computations using nonlinearmodel Boltzmann equations［J］.Journal of Computational Physics，1995，120（2）：323-339.

［26］Titarev V A，Shakhov E M.Heat transfer and evaporation from a p lane surface into a half-space upon a sudden increase in body temperature［J］.Fluid dynamics，2002，37（1）：126-137.

［27］Mieussens L.Discrete-velocitymodels and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries［J］.Journal of Computational Physics，2000，162（2）：429-466.

［28］李志辉，张涵信.稀薄流到连续流的气体运动论统一数值算法初步研究［J］.空气动力学学报，2000，18（3）：255-263. Li ZH.，Zhang H X.Study on gas-kinetic algorithm for flows from rarefied transition to continuum，Acta Aerodynamica Sinica，2000，18（3）：255-263.（in Chinese）

［29］Li Z，Zhang H.Study on gas kinetic algorithm for flows from rarefied transition to continuum［C］//Proc.of 22ndInternational Symposium on Rarefied Gas Dynamics，Sydney，Australia，Jul.9-14，2000：628-636.

［30］李志辉.从稀薄流到连续流的气体运动论统一数值算法研究［D］.绵阳：中国空气动力研究与发展中心研究生部，2001. Li Z H.Study on gas kinetic unified algorithm for flows from rarefied transition to continuum，Ph.D thesis，China Aerodynamics Research and Development Center，Jan.，2001.（in Chinese）

［31］Li ZH，Zhang H X.Study on gas kinetic unified algorithm for flows from rarefied transition to continuum［J］.Journal of Computational Physics，2004，193（2）：708-738.

［32］李志辉，张涵信.基于Boltzmann模型方程的气体运动论统一算法研究［J］.力学进展，2005，35（4）：559-576. Li Z H.，Zhang H X.Study on gas kinetic numerical algorithm using Boltzmann model equation.Advances in Mechanics 2005；35（4）：559-576.（in Chinese）

［33］Li Z，Zhang H.Gas-kinetic numericalmethod for solvingmesoscopic velocity distribution function equation［J］.Acta Mechanica Sinica，2007，23（2）：121-132.

［34］Li Z H，Zhang H X.Gas-kinetic numerical studies of threedimensional complex flows on spacecraft re-entry［J］.Journal of Computational Physics，2009，228（4）：1116-1138.

［35］李志辉，梁杰，李四新，等.箔条云跨流域整体气动特性计算研究［J］.空气动力学学报，2011，29（01）：59-67. Li Z H.，Liang J，Li SX，et al.Computing study on holistic aerodynamics of chaff cloud covering various flow regimes，Acta Aerodynamica Sinica，2011，29（1）：59-67.（in Chinese）

［36］Li Z H.Gas-kinetic unified algorithm for re-entering complex flows covering various flow regimes by solving boltzmann model equation［M］//Hall J.Advances in Spacecraft Technologies.Croatia：InTech，2011：273-332.

［37］李志辉，彭傲平，吴俊林，等.稀薄气体自由分子流到连续流跨流域气动力热绕流统一算法研究［J］.载人航天，2013，19（2）：81-91. Li ZH.，Peng A P，Wu JL，etal.Gas-kinetic unified algorithm for aerothermodynamics from free molecular regime to continuum，Manned Spaceflight，2013，19（2）：81-91.（in Chinese）

［38］Gentry A E，Smyth D N，OliverW R.The Mark IV supersonic-hypersonic arbitrary-body program，Vol.11：Program Formulation［R］.AFFDL TR-73-159，1973.

［39］Kashkovskyt A V.Computational tools for rarefied aerodynamics［C］//Rarefied Gas Dynamics：Technical Papers from the Proceedings of the Eighteenth International Symposium on Rarefied Gas Dynamics，University of British Columbia，Vancouver，British Columbia，Canada，1992：115-126.

［40］Mosanov SV，Freedlender O G，Nikiforov A P，et al.Experimental determination ofmomentum transfer coefficients in hypersonic freemolecular flow and distribution function recovery of reflected molecules［C］//Proceedings of the XIII International Symposium Rarefied Gas Dynamics，Vol.1，Plenum，New York，1985：669-676.

［41］谢砚儒，李志辉.尖钝锥诱饵再入阻力特性工程计算研究［R］.绵阳：中国空气动力研究与发展中心，中国国防科技报告，1992. Xie Y R.，Li ZH.，Study ofengineering calculation of re-entry drag characteristics for sharp-cone body［R］.Technical Report，Hypersonic Aerodynamics Institute，China Aerodynamic Research and DevelopmentCentre，1992.（in Chinese）

［42］李志辉.箔条及箔条云整体气动特性建模方法研究［R］.绵阳：中国空气动力研究与发展中心，中国国防科技报告，2005. Li Z H.，Modeling study on holistic aerodynamic characteristic of foil and foil-cloud cluster［R］.Technical Report，Hypersonic Aerodynamics Institute，China Aerodynamic Research and Development Centre，2005.（in Chinese）

［43］谢砚儒，刘涛.弹头稀薄气动特性的工程计算方法研究［R］.绵阳：中国空气动力研究与发展中心，中国国防科技报告，1990.Xie Y R.，Liu T.，Study ofengineering calculation method of rarefied aerodynamics for bullet［R］.Technical Report，Hypersonic Aerodynamics Institute，China Aerodynamic Research and Development Centre，1990.（in Chinese）

［44］粱杰，李志辉.目标稀薄气体特性工程计算［R］.绵阳：中国空气动力研究与发展中心，国防科技报告，1999. Li Z H.，Liang JEngeering calculation of target rarefied aerodynamic characteristics［R］.Technical Report，Hypersonic Aerodynamics Institute，China Aerodynamic Research and Development Centre，1999.（in Chinese）

Unified Modeling and Calculation of Aerodynam ics Characteristics during Low-O rbit Flying Control of the TG Vehicle

LIZhihui1，2，WU Junlin1，PENG Aoping1，TANG Geshi3
（1.Hypervelocity Aerodynamics Institute，China Aerodynamics Research and Development Center，Mianyang 621000，China；2.National Laboratory for Computational Fluid Dynamics，Beijing 100191，China；3.Science&Technology on Aerospace Flight Dynamics Laboratory，Beijing Aerospace Control Center，Beijing100094，China）

The surface analyticalmethod，considering surface shielding effect of the complex structure and the modified Boettcher/Legge nonsymmetric bridge correction function，was proposed to computationally model aerodynamic characteristics of the large-scale Tiangong spacecraft of irregular plate-capsule assembly.The complex configuration processing and computing rules of surface element aerodynamic coefficientswere setup by developing a general triangle elementapproximation for complex shapes.A unified fast algorithm and computing software for aerodynamic characteristics of large-scale complex spacecraft structures were developed，by computing correction of local correlation parameters during Tiangong spacecraft’s low-earth orbit flight control，in which the DSMC method and the gas-kinetic unified algorithm solving the Boltzmann model equation were applied.Itis indicated that，by computing the aerodynamics of the Tiangong spacecraft in 300～200 km altitude-orbit flight process，with various flying heights，various angles of attack and various angles between the panel and the principal axis，the aerodynamic coefficients varied remarkably，with a range of 8%～50%，with the altitude change，and the atmospheric drag in low-earth orbit flight control was the key factor of the spacecraft orbit prediction accuracy.It is validated that，if the fixed and constant unity aerodynamic coefficient is used for long-term orbit flight of a large spacecraft，the error accumulation will be huge，and protectivemeasures are needed.

TG-1 target spacecraft；aerodynamic characteristic during low-orbit control；local bridge function；correlation parameter evaluation；engineering calculation；DSMCmethod；Gas-Kinetic U-nified Algorithm（GKUA）

V411.4

A

1674-5825（2015）02-0106-09

2014-09-12；

2015-02-26

国家973计划（2014CB744100）；国家自然科学基金（11325212、91016027）；国防基础科研基金（51313030104）作者简介：李志辉（1968-），男，博士，研究员，研究方向为跨流域空气动力学。E-mail：zhli0097@x263.net