液体火箭纵向耦合振动建模及其动态特性分析

郝 雨，徐得元，杨琼梁，刘锦凡，柳征勇，唐国安

（1.复旦大学力学与工程科学系，上海 200433；2.上海宇航系统工程研究所，上海 201108）

液体火箭纵向耦合振动建模及其动态特性分析

郝 雨1，徐得元1，杨琼梁2，刘锦凡2，柳征勇2，唐国安1

（1.复旦大学力学与工程科学系，上海 200433；2.上海宇航系统工程研究所，上海 201108）

针对液体火箭结构－推进系统的纵向耦合振动问题，建立推进系统各个组件的时域动力学模型，推导出结构－推进耦合系统二阶线性微分方程形式的控制方程。该方法能快速求解得到耦合系统的全部特征值，而且还能在此基础上分析特征值关于蓄压器和泵的物理参数的灵敏度，为液体火箭推进系统的参数优化提供技术手段。

液体推进剂火箭；耦合；振动；灵敏度

大型液体运载火箭在飞行的过程中，飞行器结构系统和推进系统会发生相互作用，使液体火箭产生不稳定的纵向耦合振动，从而在响应的时间历程曲线上出现一个鼓包，对运载火箭造成诸如仪器损坏和发动机早期关机等严重后果［1－3］。这种纵向耦合振动的稳定性与推进系统的蓄压器和泵的物理参数有关，为了更好地调节这些参数，分析耦合系统特征值对推进系统结构参数的灵敏性具有重要意义。

液体火箭纵向耦合振动的频域分析较多使用传递函数法［1，4］，这些方法需要建立推进系统从贮箱出口到发动机燃烧室出口的传递函数，进而与火箭结构系统进行耦合，但是这种方法最终的控制方程是一个复杂的超越函数，求解时容易漏根。徐得元等［5］在矩阵法［6］的基础上，利用有理多项式拟合，能够快速计算耦合系统的特征值，但与矩阵法相同，系统参数仍然隐含在算式中，给耦合振动的灵敏度分析造成困难。Oppenheim等［7］建立了液体运载火箭推进系统各组件的有限元模型，该模型将参数显式体现在系统微分方程中，为灵敏度计算创造了条件。但是Oppenheim等的模型采用流量作为推进系统的自变量，这样就需要引入多通管的超单元，不能采用标准的有限元过程进行组装。Zhao等［8］应用这一模型，对一个简单少自由度系统模型进行了参数分析，但这个模型过于简单，不能满足实际工程的需要。

在本文中，对以动量、质量守恒的偏微分方程控制的直管等组件内的推进剂流动直接建立时域微分方程，对于以频域特性描述的蓄压器等部分组件将传递函数进行等效变换，得到其时域控制方程，进而将其与火箭结构系统进行耦合，建立了以二阶线性常微分方程组形式表示的结构与推进系统纵向耦合振动的时域动力学模型，并求解了系统特征值与对部分推进系统参数的灵敏度。

1 推进系统部件控制方程

火箭的推进系统主要由输液管道、蓄压器、泵及泵后管、发动机等几部分组成，其示意图如图1所示。每一个部件中的推进剂压力和流量的控制方程将统一表示成二阶线性常微分方程组的形式，部分部件中推进剂对结构的扰动力也将统一表示成结构振动加速度和推进剂压力的线性组合，下面列出推进系统每个部件的控制方程。

图1 推进系统简化模型Fig.1 Simplified model for propulsion system

在以下推导过程中，如无特殊说明，将用变量p和q（包括其向量形式p和q）代表压力和流量的扰动量，大写的P（或P）和Q（或Q）表示其Laplace变换后的值。再用大写的拉丁字母下标表示推进系统的部件，如用T表示贮箱，D代表管段，A代表蓄压器，P代表泵，C代表发动机燃烧室。

（1）贮箱出口模型

记贮箱出口处扰动压力和流出流量分别为pout和qout，贮箱柱段面积为A，贮箱出口面积为Aout，推进剂液面高度h，密度为ρ0，文献［2］给出贮箱出口的传递函数关系式

将式（1）作简单的改写即可得到贮箱出口端推进剂的控制方程。

其中：

这里p1＝pout，而q1＝－qout为流入流量。

