黄夕轩,黄敏超,胡小平

(国防科技大学 空天科学学院,湖南 长沙 421000)

0 引言

液体火箭发动机推力可变技术是航天液体动力的重要发展范畴。采用变推力技术的液体火箭发动机可以实现航天运输系统的最佳推力控制,登月软着陆的实现离不开变推力发动机[1-4],空间姿轨控任务使用变推力发动机可以提高控制系统的灵活性。对于导弹系统,变推力发动机可以提高其突防能力。目前拥有高关注度的可重用发动机技术也要求具有推力可调性能的发动机[5]。

随着可变推力发动机的发展,追求更高推力比,实现稳定高效的变推过程成为未来的发展趋势[6-7]。对于变推发动机的稳定性能研究多以数值仿真进行。陈启智[8]等人对带有可调喷注器的双组元变推力液体火箭发动机出现的低频振荡现象进行了研究,分析了设计参数对稳定性的影响。文献[9]通过动力学模型的方法分别描述了以流量信号反馈和位移信号反馈两种方式的挤压式推进剂供应变推力液体火箭发动机的传递函数以及状态空间。在状态空间的基础上对两种反馈方式的稳定条件进行了分析,并且通过试车数据验证了计算得到的系统的稳定区域与不稳定区域,证明了采用动力学模型的方法寻找系统稳定性区间的可行性。Rodden[10]等人对可调喷注器进行了数学建模并绘制了喷注器调节系统流程图,结合频域分析,给出了喷注器的设计建议。文献[11]对流量调节器及泵压式供应系统建立了集中参数线性化动力学模型,研究了系统频率响应特性以及系统的固有稳定性,通过时域与频域的对比验证了频域分析的正确性,因此本文使用的集中参数线性化动力学模型是可以应用于频域分析的。刘景华[12]采用混合仿真方法,建立了过氧化氢发动机稳定性分析模型,得到脉冲扰动下发动机响应特性以及燃烧时滞、喷注压降等参数对稳定性影响。总体看来,对于动态特性的研究主要是通过数学方法建立状态空间的模型,进而分析频率响应特性,且大部分工作是对于挤压式系统进行的。目前对于泵压式变推力发动机的变推力过程的稳定性研究较少。

泵压式液体火箭发动机是使用涡轮泵在较低的压力下从推进剂贮箱中接收液体推进剂,并以所需的流量和喷射压力将其供给发动机燃烧室。涡轮的动力来自高压气体的膨胀,而高压气体通常是由涡轮泵输送的推进剂的燃烧产物。泵压式火箭发动机应用范围广泛:如20世纪60年代为土星五号运载火箭开发的F-1和J-2发动机; 20世纪70年代航天飞机使用的主引擎(SSME),使航天器首次实现可重复使用[13]。因此,对于泵压式液体火箭发动机频率响应方面的研究进行完善是必要的。

通过进行线性分析,可以了解系统状态空间的特征,可以反映时域仿真中不能明确体现的频率响应特性,进一步了解系统高、低频率下的幅相频特性。因此,本文在模拟了发动机启动过程及变推力过程的主要部件时域响应之外,用快速傅里叶方法提取了发动机的变推力过程推力室压力的主要频率,对于加入该频率扰动信号的推力室压力进行了线性分析,得到其频率响应特性。

1 发动机主要组件建模

文献[14]中对液氧煤油发动机起动仿真结果与试车数据进行了对比,发现2组数据的参数变化趋势与过渡时间吻合性较好,验证了模型的正确性与可行性。本文使用的物理模型与文献[14]相同,故不再进行验证。

1.1 仿真系统示意图

本文所研究燃气发生器循环发动机的系统如图1所示,该系统由燃料贮箱、氧化剂贮箱、涡轮、离心泵、燃气发生器、混合比调节器、推力室阀门、管道及节流孔板所构成。

图1 泵压式燃气发生器循环发动机系统示意图Fig.1 Schematic diagram of gas generator cycle engine system

1.2 管路模型

发动机管路采用集中参数法构建模型。集中参数法是将管路按照流阻、流感和流容3类基本元件进行有限元分割的方法。在低频动力学范围内广泛使用,因该方法的前提假设为小扰动、线性化,故不适合于中高频的震荡特性的研究[15]。

在该方法中,将主要考虑流体惯性的单元称为惯性单元(I),主要考虑局部阻力的单元称为阻性单元(R),主要考虑流体压缩性的单元设置为容性单元(C),同时考虑惯性和阻性的单元称为运动单元(IR)。

1.3 离心泵模型

建立发动机离心泵模型过程中使用了泵扬程及阻力矩的性能数据。离心泵性能数据由文献[16-17]中实验数据拟合得到。泵的特性曲线是将一维数据表拟合成关于自变量θ的函数而绘成,其中,θ为反映泵工作特性的特征角,其定义为

