王慶偉, 王小軍, 張青松, 潘 輝, 譚述君

(1. 北京宇航系統(tǒng)工程研究所， 北京 100076; 2. 中國運(yùn)載火箭研究院， 北京 100076；3. 大連理工大學(xué) 航空航天學(xué)院， 大連 116024)

POGO振動是液體火箭在飛行過程中由于結(jié)構(gòu)系統(tǒng)與推進(jìn)系統(tǒng)耦合而產(chǎn)生的一種縱向自激振動，Titan II、Diamond-B、Saturn V和我國的長征二F火箭上都出現(xiàn)過不同程度的POGO振動，其產(chǎn)生的縱向振動加速度對運(yùn)載火箭的安全性和可靠性造成非常不利的影響，特別是對宇航員造成極大的生理痛苦[1-7]。因此，抑制POGO振動已成為液體火箭研制過程中重要環(huán)節(jié)之一。Rubin等將推進(jìn)系統(tǒng)簡化為單推進(jìn)劑單發(fā)動機(jī)模型，推導(dǎo)了POGO振動產(chǎn)生的機(jī)理并提出了判斷POGO穩(wěn)定性的臨界阻尼比判別法。Zhao等也利用簡化的單推進(jìn)劑單發(fā)動機(jī)推進(jìn)系統(tǒng)模型，推導(dǎo)出五個無量綱的參數(shù)來判定POGO振動的穩(wěn)定性并提出了耦合強(qiáng)度的概念。譚述君等[8]進(jìn)一步完善了耦合強(qiáng)度理論，深入研究了系統(tǒng)阻尼比對POGO穩(wěn)定性的非線性影響關(guān)系，指出臨界阻尼比法的局限性并提出了臨界耦合強(qiáng)度判別法，通過耦合強(qiáng)度和臨界耦合強(qiáng)度的大小關(guān)系即可判斷耦合系統(tǒng)穩(wěn)定性。徐得元等[9]利用有理多項式逼近的方法提出了POGO系統(tǒng)振動特性的一種快速求解方法。郝雨等[10]得出了結(jié)構(gòu)-推進(jìn)耦合系統(tǒng)的二階線性微分方程，并提出了通過求靈敏度進(jìn)行推進(jìn)系統(tǒng)參數(shù)優(yōu)化的方法。唐冶等[11]對液體火箭推進(jìn)系統(tǒng)頻率特性進(jìn)行了靈敏度分析。Oppenheim等[12]提出了一種POGO振動模型建模的有限元法，得出了狀態(tài)空間形式的POGO振動模型。該方法廣泛地應(yīng)用于宇宙神二/半人馬座和長征二F的POGO問題研究中。Wang等[13]在Oppenheim等的基礎(chǔ)上提出了一種改進(jìn)的模型，只選取節(jié)點(diǎn)上的相對脈動質(zhì)量位移為變量，得出了維數(shù)減半的微分方程形式的POGO振動模型，不僅降低了計算量，而且可以用于時域仿真。張青松等[14]提出了基于鍵合圖理論的POGO建模方法，建立了包含氣路單元的POGO模型，通過特征值分析POGO振動的穩(wěn)定性。趙治華[15]利用多體系統(tǒng)動力學(xué)方法，分析了包含結(jié)構(gòu)縱橫扭模型的液體火箭POGO振動穩(wěn)定性。