贮箱对火箭结构的扰动力包括两部分，

其中：

（2）直管段推进剂模型

考虑一段直管内推进剂的流动，直管截面积为A、推进剂的密度为ρ0、体积压缩模量为K0（考虑进管段弹性影响）。直管内推进剂的沿定常流动方向x的动量方程以及物态方程分别为

考虑到方程的相似性，将直管内推进剂压力和流量扰动可以用弹性直杆的纵向振动进行比拟。参照直杆振动的有限元法，可以建立直管内推进剂的时域模型。将长度为l的直段推进剂划分成n－1个均匀单元，则n个节点上的扰动压力和流量满足

对于弯管和波纹管，文献［1］给出了它们的传递函数，与直管段的传递函数形式上相同，其时域表达式也将具有类似形式，可以按照上面叙述的方法进行处理。

（3）蓄压器段推进剂模型

在Laplace域上，记连接蓄压器的上游直管出口、下游直管入口的扰动压力和流量在频域中分别为Pin（s），Qin（s），Pout（s）和Qout（s），其中s是Laplace变量。蓄压器的柔度、惯性和阻力系数分别为C，L和R。文献［2］给出了蓄压器频域的传递关系

（4）泵和泵后管组合段模型

根据泵与泵后管的传递函数，可以将泵和泵后管合为一部分进行分析和建模。记泵的动力学增益为m、气蚀柔度为C、入口面积为A、惯性系数和阻力系数为L′和R′，泵的入口扰动压力和流量在频域中分别为Pin（s），Qin（s）；泵后管的惯性系数和阻力系数为L′′和R′′，出口处的扰动压力和流量为Pout（s）和Qout（s）。则组合段传递函数为

其中：Z＝sL＋R，L＝L′＋L′′，R＝R′＋R′′为组合段的总阻抗。

这里：P1＝Pin，P2＝Pout，Q1＝Qin，Q2＝Qout。

泵前后的压力动量变化会对箭体结构产生一个作用力，由此造成的泵与箭体结构安装点的受力与泵入口处的压力p1、流量q1，泵的纵向运动速度x··p以及泵之前管路的横截面积Apre、泵的稳态体积流量qs有关，将流量通过式（12）消去，可以写成

2 耦合系统控制方程

2.1 推进系统控制方程

其中：上标k可取值1到7，分别表示贮箱、波纹管/直管/多通道联通器、波纹管＋弯管＋波纹管、蓄压器、摇摆软管＋活动阀门、泵及泵后管和发动机。p（k）和q（k）分别表示各部件中推进剂的节点压力（包括引入的辅助变量）；xT，xP分别代表贮箱和泵处的纵向位移，这两处与箭体结构相连，可以视为推进系统的输入，且xq＝｛xTxP｝T；M，C，K，S，E，F，G，H是仅与时间有关的系数矩阵。

组装推进系统整体模型时，需要对每个节点压力、包括引入的辅助变量编排序号，组件连接处的公共节点共用一个相同的编号。将每个部件的M（k），C（k），K（k）和q（k）等视为结构动力学有限元方法中超单元的质量、阻尼、刚度矩阵和节点力向量，采用标准的装配程序进行组装，形成推进系统的控制方程和扰动力表达式

由于讨论的是一个多发动机模型，假设每个发动机系统均相同，组装过程中，对于出现在分支上的组件（k＝3～7），应先将M（k），C（k），K（k）和q（k）分别乘以N以后再组装。这里N是发动机数目。

在进行组装时，各部件的入口、出口均会与相邻组件的入口、出口合并。在对推进系统的节点进行编号时，若将泵及泵后管组件引入的辅助变量编为最后，从第1节中可以看出，在（19）中矩阵Mf和Sf的最后一行全部为零，而剩下部分的质量矩阵满秩。这样可将推进系统的压力向量pf和流量向量qf分解成

2.2 推进系统与箭体的耦合

根据液体火箭结构系统纵向振动的动力学特性，在推进－结构纵向耦合振动分析中，采用弹簧－质量模型将火箭结构系统简化为具有集中参数的多自由度振动系统［7］。

结构的有限元动力学方程

其中fsf是液体对结构的作用力。根据前面叙述，fsf非零分量来自ff的相应位置。