θ=π+arctan(ν/n)

(1)

式中ν和n分别为流量和转速的无量纲特征数,其定义为

(2)

扬程和转矩的无量纲参数h和β定义如下

(3)

式中：H为实际扬程;HR为参考扬程, m;ηR为泵参考效率;ρ为泵入口和出口的工质平均密度,kg/m3。

1.4 燃烧室模型

由于燃烧过程对发动机低频特性影响较小,燃烧零维模型被广泛用于发动机动态特性分析[20-21]。燃烧室模型采用基于燃烧时滞的燃烧动力学模型,且假设燃烧时滞是一个常数,燃气在燃烧室是均匀分布的,各参数不随空间坐标变化。

质量微分方程

(4)

能量微分方程

(5)

混合比微分方程

(6)

式中：ρ为燃烧室平均密度;p为燃烧室室压;kc为燃烧室混合比(O/F);γ为定压比热与定容比热之比,是关于时间的函数;qmf为燃料质量流量;qmo为氧化剂质量流量;V为燃烧室体积;hf和ho分别为燃料和氧化剂的比焓;R为燃气气体常数;T为燃气温度。

2 结果与讨论

2.1 工况说明

本文设置了2个工况：其中每个工况对应的燃料及氧化剂入口流量如表1所示,入口流量调节是通过调节氧化剂主阀与燃料主阀开度实现的。

表1 设置工况说明

本文设置了2个工况转换点,分别位于仿真时间的4 s由高工况转变为低工况和7 s由低工况转变为高工况。切换工况时推进剂主阀开度改变,同时主系统工况切换阀和副系统工况切换阀切换工作状态。

2.2 发动机启动及变工况仿真结果

图2～图4分别给出启动后经过工况转换的发动机系统主要部件主要参数随时间变化曲线。推力室室压在启动时出现超调,在高工况时稳定在4.2 MPa,低工况时稳定在1.7 MPa,在2个工况切换点均未出现超调。离心泵转速启动后在升至36 000 r/min中无超调,在转换为低工况后有极小的超调,之后稳定在37 500 r/min。涡轮功率在两个工况切换点均出现了超调,且超调量较大,反映了副系统受变工况的影响较大。

图2 推力室室压随时间变化曲线Fig.2 Thrust chamber pressure vs time

图3 离心泵转速随时间变化曲线Fig.3 Centrifugal pump speed vs time

图4 涡轮功率随时间变化曲线Fig.4 Turbine power vs time

2.3 系统特征频率与阻尼系数分析

在燃气发生器循环液体火箭发动机系统中,燃气发生器及推力室参数被作为主要被控对象。将燃气发生器与推力室压力等参数设置为观测变量,分别在2种工况达到稳定状态时考察系统的特征频率及阻尼系数。根据各特征频率下观测变量的模态振型,系统的特征频率可分为2类：对应于燃气发生器的特征频率与对应于4台并联的推力室的特征频率。由表2可以看出,4台并联的推力室特征频率非常接近。

表2 高工况系统阻尼系数及频率

Tab.2 Eigenvalues and frequency of system at high operating condition

特征频率/Hz阻尼比阻尼系数实部阻尼系数虚部6.476 01.00-40.691 40.0056.860 21.00-357.264 00.0056.861 01.00-357.268 00.0056.861 51.00-357.271 00.0056.861 61.00-357.272 00.00110.777 01.00-696.036 00.00110.778 01.00-696.042 00.00110.779 01.00-696.045 00.00110.779 01.00-696.046 00.00

表2所示的高工况9个特征频率中,对应于燃气发生器压力的特征频率是6.476 Hz。4个在56.86 Hz附近的和4个在110.78 Hz附近的特征频率对应4台推力室的特征频率。高工况下系统所有阻尼系数均为实数,且实部为负数。这表明系统是稳定的。

表3中所示的低工况9个特征频率中,燃气发生器特征频率为17.997 Hz。其中4个为46.563 1 Hz和4个在94.896 4 Hz附近的特征频率为4台推力室的特征频率。相对于高工况,低工况下系统所有特征频率及阻尼系数均发生了改变,且阻尼系数全部为负实数。这是由于变工况时,主阀及主、副系统工况切换阀工作状态改变,导致通流面积改变,因而引起管道局部阻力的变化,而流阻的变化必将改变系统的动力学方程参数,进而影响系统状态空间,最终表现为特征频率与阻尼系数的改变。

表3 低工况系统阻尼系数及频率

Tab.3 Eigenvalues and frequency of system atlow operating condition

