杨盛庆，贾艳胜，崔 佳，杜耀珂，王文妍，吴敬玉

(1.上海航天控制技术研究所，上海 201109； 2.上海市空间智能控制技术重点实验室，上海 201109)

0 引 言

航天器的姿轨控主要采用脉冲式输出的落压推力器。推力器贮箱通常由气垫将燃料和气体隔离，在平衡状态下贮箱压力为气垫内气体压力。当推力器工作时，由于贮箱与外界真空环境存在压差，燃料流经管路在喷嘴处产生推力，贮箱内燃料体积和气垫体积随之产生变化。根据热力学的理想气体方程，气垫内气体分子量恒定，气体压力和体积的乘积为常数。为提高燃料使用效率，推力器通常携带增压气瓶[1-2]。当增压气瓶工作时，贮箱气垫内的气体分子量增大，推力器推力发生明显变化；随着推力器贮箱压力逐步减小，推力器寿命初期、末期的推力性能存在明显差别。因此，有必要对携气瓶推力器的性能进行仿真分析。

卫星在轨自主编队通常采用开环控制模式[3]，由星载计算机(AOCC)根据控制策略生成的速度增量和推进系统工作状态来计算喷气时长，从而控制推力器阀门，执行编队控制。目前，国内外卫星编队飞行多采用自主编队保持模式，控制频率基本为每天1组多脉冲控制[3-4]。随着推力器燃料消耗和气体压力发生变化，实际的推力输出会产生相应变化。对于单次控制需要较大速度增量或1组多次喷气控制的情况，采取单点恒定推力计算喷气时长的方法，其控制精度有待提高。在传统工程应用中，开环控制形式的轨道/编队控制任务需求较少，对推力器使用的精度要求不高。目前，推力器的喷气时长优化通常依靠推力器的在轨标定，但推力器标定需要依据控制目标与实际控制结果的偏差进行数值解算，耦合了多项因素，容易造成不必要的标定系数偏差，以及推进系统不必要的燃料消耗。本文通过推力器的状态传递和推力预测，构造了以喷气时长T为变量的代数方程，根据实时的推力状态变化，有效减小喷气时长的计算误差，进而提高控制精度、减小燃料消耗。

1 推力器性能分析

1.1 推力器的推力特性

推力器工作的基本原理为热力学的理想气体方程[5-6]，其表达式为

PV/T=C

(1)

式中：P为气体压力；V为气体体积；T为气体热力学温度；C为理想气体常数。

假设贮箱初始热力学温度为Tc0，则根据贮箱初始热力学温度可计算初始贮箱内燃料密度ρ0[7-8]，即

ρ0= 1 025.817-0.874 2(Tc0-273.15)-

0.000 5(Tc0-273.15)2

(2)

贮箱内初始气垫体积为

(3)

式中：Vc为推力器贮箱容积；Mc0为燃料初始质量。

推力器工作时，根据贮箱中压力传感器测量的压力值Pc，可给出推力的二阶近似估计值[9-10]，其公式为

(4)

式中：Pc为贮箱压力。每秒气体流量

(5)

式(4)、(5)中：ai,bi(i=0,1,2)为产品工程参数，其由推进系统的热试车数据拟合得到。

1.2 气瓶工作原理

气瓶和贮箱由1个单向阀连接，单向阀开关触发的压力阈值为ΔP。当气瓶与贮箱的压力差大于ΔP时，气瓶中的加压气体进入贮箱的气垫中。贮箱气垫和气瓶气体混合平衡后，贮箱压力增大，推力产生脉冲式跳变。

气瓶中气体进入贮箱气垫后，贮箱气垫内气体分子量增大。测量贮箱压力、气瓶压力和各初始状态量，理想气体满足方程为

Vq=Tc×

(6)

式中：Pc0为压力传感器初态测量值；Vc为贮箱体积固定值；Tc为贮箱热力学温度测量值；Tc0为贮箱热力学温度初态测量值；Vq为贮箱气垫体积计算值；Vq0为贮箱气垫体积初态计算值；Pg为气瓶压力测量值；Pg0为气瓶压力初态测量值；Vg为气瓶体积固定值；Tg为气瓶热力学温度测量值；Tg0为气瓶热力学温度初态测量值。

1.3 剩余燃料计算

根据测量的贮箱热力学温度Tc，计算贮箱内剩余燃料量[11]为

Mc=ρ(Vc-Vq)

(7)

式中：Vq为贮箱气垫体积；ρ为贮箱内燃料密度，ρ=1 025.817 -0.874 2(Tc-273.15)-0.000 5(Tc-273.15)2。卫星当前质量

Ms=M0+Mc

(8)

式中：M0为卫星净质量；Mc为剩余燃料量。

1.4 寿命期间的性能分析

