CN102169047A - 一种10n推力器羽流场热效应及动力学效应确定方法 - Google Patents

Abstract

一种10N推力器羽流场热效应及动力学效应确定方法，(1)建立卫星几何尺寸及卫星坐标系、10N推力器本体坐标系下羽流场位置及指向、太阳翼本体坐标系下其几何尺寸及对日定向角度，以及通信天线本体坐标系下其几何尺寸与形面法矢参数；(2)将步骤(1)中涉及的推力器、太阳翼及通信天线的参数转换到卫星坐标系下；(3)对太阳翼及通信天线进行网格划分成单元，计算所有网格的羽流场原点至太阳翼、通信天线表面第i单元格的矢量

、矢量

与太阳翼、通信天线表面第i单元格的法向矢量间夹角β、矢量与推力器轴线间夹角κ；(4)根据上述计算的每个单元格的矢量、β、κ，计算推力器的羽流场干扰力和热流量分布。

Description

一种10N推力器羽流场热效应及动力学效应确定方法

技术领域

本发明涉及卫星总体布局设计及其优化，并同时可供卫星热设计、卫星控制设计进行参考。

背景技术

羽流是指卫星发动机工作在高真空下时喷射出来的气流，它在喷嘴处形成的形状像羽毛一样，所以称之为羽流。羽流中的成分复杂，包含分子污染物和粒子污染物，如未完全燃烧的燃料液滴，燃烧室壁面颗粒物等。在空间环境下，喷射出的气体将迅速膨胀，变得越来越稀薄，并且可以延伸得很远，形成羽流，有一部分气体会形成回流并撞击在航天器自身表面上，直接影响航天器的工作性能和寿命，特别是对热控涂层以及太阳能电池阵表面产生有害影响。

对于羽流场分析问题，曾经广泛使用基于连续介质假设的Simons方法和特征线法(MOC)。Simons方法是在1972年提出的一种羽流点源模型，这是一种工程估算方法。MOC方法是在连续介质流体力学理论的基础上发展起来的，其同样是将N-S偏微分控制方程组将退化为常微分方程组的工程估算方法。基于连续介质假设的Simons模型和MOC方法虽然处理速度相对较快，但对于众多参数的依赖性较大，较难在工程应用早期对羽流场效应做出准确有效的判断。1990年代后期，计算机技术得到了迅速发展，一般不再使用基于连续介质假设的Simons模型和MOC方法。

由于羽流场包含很大部分的过渡流区和自由分子流区，连续介质方法不能准确求解这些流区，必须以分子的观点来描述。迄今为止，直接模拟蒙特卡罗(Direct Simulation Monte Carlo，DSMC)方法是应用最广泛的分子型气体流动模拟方法。1963年，澳大利亚学者G.A.Bird提出了DSMC方法，该方法通过模拟大量仿真分子的运动和相互碰撞过程来模拟真实气体流动。但限于当时计算机水平的约束，DSMC方法一直未能得到广泛的应用与信任。随着DSMC方法的逐步发展，以及计算机技术的发展，伴随着1992年，Penko和Chung分别通过实验数据验证了DSMC羽流场仿真结果，DSMC终于得到广泛的接受。DSMC方法有一个致命的缺点就是由于大量的计算而造成对于计算机能力的需求增加。目前，通常使用计算机并行计算克服DSMC方法计算量大的难题，但对于卫星工程应用来说，这又是不太现实的，尤其是不能够在卫星早期设计阶段给出高效、明确的分析结论。

发明内容

本发明的技术解决问题是：对于羽流场效应的分析主要涉及热流影响和动力学影响两个方面。以往的分析方法中，要么是模型本身具有较大的局限性(例如Simons方法)，要么是对于计算机的计算和存储能力具有很高的要求(例如DSMC方法)，一时很难较快较容易的直观得到羽流场效应的工程评价结果。通过本发明方法能够在卫星设计早期进行快速的羽流影响效应评估，后续根据评估结果采取必要防护措施，相对于在制造后期再引入低效的防护手段可以显著地节约成本，提高效益。

本发明的技术解决方案是：一种10N推力器羽流场热效应及动力学效应确定方法，步骤如下：

(1)首先建立卫星几何尺寸及卫星坐标系；同时分别建立10N推力器本体坐标系下羽流场位置及指向、太阳翼本体坐标系下其几何尺寸及对日定向角度，以及通信天线本体坐标系下其几何尺寸与形面法矢参数；

