类天宫航天器迎风面积建模及变化特性分析

2020-09-02李志辉

载人航天 2020年4期

杨成，李勰，孙军，李志辉

（1.航天飞行动力学技术重点实验室，北京100094；2.北京航天飞行控制中心，北京100094；3.中国空气动力研究与发展中心超高速空气动力研究所，绵阳621000；4.空气动力前沿技术研究中心国家计算流力学实验室，北京100191）

1 引言

大气阻力是航天器在轨运行的一种重要非保守力，对轨道处理特别是中长期轨道预报精度具有重大影响［1-3］。航天器迎风面积是计算大气阻力的重要因素，定义为沿航天器相对于当地大气速度方向航天器几何结构的投影面积。

大型航天器一般由多功能舱段组成，且具有较大的太阳帆板，构成一个复杂的几何外形。航天器在轨运行中姿态的变化导致迎风面也在不断变化，同时太阳帆板会根据太阳方位不断调整角度，对航天器的外形产生影响。因此，在航天器运行过程中，其迎风面积将发生非常大的变化。

工程计算中通常将迎风面积作为待估计参数，使用大气阻力系数与面质比乘积的等效弹道系数法［3-4］，结合历史轨道数据，选取一定时间范围（一周或更长时间），根据该时段大型航天器轨道半长轴衰减情况，估算航天器等效弹道系数，进行定轨与轨道外推寿命预报。根据前期数据将迎风面积作为参数进行拟合，在航天器飞行预报阶段会引入较大误差，从而影响轨道预报精度。因此，根据航天器的实际几何外形、运行姿态等信息，实时地计算航天器迎风面积，对提高轨道处理精度有重要意义。

航天器一般包括本体和太阳帆板两部分，为便于处理，本体简化为圆柱形，太阳帆板简化为长方体。唐歌实等［3］提出了一种航天器俯仰滚动偏航角均为0°的正飞姿态情况下的等效迎风面积计算方法。首先计算太阳帆板在绕旋转轴一周的面积变化，求取平均值，加上本体的截面积即为迎风面积，实际所得是航天器正飞状态的平均迎风面积。这种方法是不同状态下的数值平均，并不能代替实际情况下的迎风面积变化；另一类为近似外形几何投影方法［5-6］。首先将航天器外形简化，分解为方形、柱形等基础形状的组合。分别对这些简单形状在航天器速度方向进行投影处理，获得各自在投影面上的投影区域，得到航天器的迎风面积。这种方法采用简单几何体来对航天器外形进行简化，在计算精度上难以满足精密定轨的要求，并且计算过程较为复杂，缺乏通用性。

杨成等［7］提出了一种应用阴影图计算迎风面积的方法，通过对航天器几何外形进行精确网格建模，同时考虑航天器姿态变化和帆板转动情况，利用几何投影方式对航天器进行投影；根据投影面积，确定航天器的迎风面积，并以天宫一号航天器为例，对其理论上的迎风面积进行了计算。朱战霞等［8］基于微元划分思想提出了射击线扫描法以求解有效迎风面积。这类方法在外形描述和投影计算方面具有较好的通用性，但计算过程需要的航天器状态参数较多，在实际工程应用中受测控条件影响［9］，在测控弧段外航天器状态数据缺失，此时如何建立合适的计算模型并解决模型参数获取是个问题。

本文提出了一种类天宫航天器迎风面积建模方法。该模型按照地球经纬度定义方式，在航天器本体系建立经纬度网格，将计算迎风面积的投影方向转换为经度和纬度，与太阳帆板的转动角度一起构成迎风面积模型的参数；在模型参数取值范围内按照一定的步长，利用几何投影［7］计算所有状态的迎风面积数据，形成航天器的迎风面积模型；在轨道处理过程中，按照航天器长期飞行期间的控制模式，计算航天器各时刻的姿态和帆板转动角度，确定对应的模型参数，通过插值获得迎风面积。

2 迎风面积几何投影算法