由于燃料輸送管較短，燃料輸送系統(tǒng)的頻率遠(yuǎn)高于結(jié)構(gòu)系統(tǒng)的低階頻率。因此，將推進(jìn)系統(tǒng)簡化為單推進(jìn)劑-單發(fā)動機(jī)模型，不僅可降低模型維數(shù)，且更便于研究POGO振動產(chǎn)生的機(jī)理及開展參數(shù)靈敏度分析。然而目前對單推進(jìn)劑-單發(fā)動機(jī)推進(jìn)系統(tǒng)模型的研究通常將輸送管的方程簡化為一個物理方程，或者直接忽略管路的柔性，作為不可壓縮管路建模，只保留蓄壓器的柔性和泵的氣蝕柔性。該方法僅能體現(xiàn)出推進(jìn)系統(tǒng)一階頻率與結(jié)構(gòu)系統(tǒng)的耦合，無法體現(xiàn)推進(jìn)系統(tǒng)高階頻率與結(jié)構(gòu)系統(tǒng)的耦合。然而隨著液體火箭的規(guī)模增大，液體輸送管路長度不斷增加，氧化劑輸送管路的柔性體降低了推進(jìn)系統(tǒng)二階等低階頻率，且火箭規(guī)模的增大也導(dǎo)致結(jié)構(gòu)系統(tǒng)的低頻密集分布，推進(jìn)系統(tǒng)的二階及以上模態(tài)仍可能與結(jié)構(gòu)系統(tǒng)不同模態(tài)產(chǎn)生耦合導(dǎo)致不穩(wěn)定。然而由于泵、三通和蓄壓器部件的存在，使得文獻(xiàn)[13]中推進(jìn)系統(tǒng)模型系數(shù)質(zhì)量矩陣雖為非奇異陣，但是其剛度陣為非對稱陣，阻尼陣為非比例阻尼或非瑞利阻尼陣，因此無法以傳統(tǒng)的振型分解法將推進(jìn)系統(tǒng)模型振型分解。進(jìn)一步分析文獻(xiàn)[13]中的模型可得出，輸送管路單元的變量占了模型總變量的大多數(shù)。而且輸送管路的不同單元具有相同的表現(xiàn)形式，因此可以將管路模型的質(zhì)量陣、剛度陣和阻尼陣解耦。

為此，本文提出了一種推進(jìn)系統(tǒng)模型縮聚的方法，以貯箱出口和蓄壓器入口作為管路的邊界條件，利用廣義逆迭代法將所有輸送管路的單元的質(zhì)量陣、剛度陣和阻尼陣模態(tài)分解。選取重點(diǎn)關(guān)注階次的推進(jìn)系統(tǒng)管路的模態(tài)方程，與其余部件三通、蓄壓器、泵和推力室組成縮聚后推進(jìn)系統(tǒng)的模型。該模型既具有單推進(jìn)劑-單發(fā)動機(jī)模型的維數(shù)低的優(yōu)點(diǎn)，又包括推進(jìn)系統(tǒng)不同階次模態(tài)。因此可以方便地用于研究推進(jìn)系統(tǒng)與結(jié)構(gòu)系統(tǒng)不同階次模態(tài)之間的耦合穩(wěn)定性。

1 推進(jìn)系統(tǒng)的模型

隨著推進(jìn)劑的消耗，結(jié)構(gòu)系統(tǒng)的頻率增大，推進(jìn)系統(tǒng)的一階頻率往往與結(jié)構(gòu)系統(tǒng)的某一階頻率存在交叉區(qū)，存在POGO失穩(wěn)的風(fēng)險。由于燃料輸送系統(tǒng)的頻率遠(yuǎn)高于結(jié)構(gòu)系統(tǒng)的頻率，因此本文建立的推進(jìn)系統(tǒng)中只考慮氧化劑輸送系統(tǒng)，但建模方法同樣適用于燃料輸送系統(tǒng)。氧化劑輸送系統(tǒng)的示意圖如圖1所示，輸送管路劃分為n段，根據(jù)文獻(xiàn)[12]，每段管路單元的最大長度L滿足

(1)

式中：a為有效當(dāng)?shù)芈曀伲琭u是關(guān)心的最高頻率。

圖1 推進(jìn)系統(tǒng)示意圖

文獻(xiàn)[12]中給出了每種單元動力學(xué)方程的詳細(xì)推導(dǎo)過程，在此不再贅述。圖1所示模型中六類單元的方程形式如下所述。貯箱出口單元的方程為

(2)

氧化劑輸送管路劃分為n段，每一段的方程為

(3)

(4)

式中：Mp、Rp和Kp分別為質(zhì)量陣、阻尼陣和剛度陣，由所有管路單元的剛度項、阻尼項和剛度項組成；U矩陣由各個管路單元方程中結(jié)構(gòu)系統(tǒng)對推進(jìn)系統(tǒng)的作用項組成；p1和ps分別為管路入口節(jié)點(diǎn)和出口節(jié)點(diǎn)的脈動壓力。

三通是連接管路末端出口、蓄壓器和泵入口的元件，其慣性和阻尼特性在與其連接的管路單元中考慮，其動力學(xué)方程為

ps=pa=pp

ws=wa+wp

(5)

式中：pa和pp分別是蓄壓器入口和泵入口處的脈動壓力；ws、wa和wp分別是管路末端出口、蓄壓器入口和泵入口處相對脈動質(zhì)量位移。

