无拖曳卫星氮气微推进系统填充与开机工作特性的仿真研究

汪旭东，李国岫，陈 君，李洪萌，虞育松

(1. 北京交通大学机械与电子控制工程学院，北京 100044；2. 北京控制工程研究所，北京 100190；3. 北京市高效能及绿色宇航推进工程技术研究中心，北京 100190)

0 引 言

无拖曳卫星在引力波测量[1]、地球重力场探测[2]、卫星导航、空间谱分辨率探测[3]等科研与军事领域的研究中具有非常重要的应用价值。如欧空局在2009年发射的GOCE卫星主要用于地球重力场测量[4-5]，可对地球重力场和海洋流动进行高精度测量。该卫星安装有2台离子推进系统和1套冷气推进系统。另外，LISA Pathfinder卫星则主要用于空间基础科学研究，进行空间引力波探测[6]。

无拖曳卫星在工作时需要使用微推进系统来实施连续、精确的主动控制，以抵消卫星在低轨道飞行时受到的大气阻力和多种环境干扰力(如地球磁场、太阳辐射导致的光压等)，从而使卫星处于完全的自由落体状态[7-9]。为实现卫星的无拖曳飞行状态，卫星需要将残余扰动加速度控制在10-9m/s2范围内[10]，从而实现卫星的超稳定控制[11]。这对卫星控制及推进系统提出了非常高的要求。

氮气微推进系统是冷气推进系统大类中的一种。这类型推进系统以氮气作为工作介质，具有推力控制精度高、系统干质量小、结构简单、耗能低、无污染、性能稳定等优点[12]。特别适用于无拖曳卫星微小推力控制和需求。如NASA制定的无拖曳CubeSat计划，其卫星中采用了VACCO公司研发的氮气微推进系统MiPS[13]。

氮气微推进系统的主要部件一般包括气瓶、压力调节阀、推力器组件、传感器等[14]。推力由高压氮气流经喷管产生。国外对氮气微推进系统的研究较多。2017年Ranjan等[15]采用试验与仿真相结合的方法对氮气推进系统的推力性能进行了分析。2018年Ranjan等[16]对CubeSats上的氮气微推进系统进行了试验研究，分析了真空与大气条件下推力器工作时的推力、比冲结果。另外Samuel等[17]也开展了有关氮气微推进系统方面的分析。

国内对氮气微推进系统的试验与仿真研究相当匮乏。2010年周伟勇等[18]提出了一种面向目标任务的小卫星冷气推进系统设计方法，并设计了一套20 mN氮气推进系统。2014年北京交通大学的王璐等[19]简单开展了微小超高压冷气推进系统动态工作性能的仿真研究。但上述系统的推力较大(200 mN)，且推力阀组件为电磁驱动，响应较慢。

本文以应用在某型无拖曳卫星上的氮气微推进系统为研究对象。该推进系统的推力组件为压电方式驱动，且推力较小(100 mN级)。目前国内尚无针对压电驱动方式、推力为100 mN级的氮气微推进系统工作性能方面研究的相关报道。本文的研究工作如下：采用AMESim仿真软件建立上述氮气微推进系统的模型，并对氮气微推进系统的氮气填充和开机过程的系统工作性能进行深入分析，从而为无拖曳卫星氮气微推进系统的研发提供必要依据。

1 数学模型

1.1 流动控制方程

1) 氮气瓶状态热力学模型

高压氮气瓶存贮在两个高压气瓶中。假定气瓶为均匀压力和温度分布的单开口绝热热力系统。

气瓶内氮气的质量守恒方程：

(1)

气瓶内氮气的能量守恒方程：

(2)

式中：V为气瓶容积；ρ，p和T分别为氮气密度、压力和温度；h为气体焓值，cv为氮气定容比热。

2) 气体管路模型

冷气推进系统中包含总长约为0.5 m的气体管路，连接了贮箱、自锁阀、减压阀和整合了喷管的压电比例阀。AMESim将管路内气体流动近似为一维瞬态流动。管路内的气体流动控制方程包含质量、动量和能量守恒方程[20]。

(3)

(4)

(5)

式中：u为x方向速度分量，τ为x方向的应力张量分量，p为气体压力，e为比内能，Sτ和St分别为气体黏性力做功、管路与环境的换热。

3) 气体状态方程

由于氮气在气瓶内的存贮压力达到1.5×107Pa，为了准确求解气体热力状态。本次研究采用RK(Redlich-Kwong)方程[21]描述气体压力、密度与温度之间的关系。