由于航天器在轨飞行过程中姿态变化，帆板等大尺寸部件发生相对运动，迎风面积将发生较大的变化。精确迎风面积计算方法［7］通过构建航天器精确的几何外形，根据实际航天器姿态、活动部件运动变形情况，在相对速度方向进行几何投影，利用产生的投影面积即可确定对应迎风面积的大小。基本原理包括模型表示和几何投影两部分，如图1所示。利用非结构化三角网格表征航天器的几何外形，使用树结构表征航天器部件之间的关系，根据航天器姿态、帆板转动角度，在航天器相对大气速度方向进行几何投影，获得航天器该状态下的迎风面积。

2.1 模型表示

使用树结构表示航天器活动部件之间的层次连接关系和运动状态。树结构的节点对应一个航天器部件，使用非结构三角网格表征该部件实体的几何外形信息，可逼近任何复杂外观和结构的航天器，提高模型表示精度［10-11］，同时降低网格数量。利用非结构网格表示航天器如图2所示，几何形状信息里包括顶点和面片两部分信息，其中顶点信息包括所有顶点的三维坐标，面片信息包括组成几何形状外表面所有面片的顶点序号。按照航天器的结构，对不同部分使用单独的三角网格进行表示，平缓表面网格稀疏，弯曲表面使用更稠密网格拟合。

图1 类天宫航天器迎风面积计算方法Fig.1 Calculation method of cross section area of TG-like spacecraft

图2 航天器部件几何外形的非结构网格Fig.2 Unstructured triangular surface elementmesh of spacecraft com ponent

2.2 几何投影

航天器迎风面积几何投影原理如图3所示。以航天器相对大气的速度矢量方向为投影方向，利用计算机图形学方法［11］，对航天器几何网格进行投影，在垂直于投影方向的投影面产生航天器投影的影像。通过计算影像所占面积大小，即可确定航天器的投影面积。

3 迎风面积建模方法

轨道计算一般包括定轨和预报两部分，需要对一定弧段内若干时刻的迎风面积进行计算。在轨道处理过程中，直接采用2.1节的方法，若频繁调用迎风面积计算模块，以图形学几何投影原理来获取某个时刻的迎风面积数值，将影响轨道处理速度。

图3 航天器几何投影示意图Fig.3 Geometric projection diagram of spacecraft

根据类天宫航天器外形和在轨飞行特点，本文使用投影经度、投影纬度、帆板角度3个参数对计算状态进行表示，根据各参数取值范围，以一定步长计算状态空间中所有可能状态的迎风面积，离线建立迎风面积模型。在轨道处理过程中，根据飞行状态和飞行模式，计算航天器相对大气的速度、飞行姿态［12］、帆板转角，并转化为投影经度、投影纬度、帆板角度3个模型参数，使用线性插值的方式从模型中获得迎风面积，如图4所示。

图4 类天宫航天器迎风面积建模方法Fig.4 M odeling on cross section area for Tiangongtype spacecraft

3.1 离线建模

类天宫航天器的主要结构包括本体和太阳帆板，本体近似圆柱体，太阳帆板连接在本体两侧，其迎风面积可认为与投影方向和帆板转角有关［12］。引入经纬度球的概念，将投影方向矢量从三维转化为经度、纬度二维，降低模型参数个数。参照地球地理信息系统中的经纬度定义［3］，在航天器本体坐标系建立经纬度球，其中：x轴和z轴构成的平面指向x轴正方向对应经度为0°，y轴和z轴构成的平面指向y轴正方向对应经度为90°，经度范围为［0°，360°）；x轴和y轴构成的平面对应纬度为0°，沿z轴正方向纬度增加，北极为90°，沿z轴负方向纬度降低，南极为-90°，纬度范围为［-90°，90°］，如图5所示。

图5 航天器本体系下的相对速度Fig.5 Relative velocity of spacecraft in body fixed coordinate system

航天器相对大气的速度归一化为单位向量［Vx，Vy，Vz］，即为计算迎风面积的投影方向，对应航天器经纬度球面的一个点，可以用该点的经度Pj和纬度Pw进行表示。［Vx，Vy，Vz］和［Pj，Pw］之间的转化关系如下：

已知投影方向对应的经度、纬度和帆板转角即可确定迎风面积计算状态，本文的迎风面积模型使用投影经度、投影纬度、帆板转角作为模型参数，记为［Pj，Pw，Rf］，其中Pj∈ ［0，360），Pw∈［-90，90］，Rf∈ ［0，180）。