蓄壓器作為調(diào)節(jié)推進(jìn)系統(tǒng)頻率特性的主要元件，其動力學(xué)方程為

(6)

式中：Ia、Ra和Ka分別為蓄壓器的慣性、阻尼和剛度。當(dāng)只將蓄壓器作為純?nèi)嵝栽r，方程中只考慮蓄壓器的剛度項，其方程為

pa=Kawa

(7)

泵單元的方程為

pp=Kp(wp-wc)

(8)

式中：Ip、Rp和Kp分別是泵的慣性、阻尼和剛度，m+1為泵的增益。將泵后管的慣性和阻尼添加到泵的慣性和阻尼項中，泵出口的脈動壓力pc等同于推力室中的脈動壓力。

推力室的方程為

(9)

式中：Rc為推力室的阻尼。

1.1 忽略蓄壓器慣性和阻尼時的推進(jìn)系統(tǒng)模型

因每個管路單元的方程形式相同，所以Mp和Kp為對稱陣，Rp為瑞利阻尼陣，因此可以將管路單元的模型進(jìn)行嚴(yán)格解耦。廣義逆迭代法[16-17]，避免了降階擴(kuò)維的復(fù)雜運(yùn)算，僅通過簡單的迭代即可將管路系數(shù)矩陣模態(tài)分解，得到管路的前m階模態(tài)方程。

將蓄壓器作為純?nèi)嵝栽r，綜合式(5)～式(8)，可得

pa=K*(ws-wc)

(10)

式中：K*=KaKp/(Ka+Kp)為等效剛度。根據(jù)式(5)，將ps的表達(dá)式代入式(4)中，修正Kp矩陣的最后一個元素。以(Mp,Rp,Kp)得出的振型矩陣，對式(4)模態(tài)分解，令w=Φpz，得到

(11)

(m+1)K*Φp(n)z=0

(12)

式(11)和(12)即推進(jìn)系統(tǒng)的縮聚模型，可以看出，當(dāng)只取管路的某一階模態(tài)時，推進(jìn)系統(tǒng)僅由兩個方程組成。

1.2 考慮蓄壓器慣性和阻尼的推進(jìn)系統(tǒng)模型

當(dāng)考慮蓄壓器慣性和阻尼時，三通出口端變量wp是一個獨(dú)立變量，將式(6)代入式(4)中，令w=Φpz，對式(4)模態(tài)分解得出

(13)

綜合式(6)和式(8)，可得

Ka(Φp(n)z-wp)-Kp(wp-wc)=0

(14)

(15)

2 結(jié)構(gòu)系統(tǒng)的模型

推進(jìn)系統(tǒng)對結(jié)構(gòu)的作用力包括管路單元、泵和推力室脈動壓力對結(jié)構(gòu)的作用力，因此在推進(jìn)系統(tǒng)作用力下的結(jié)構(gòu)系統(tǒng)的模態(tài)方程為

(16)

式中：Ms、Zs和Ωs分別為結(jié)構(gòu)系統(tǒng)的模態(tài)質(zhì)量陣、阻尼陣和剛度陣，q為結(jié)構(gòu)系統(tǒng)的模態(tài)坐標(biāo)，Φs為結(jié)構(gòu)系統(tǒng)的振型矩陣。fd、fp和fc分別為管路單元、泵單元和推力室單元給結(jié)構(gòu)的作用力。

每一段管路單元對結(jié)構(gòu)系統(tǒng)的作用力為

(17)

式中：Ai和Aj分別為管單元入口和出口的截面積，Ni和Nj分別為管路單元入口和出口處的方向向量，g為重力加速度。

泵中流體對結(jié)構(gòu)作用力為

fp=KpAp(wp-wc)N

(18)

式中：Ap為泵入口的截面積，N為泵入口處的方向向量。

推力室對結(jié)構(gòu)的作用力為

fc=-NAtCfpc

(19)

式中：At為推力室喉部的截面積，Cf為推力室的推力系數(shù)。

將式(17)～(19)及推進(jìn)系統(tǒng)節(jié)點(diǎn)上的脈動壓力表達(dá)式代入式(16)中，可以得出將蓄壓器作為純?nèi)嵝栽r結(jié)構(gòu)系統(tǒng)的方程為

(20)