(2)建立推力器、太阳翼及通信天线间的坐标系转换关系，将步骤(1)中涉及的推力器、太阳翼及通信天线的参数转换到卫星坐标系下；

(3)对太阳翼及通信天线进行网格划分成单元，从第一个单元格i＝1开始执行步骤(4)～(6)；

(4)计算羽流场原点至太阳翼、通信天线表面第i单元格的矢量

(5)计算矢量与太阳翼、通信天线表面第i单元格的法向矢量间夹角β；

(6)计算矢量与推力器轴线间夹角κ；

(7)令i＝i+1，重复执行步骤(4)～(6)，直至所有单元格计算完毕，执行步骤(8)；

(8)根据上述计算的每个单元格的矢量β、κ，计算推力器的羽流场干扰力和热流量分布，羽流场干扰力密度分布函数F(r，κ，β)和热流密度分布函数H(r，κ，β)分别如下：

$F(r,\mathrm{\κ},\mathrm{\β})=F_1\left(r_0\right)\·{\left(\frac{r_0}{r}\right)}^2\·F_2\left(\mathrm{\κ}\right)\·{\cos }^2\mathrm{\β}$

$H(r,\mathrm{\κ},\mathrm{\β})=H_1\left(r_0\right)\·{\left(\frac{r_0}{r}\right)}^2\·H_2\left(\mathrm{\κ}\right)\·\cos \mathrm{\β}$

式中，r₀为实验中所取的距离，F₁(r₀)、H₁(r₀)分别为对应r₀的干扰力密度、热流密度试验结果；r为矢量

的标量值；F₂(κ)、H₂(κ)分别为仅与κ相关的推力器羽流测试的归一化试验结果干扰力和热流量的拟合函数。

所述的仅与κ相关的推力器羽流测试的归一化试验结果干扰力密度的拟合函数F₂(κ)为：

J₀(κ)为κ的0阶贝塞尔函数，κ以“度”表示。

所述的推力器羽流测试的归一化试验结果热流密度的拟合函数H₂(κ)为：

J₀(κ)为κ的0阶贝塞尔函数，κ以“度”表示。

本发明与现有技术相比有益效果为：

(1)本发明通过卫星布局构型阶段引入10N推力器羽流场效应的工程评价方法，具有积极的效果：结合星表太阳翼、通信天线的构型和布局，能够快速有效地判断10N推力器羽流场对其所产生的热效应及动力学效应。其并不依赖于10N推力器本身的设计特征，且计算过程对计算机的配置无特殊要求，使得计算过程效率大大提高。通过与国外同类卫星10N推力器羽流场计算软件结果进行比对，具有相似一致性，从而验证本方法的正确性，及工程实用性。

(2)本发明中对羽流场分布的描述基于测试结果。将测试结果首先进行数值插值处理，并将数值插值结果使用解析函数进行表达，从而进一步提高其工程适用范围，提高计算速度与效率。

附图说明

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

图2为本发明羽流场分析模型示意图；

图3为本发明模拟的卫星太阳翼布局示意图；

图4为本发明模拟的通信天线布局示意图；

图5为本发明对羽流场试验值的三次样条函数插值结果；

图6为本发明模拟出的羽流热流密度场分布示意图。

具体实施方式

下面结合附图及实施例详细介绍本发明的实现过程。

本发明进行羽流场分布的描述，其最终目的是为了分析羽流场对卫星表面太阳翼、通信天线等敏感器件的热效应及动力学效应影响。如图1所示，本发明具体步骤如下：

(1)首先建立卫星几何尺寸及卫星坐标系；同时分别建立10N推力器本体坐标系下羽流场位置及指向、太阳翼本体坐标系下其几何尺寸及对日定向角度，以及通信天线本体坐标系下其几何尺寸与形面法矢参数；

对于通信卫星而言，其参考坐标系(即卫星坐标系)原点定义为与运载对接框的理论圆心，X轴向指向卫星同步轨道的飞行方向，Z轴向指向卫星同步轨道的星下点方向，X、Y、Z轴呈右手螺旋分布。羽流场在本体坐标系下分布如表1所示。太阳翼对日定向角度一般为[-180°，180°]。其余参数根据卫星设计的总体性能参数确定，会有一定的变化。例如某型号卫星本体长、宽、高尺寸分别为2360mm，2100mm和3600mm；其上安装14台10N推力器，推力器口面中心与卫星本体间的坐标转换关系如表所示；太阳翼几何尺寸为2360mm×9900mm，安装于卫星的南北板；某通信卫星安装4副通信天线，其基本几何尺寸如表所示。