按照模型参数Pj、Pw、Rf的取值范围，对每个参数以一定步长取值，产生所有参数对，根据公式（1）计算投影方向，然后利用投影方向和帆板转角，应用几何投影法计算该参数对应状态的迎风面积，形成迎风面积模型。

3.2 实时计算

在对轨道数据进行定轨和利用轨道根数进行轨道预报过程中，都需要轨道处理弧段内的航天器迎风面积数据。航天器在轨长期飞行期间，为保证太阳帆板获取太阳能的效率，姿态控制模式在正飞和偏航2种模式之间切换，如图6所示。定义地心到太阳的方向为太阳矢量方向，当太阳矢量与航天器轨道面之间夹角在阈值以内时，采取正飞模式，偏航、俯仰、滚动3个姿态角为0°，帆板通过调整角度，使得太阳与航天器连线位于帆板转轴与帆板法线组成的平面内，以便在保持正飞姿态前提下太阳帆板获得最大太阳能量；当太阳矢量与航天器轨道面之间夹角大于阈值时，采取偏航模式，俯仰、滚动保持0°，同时调整偏航角度和帆板转动角度，使得太阳与航天器连线垂直于帆板平面，太阳帆板获得最大太阳能量。

图6 航天器飞行模式与太阳方位关系示意图Fig.6 Diagram of relationship between spacecraftflightmode and solar orientation

在地球惯性系下，记航天器位置为P=［x，y，z］T，速度为下面给出按照飞行模式确定迎风面积模型参数的方法。

3.2.1 正飞模式

正飞模式下，航天器x轴指向速度V方向，z轴指向地心方向，与x轴垂直，y轴与z、x轴形成右手坐标系。航天器相对大气速度转换到航天器本体坐标系下，按照公式（2）可计算出投影经纬度［Pj，Pw］；将太阳方向矢量投影到航天器z、x轴构成的平面上，并归一化为单位向量，在航天器本体系下的坐标记为S=［Sx，0，Sz］T，则太阳帆板转动角度Rf计算如式（3）所示：

3.2.2 偏航模式

偏航模式下，航天器将根据太阳方向进行偏航，使得航天器的x轴沿z轴转动一个偏航角到达x'位置，与太阳方向矢量S和z轴共面，从而帆板转动一定角度就能使得帆板法线与太阳矢量S重合，获得最大太阳能转化效率，如图7所示。x'为太阳方向矢量S在x-y平面的投影，通过S与X矢量的内积LSX=S·X和S与Y矢量的内积LSY=S·Y，可计算出偏航角度ψ为式（4）所示：

图7 偏航模式下偏航角与太阳方位示意图Fig.7 Diagram of yaw angle and sun orientation in yaw mode

航天器偏航ψ角度后，其本体系坐标轴转动到x'、y'和z位置，y'为帆板转动轴方向，由于太阳方向矢量S在x'和z轴构成的平面内，y'与S是垂直关系，通过转动帆板即可使得帆板法向量与太阳方向矢量S重合。同理，根据公式（3）可计算太阳帆板转动角度。

利用航天器飞行模式确定该时刻迎风面积模型的参数［Pj，Pw，Rf］，通过线性插值，得到对应迎风面积。

4 模型方法验证

天宫一号为中国目前最大的复杂航天器，主要由本体、太阳帆板组成，其中太阳帆板可绕本体系的y轴旋转。本体由资源舱和实验舱组成，如图8所示，总长约为10 m，直径约为3 m，太阳帆板展开后长约为20 m。由于该航天器尺寸较大，大气阻力对轨道的影响较为显著。本文以天宫一号为研究对象，验证类天宫航天器迎风面积模型，分析其在轨飞行期间迎风面积的变化特性。

图8 类天宫航天器几何外形Fig.8 Geometric surface of Tiangong-type spacecraft

4.1 模型建立