式中：方程最右端項的出現(xiàn)是因為將相對脈動質(zhì)量位移替換作用力表達(dá)式中脈動壓力時得出的項，體現(xiàn)了推進(jìn)系統(tǒng)與結(jié)構(gòu)系統(tǒng)耦合作用對結(jié)構(gòu)慣性項的影響，是引起結(jié)構(gòu)耦合頻率變化的主要因素。特別地，將式(20)右端項中的K*替換為Kp即為考慮蓄壓器慣性和阻尼時的結(jié)構(gòu)系統(tǒng)方程。

3 POGO振動模型

3.1 忽略蓄壓器慣性和阻尼時POGO振動模型

將蓄壓器當(dāng)作純?nèi)嵝栽r，以x=[z,wc,q]為狀態(tài)變量，綜合式(11)，(12)和式(20)，可得到縮聚后的POGO振動模型為

(21)

式中：

(22)

式(21)即模態(tài)縮聚后的POGO振動模型，可以看出，耦合模型的變量由管路單元的模態(tài)位移、推力室的相對脈動質(zhì)量位移和結(jié)構(gòu)的模態(tài)位移組成。耦合模型的維數(shù)N與管路模態(tài)變量個數(shù)N(z)、結(jié)構(gòu)模態(tài)位移的個數(shù)N(q)的關(guān)系滿足N=N(z)+N(q)+1。其中，N(z)和N(q)可以根據(jù)實際情況選擇不同的維數(shù)。特別地，當(dāng)N(z)與和N(q)分別取1時，耦合系統(tǒng)模型的維數(shù)為3，與改進(jìn)的Rubin模型相比，模型的維數(shù)和計算量均大幅降低，且M矩陣、K矩陣也為非奇異矩陣，可以直接通過QR法計算耦合模型的特征值以判斷系統(tǒng)穩(wěn)定性。利用該縮模型，可以更有針對性地研究推進(jìn)系統(tǒng)重要參數(shù)對POGO穩(wěn)定性的影響。

3.2 考慮蓄壓器慣性和阻尼時的POGO振動模型

綜合式(13)～(15)以及式(20)(Kp替換K*時的結(jié)構(gòu)系統(tǒng)方程)，以x=[z,wp,wc,q]為狀態(tài)量，即得考慮蓄壓器慣性和阻尼項時的與式具有相同形式的POGO振動模型，其中M、R和K分別為

(23)

可以看出，M陣仍然為非奇異陣，也可以直接用QR法求解系統(tǒng)的特征值。模型的維度比式(22)的模型多了一個維度，即N=N(z)+N(q)+2。

4 算 例

用本文推導(dǎo)出的縮聚模型與改進(jìn)的Rubin模型分別計算了某型號液體火箭在某個飛行時間段POGO穩(wěn)定性，其中橫坐標(biāo)除以最大飛行時間進(jìn)行了無量綱化，蓄壓器設(shè)計為純?nèi)嵝栽?，采用式所示的模型計算。圖2和圖3給出了不同飛行秒點(diǎn)結(jié)構(gòu)系統(tǒng)與推進(jìn)系統(tǒng)耦合頻率與阻尼比的計算曲線。可以看出，隨著飛行時間的增大，結(jié)構(gòu)系統(tǒng)與推進(jìn)系統(tǒng)的頻率在0.9 s(s僅代表無量綱時刻)附近出現(xiàn)了接近和穿越，從圖3中可以看出結(jié)構(gòu)耦合阻尼比在0.9 s附近大幅降低，表明結(jié)構(gòu)系統(tǒng)的穩(wěn)定性降低，而推進(jìn)系統(tǒng)耦合阻尼比大幅增大。以改進(jìn)的Rubin模型的計算結(jié)果為基準(zhǔn)值，縮聚模型的耦合頻率與阻尼比的最大相對偏差如表1所示，結(jié)構(gòu)系統(tǒng)耦合頻率和阻尼比的最大相對偏差為0.032%和1.9%，推進(jìn)系統(tǒng)耦合頻率和阻尼比的最大相對偏差為0.55%和0.20%，因此本文推導(dǎo)出POGO振動縮聚模型不僅維數(shù)大幅降低，而且具有相當(dāng)高的精度。

圖2 不同飛行秒點(diǎn)結(jié)構(gòu)系統(tǒng)與推進(jìn)系統(tǒng)耦合頻率

圖3 不同飛行秒點(diǎn)結(jié)構(gòu)系統(tǒng)與推進(jìn)系統(tǒng)耦合阻尼比