1.2 喷管推力模型

喷管推力计算采用Laval喷管理论解析解进行近似计算。假定整个流动过程为一维等熵流动状态。喷管出口马赫数Me与扩展比ε存在以下关系：

(6)

式中：ε=A*/Ae，喷管喉口面积为A*，喷管出口面积为A。

喷管喷射质量流量：

(7)

喷管出射平均速度：

(8)

推力计算：

(9)

其中，γ为工质比热比，R为氮气气体常数，pe和pt分别为喷管出口与入口的气体压力，Tt为喷管入口气体温度，下标e表示喷管出口位置。

喷管流通状态由针阀升程控制(见图1)。为了考虑针阀运动对喷管流通与推力的影响，建模时将喷管喉口有效流通面积A*表示为针阀升程S和针阀针尖锥角的函数。

图1 喷管喉口有效流通截面Fig.1 Schematic diagram of the effective flow area ofthe nozzle

喉口入口端与针阀锥面的距离L：

L=Ssinα

(10)

喷管喉口有效流通面积A*：

(11)

式中：r为喷管喉口半径，α为针阀针尖半锥角。

2 氮气微推进系统建模

本次研究的氮气微推进系统主要包含：2个高压氮气瓶、防爆阀、自锁阀、减压阀、压电比例阀、压力和温度传感器等组成。在实际推进系统中会涉及多个相同压电比例阀组成的推力器组(以满足实际工作任务需求)，但本次研究只分析单个压电比例阀的工作特性。由于防爆阀在正常工作时不工作，而只在异常情况下启动，因此本文在推进系统建模时不予考虑(见图2)。

图2 氮气微推进系统组成Fig.2 Composition of nitrogen micro propulsion system

氮气微推进系统的工作过程分为3个阶段。1)关机阶段。自锁阀和压电比例阀处于关闭，减压阀处于开启状态。2)填充阶段。自锁阀开启，气瓶内的高压氮气通过自锁阀和减压阀，系统管路和阀门组件内逐渐被氮气填充，压力逐渐上升至目标压力。减压阀至压电比例阀之间的氮气压力由减压阀进行控制，以保证喷管上游压力腔内的压力处于稳定状态，实现推力的精确控制。3)推力器工作阶段。压电比例阀开启，高压气体喷出喷管形成推力。

减压阀是调节下游压电比例阀入口压力的关键部件。减压阀的工作原理：当反馈腔内压力低于目标压力值时，由于弹簧力大于氮气压力，球阀运动并打开阀孔。上游高压氮气会进入下游管路和压电比例阀腔。当反馈腔内的氮气压力达到目标值时，弹簧会推动连杆运动，控制球阀芯关闭阀孔，截断上游来流，实现维持压力的作用。

图3 压电比例阀示意图Fig.3 Schematic diagram of cross section of thepiezoelectric valve

图4 采用AMESim软件建立的氮气微推进系统模型Fig.4 AMESim models for the nitrogen gas propulsion system

压电比例阀是氮气微推进系统的核心部件，主要包括压电驱动电路、电压阀芯、阀腔及其氮气流道、针阀和喷管等几个部分(见图3)。压电比例阀的工作性能决定了系统的推力响应速率和推力控制精度。压电比例阀腔及其氮气流道的AMESim模型搭建考虑了阀芯下侧的4个节流孔、阀芯内部流道和上下侧的容积腔(见图4)。

喷管推力器元件的构建采用AMESET模块根据第1.2节的喷管推力模型进行构建。喷管模型在创建时采用了AMESim中的Supercomponent功能进行封装。推进系统的管路总长为0.5 m。氮气微推进系统整体的AMESim模型见图4所示。主要系统参数如表1所示。计算初始时刻，氮气瓶内压力为1.5×107Pa，而系统所有管路、阀内腔和环境压力均为1.0×10-5Pa。计算时间步长为1.0×10-5s。压电比例阀的驱动电压为80 V。

3 结果与讨论

本次研究分别针对氮气填充、开机两个阶段的工作过程。氮气填充过程主要分析氮气瓶和减压阀的工作特性。开机过程主要分析压电比例阀的工作特性，探讨压电比例阀响应特性、驱动电压对该组件响应与推力性能的影响。

表1 推进系统主要参数Table 1 Main parameters of the propulsion system

3.1 填充过程中的气瓶与减压阀工作特性

