液體火箭縱向耦合振動建模及其動態特性分析

郝 雨，徐得元，楊瓊梁，劉錦凡，柳征勇，唐國安

（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月生