基于特征空間變換的運(yùn)載火箭Pogo模型降階方法

譚述君，高 強(qiáng)，趙 旺，劉錦凡

(1.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024；2.大連理工大學(xué)遼寧省空天飛行器前沿技術(shù)重點(diǎn)實(shí)驗(yàn)室，大連 116024；3.上海宇航系統(tǒng)工程研究所，上海 201108)

0 引 言

Pogo是大型液體運(yùn)載火箭在飛行過程中由推進(jìn)系統(tǒng)液路脈動(dòng)與結(jié)構(gòu)振動(dòng)相互耦合而產(chǎn)生的一種自激振動(dòng)現(xiàn)象[1]。Pogo會(huì)導(dǎo)致運(yùn)載火箭結(jié)構(gòu)破壞、影響飛行穩(wěn)定性，甚至威脅到航天員的生命等[2-4]，而且Pogo在不同型號(hào)火箭、同一型號(hào)火箭的不同任務(wù)飛行中都曾出現(xiàn)過，具有任務(wù)依賴性，發(fā)生機(jī)理復(fù)雜[5]。因此，開展Pogo建模方法研究對(duì)揭示Pogo發(fā)生機(jī)理和抑制Pogo振動(dòng)具有重要意義。

國(guó)內(nèi)外學(xué)者對(duì)運(yùn)載火箭Pogo建模開展了大量研究，其中傳遞矩陣法和閉環(huán)傳遞函數(shù)法因其簡(jiǎn)潔得到了廣泛應(yīng)用[6-9]，然而這些方法難以處理捆綁火箭的助推與芯級(jí)耦合、多模態(tài)耦合、多作用力耦合等問題[3]。王建民等[10]提出對(duì)捆綁火箭的Pogo問題必須基于空間模態(tài)進(jìn)行設(shè)計(jì)。Oppenheim與Rubin[11]面向工程中的復(fù)雜液路系統(tǒng)提出了一種基于有限元描述的狀態(tài)方程法，可以推廣到復(fù)雜三維管路的建模，已成功應(yīng)用于中國(guó)CZ-2F火箭Pogo問題的研究[3]。在此基礎(chǔ)上，Tan等[12]和Wang等[13]基于獨(dú)立流量位移提出改進(jìn)的狀態(tài)方程法，解決了Rubin方法的奇異性問題。文獻(xiàn)[14-15]則通過對(duì)阻抗函數(shù)矩陣進(jìn)行有理多項(xiàng)式逼近，提出了快速求解Pogo振動(dòng)特征值和穩(wěn)定性分析的方法。牛澤雄等[16]采用分段三次Hermite插值得到了管路系統(tǒng)描述的新形式。Park等[17]改進(jìn)了結(jié)構(gòu)振動(dòng)的三維和一維單元混合的有限元建模方法以更精確地進(jìn)行Pogo分析和預(yù)測(cè)。趙治華[18]則提出了Pogo振動(dòng)的多體動(dòng)力學(xué)建模方法，分析了包含結(jié)構(gòu)縱橫扭模態(tài)的液體火箭Pogo振動(dòng)穩(wěn)定性。上述建模方法更真實(shí)地描述了Pogo特性，同時(shí)也導(dǎo)致狀態(tài)空間模型的維數(shù)過高和奇異性問題。對(duì)此，王慶偉等[19]提出了一種Pogo模型縮聚方法，通過對(duì)液體輸送管路的局部降階實(shí)現(xiàn)對(duì)整體Pogo模型的降階，驗(yàn)證了降階的正確性和可行性。然而該方法僅限于對(duì)液體輸送管路的降階，不能處理眾多推進(jìn)系統(tǒng)部件、支路連接產(chǎn)生的高維問題。因此，需要進(jìn)一步發(fā)展適用于不同型號(hào)推進(jìn)系統(tǒng)的運(yùn)載火箭Pogo模型的一般性方法。

本文基于特征空間變換理論提出了一種面向運(yùn)載火箭Pogo狀態(tài)空間模型的整體降階方法。該方法具有很好的通用性，并通過對(duì)兩種不同型號(hào)推進(jìn)系統(tǒng)的運(yùn)載火箭Pogo問題的分析和仿真，校驗(yàn)了所提出的降階方法的正確性和通用性。

1 Pogo狀態(tài)空間模型