特征频率/Hz阻尼比阻尼系数实部阻尼系数虚部17.997 01.00-113.078 30.0046.563 11.00-292.565 00.0046.563 11.00-292.565 00.0046.563 11.00-292.565 00.0046.563 11.00-292.565 00.0094.896 41.00-596.252 00.0094.896 41.00-596.252 00.0094.896 41.00-596.252 00.0094.896 41.00-596.252 00.00

总体看来,系统一阶特征频率对应的阻尼系数最小,表明系统对该频率较敏感。该特征频率是燃气发生器压力的特征频率,说明由于副系统对低频信号更敏感。

由阻尼系数变化可知系统由高工况变为低工况后系统对低频敏感性减小,对高频敏感性增加,但影响系统稳定性的主要信号仍然是低频信号。

2.4 频率响应分析

本文对第一个变工况点后的变推力段推力室室压进行快速傅里叶变换,如图5所示,得到其最主要的频率成分在10 Hz以内,在5.34 Hz达到了峰值,验证了上文对特征频率的分析。为了研究该频率下的入口压力震荡对系统稳定性的影响,将扰动设置为一正弦波,因此推进剂入口压力波方程为

(7)

式中：A=10%PIn0；f=5.34 Hz；PIn0为推进剂入口压力；φ为相位差；PIn0为0.49 MPa；φ=0°。本文选取4.01 s,4.02 s,4.03 s,4.04 s,4.05 s 5个时间点,将推力室压力设为输出信号,绘制以入口压力为输入信号的Nyquist图,结果如图6～图10所示。

图5 变推力过程室压FFT分析 Fig.5 FFT analysis of chamber pressure of thrustvarying process

图6 4.01 s推力室压力的Nyquist曲线Fig.6 Nyquist curve of thrust chamber pressure at 4.01 s

图7 4.02 s推力室压力的Nyquist图Fig.7 Nyquist curve of thrust chamber pressure at 4.02 s

图8 4.03 s推力室压力的Nyquist图Fig.8 Nyquist curve of thrust chamber pressure at 4.03 s

图9 4.04 s推力室压力的Nyquist图Fig.9 Nyquist curve of thrust chamber pressure at 4.04 s

图10 4.05 s推力室压力的Nyquist图Fig.10 Nyquist curve of thrust chamber pressure at 4.05 s

Nyquist图反映了相对于控制信号(入口压力),输出信号(推力室压力)的相位差以及幅值比,可以观察到在高频区曲线的行为趋向于二阶震荡环节,但由于幅值较小没有影响系统稳定性。

根据Nyquist判据,在开环系统稳定时,Nyquist曲线不包围(-1,j0)点,闭环系统也稳定。图6～图10中Nyquist曲线均未包围(-1,j0)点。因此将该系统在所设置扰动下变工况过程稳定。

系统的幅相频特性曲线的最大幅值小于1时可知系统的相位不会改变系统的稳定性,系统的幅值裕度可以由曲线与负实轴交界处幅值求得。各测点的稳定裕度如表4所示。在系统变工况过程中, 4.03 s低频段出现了幅值过大的现象,且在相角滞后大于30.856°后会出现系统不稳定现象。4.05 s低频段时出现了较低的稳定裕度,目前系统的开环增益增大3.404 4倍则系统失稳。由此可知变工况的过程中低频扰动信号对于稳定性影响不容忽略。

表4 各观测点稳定裕度

3 结论

本文对泵压式双组元燃气发生器循环液体火箭发动机系统变工况频率响应特性进行了研究,可以得到以下一些结论：

1)本文所研究的发动机系统在高工况与低工况下的特征频率存在明显差别：系统的一阶特征频率(6.476 Hz和17.996 Hz)均为燃气发生器压力对应的特征频率,属于低频范围,是影响系统稳定性的主要频率。并联的4台推力室特征频率十分接近,特征频率相对较高,且阻尼系数较大,对系统影响相对较小。

2)由高工况下转变为低工况过程中,该发动机系统阻尼系数改变。由阻尼系数变化可以得出：系统由高工况变为低工况后系统对低频敏感性减小,对高频敏感性增加,但影响系统稳定性的主要信号仍然是低频信号。因此对于该发动机系统推力较大的工况,应当特别注意频率为10 Hz以内的信号的影响。

3)变工况过程中, Nyquist曲线反映了高频段系统偏向于二阶震荡环节但幅值小不影响系统稳定性。低频段曲线反映出的加入低频压力扰动信号可以造成推力室压力的稳定裕度变化,出现了较小的幅值裕度(3.404 4)。由于出现相角裕度小于40°的点,说明系统若存在相角延迟容易造成系统失稳的现象,进一步说明变工况过程容易受到低频压力入口扰动影响而失稳。系统在工作时,为了维持更高的系统稳定性裕度应该尽可能避免低频入口压力扰动信号的出现。