1 引言
肼类单组元推进剂有良好的处理特性、易于贮存、分解产物清洁，故肼类单组元动力系统广泛应用于各类空间飞行器的姿态控制、轨道转移以及末速修正等任务中，其性能以及稳定性直接影响着航天器的控制精度、寿命和可靠性[1-3]。集合了多台发动机的动力系统在大部分工作时间都处在脉冲工作状态，系统内多台发动机频繁地起动与关机，发动机之间的相互影响不可避免。这种影响会导致响应延迟，进而影响飞行器控制[4]，因此需要对动力系统的动态过程开展仿真预示研究。国外Moore[5]对某推进系统的气路增压模块开展了试验与仿真研究工作，分析了初始气瓶压力、贮箱压力、贮箱体积、等因素对增压时间的影响，但研究中未考虑液路系统及推力室组件。国内孙得川[6]采用CFD方法对单台HAN基单组元发动机的起动过程和稳态工作过程进行了研究，其他学者[7-8]对双组元姿控动力系统进行了动态仿真平台的搭建与研究，但缺少单组元发动机的系统仿真模块；并且动力系统动态过程的研究以相同推力量级发动机之间的相互影响为主[9]，未见研究不同推力量级发动机动态过程的相互影响。
因此，本文在借鉴已有研究成果的基础上，基于MWorks软件平台开发含有不同推力量级单组元发动机的动力系统仿真模型，计算分析不同推力量级发动机动态过程之间的相互影响规律，旨在为动力系统设计和飞行器姿态控制策略提供参考。
2 典型组件的动力学模型
单组元姿控动力系统是由气路增压系统、液路推进剂供应系统以及推力室组成。对于恒压挤压式姿控动力系统，因贮箱增压压力基本保持不变[8]，可忽略气路增压系统特性对姿控动力系统发动机的起动、关机等工作过程影响，故仅建立液路推进剂管路、阀门及推力室等典型组件的动力学模型。
2.1 管路模型
在建立液体管路模型时不考虑管道变形，认为管道为刚性；管道横截面上的速度分布是均匀的，即液体流动为一维。根据所研究的动力学范围，采用分段集中参数模型。根据经验，采用的分段原则是将分段管长取为波长的4%以内，即：Lmax≤0.04a/fmax。
对于一段长度为L、截面积A的管道，考虑惯性、黏性时的流动方程为：
(1)
式中，qm为管路内推进剂流量；p1、p2分别为管路入口和出口压力； ρ为推进剂密度； ξ为总的流阻系数，由管路的沿程损失和局部损失组成。
考虑液体压缩性的方程为：
(2)
其中z=V/a2表示管路的流容系数，表征流体的压缩性。
2.2 电磁阀模型
在液体姿控动力系统中一般采用常闭式电磁阀，主要用于实现推进剂流路的打开和关闭，因此电磁阀的动态特性直接影响动力系统的动态过程。
对于系统仿真中所关注的电磁阀响应时间和阀门开度进行建模，
(3)
Γ是电磁阀开度与流量系数的乘积(即Γ=τCv)。它主要由实验来确定，将实验结果用列表法或曲线拟合给出Γ与时间的对应关系。当缺乏实验数据时，通常采用下述方法来确定。
当阀门打开时：
(4)
当阀门关闭时：
(5)
式中，te，k和te，g为电磁阀阀芯开启与关闭的动作时间；tc，k和tc，g为电磁阀开启与关闭的响应时间；m为试验曲线拟合指数，主要由Γ-t的曲线形状而定，一般情况下m=1～3[10]。
2.3 催化床模型
考虑到催化床内复杂的物理化学过程，引入燃烧室时滞模型，假设推进剂在一定转化时间之后分解为燃气(此时间包括了推进剂蒸发、流动、分解时间的总和)。假设在该转化时间内推进剂的转化速率均匀，并且分解后的气体在任何瞬时也是均匀分布的，即氨的解离过程与DT-3的分解过程同时进行，则液体推进剂及气体分解产物在催化床中的质量变化分别为：
(6)
(7)
式中，
为催化床出口的气体质量流量，
为喷管的喉部截面积；ml、mg分别为催化床内积存的液体推进剂质量和气体质量；
分别为流入催化床内的液体推进剂质量流量和流出推力室的气体质量流量；τ为推进剂的转化时间，根据推进剂的分解延迟期试验得到约为10 ms。
由于催化床内装载着紧密的催化剂，因此推进剂流过催化床时需要克服阻力(床流阻)。假设催化床的流阻都通过一个压降面来实现，故催化床内的压力平衡方程为：
pv=pb+pc
(8)
式中，pv、pb、pc分别为催化床入口压力、床流阻和推力室室压。催化剂的装填工艺对于床流阻有较大的影响，考虑到催化床内复杂的物理化学过程，床流阻常采用文献中[1]的经验公式计算。
假设催化床内气体服从完全气体状态方程，则催化床内压力的变化为[7]：
(9)
在液体推进剂输送系统中，由于多个发动机共用一套推进剂输送系统，其中存在多个推进剂管路的分支；通常忽略分支处的压力损失，而采用分支点处各支路压力相等、流量平衡来计算。
3 仿真模型建立与校验
以上述单组元发动机典型组件的数学模型为基础，将各组件定义为具有一定输入输出关系的模块，利用Modelica语言[]的非因果建模和面向对象功能，在多领域统一建模与仿真分析软件平台MWorks[12](MWorks由华中科技大学开发，支持Modelica语言，适用于复杂工程系统的建模、仿真、分析与优化。)上编制开发了发动机各典型组件的参数化、通用化的仿真模块。
为了验证模型的正确性，选取如图1所示的某单组元姿控动力系统为对象，以某次热试车试验结果作为验证依据。
图1 某单组元动力系统原理示意图
Fig.1 Schematic of the monopropellant propulsion system
该动力系统中配置了多台100 N量级的轨控发动机与10 N量级的姿控发动机，其中A类为轨控发动机，B类为姿控发动机。热试车过程中，为了考验发动机之间的工作协调性，在多台发动机处于工作状态时，向处于关机状态的A1、B1发动机发出开机指令，同时向处于工作状态的A2、B2发动机发出关机指令。仿真采用定步长的四阶Runge-Kutta法，计算时间步长取10-6 s。
图2为A1发动机起动、A2发动机关机的计算结果和试车结果对比图。
图3为B1发动机起动、B2发动机关机的计算结果和试车结果对比图。
图2 轨控发动机动态过程计算值与测量值曲线
Fig.2 Comparison between the calculated and test results of the transient process of the orbit control engines
图3 姿控发动机动态过程计算值与测量值曲线
Fig.3 Comparison between the calculated and test results of the transient process of the attitude control engines
从图2可以看到，与A1单机起动过程相比，多机动作时A1发动机起动过程的压力波动幅度更大；2种状态下A2发动机的关机过程差异则更为显著，多机动作时室压曲线存在一个明显的凹坑，而A2单机关机时则不存在此现象。这是由于在两台较大推力的轨控发动机同时动作段，A1起动和A2关机引起供应系统内供应压力变化相互影响所致。而从图3中发现，在两台较小推力姿控发动机单独动作和同时动作状态下，B1起动和B2关机的室压变化都较为一致，无明显差异。因此，在单组元液体姿控动力系统中，较大推力的轨控发动机动作会影响整个供应系统，进而对同时动作的其他发动机产生影响。
另外，在多机同步动作状态下的计算结果与试验结果吻合度较好，证实了所建立的仿真模型较为准确。
4 多机耦连的影响分析
尽管上节的结果表明，在单组元姿控动力系统中应尽量避免较大推力量级发动机动作的同时其余位置发动机动作，但这在飞行器需要姿控调整和轨道保持时难以避免。并且姿控动力系统中不同位置发动机的阀门实际动作时间存在一定的散差，不同发动机的动作难以达到完全同步，会存在一定的时间差。因此，下文将计算分析连续工作和脉冲工作2种状态下，2种推力发动机动作存在的时间差对姿控发动机动态过程的影响。考虑到仿真计算的实际设置问题，此时间差可以通过调整对不同位置发动机发出控制指令的时间差(即指令时差)来实现。
4.1 指令时差对连续工作的影响
图4给出了A1和B发动机在不同指令时差条件下，单台A1发动机起动分别对B1发动机起动和B2发动机关机室压曲线的影响。
可以看到，当A1发动机起动提前于B1和B2发动机动作时(即图4中Δt=-20 ms)，B1发动机的起动过程不受A发动机动作所影响，而B2发动机的关机过程将受到较大的影响，会在关机前产生明显的压力波动。而在A1发动机滞后于B1和B2发动机动作时(即图4中Δt=20 ms)，B1发动机的起动室压波动幅值和频率有明显的增加，B2发动机的关机过程则不会受到A1发动机的滞后动作影响。在A1发动机和B1、B2发动机同时动作(即图4中Δt=0 ms)时，B1发动机的起动过程受到的影响较小，而B2发动机的室压曲线出现了明显波动。
产生这种情况的原因是，A1发动机的提前起动会使动力系统产生较大压力波动，使得B2发动机关机产生较为明显的波动；而在A1发动机滞后起动时，B2发动机的关机过程已结束，所以没有受到明显影响。B1发动机的起动受到A类发动机滞后起动的影响，这是由于A1发动机滞后起动时，B1发动机的室压已上升接近稳态值，A1发动机此时起动使得系统压力产生一定波动，因此B1发动机受到了较为明显的影响。
图4 A1起动与B类发动机动作指令时差的影响曲线
Fig.4 Effect of command time lag between engine-A1 starting and engine-B action
图5给出了A1和B发动机在不同指令时差条件下，单台A1发动机关机分别对B1发动机起动和B2发动机关机室压曲线的影响。
图5 A1发动机关机与B类发动机动作指令时差的影响曲线
Fig.5 Effect of command time lag between engine-A1
shutdown and engine-B action
可以看到，当A1发动机关机提前于B1和B2发动机动作时(即图5中Δt=-20 ms)，B1发动机的起动过程不受A1发动机动作影响，而B2发动机的关机过程将受到较大的影响，在关机前产生较大的波动。而当A1发动机滞后于B1和B2发动机动作时(即图5中Δt=20 ms)，B1发动机的起动室压波动幅值增大，B2发动机的关机过程则不受A1发动机的滞后影响。
上述结果产生的原因是，由于A1发动机提前关机会使动力系统产生较大压力波动，使得B2发动机关机时产生较为明显的波动，而在A1发动机滞后关机时，B2发动机的关机过程已结束，所以没有受到明显影响。B1发动机的起动受到A1发动机滞后关机的影响，这是由于A1发动机滞后起动时，关机引起的压力波动较晚传播到B1发动机前，故B1发动机的室压上升至接近稳态值时有一定的波动增大现象。
4.2 指令时差对脉冲工作的影响
图6为B1姿控发动机在单独工作和轨控A1发动机脉冲工作影响下的一段脉冲工作过程室压曲线，其中实线表示单独工作，虚线表示存在轨控发动机脉冲工作的情况。
图6 不同状态下B1发动机脉冲工作过程曲线
Fig.6 Pulse working of engine-B1 under different states
可以看到，前几个脉冲的关机后效冲量较大，尤其是第1个脉冲。而在A1发动机影响的情况下，各脉冲的关机压力拖尾更小一些。这是受到阀门关闭特性的影响，当A1发动机工作时，供应系统的压降更高，使小发动机的入口压力有所降低，阀门的关阀响应更快。
发动机的脉冲工作过程是以单个脉冲循环中的冲量为主要特性参数的。图7表示了不同指令时序差情况下，B1发动机前5次脉冲工作的冲量变化。
从图7观察到，随着指令时差的增加，第1个脉冲冲量出现了明显的降低，但在第2个脉冲又迅速增加。脉冲冲量随连续脉冲数增加而逐渐减小，当轨控和姿控发动机脉冲工作指令同时发出时，姿控(B1)发动机的脉冲冲量整体较小，这与稳态工作的影响结果相似。较大推力的轨控发动机脉冲工作指令先于较小推力的姿控发动机脉冲工作指令时，对动力系统耦合工作的影响较小。
4.3 工作占空对脉冲工作的影响
工作占空np定义为发动机脉冲工作中开机时间与关机时间之比。对于卫星飞行器轨道调整而言，飞行器所需脉冲总冲量为一固定值。因此，脉冲工作中的不同工作占空由脉冲过程的关机时间长短决定，而开机时间相同。
图7 不同指令时差下B1发动机脉冲工作冲量曲线
Fig.7 Pulse impulse of engine-B1 under different
command time lag
在A1发动机脉冲工作情况下，工作占空分别为1和0.5的B1发动机脉冲工作过程如图8所示。实线表示工作占空为1，虚线表示工作占空为0.5。
图8 不同工作占空下姿控发动机工作过程曲线
Fig.8 Engine pulse working under different duty cycle
可以看到，不同工作占空时，姿控发动机工作段内的压力上升过程和平台压力近似相同。但当工作占空为0.5时，由于阀门关闭时间的增长，姿控发动机脉冲关机压力拖尾明显减小。这会对整个姿控发动机脉冲工作产生明显的影响。
在飞行器和卫星调整姿态时，需要动力系统中各类发动机都提供相对稳定的总冲量，一般为额定冲量的±10%。为此计算了不同工作占空下，姿控B1发动机的脉冲总冲量，结果如图9所示。图中无量纲冲量为发动机十次脉冲的平均无量纲冲量。
可以看到，随着工作占空的逐渐减小，即在一个脉冲中关机时间逐渐增加，冲量逐渐增加，当工作占空为0.25时，冲量已经达到稳定状态。
在所计算的算例中，当工作占空大于0.5时，冲量小于额定冲量的90%，该状态下无法满足总体的稳定工作冲量要求。这表明在轨姿控动力系统脉冲工作重合的情况下，姿控发动机会存在一个临界稳定工作占空，当工作占空大于此值时，脉冲总冲量会减小至相对稳定冲量以下。因此，在轨姿控动力系统的设计和飞行器总体姿态动力学分析中需要考虑该因素的影响。
图9 工作占空对脉冲总冲量的影响曲线
Fig.9 Effect of duty cycle on total pulse impulse
5 结论
1) 轨控发动机的瞬态动作会对姿控发动机动态过程产生明显影响，在姿控发动机起动、关机等动态过程中应尽量避免轨控发动机瞬态动作。
2) 在轨控发动机与姿控发动机脉冲工作重合时间内，姿控发动机的工作占空如果大于临界值，姿控发动机将不能提供所需的额定冲量。因此姿控发动机工作占空应小于临界值，才能保证轨姿控发动机脉冲重合工作段的脉冲冲量。
[1] 周汉申.单组元液体火箭发动机设计与研究[M].北京:中国宇航出版社，2009.
[2] 潘海林，丁风林，李永,等.空间推进[M].西安:西北工业大学出版社，2016.
[3] 肖明杰.卫星姿态控制和反作用控制用单组元推进剂供应系统的落压过程和水击过程[J].火箭推进，1998(04):52-59.
Xiao M J.Blowdown process and water-hammer process of monopropellant feed system for satellite attitude control and reaction control[J].Journal of Rocket Propulsion，1998(04):52-59.
[4] 臧月进，李仁俊，周莎.动能拦截器姿控地面验证试验误差建模与仿真[J].兵器装备工程学报，2018，39 (04):151-154.
Zang Y J，Li R J，Zhou L S.Error modeling and simulation of kill vechicle attitude control system ground test[J].Journal of Ordnance Equipment Engineering,2018,39(04):151-154．
[5] Moore J D.Activation time analysis in a propulsion system[J].AIAA 2017-4763.
[6] 孙得川，金东洙，于泽游.硝酸羟胺基单组元发动机起动过程数值模拟[J].兵器装备工程学报，2018，39(05):5-10.
Sun D C，Kim D S，Yu Z Y.Numerical simulation of start-up process for HAN-based monopropellant rocket engine[J].Journal of Ordnance Equipment Engineering，2018，39(05):5-10．
[7] 林西强，程谋森，刘昆,等.液路耦连多发动机系统开关机动态特性研究[J].推进技术，1999，20(04):22-25.
Lin X Q，Cheng M S，Liu K，et al.Dynamic Characteristics during startup and shutdown operation of liquid feed line coupled multi-engine system[J].Journal of Propulsion Technology，1999，20(04):22-25.
[8] 张黎辉，李家文，张雪梅,等.航天器推进系统发动机动态特性研究[J].航空动力学报，2004，19(04):546-549.
Zhang L H，Li J W，Zhang X M，et al.Dynamic characteristics study of spacecraft propulsion system engine[J].Journal of Aerospace Power，2004，19(04):546-549.
[9] 王申.航天器推进系统动态特性数值仿真与分析[D].长沙:国防科技大学，2007.
[10] 陈新华，田希晖，苏凌宇,等.航天器推进理论[M].北京:国防工业出版社，2013.
[11] Tiller M M.Introduction to physical modeling with modelica[M].Holland:Kluwer Academic Publishers，2001.
[12] 陈立平，周凡利，丁建完,等.多领域物理统一建模语言MODELICA与MWORKS系统建模[M].武汉:华中科技大学出版社，2017.