運(yùn)載火箭的Pogo振動(dòng)方程由結(jié)構(gòu)系統(tǒng)的振動(dòng)方程和推進(jìn)系統(tǒng)的脈動(dòng)方程相互耦合生成，其中推進(jìn)系統(tǒng)由眾多不同特性的物理部件組成，動(dòng)力學(xué)方程復(fù)雜，因此推進(jìn)系統(tǒng)脈動(dòng)方程的建立是其關(guān)鍵。Oppenheim等[11]提出的狀態(tài)方程法采用了有限元描述推進(jìn)系統(tǒng)各物理部件的系統(tǒng)方程，可以嚴(yán)格處理流體容器的任意運(yùn)動(dòng)以及作用于結(jié)構(gòu)系統(tǒng)的所有力，從而建立了一套系統(tǒng)的分析方法。對(duì)于典型的運(yùn)載火箭推進(jìn)系統(tǒng)所包含的貯箱、燃料輸送管路、蓄壓器、泵、推力室等物理部件，該方法的優(yōu)點(diǎn)是可將任一部件統(tǒng)一描述為不超過二階微分的動(dòng)力學(xué)方程。

以可壓縮直管為例[11]，如果將管路劃分為Nd段，只要每段管路單元的長(zhǎng)度ld滿足

(1)

式中：a為液路壓力傳播的有效聲速，fu是系統(tǒng)分析所關(guān)心的最高頻率,那么每一段的動(dòng)力學(xué)方程可表示為

(2)

(3)

對(duì)于結(jié)構(gòu)系統(tǒng)，采用二階微分的模態(tài)振動(dòng)方程來(lái)描述。如第n階模態(tài)振動(dòng)方程為

(4)

式中：qn為結(jié)構(gòu)振動(dòng)的模態(tài)位移，Mn，ωn和ξn分別為第n階模態(tài)質(zhì)量、模態(tài)頻率和模態(tài)阻尼比，fi是第i個(gè)推進(jìn)系統(tǒng)單元作用于結(jié)構(gòu)的力向量，φin則為第i個(gè)推進(jìn)系統(tǒng)單元處的結(jié)構(gòu)第n階模態(tài)陣型向量，右端項(xiàng)體現(xiàn)了推進(jìn)系統(tǒng)對(duì)結(jié)構(gòu)振動(dòng)的影響。

由于結(jié)構(gòu)系統(tǒng)和推進(jìn)系統(tǒng)各部件的動(dòng)力學(xué)方程都是不超過二階微分的動(dòng)力學(xué)方程，因此引入狀態(tài)變量v，包括推進(jìn)系統(tǒng)的脈動(dòng)重量位移w=[…,wi, …]T和脈動(dòng)壓力p=[…,pi, …]T，以及結(jié)構(gòu)系統(tǒng)的模態(tài)位移q=[…,qn, …]T，那么結(jié)構(gòu)系統(tǒng)和推進(jìn)系統(tǒng)耦合的Pogo方程可以寫成標(biāo)準(zhǔn)的二階微分矩陣方程形式[11]，

(5)

(6)

其中,

(7)

基于狀態(tài)空間方程式，通過求解廣義特征值問題,可以進(jìn)行Pogo穩(wěn)定性分析。針對(duì)式(6)中E矩陣奇異的問題，文獻(xiàn)[12-13]提出了改進(jìn)方法，導(dǎo)出了式(6)形式的Pogo非奇異狀態(tài)空間方程，可同時(shí)用于時(shí)域仿真和穩(wěn)定性分析。

2 Pogo模型的降階方法

狀態(tài)方程法為Pogo分析提供了一套系統(tǒng)方法，在實(shí)際工程中得到了成功的應(yīng)用。然而，正如式(1)所述，為了提高Pogo分析的精度，需要將輸送管路劃分為更多的段數(shù)，這將導(dǎo)致狀態(tài)方程式(6)的維數(shù)更高，增加計(jì)算困難。另一方面，文獻(xiàn)[12-13]解決了液路推進(jìn)系統(tǒng)狀態(tài)空間模型的奇異性問題，但對(duì)于含氣路推進(jìn)系統(tǒng)的Pogo建模，仍然存在奇異問題。這給時(shí)域仿真和控制系統(tǒng)設(shè)計(jì)帶來(lái)了極大困難。因此，本節(jié)給出Pogo狀態(tài)空間模型的降階方法。

不失一般性，將建立的Pogo狀態(tài)空間方程描述為

(8)

式中：E和A是系統(tǒng)矩陣，f是外部干擾或控制力。對(duì)矩陣E的奇異性不做要求。

定理1. 矩陣(E,A)和矩陣(ET,AT)具有相同的廣義特征值。