天宫一号本体部分的几何外形使用8342个三角面片表示，帆板部分几何构型比较规则，包含648个三角面片。迎风面积模型使用投影经度、投影纬度、帆板转角作为模型参数，记为［Pj，Pw，Rf］，取值范围Pj∈［0，360），Pw∈［-90，90］，Rf∈ ［0，180）。参数步长设置为1°，这样Pj存在360种取值可能性，Pw存在181种取值可能性，Rf存在180种取值可能性，组合起来对应360×181×180种状态。几何投影功能基于OpenGL图形库实现［10］，处理电脑配置为CPU i3、内存4G、显存1G，计算所有状态的迎风面积耗时2 h。

4.2 在轨运行期间迎风面积变化情况

利用2015年天宫一号在轨飞行数据，验证本文迎风面积计算方法。数据为153天的连续观测弧段，期间航天器的飞行姿态模式分为正飞和偏航两种。整个实验过程中，随着太阳矢量方向与轨道面夹角的变化，天宫一号进行了多次正飞和偏航模式的切换。

首先根据航天器飞行模式，对其偏航角度和帆板转动角度的计算方法进行实验。受航天器机械控制系统的影响，实际角度与本文计算角度会存在一定误差。2015年11月1日天宫一号处于偏航模式，在16时左右存在两段约8 min的测控弧段，偏航角度和帆板转动角度数据通过遥测进行下传，使用该数据与本文方法计算的偏航角度和帆板转动角度数据进行对比验证。

图9为帆板转动角度对比情况。天宫一号帆板安装在资源舱左右两侧，将左右帆板转动当作相同的，实际由于受2套独立机构驱动，左右帆板的实际转动角度会有不同。图9中虚线为本文计算的帆板转动角度，实线为左右帆板通过遥测下传的实际转动角度。经统计，帆板转动角度的平均误差为2.65°。图10为偏航计算角度和遥测下传实际角度的对比，偏航角度平均误差为0.53°。从图9、图10可以看出，按照本文方法计算角度和遥测下传实际角度在角度变化上趋势一致，角度值误差较小。

图9 迎风面积方法计算的帆板转动角度和遥测实际数据比较Fig.9 Comparison between telemetry data and calculation of the sailboard rotation angle

图10 迎风面积模型偏航角度和遥测实际数据比较Fig.10 Com parison between telemetry data and calcu lation of the yaw angle

每种模式飞行期间内，可以根据该模式计算该时刻航天器的姿态和帆板转动角度，得到模型参数，然后采用插值方式计算迎风面积。在153天实验弧段内，每隔30 s进行一次时间采样，共计算440 640个时刻的迎风面积，结果如图11所示。图中可以看出，正飞过程迎风面积变化范围更大，可以达到理论上的极值，该模式下近似为本体前端面和太阳帆板产生迎风面积，本体前端面为圆柱形的截面积保持不变，太阳帆板循环转动360°，当帆板转动到水平和垂直状态，迎风面积就会接近最小值和最大值；偏航模式中，由于本体存在偏航转动，本体的前端面和侧面都会产生迎风面积，所以迎风面积的最小值比正飞模式大，由于偏航时本体和太阳帆板之间的遮挡抵消，迎风面积的最大值比正飞模式小。

图11 在轨飞行过程中迎风面积变化情况Fig.11 Variation of cross section area during orbit flight

图12为每隔3 h统计一次的迎风面积平均值，由于天宫一号的轨道周期约1.5 h，该平均值反映的是轨道周期时间尺度上迎风面积的变化情况。可以看出，正飞模式下，航天器姿态固定，在每个轨道周期内太阳帆板旋转360°，以轨道周期时间范围进行平均后，迎风面积平均值约为36.5～37.5 m2，变化幅度很小；偏航模式下，轨道面与太阳位置关系确定了航天器的偏航角度，迎风面积平均值约为32.6～40.3 m2，在一个偏航模式范围内先变小再变大，变化范围比正飞模式大。分析可以看出，正飞模式下天宫一号迎风面积受帆板转动影响变化幅度很大，但轨道周期尺度的均值几乎恒定；偏航模式下迎风面积的下限更高，轨道周期尺度的均值具有较大的变化幅度。因此，天宫一号迎风面积在不同飞行模式下具有不同的变化情况，轨道计算中使用固定等效面积或者弹道系数的方法不能反映迎风面积的这种变化，带来了较大误差，需要根据实际状态实时计算。