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

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

0 引 言

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

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

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

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

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

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

(1)

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

(2)

(3)

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

(4)

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

由于結(jié)構(gòu)系統(tǒng)和推進(jìn)系統(tǒng)各部件的動力學(xué)方程都是不超過二階微分的動力學(xué)方程，因此引入狀態(tài)變量v，包括推進(jìn)系統(tǒng)的脈動重量位移w=[…,wi, …]T和脈動壓力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)定性分析。針對式(6)中E矩陣奇異的問題，文獻(xiàn)[12-13]提出了改進(jìn)方法，導(dǎo)出了式(6)形式的Pogo非奇異狀態(tài)空間方程，可同時用于時域仿真和穩(wěn)定性分析。

2 Pogo模型的降階方法

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

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

(8)

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

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

證.設(shè)λ是矩陣(E,A)的廣義特征值，那么廣義特征值矩陣(λE-A)的行列式為零，又因為其轉(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)的廣義特征值。因為λi是矩陣(E,A)的廣義特征值，所以存在對應(yīng)的特征向量φi滿足

λiEφi=Aφi

(10)

因為λj是矩陣(ET,AT)的廣義特征值，所以存在對應(yīng)的特征向量φt,j滿足

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

(11)

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

(12)

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

(13)

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

(14)

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

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

EΦΛ=AΦ

(15)

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

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

(16)

(17)

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

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

ETΦtΛ=ATΦt

(18)

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

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

(19)

利用特征向量矩陣Φ對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ù)高的問題，可以方便地開展時域仿真、控制系統(tǒng)設(shè)計等工作。求解得到η之后，利用變換式(24)可以求得原狀態(tài)x的解。

幾點說明：

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

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

3)η初值η(t0)的確定?；跔顟B(tài)微分方程(26)求解η(t)時，需要確定初值η(t0)。因為原狀態(tài)初值x(t0)已知，因此可以根據(jù)變換(24)的反變換式得到

(27)

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

(28)

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

(29)

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

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

(30)

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

(31)

3 仿真算例

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

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

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

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

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

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

圖3 推力室脈動壓力Fig.3 Perturbed pressures of the chamber

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

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

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

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

某型號火箭采用新型液氧-煤油富氧補(bǔ)燃循環(huán)發(fā)動機(jī)推進(jìn)系統(tǒng)[21]，除了傳統(tǒng)的液路還包含補(bǔ)燃的氣路，屬于氣-液路耦合的推進(jìn)系統(tǒng)，結(jié)構(gòu)系統(tǒng)取兩階縱向模態(tài)振動方程。由于氣路的存在，導(dǎo)出的如式(8)所示Pogo狀態(tài)空間模型是奇異的，維數(shù)是60階?；谠摖顟B(tài)空間模型進(jìn)行穩(wěn)定性分析，圖5給出了利用該模型分析得到的不同秒點的Pogo系統(tǒng)結(jié)構(gòu)阻尼比?？梢钥闯鲈揚ogo系統(tǒng)在0.2時間點附近第一階結(jié)構(gòu)阻尼比小于零，振動失穩(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 推力室脈動壓力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時間點附近的區(qū)段明顯發(fā)散，隨后收斂，這與圖5中的結(jié)構(gòu)第一階阻尼比曲線是一致的。同時結(jié)構(gòu)第二階模態(tài)位移在該時間段也出現(xiàn)一定程度的發(fā)散，而圖5中結(jié)構(gòu)第二階阻尼比一直大于零，說明結(jié)構(gòu)第二階模態(tài)雖然并未失穩(wěn)，但由于結(jié)構(gòu)第一階模態(tài)位移發(fā)散引起激勵項的顯著增大，從而導(dǎo)致第二階模態(tài)位移響應(yīng)也出現(xiàn)發(fā)散。圖7給出了推力室脈動壓力響應(yīng)，可以看出推力室脈動壓力尚未出現(xiàn)明顯的發(fā)散情況。從上面的分析來看，采用降階模型的仿真結(jié)果是合理的，校驗了降階方法的正確性。同時，也說明時域仿真與穩(wěn)定性分析相結(jié)合可以更好地揭示Pogo特征和機(jī)理。

4 結(jié) 論

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