耦合計算結(jié)果最大相對偏差結(jié)構(gòu)系統(tǒng)頻率0.032%阻尼比1.9%推進(jìn)系統(tǒng)頻率0.55%阻尼比0.20%

蓄壓器能量值對不同秒點(diǎn)模型耦合穩(wěn)定性的影響如圖4和圖5所示，其中蓄壓器能量值除以設(shè)計工況的額定值進(jìn)行了無量綱化。由圖4可以看出，蓄壓器的能量值對不同秒點(diǎn)結(jié)構(gòu)系統(tǒng)耦合阻尼比影響是不同的。隨著蓄壓器能量值增大，0.7 s結(jié)構(gòu)系統(tǒng)耦合阻尼比變化很小，略有增大，0.93 s耦合模型的結(jié)構(gòu)系統(tǒng)耦合阻尼比先大幅減小又大幅增大。說明0.93 s耦合模型的穩(wěn)定性對蓄壓器的靈敏度較高。與圖3對比可以看出，0.7 s耦合模型的結(jié)構(gòu)系統(tǒng)與推進(jìn)系統(tǒng)的頻率差別較大，0.93 s耦合模型的結(jié)構(gòu)系統(tǒng)與推進(jìn)系統(tǒng)頻率差別較小，兩者頻率的接近造成了耦合作用。

圖4 蓄壓器能量值對不同秒點(diǎn)結(jié)構(gòu)系統(tǒng)耦合阻尼比的影響

圖5給出了0.7 s和0.9 s推進(jìn)系統(tǒng)的阻尼比隨蓄壓器能量值增大的變化情況，可以看出，隨著蓄壓器能量值的增大，兩個秒點(diǎn)的耦合模型中的推進(jìn)系統(tǒng)阻尼比都大幅降低。但0.9 s耦合模型推進(jìn)系統(tǒng)的阻尼比在0.6 s～1.5 s區(qū)間有一向上的“鼓包”，通過對比圖4可以得出是推進(jìn)系統(tǒng)與結(jié)構(gòu)系統(tǒng)耦合作用的結(jié)果。

圖5 蓄壓器能量值對不同秒點(diǎn)推進(jìn)系統(tǒng)耦合阻尼比的影響

當(dāng)考慮蓄壓器的慣性和阻尼時，可以用式(23)所示的模型研究蓄壓器慣性和阻尼等參數(shù)對POGO穩(wěn)定性的影響。圖6給出了蓄壓器慣性對0.7 s和0.9 s耦合模型中結(jié)構(gòu)系統(tǒng)耦合阻尼比的影響曲線?？梢钥闯觯顗浩鞯膽T性對0.9 s模型中的結(jié)構(gòu)耦合阻尼比影響較大，隨著蓄壓器的慣性的增大，且這種影響也是非單調(diào)的。

圖6 蓄壓器慣性對不同秒點(diǎn)結(jié)構(gòu)系統(tǒng)耦合阻尼比的影響

5 結(jié) 論

本文提出了建立縮聚POGO振動模型的方法，分別推導(dǎo)了將蓄壓器作為純?nèi)嵝栽涂紤]蓄壓器慣性和阻尼時的POGO振動縮聚模型。與改進(jìn)的Rubin模型相比，縮聚模型的維度大幅降低，從而大幅降低計算量。算例表明縮聚后的POGO振動模型具有很高的精度?；赑OGO縮聚模型分別計算了某型號液體火箭蓄壓器能量值、慣性對0.7 s和0.93 s耦合模型POGO穩(wěn)定性的影響。研究得出，蓄壓器的能量值和慣性對不同時間點(diǎn)結(jié)構(gòu)耦合阻尼比具有不同的影響，且這種影響是非單調(diào)的。當(dāng)結(jié)構(gòu)系統(tǒng)與推進(jìn)系統(tǒng)頻率較接近而導(dǎo)致結(jié)構(gòu)耦合阻尼比大幅降低時，若蓄壓器設(shè)計為純?nèi)嵝栽侠碚{(diào)節(jié)蓄壓器的能量值可以增大結(jié)構(gòu)耦合阻尼比，當(dāng)考慮蓄壓器的慣性和阻尼時，合理的增大蓄壓器的慣性也可以增大結(jié)構(gòu)系統(tǒng)的耦合阻尼比，增強(qiáng)結(jié)構(gòu)系統(tǒng)的穩(wěn)定性。