对1 N携气瓶推力器的性能进行仿真分析，推力器性能如图1所示。图中可见：推力在寿命初期和末期变化明显。图1(a)中，气瓶工作时推力器的推力存在有限增幅的脉冲式跳变，并随着贮箱压力的减小直至下一次气瓶工作。图1(b)中，在第1个10 kg燃料消耗区间内，气瓶工作次数达15次，随后迅速减小，在第3个10 kg燃料消耗区间内工作次数减小到4次。

图1 推力器寿命期间性能Fig.1 Capability of thruster during operating life

2 推力器喷气时长计算方法

考虑到AOCC计算能力有限，推力器动力学模型中使用的速度增量关机方式不再适用于喷气时长计算。由于传统的基于单点测量的推力器喷气时长计算方法不能有效体现推力变化和燃料消耗，设计了一种基于推力预测的喷气时长计算方法。该方法通过估计单位时间内管路的过肼量，给出当前状态下推力的变化率，构造以喷气时长T为未知变量的代数方程；针对一组多脉冲控制的情况，设计了基于状态传递和推力预测的喷气时长计算方法，给出了各脉冲喷气时长的计算公式。

2.1 推力的估计公式

当推力器工作时，可根据贮箱压力变化，估计推力值。由式(5)可得每秒气垫体积变化，即

(9)

每秒气垫体积变化导致的喷气1 s后贮箱压力值为

(10)

(11)

2.2 一次喷气的时长计算

(12)

式中：θ为安装角度；N为推力器工作时的喷嘴数量。当多个推力器同时工作时，不同喷嘴的过肼量一致性(推力一致性)由管路设计保证，并已通过推力器热试车验证。因此，构造关于喷气时长T的代数方程，其表达式为

(13)

2.3 一组多次的喷气时长计算

当卫星编队在轨运行时，通常需要生成1组多脉冲的编队保持控制策略。以1组3脉冲喷气为例，给出了其喷气时长的计算公式，见表1。对于1组更多脉冲数的喷气控制，同样可以依据状态传递原理给出喷气时长计算公式。具体的计算步骤为：组内的第N+1次喷气，需要根据计算得到的组内前N次的喷气时长，估算当前贮箱中气垫体积，并根据理想气体方程得到当前的气体压力，进而得到推力器当前的推力和推力变化率。

表1 1组3脉冲喷气的喷气时长计算公式

3 推力器喷气时长计算方法的仿真验证

3.1 喷气时长的动力学理论值

设计数值试验对推力器喷气时长计算方法的有效性进行仿真验证。其中，喷气时长的动力学理论值计算，依据速度增量关机的方式搭建仿真模型[12]。每个仿真步长dt内，需要计算推力f(t)、卫星总质量m(t)，速度增量的积分形式为

(14)

该喷气仿真模型中，当卫星积分得到的速度增量达到要求的速度增量时，喷气停止。仿真步长dt足够小时，由该喷气仿真模型计算得到的喷气时长Tx为喷气时长的动力学理论值。

3.2 喷气时长计算方法对照试验

构建以喷气启控时刻、速度增量为仿真输入的精度分析模型，模型如图2所示。

图2 喷气时长计算方法的精度分析模型Fig.2 Analysis model for different fire-time algorithms

图3 不同速度增量的喷气时长计算方法比对Fig.3 Comparison of different fire-time algorithms for multiple delta-v

基于精度分析模型，本文采取本文算法和传统算法对以下2种工况进行比对试验。

1) 不同速度增量的单次喷气对比。等间隔覆盖图3中的速度增量，不同的速度增量单次调用精度分析模型，即可得到喷气时长的动力学理论值、原有算法和本文算法计算的喷气时长值。不同速度增量的单次喷气对比结果如图3所示。由图可知，本文方法计算的喷气时长与理论喷气时长的偏差小于传统的喷气时长计算方法，且速度增量越大，本文方法的优越性越明显。

2) 1组3脉冲喷气的对比试验。仿真输入1组3个喷气启控时刻和3个速度增量。1组3脉冲喷气的比对结果见表2。由表可见，本文方法计算的各脉冲喷气时长与理论喷气时长更接近，优于传统的喷气时长计算方法。

表2 1组3脉冲的推力器喷气时长计算方法的比对结果

4 结论

卫星自主编队保持属于一类开环控制，高精度的控制策略计算对推力器的精确使用尤为重要。本文针对携气瓶推力器进行了性能分析，介绍了推力器寿命期间气瓶工作、推力变化情况，提出了一种针对时变推力的AOCC喷气时长计算方法，通过推力预测，构造了关于待求喷气时长的代数方程，适用于单次喷气、一组多脉冲喷气等情况。本文方法能有效减小喷气时长的计算误差，提高卫星姿轨控系统的控制精度。