表110N推力器口面中心与卫星本体坐标转换关系

表2通信天线几何尺寸

(2)建立推力器、太阳翼及通信天线间的坐标系转换关系，将步骤(1)中涉及的推力器、太阳翼及通信天线的参数转换到卫星坐标系下；

推力器与卫星本体的坐标转换关系矩阵如表所示，南北太阳翼与卫星本体的坐标转换关系矩阵如表3所示，通信天线与卫星本体的坐标转换关系矩阵如表4所示。模拟的太阳翼布局与通信天线布局如错误！未找到引用源。及错误！未找到引用源。所示。

表3南北太阳翼与卫星本体坐标系转换关系矩阵

表4通信天线与卫星本体坐标系转换关系矩阵

(3)对太阳翼及通信天线进行网格划分，从第一个单元格i＝1开始执行步骤(4)～(6)；

(4)计算羽流场原点至太阳翼、通信天线表面第i单元格的矢量

$\stackrel{\→}{r}=\left(\begin{array}{ccc}{x}_t-x_p& y_t-y_p& z_t-z_p\end{array}\right)\·{\left(\begin{array}{ccc}\hat{x}& \hat{y}& \hat{z}\end{array}\right)}^T$

其中x_t、y_t、z_t分别表示太阳翼、通信天线第i单元格中心在卫星坐标系下的三个坐标，x_p、y_p、z_p分别表示羽流场原点在卫星坐标系下的三个坐标，

分别表示卫星坐标系下三个方向上的单位矢量。

(5)计算矢量与太阳翼、通信天线第i单元格的法向矢量间夹角β；

$\mathrm{\β}=\mathrm{\π}-\arccos (\frac{\stackrel{\→}{r}}{\left|\stackrel{\→}{r}\right|}\·\hat{n})$

其中为卫星坐标系下太阳翼、通信天线第i单元格中心的法向单位矢量。

(6)计算矢量与推力器轴线间夹角κ；

$\mathrm{\κ}=\arccos (\frac{\stackrel{\→}{r}}{\left|\stackrel{\→}{r}\right|}\·{\hat{P}}_0)$

其中为卫星坐标系下羽流场方向的单位矢量。

(7)令i＝i+1，重复执行步骤(4)～(6)，直至所有单元格计算完毕，执行步骤(8)；

(8)根据上述计算的每个单元格的矢量、β、κ，计算推力器的羽流场干扰力和热流量分布，下面重点详细的介绍此步骤的实现。

本发明的实验模型如图2所示，由图中可见，可能对卫星表面敏感器件产生动力学干扰和热流影响的分布与以下三个参数有关：从推力器出口到敏感器件表面的距离；羽流场的分布角度κ，以及羽流场与敏感器件表面法线之间的夹角β。

于是，敏感器受羽流效应产生的干扰力密度或热流场密度的分布可以使用函数F(·)或H(·)表示，并且有

$\begin{array}{c}F(\·)=F(r,\mathrm{\κ},\mathrm{\β})\\ H(\·)=H(r,\mathrm{\κ},\mathrm{\β})\end{array}---\left(1\right)$

进一步地，由于

κ，β相互独立，方程可以经分离变量法进一步表示为(以下仅以H(·)为例讨论)：

H(r，κ，β)＝H₁(r)·H₂(κ)·H₃(β)               (2)

假设羽流到达自由分子状态前总通量保持不变，则由场论可知热流场与传播距离之间的关系为：

$H(r,\mathrm{\κ},\mathrm{\β})\∝\frac{1}{r^2}---\left(3\right)$

即：

H₁(r)＝H₁(r₀)·(r₀/r)²           (4)

表5羽流试验的归一化结果

表5表示对10N推力器羽流测试的归一化试验结果，其中除κ外，r与β保持不变。对以上数据进行三次样条函数插值，所得曲线如图5所示。

插值结果为数值结果，为在后续计算中降低计算难度并进一步提高计算效率，将其进行解析表达：

其中J₀(κ)为0阶贝塞尔函数，κ以“度”表示。

羽流冲击力密度及热流密度随撞击角的函数β的关系为：