證.設(shè)λ是矩陣(E,A)的廣義特征值，那么廣義特征值矩陣(λE-A)的行列式為零，又因?yàn)槠滢D(zhuǎn)置矩陣具有相同的行列式[20]，即

0=|λE-A|=|(λE-A)T|=|λET-AT|

(9)

所以λ也是矩陣(ET,AT)的廣義特征值。即矩陣(E,A)和矩陣(ET,AT)具有相同的廣義特征值。證畢。

證. 設(shè)λi和λj分別是矩陣(E,A)的第i階和第j階廣義特征值，根據(jù)定理1，λi和λj也是矩陣(ET,AT)的廣義特征值。因?yàn)棣薸是矩陣(E,A)的廣義特征值，所以存在對(duì)應(yīng)的特征向量φi滿足

λiEφi=Aφi

(10)

因?yàn)棣薺是矩陣(ET,AT)的廣義特征值，所以存在對(duì)應(yīng)的特征向量φt,j滿足

λjETφt,j=ATφt,j

(11)

將式(11)兩邊轉(zhuǎn)置，并同時(shí)右乘向量φi，得到

(12)

將式(10)代入式(12)，整理得到

(13)

因?yàn)榫仃?E,A)是Pogo狀態(tài)空間方程的系統(tǒng)矩陣，每一個(gè)特征值對(duì)應(yīng)系統(tǒng)的某一階模態(tài)，包括結(jié)構(gòu)振動(dòng)模態(tài)與推進(jìn)系統(tǒng)的液路脈動(dòng)模態(tài)，所以各階特征值不同，故當(dāng)i≠j時(shí)，λi≠λj。因此由式(13)可以導(dǎo)出

(14)

接下來(lái)，本節(jié)將基于上述定理，通過特征空間變換理論進(jìn)行Pogo狀態(tài)空間方程(8)的降階研究。

對(duì)Pogo狀態(tài)空間方程的系統(tǒng)矩陣(E,A)進(jìn)行廣義特征值分析，即

EΦΛ=AΦ

(15)

式中：Λ和Φ分別是對(duì)角的特征值矩陣和特征向量矩陣。將Λ按特征值的模從小到大排列，特征向量Φ也對(duì)應(yīng)地排列，如下所示

Λ=diag(λ1，λ2，…，λN)

(16)

(17)

式中：λ1,λ2,…,λN是特征值，且|λ1|≤|λ2|≤…≤|λN|。

對(duì)矩陣(ET,AT)進(jìn)行廣義特征值分析。由定理1可知，式(15)和式(16)中的Λ也是矩陣(ET,AT)的廣義特征值，即

ETΦtΛ=ATΦt

(18)

式中：特征向量矩陣Φt也按式(16)中的特征值對(duì)應(yīng)排列，如下所示

Φt=[φt,1, φt,2, …, φt,N]

(19)

利用特征向量矩陣Φ對(duì)Pogo狀態(tài)空間方程式(8)的原狀態(tài)變量x進(jìn)行變換如下

x=Φη

(20)

將式(20)代入式(8)，得到

(21)

將式(15)代入式(21)，得到

(22)

(23)

(24)

(25)

(26)

式(26)就是用狀態(tài)η描述的非奇異Pogo降階模型，解決了原Pogo模型的奇異性和維數(shù)高的問題，可以方便地開展時(shí)域仿真、控制系統(tǒng)設(shè)計(jì)等工作。求解得到η之后，利用變換式(24)可以求得原狀態(tài)x的解。

幾點(diǎn)說明：

1)一般來(lái)說，振動(dòng)系統(tǒng)的響應(yīng)通過截?cái)嗌贁?shù)低階特征值及其特征向量就可以實(shí)現(xiàn)高精度的近似，因此對(duì)于降階模型，感興趣的特征值通常取式(16)中的低階特征值。

2)即使矩陣E是奇異的，選取前面低階特征值及其特征向量導(dǎo)出的降階系統(tǒng)(26)也是非奇異的。因此該方法也可以實(shí)現(xiàn)原狀態(tài)方程式(8)的去奇異化處理。

3)η初值η(t0)的確定。基于狀態(tài)微分方程(26)求解η(t)時(shí)，需要確定初值η(t0)。因?yàn)樵瓲顟B(tài)初值x(t0)已知，因此可以根據(jù)變換(24)的反變換式得到

(27)