用Lps表征管路系统对箭体结构作用位置的布尔矩阵，有：

xT和xP均来自结构节点位移向量xs的相应元素。用Lqs表征箭体对推进系统的作用位置的布尔矩阵，则：

由于推进系统中的结点变量是压力，而在箭体结构中其变量是位移，二者量纲不一致，其矩阵元素可能相差几个数量级，可能导致严重的数值问题。因此在耦合时首先将推进系统节点压力转化成位移量纲

2.3 灵敏度分析

上面得到了耦合系统状态空间的矩阵A，B，它们都是非对称的，其特征值问题可以表述为［9］

其中：λi是系统的第i阶特征值，xi和yi分别是对应的左、右特征向量。分析这个系统的特征值λi对设计变量的灵敏度。

假设设计变量为v，系统特征值随设计变量的偏导数

3 耦合系统的动态特性及灵敏度分析

3.1 耦合系统的特征值

矩阵法需要先计算出耦合系统的闭环或者开环传递函数，然后通过非线性方程求解得到复数特征根或者绘制Bode和Nyquist图等手段来进行稳定性分析。这一过程比较繁琐，也容易漏根。本文提出的方法通过矩阵特征值求解可以直接得到全部特征值。对于中小规模并且矩阵不出现病态的问题（一般工程问题均满足此条件），只要选取合适的特征值数值求解方法，就不会发生漏根现象。

作为应用，根据某液体火箭在100 s时的耦合系统参数，计算得到前10阶特征值，如表1所示，其中“矩阵法”系指利用传递函数在复平面中找根的方法，“有限元方法”指本文提出的方法。从表中可以看出，除少数几阶特征值的实部外，耦合系统的前几阶特征值均与矩阵法计算结果相差很小，验证了本文方法的正确性。

表1 不同方法计算的耦合系统特征值Tab.1 Eigenvalues of coup led system calculated by differentmethods

从计算结果中可以发现，耦合系统在第7阶特征值（38 Hz左右）出现正实部，表明耦合系统受扰动后的响应有可能会不断增加、出现失稳。由于箭体结构的时变性和非线性效应，特征值出现正实部并不意味着系统失稳，而是系统发生较大幅度的动态响应［7］。为了降低系统的动态响应、改善有效载荷在运载火箭主动段飞行时的动力学环境，调节推进系统管路中的蓄压器柔度等参数是较为常用的技术手段。在此可以通过分析管路中某些参数对耦合系统特征值的灵敏度，优化系统参数、减小特征值的正实部的绝对值，进而提高系统稳定性或降低系统的动态响应。

3.2 系统特征值对参数的灵敏度

工程设计中，蓄压器的柔度和液路惯性可以近似表示为其中：P0是蓄压器初始工作压力，V0是蓄压器初始工作容积，ν比热比，Pt是蓄压器工作压力，ρ是氧化剂密度，L是蓄压器通道长度，S是蓄压器通道截面积。

针对现役型号，推进系统设计已经定型，泵的特性无法改变，蓄压器的液路惯性与蓄压器结构相关参数也无法调整，从蓄压器的柔度近似公式（式（29）），可以通过改变蓄压器的初始工作压力和初始工作容积来实现改变蓄压器的柔度。针对新设计的蓄压器、泵，各个参数都是可以调节的。

且除了泵所对应的三个自由度的元素外，其他部分导数均为0。式中Lp是泵的惯性，其他参数参见第1节。

表2 灵敏度与差分近似值的对比Tab.2 Com paration of sensitivity and finite difference approximation

观察到这一阶特征值的实部随上述四个参数都有较明显的变化，如果将CA在原基础上减小1%，则第7阶特征值的实部将减小0.018 7。同理调整泵的柔度CP等参数也可以达到抑制发散的目的。

从表2中可以发现，以上4个参数中，系统特征值的正实部随蓄压器惯性变化最明显，在其他参数不变的情况下，只需将蓄压器惯性减小一半，系统第7阶特征值的实部将由正值变为负值，进而提高系统的稳定性。

4 结 论

上面建立了推进系统的时域动力学模型，并与箭体结构系统进行耦合，进而得到耦合系统的时域模型。该模型具有如下特点：

（1）耦合系统的控制方程与结构动力学形式上类似，便与应用标准结构动力学方法进行求解；

