液體火箭POGO振動縮聚模型研究

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

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

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

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

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

1 推進系統的模型

隨著推進劑的消耗，結構系統的頻率增大，推進系統的一階頻率往往與結構系統的某一階頻率存在交叉區，存在POGO失穩的風險。由于燃料輸送系統的頻率遠高于結構系統的頻率，因此本文建立的推進系統中只考慮氧化劑輸送系統，但建模方法同樣適用于燃料輸送系統。氧化劑輸送系統的示意圖如圖1所示，輸送管路劃分為n段，根據文獻[12]，每段管路單元的最大長度L滿足

(1)

式中：a為有效當地聲速，fu是關心的最高頻率。

圖1 推進系統示意圖

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

(2)

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

(3)

(4)

式中：Mp、Rp和Kp分別為質量陣、阻尼陣和剛度陣，由所有管路單元的剛度項、阻尼項和剛度項組成；U矩陣由各個管路單元方程中結構系統對推進系統的作用項組成；p1和ps分別為管路入口節點和出口節點的脈動壓力。

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

ps=pa=pp

ws=wa+wp

(5)

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

蓄壓器作為調節推進系統頻率特性的主要元件，其動力學方程為

(6)

式中：Ia、Ra和Ka分別為蓄壓器的慣性、阻尼和剛度。當只將蓄壓器作為純柔性元件時，方程中只考慮蓄壓器的剛度項，其方程為

pa=Kawa

(7)

泵單元的方程為

pp=Kp(wp-wc)

(8)

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

推力室的方程為

(9)

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

1.1 忽略蓄壓器慣性和阻尼時的推進系統模型

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

將蓄壓器作為純柔性元件時，綜合式(5)～式(8)，可得

pa=K*(ws-wc)

(10)

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

(11)

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

(12)

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

1.2 考慮蓄壓器慣性和阻尼的推進系統模型

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

(13)

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

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

(14)

(15)

2 結構系統的模型

推進系統對結構的作用力包括管路單元、泵和推力室脈動壓力對結構的作用力，因此在推進系統作用力下的結構系統的模態方程為

(16)

式中：Ms、Zs和Ωs分別為結構系統的模態質量陣、阻尼陣和剛度陣，q為結構系統的模態坐標，Φs為結構系統的振型矩陣。fd、fp和fc分別為管路單元、泵單元和推力室單元給結構的作用力。

每一段管路單元對結構系統的作用力為

(17)

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

泵中流體對結構作用力為

fp=KpAp(wp-wc)N

(18)

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

推力室對結構的作用力為

fc=-NAtCfpc

(19)

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

將式(17)～(19)及推進系統節點上的脈動壓力表達式代入式(16)中，可以得出將蓄壓器作為純柔性元件時結構系統的方程為

(20)

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

3 POGO振動模型

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

將蓄壓器當作純柔性元件時，以x=[z,wc,q]為狀態變量，綜合式(11)，(12)和式(20)，可得到縮聚后的POGO振動模型為

(21)

式中：

(22)

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

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

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

(23)

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

4 算 例

用本文推導出的縮聚模型與改進的Rubin模型分別計算了某型號液體火箭在某個飛行時間段POGO穩定性，其中橫坐標除以最大飛行時間進行了無量綱化，蓄壓器設計為純柔性原件，采用式所示的模型計算。圖2和圖3給出了不同飛行秒點結構系統與推進系統耦合頻率與阻尼比的計算曲線?？梢钥闯?，隨著飛行時間的增大，結構系統與推進系統的頻率在0.9 s(s僅代表無量綱時刻)附近出現了接近和穿越，從圖3中可以看出結構耦合阻尼比在0.9 s附近大幅降低，表明結構系統的穩定性降低，而推進系統耦合阻尼比大幅增大。以改進的Rubin模型的計算結果為基準值，縮聚模型的耦合頻率與阻尼比的最大相對偏差如表1所示，結構系統耦合頻率和阻尼比的最大相對偏差為0.032%和1.9%，推進系統耦合頻率和阻尼比的最大相對偏差為0.55%和0.20%，因此本文推導出POGO振動縮聚模型不僅維數大幅降低，而且具有相當高的精度。

圖2 不同飛行秒點結構系統與推進系統耦合頻率

圖3 不同飛行秒點結構系統與推進系統耦合阻尼比

耦合計算結果最大相對偏差結構系統頻率0.032%阻尼比1.9%推進系統頻率0.55%阻尼比0.20%

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

圖4 蓄壓器能量值對不同秒點結構系統耦合阻尼比的影響

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

圖5 蓄壓器能量值對不同秒點推進系統耦合阻尼比的影響

當考慮蓄壓器的慣性和阻尼時，可以用式(23)所示的模型研究蓄壓器慣性和阻尼等參數對POGO穩定性的影響。圖6給出了蓄壓器慣性對0.7 s和0.9 s耦合模型中結構系統耦合阻尼比的影響曲線。可以看出，蓄壓器的慣性對0.9 s模型中的結構耦合阻尼比影響較大，隨著蓄壓器的慣性的增大，且這種影響也是非單調的。

圖6 蓄壓器慣性對不同秒點結構系統耦合阻尼比的影響

5 結 論

本文提出了建立縮聚POGO振動模型的方法，分別推導了將蓄壓器作為純柔性元件和考慮蓄壓器慣性和阻尼時的POGO振動縮聚模型。與改進的Rubin模型相比，縮聚模型的維度大幅降低，從而大幅降低計算量。算例表明縮聚后的POGO振動模型具有很高的精度?；赑OGO縮聚模型分別計算了某型號液體火箭蓄壓器能量值、慣性對0.7 s和0.93 s耦合模型POGO穩定性的影響。研究得出，蓄壓器的能量值和慣性對不同時間點結構耦合阻尼比具有不同的影響，且這種影響是非單調的。當結構系統與推進系統頻率較接近而導致結構耦合阻尼比大幅降低時，若蓄壓器設計為純柔性元件，合理調節蓄壓器的能量值可以增大結構耦合阻尼比，當考慮蓄壓器的慣性和阻尼時，合理的增大蓄壓器的慣性也可以增大結構系統的耦合阻尼比，增強結構系統的穩定性。