4)運(yùn)算過程的實(shí)數(shù)化。由于E和A是非對(duì)稱實(shí)矩陣，因此特征值Λ和特征向量矩陣Φ可能是復(fù)數(shù)，而如果是復(fù)數(shù)則必然是復(fù)共軛成對(duì)出現(xiàn)的[20]，與結(jié)構(gòu)系統(tǒng)或推進(jìn)系統(tǒng)的二階振動(dòng)方程相對(duì)應(yīng)。因此矩陣的特征值和特征向量可以描述為

(28)

Φ=[Φr+iΦi,Φr-iΦi]

(29)

式中：Λr和Λi分別是特征值Λ的實(shí)部和虛部，Φr和Φi分別是特征向量矩陣Φ的實(shí)部和虛部。

將式(28)和式(29)代入式(15)，分離實(shí)部和虛部，得到

(30)

因此，前述方法中的特征值和特征向量可以用下面實(shí)數(shù)矩陣描述，即

(31)

3 仿真算例

本節(jié)分別利用液路推進(jìn)系統(tǒng)的氣-液路耦合推進(jìn)系統(tǒng)的兩種不同型號(hào)火箭的Pogo分析與仿真，來(lái)校驗(yàn)本文降階方法和去奇異化方法的正確性。需要說明的是，在不影響對(duì)比驗(yàn)證的前提下，根據(jù)保密要求，本文對(duì)仿真結(jié)果的數(shù)據(jù)進(jìn)行了適當(dāng)歸一化處理。

算例1：液路推進(jìn)系統(tǒng)的某型號(hào)火箭

某型號(hào)火箭采用液氫-液氧發(fā)動(dòng)機(jī)，推進(jìn)系統(tǒng)完全由液體輸送管路連接構(gòu)成，結(jié)構(gòu)系統(tǒng)采用兩階縱向模態(tài)振動(dòng)方程。對(duì)此，可采用文獻(xiàn)[12-13]提出的改進(jìn)的Rubin狀態(tài)空間建模方法，得到式(8)所示的Pogo非奇異狀態(tài)空間模型，直接應(yīng)用于時(shí)域仿真。然而，由于輸送管路需要?jiǎng)澐趾芏喽螖?shù)，導(dǎo)出的狀態(tài)空間模型維數(shù)高達(dá)45階(如果采用原Rubin方法建模[11]，導(dǎo)出的狀態(tài)空間模型奇異，且維數(shù)增加近一倍)。

圖1 預(yù)壓泵入口隨機(jī)脈動(dòng)干擾Fig.1 Random perturbation at the inlet of prepressure pump

圖2和圖3分別給出了隨機(jī)脈動(dòng)干擾下的結(jié)構(gòu)模態(tài)位移和推力室脈動(dòng)壓力的響應(yīng)曲線。圖4給出了降階模型與原模型仿真響應(yīng)相對(duì)誤差曲線。可以看出相對(duì)誤差小于10-1量級(jí)，大部分在10-2量級(jí)。相對(duì)于求解器精度10-3來(lái)說，降階模型的維數(shù)雖然僅為原模型的1/4，但是仿真結(jié)果與原模型高度吻合，校驗(yàn)了降階模型的正確性。這也說明降階模型保留了原模型的主要特性，由于管路劃分產(chǎn)生的高階模態(tài)是非主要模態(tài)，是可以忽略的，從而為Pogo主動(dòng)抑制設(shè)計(jì)提供了更合適的控制系統(tǒng)模型。

圖2 結(jié)構(gòu)模態(tài)位移Fig.2 Modal displacements of structure systems

圖3 推力室脈動(dòng)壓力Fig.3 Perturbed pressures of the chamber

圖4 降階模型與原模型仿真響應(yīng)相對(duì)誤差Fig.4 Relative errors between the reduced model and the original model

表1進(jìn)一步給出了MATLAB的ode45(·)函數(shù)在不同求解精度下，采用降階模型和原模型仿真耗費(fèi)的時(shí)間(程序在Intel(R) Core(TM) i7-8550U CPU處理器的筆記本電腦上運(yùn)行)。表1中數(shù)據(jù)顯示當(dāng)求解精度較低時(shí)(10-3)，利用降階模型仿真的效率是原模型的14倍；當(dāng)求解精度較高時(shí)(10-6)，利用降階模型仿真的效率是原模型的50倍。因此，利用降階模型仿真的計(jì)算效率大大提高，提高了一個(gè)數(shù)量級(jí)以上。從表1還可以看出，當(dāng)求解精度從10-3提高到10-6時(shí)，原模型由于維數(shù)較高(其中含有大量不起主要作用的高階模態(tài))，計(jì)算時(shí)間增加了298%；而降階模型由于顯著降低了維數(shù)，計(jì)算時(shí)間僅增加了11%，計(jì)算效率優(yōu)勢(shì)更為明顯。