（2）推进系统各个部件能够用标准有限元程序进行组装，可以更好地适应推进系统结构的改变，通用性较强；

（3）控制方程中推进系统的每个部件在矩阵中都是相对分离的，并且都能够用显式表示，容易分析部件矩阵关于推进系统特性参数的关系。

应用此模型，根据某型火箭的推进系统管路参数，计算了耦合系统的特征值，并与传递函数法进行了比较，验证了本文方法的正确性。分析了耦合系统特征值对蓄压器、泵等推进系统参数的灵敏度，为如何设计推进系统管路参数，减小推进系统与火箭结构的耦合效应，提供了理论依据。

［1］马道远，王其政，荣克林.液体捆绑火箭POGO稳定性分析的闭环传递函数法［J］.强度与环境，2010，37（1）：1－7.

MA Dao-yuan，WANG Qi-zheng，RONG Ke-lin.Close-loop transfer function of POGO stability analysisfor bindingliquidpropellant rocket［J］.Structure＆Environment Engineering，2010，37（1）：1－7.

［2］王其政.结构耦合动力学［M］.北京：宇航出版社，1999：251－280.

［3］杨明，杨智春，史淑娟.液体火箭POGO振动时域分析方法研究［J］.强度与环境，2010，37（2）：29－36.

YANG Ming，YANG Zhi-chun，SHI Shu-juan.Method research of time domain analysis of POGO vibration of the liquid rocket［J］.Structure＆Environment Engineering，2010，37（2）：29－36.

［4］Rubin S.Longitudinal instability of liquid rockets due to propulsion feedback（POGO）［J］.Journal of Spacecraft and Rockets，1966，3（8）：1188－1195.

［5］徐得元，郝雨，杨琼梁，等.液体火箭结构－推进纵向耦合振动特性的快速求解方法［J］.宇航学报，2013.34（12）：1－7.

XU De-yuan，HAO Yu，YANG Qiong-liang，et al，Fast analysis method for structure-propulsion coupled vibration characteristics of liquid rocket［J］.J.Astronautics，2013，34（12）：1－7.

［6］Lock MH，Rubin S.Passive suppression of POGO on the space shuttle［R］.NASA CR－132452，1974.

［7］Oppenheim B W，Rubin S.Advanced pogo stability analysis for liquid rockets［J］，J.Spacecraft Rockets，1993，30（3）：360－373.

［8］Zhao Z H，Ren G X，Parameter study on pogo stability of liquid rockets［J］.J.Spacecraft Rockets，2011，48（3）：537－541.

［9］张义民，刘巧玲，闻邦椿.一般实矩阵特征值问题的多维灵敏度分析［J］.应用力学学报，2003，20（1）：107－109.

ZHANG Yi-min，LIU Qiao-ling，WEN Bang-chun.Multidimensional sensitivity analysis for eigenvalue problem of generalmatrices［J］.Chinese Journal of Applied Mechanics，2003，20（1）：107－109.

Modeling and dynamic characteristic analysis for longitudinal coup led vibration of a liquid-propulsion rocket

HAO Yu1，XU De-yuan1，YANGQiong-liang2，LIU Jin-fan2，LIU Zheng-yong2，TANGGuo-an1
（1.Department ofMechanics and Engineering Science，Fudan University，Shanghai200433，China；
2.Shanghai Academy of Spaceflight Technology，Shanghai201108，China）

For the problem of structure-propulsion system coupled longitudinal vibration of a liquid-propulsion rocket，dynamic model of each component of the propulsion system was established in time domain.The governing equation of the structure-propulsion coupled system was deduced in the form of second-order linear differential equations.With thismethod，all eigenvalues of the coupled system were obtained quickly，and their sensitivities with respect to physical parameters of the pressure accumulator and pump were calculated easily.The study results provided a technique for parametric optimization of rocketswith liquid-propulsion.

liquid-propulsion rocket；propulsion system-structure coupled vibration，sensitivity analysis

V475.1

A

10.13465/j.cnki.jvs.2014.24.012

国家自然科学基金（11202052）；973项目（613133）

2013－10－21 修改稿收到日期：2014－01－02

郝雨男，博士生，1988年12月生

唐国安男，教授，博士生导师，1962年10月生