填充过程为自锁阀打开后，气瓶中的高压氮气经由减压阀填充整个系统管路并达到准稳定状态的过程。减压阀的关闭动作由反馈腔内的氮气压力控制。减压阀将反馈腔内的氮气压力控制在目标值内。图5给出了其中一个氮气瓶在填充期间的压力和质量流量变化情况。两个氮气瓶对称布置，其工作状态相同，因此本文只给出其中一个氮气瓶的状态过程。计算结果显示，填充过程中气瓶的压力出现微弱下降，从初始的1.5×107Pa降低到0.3 s时的1.49947×107Pa。气瓶氮气质量流量出现梯形变化过程。填充开始瞬间气瓶质量流量迅速降低(由于管路内初始压力为1.0×10-5Pa。气瓶与管路的大压差导致了显著的初始质量流量)，随后维持约1.5 s的高质量流量稳定阶段，当减压阀腔内压力达到稳定压力后，该质量流量逐渐减小，在填充开始后0.3 s后接近0 g/s，并存在微弱的波动。该波动是由减压阀球阀芯的开关动作引起的。

图5 氮气瓶的压力和质量流量Fig.5 Pressure and mass flow rate of the nitrogen tank

图6给出了氮气填充过程中的减压阀入口及反馈腔的压力变化过程。减压阀反馈腔内的压力在填充开始后迅速增加并到达稳定压力值。计算得到的实际准稳定压力值为1.98×105Pa，与控制的目标值偏差为1%，满足2×105Pa的设计要求。模型设定的参数，如阀孔实际形状及其流量系数(当前设置为0.72)可能引起仿真偏差。减压阀入口端管路的氮气压力稳定时间相对较长(约0.3 s)。高压氮气瓶中的气体压力较高，导致减压阀至氮气瓶之间管路的填充稳定时间较长。

图6 减压阀出入口端和反馈腔内的氮气压力Fig.6 Nitrogen pressure in exit, entrance and feedbackcavity of the pressure reducing valve

图7 减压阀反馈腔的氮气压力变动率Fig.7 Nitrogen pressure fluctuations in the feedback cavityof the pressure reducing valve

为保持反馈腔压力稳定，减压阀内球芯会在气压与弹簧力作用下开关阀门(见图8，纵坐标为log10显示)，控制流通特性。球阀芯的高频运动造成减压阀反馈腔压力的波动。计算结果显示，整个填充过程，减压阀的阀芯并未完全关闭。填充开始后约0.15 s，阀芯就开始出现快速振荡，振幅随着时间逐渐减小(见图8)。

图8 减压阀阀芯的瞬态位移Fig.8 Transient displacement of the relief valve core

图9给出了减压阀出口的瞬态质量流量结果。自锁阀开通后，高压氮气迅速进入减压阀，导致减压阀出口质量流量迅速增加，在0.02 s附近出现瞬时峰值。随后质量流量迅速降低。当阀芯接近关闭后，减压阀出口质量流量减至0附近，并出现小幅高频振荡。

图9 减压阀出口的质量流量Fig.9 Mass flow rate in the exit of the pressurereducing valve

3.2 开机过程的压电比例阀工作特性

氮气在系统管路内填充完毕后(约0.15 s)，压电比例阀即可进入工作状态。本次研究假定在0.5 s时压电比例阀开启(氮气填充过程在0.3 s时完成)，推进系统进入工作状态。从开机到稳定过程总时为2.5 s。图10给出了压电比例阀内腔内(节流孔上游腔体)的氮气压力。填充完成后，压电比例阀内腔体的氮气压力已经处于稳定状态。一方面压电比例阀内部腔体容积设计较大，具有稳压的效果；另一方面，针阀开度较小，对上游腔体氮气的消耗很慢。上述两方面导致压电比例阀内腔氮气压力在开机过程中非常稳定。

图10 压电比例阀内腔体的氮气压力Fig.10 Nitrogen pressure in the cavity of thepiezoelectric valve

压电阀芯的运动特性是决定氮气微推进系统推力和响应性能的主要因素。图11给出了开机过程中的压电比例阀的阀芯位移。结果显示，阀芯在驱动电路通电后迅速开启。从开机时刻到阀芯开启至稳定位移的响应时间约为0.64 ms。阀芯开启的稳定位移为3.67 μm。试验测量当前驱动电压(80 V)的阀芯位移为3.85 μm，响应为0.69 ms。阀芯位移和响应的仿真预测误差分别为4.90%和7.81%。说明当前建立的压电比例阀模型具有良好的预测精度，可准确计算阀芯运动。