表1 降階模型和原模型仿真的計(jì)算時(shí)間Table 1 Computational time for the reduced model and the original model

算例2：氣-液路耦合推進(jìn)系統(tǒng)的某型號(hào)火箭

某型號(hào)火箭采用新型液氧-煤油富氧補(bǔ)燃循環(huán)發(fā)動(dòng)機(jī)推進(jìn)系統(tǒng)[21]，除了傳統(tǒng)的液路還包含補(bǔ)燃的氣路，屬于氣-液路耦合的推進(jìn)系統(tǒng)，結(jié)構(gòu)系統(tǒng)取兩階縱向模態(tài)振動(dòng)方程。由于氣路的存在，導(dǎo)出的如式(8)所示Pogo狀態(tài)空間模型是奇異的，維數(shù)是60階?；谠摖顟B(tài)空間模型進(jìn)行穩(wěn)定性分析，圖5給出了利用該模型分析得到的不同秒點(diǎn)的Pogo系統(tǒng)結(jié)構(gòu)阻尼比?？梢钥闯鲈揚(yáng)ogo系統(tǒng)在0.2時(shí)間點(diǎn)附近第一階結(jié)構(gòu)阻尼比小于零，振動(dòng)失穩(wěn)。

圖5 Pogo系統(tǒng)的結(jié)構(gòu)阻尼比Fig.5 Structural damp ratios of the Pogo system

圖6 結(jié)構(gòu)模態(tài)位移Fig.6 Modal displacements of structure systems

圖7 推力室脈動(dòng)壓力Fig.7 Perturbed pressures of the chamber

圖6給出了Pogo系統(tǒng)的結(jié)構(gòu)模態(tài)位移響應(yīng)曲線，從中可以看出結(jié)構(gòu)第一階模態(tài)位移在0.2時(shí)間點(diǎn)附近的區(qū)段明顯發(fā)散，隨后收斂，這與圖5中的結(jié)構(gòu)第一階阻尼比曲線是一致的。同時(shí)結(jié)構(gòu)第二階模態(tài)位移在該時(shí)間段也出現(xiàn)一定程度的發(fā)散，而圖5中結(jié)構(gòu)第二階阻尼比一直大于零，說明結(jié)構(gòu)第二階模態(tài)雖然并未失穩(wěn)，但由于結(jié)構(gòu)第一階模態(tài)位移發(fā)散引起激勵(lì)項(xiàng)的顯著增大，從而導(dǎo)致第二階模態(tài)位移響應(yīng)也出現(xiàn)發(fā)散。圖7給出了推力室脈動(dòng)壓力響應(yīng)，可以看出推力室脈動(dòng)壓力尚未出現(xiàn)明顯的發(fā)散情況。從上面的分析來(lái)看，采用降階模型的仿真結(jié)果是合理的，校驗(yàn)了降階方法的正確性。同時(shí)，也說明時(shí)域仿真與穩(wěn)定性分析相結(jié)合可以更好地揭示Pogo特征和機(jī)理。

4 結(jié) 論

狀態(tài)方程法是當(dāng)前運(yùn)載火箭Pogo振動(dòng)研究中非常有優(yōu)勢(shì)的一種方法，同時(shí)也存在系統(tǒng)維數(shù)過高、存在奇異性等問題。本文采用特征空間變換理論開展了Pogo狀態(tài)空間模型降階方法的研究，并在兩種不同類型推進(jìn)系統(tǒng)的火箭Pogo問題中進(jìn)行了仿真驗(yàn)證。研究結(jié)果表明，本文所提的降階方法具有一般性，對(duì)奇異、非奇異Pogo狀態(tài)空間模型都能得到正確的降階模型，而且降階效果明顯，顯著提高了計(jì)算效率。同時(shí)，研究顯示基于降階模型可以方便地進(jìn)行Pogo時(shí)域仿真，而時(shí)域仿真結(jié)果能為Pogo特征描述和機(jī)理揭示提供重要的信息。進(jìn)一步，降階模型也更有利于Pogo主動(dòng)抑制器的設(shè)計(jì)，為正引起關(guān)注的Pogo主動(dòng)抑制研究[22]打下良好的基礎(chǔ)。