$\left\{\begin{array}{c}{F}_3\left(\mathrm{\β}\right)={\cos }^2\mathrm{\β}\\ H_3\left(\mathrm{\β}\right)=\cos \mathrm{\β}\end{array}\right.---\left(6\right)$

至此，可以看出10N双组元推力器羽流场干扰力和热流量的函数分布关系可以表示为：

$\begin{array}{c}F(r,\mathrm{\κ},\mathrm{\β})=F_1\left(r_0\right)\·{\left(\frac{r_0}{r}\right)}^2\·F_2\left(\mathrm{\κ}\right)\·{\cos }^2\mathrm{\β}\\ H(r,\mathrm{\κ},\mathrm{\β})=H_1\left(r_0\right)\·{\left(\frac{r_0}{r}\right)}^2\·H_2\left(\mathrm{\κ}\right)\·\cos \mathrm{\β}\end{array}---\left(7\right)$

同理，可以推出F₂(κ)为：

上面公式中出现的字母含义，r₀为实验中所取的距离，F₁(r₀)、H₁(r₀)分别为对应r₀的干扰力密度、热流密度试验结果；r为矢量的标量值；F₂(κ)、H₂(κ)分别为推力器羽流测试的归一化试验结果干扰力密度和热流密度的拟合函数。

测试条件：H₁(r₀)＝12.02W/cm²，r₀＝550mm，按照这一数据，根据上面提及的本发明步骤模拟出的羽流热流密度场分布如图6所示，将该图所示的只有在推力器存在的时候，其羽流场分布情况与背景技术中提及的采用直接模拟蒙特卡罗(Direct Simulation Monte Carlo，DSMC)方法的结果进行比较，二者吻合情况较好，结果表明本发明方法简单、可行。

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

Claims (3)

1.一种10N推力器羽流场热效应及动力学效应确定方法，其特征在于步骤如下：

(1)首先建立卫星几何尺寸及卫星坐标系；同时分别建立10N推力器本体坐标系下羽流场位置及指向、太阳翼本体坐标系下其几何尺寸及对日定向角度，以及通信天线本体坐标系下其几何尺寸与形面法矢参数；

(2)建立推力器、太阳翼及通信天线间的坐标系转换关系，将步骤(1)中涉及的推力器、太阳翼及通信天线的参数转换到卫星坐标系下；

(3)对太阳翼及通信天线进行网格划分成单元，从第一个单元格i＝1开始执行步骤(4)～(6)；

(4)计算羽流场原点至太阳翼、通信天线表面第i单元格的矢量

(5)计算矢量

与太阳翼、通信天线表面第i单元格的法向矢量间夹角β；

(6)计算矢量与推力器轴线间夹角κ；

(7)令i＝i+1，重复执行步骤(4)～(6)，直至所有单元格计算完毕，执行步骤(8)；

(8)根据上述计算的每个单元格的矢量

β、κ，计算推力器的羽流场干扰力和热流量分布，羽流场干扰力密度分布函数F(r，κ，β)和热流密度分布函数H(r，κ，β)分别如下：

$F(r,\mathrm{\κ},\mathrm{\β})=F_1\left(r_0\right)\·{\left(\frac{r_0}{r}\right)}^2\·F_2\left(\mathrm{\κ}\right)\·{\cos }^2\mathrm{\β}$

$H(r,\mathrm{\κ},\mathrm{\β})=H_1\left(r_0\right)\·{\left(\frac{r_0}{r}\right)}^2\·H_2\left(\mathrm{\κ}\right)\·\cos \mathrm{\β}$

式中，r₀为实验中所取的距离，F₁(r₀)、H₁(r₀)分别为对应r₀的干扰力密度、热流密度试验结果；r为矢量

的标量值；F₂(κ)、H₂(κ)分别为仅与κ相关的推力器羽流测试的归一化试验结果干扰力和热流量的拟合函数。

2.根据权利要求1所述的一种10N推力器羽流场热效应及动力学效应确定方法，其特征在于：所述的仅与κ相关的推力器羽流测试的归一化试验结果干扰力密度的拟合函数F₂(κ)为：

J₀(κ)为κ的0阶贝塞尔函数，κ以“度”表示。

3.根据权利要求1所述的一种10N推力器羽流场热效应及动力学效应确定方法，其特征在于：所述的推力器羽流测试的归一化试验结果热流密度的拟合函数H₂(κ)为：

J₀(κ)为κ的0阶贝塞尔函数，κ以“度”表示。

Images