图11 压电比例阀阀芯位移Fig.11 Displacement of the piezoelectric valve core

图12给出了开机过程中，压电比例阀的氮气质量流量(以喷管喉口截面统计)结果。开机后7.6 ms，氮气通过喉口的质量流量达到稳定值6.79×10-4g/s。

图12 喷管喉口截面的氮气质量流量Fig.12 Nitrogen mass flow rate passing throughthe nozzle throat

图13为开机过程的喷管推力结果。喷管推力的建立过程略滞后于针阀开启过程。针阀开启后，氮气从阀体内腔流经喷管并建立推力需要一定的迟滞时间。开机后约8 ms时推力达到稳定值0.588 mN。

图13 开机过程的喷管推力Fig.13 Thrust force of the nozzle during starting processes

3.3 驱动电压对压电阀工作特性的影响

提高驱动电压可增加陶瓷压电片的形变量，增加压电阀针阀位移，对喷管推力具有显著影响。本次计算了驱动电压为60 V、80 V、100 V和120 V条件下压电比例阀的工作性能。图14为不同驱动电压下压电阀内针阀的位移结果。开机初期，针阀位移均出现振荡并逐渐趋于稳定的过程。随着驱动电压的增加，针阀稳态位移逐渐增加，但响应时间差异较小。当驱动电压分别为60 V、80 V、100 V和120 V时，针阀稳态位移分别为3.17 μm、3.67 μm、4.10 μm和4.49 μm。

图14 不同驱动电压下的针阀位移Fig.14 Needle displacements at the different driving voltages

图15为驱动电压为60 V、80 V、100 V和120 V条件下喷管推力的变化过程。由于推力在开机后约10 ms便进入稳定状态，因此图中仅给出计算时间为0.5 s(开机时刻)至0.7 s之间的推力值。结果显示，随着驱动电压的增加，喷管推力显著增加。增大驱动电压引起针阀开度增加，增加了喷管流通面积，增加了氮气喷射质量流量，从而增加了推力。当驱动电压分别为60 V、80 V、100 V和120 V时，喷管稳态推力分别为0.404 mN、0.588 mN、0.807 mN和1.021 mN。

图15 不同驱动电压下的开机过程喷管推力Fig.15 Thrust force at the different driving voltagesduring starting processes

图16对不同驱动电压下的针阀位移和喷管推力特性进行了统计对比。经过线性拟合，驱动电压U(V)与喷管稳态推力F(mN)之间的关系式为：

F=-0.225+0.0103U

(12)

拟合式(12)的R2为0.998，说明驱动电压与喷管稳态推力之间存在良好的线性关系。

图16给出了不同驱动电压下的针阀稳态位移、推力、针阀响应时间和推力响应时间结果。当驱动电压增加时，针阀开启和喷管推力响应时间均有所增大。驱动电压从60 V增加至120 V时，针阀开启的响应时间从0.62 ms增加至0.79 ms(增幅为27.4%)，喷管推力响应时间从7.1 ms增加至10.3 ms(增幅为45.1%)。喷管推力响应时间的增幅大于针阀开启响应时间的增幅。计算结果说明，当前系统的喷管推力可通过改变驱动电压进行mN级的线性控制。

图16 驱动电压对针阀位移、喷管推力及其响应时间的影响Fig.16 Effects of driving voltages on the needle displacement,thrust force and corresponding response time

4 结 论

本文采用数值模拟方法对无拖曳卫星上氮气微推进系统进行建模与仿真研究，主要结论如下：

1)在填充开始后0.09 s，减压阀反馈腔的氮气压力达到目标值(压力变动率小于1%)。填充开始后约0.15 s，减压阀阀芯就开始出现快速振荡，振幅随着时间逐渐减小。阀芯的高频开关维持了减压阀反馈腔内压力的稳定。

2)当驱动电压为80 V时，从开机时刻到阀芯开启至稳定位移的响应时间约为0.64 ms，阀芯开启的稳定位移为3.67 μm。开机后约8 ms时，喷管推力达到稳定值(0.588 mN)。

3)当驱动电压分别为60 V、80 V、100 V和120 V时，针阀稳态位移分别为3.17 μm、3.67 μm、4.10 μm和4.49 μm，而喷管稳态推力分别为0.404 mN、0.588 mN、0.807 mN和1.021 mN。驱动电压与喷管稳态推力之间存在良好的线性关系。