確定貯箱流固耦合縱向振子模型參數的有限元方法

韋 笑，楊瓊梁，張美艷，唐國安

（1.復旦大學 力學與工程科學系，上海 200433；2.上海宇航系統工程研究所，上海 201109）

0 引言

液體運載火箭的結構振動與推力脈動存在耦合，這種耦合可能使運載火箭在主動飛行段產生較大的動態響應，甚至使其發生縱向不穩定振動。因此，對運載火箭進行動力穩定性的分析對火箭飛行安全有重要意義。分析穩定性時，需建立貯箱及推進劑的縱向振動模型，其中彈簧振子模型因其簡便、高效而在工程中被廣為采用［1］。1961年，文獻［2］提出了貯箱液固耦合的縱向質量－彈簧模型。文獻［3］對Wood模型中的貯箱箱底模型進行了進一步研究，修正了貯箱和推進劑彈簧振子模型中的3個剛度系數。文獻［4］將貯箱筒壁及箱底置于同一分析模型，用實驗進行了對比分析。文獻［5］用經典板殼理論中的無矩板理論得到貯箱筒壁的動力學方程，對Wood模型中貯箱筒段模型部分進行了細化與改進。文獻［6］提出了一種貯箱流體多自由度的質量－彈簧模型。貯箱和推進劑的彈簧振子建模方法及剛度系數計算公式至今仍被應用，但為能用殼體理論推導出剛度系數的計算公式，必須對貯箱的變形和受力引入多種假設，如筒段和箱底被認為是等厚度的各向同性均質薄殼，且殼體內只有面內的張力和剪力，而無彎矩和扭矩。對實際工程中通過加筋等工藝提高強度和剛度的薄殼容器，必須經壁厚等效等處理后才能套用經典的殼體理論。此外，新型復合材料殼體具明顯的各向異性，這也將增大經典殼體理論分析的難度。

有限元方法是一種通用的方法，不受各向異性、非均質／變厚度殼體等限制。用有限元方法計算貯箱在液體壓力及上面級箭體載荷作用下的箱體變形和液體質心位置變化可獲得較高精度。有限元與縱向振動彈簧振子的建模方法結合，能克服經典殼體理論的局限，適應更復雜箱體結構的貯箱和推進劑建模。本文用有限元方法對貯箱流固耦合縱向振子模型參數的確定進行了研究。

1 貯箱和推進劑縱向振動的彈簧振子模型

液體火箭推進劑貯箱常為圓柱形箱體，且其結構部分一般由前底、筒段、后底三部分組成。箱底的底形可分為半球形底、半橢球形底、錐形底和三心底等，筒段則多為圓柱殼。半橢球底柱形貯箱如圖1（a）所示。在縱向振動簡化分析時，箱內液體以及與液體接觸部分的殼體（貯液段）可等效為由彈簧和質量組成的振子模型，如圖1（b）所示。通常，振子模型包含1個質量參數Mp和3個剛度參數k1，k2，k3，其中Mp為貯箱內液體質量。

圖1 貯箱箱體與貯液縱向耦合振動的振子模型Fig.1 Structure and deformation of tank

為確定振子k1，k2，k3，考慮以勻加速度上升的貯箱段箱體，在液體壓力p（x）和上面級箭體載荷F的作用下，貯液段箱體將發生變形，設筒段上下兩端的位移分別為x1，x0，并約定位移方向向上為正。箱內液體的質心也將發生改變，位移量為xp。用振子等效該段箱體和箱內液體縱向耦合振動模型的原則為：當圖1（b）所示的振子上端受F、質量塊受慣性力作用時，振子的上下端點及質量塊分別出現與x1，x0，xp相同的位移量。據此原則，確定振子剛度參數的關鍵是計算在F，p（x）作用下貯液段箱體的變形。

用式（1）、（2）消去變量F，并由關系x1－x0＝（x1－xp）＋（xp－x0）可得液體質心位移控制方程

對圖1（b）的彈簧－質量振子，易得液體質心位移控制方程為

由式（3）、（4），可確定振子模型中的二個參數

振子模型中3個彈簧串－并聯后的剛度就是貯箱貯液段的縱向剛度K，則第三個剛度參數

在微小動態擾動前提下，無論結構是否均質或各向異性，位移都可用柔度系數表示為外載荷的線性函數

2.4 影響新生兒黃疸嚴重程度相關危險因素的Logistic回歸分析 以新生兒黃疸病情輕重程度作為因變量，將胎齡、胎兒出生時體質量、頭顱血腫、喂養方式、開奶時間、母嬰血型不合6因素作為回歸性分析的自變量，建立Logistic回歸模型。結果顯示:開奶時間、母嬰血型不合、喂養方式是新生兒黃疸病情嚴重程度的主要危險因素(P<0.05)。見表3。

根據位移互等定理和柔度矩陣的正定性，有c21＝c12＞0且c11c22＞c12c21［7］。從式（7）中消去變量F，可得

比較式（4）、（8），可確定振子模型中的第一、二個彈簧的剛度系數

與式（6）類似，第三個彈簧的剛度系數

對柔度系數，可用通用程序用有限元方法計算。設計兩組計算工況：

則，可得式（7）中的柔度系數為

2 液體質心在貯箱變形后位置計算

建立貯箱貯液部分的有限元模型如圖2（a）所示。此處僅示意了無加強筋的箱體結構，但由有限元方法的特點可知，對更復雜的箱體結構變形也能進行建模計算。有限元計算能給出模型中每個節點的位移分量，且可認為模型中每個單元在變形后的位置是確定的。

定義液體的質心坐標為

其中積分區域為變形后液體所占據部分的三維空間。

圖2 貯箱有限元模型Fig.2 Finite element mesh of propellant tank

根據高斯公式，有

式中：?Ω為貯箱內液體的整個表面，包括自由面以及與箱壁的接觸面；Γ為對應的表面積；l為液體表面外法線方向余弦［lmn］的x軸分量；xf為變形后液體自由面高度［8］。

在自由面上x＝xf，故只須在液體與箱體接觸面上積分，且該積分可通過對每個四邊形殼體單元進行數值積分后求和而得，即

定義式（12）中的分母也可用相同方式計算積分

變形后液體的自由面高度xf與箱體變形有關，須根據液體的不可壓縮性確定。將積分式（15）理解為平面x＝xf與箱體曲面圍成的三維單連通區域的體積，該體積為變量xf的函數。因貯箱內推進劑體積恒為Mp／ρp，故可求解非線性方程

確定自由面位置。此處：ρp為貯箱內推進劑的密度。式（16）可用二分法求得其數值解［9］。

箱體變形、液面下降后，液體自由面通常不再與箱體有限元模型的單元邊界重合，如圖2（b）所示。為此，可用計算機圖形學中的多邊形裁剪算法，先對自由面附近的箱體單元進行裁剪，再對式（14）、（15）進行積分［10］。

3 運載火箭縱向振動的貯箱建模應用

為驗證本文方法，以橢球底柱形貯箱為例進行數值計算。該貯箱貯液部分表面分布節點6 041個，共被劃分為6 000個有限元網格。取貯箱數據為：筒段內貯液高度h＝8m；筒段半徑R＝1.675m；后底短軸r＝1.047m；短長軸之比n＝0.625；貯箱筒段厚度t＝4.59mm；后底厚度tBH＝2.00mm；E＝6.8×1010Pa；v＝0.33；貯箱材料密度ρ＝2.8×103kg／m3；貯液密度ρ0＝1.0×103kg／m3。算例中貯箱內加注的液體為水。計算柔度時，取Mp＝76 575.6kg，x··p＝9.8m／s2；F0＝750 441N。有限元計算程序選用MSC.Nastran。用本文有限元計算方法計算各項數值結果為：c11＝2.435 6×10－9N／m，c21＝7.346 9×10－10N／m，c22＝4.759 6×10－9N／m。代入式（9）、（10）所得振子模型的3個縱向等效剛度系數見表1，與用文獻［2］無矩殼方法所得相應的剛度系數比較，兩種方法的差異見表1。

表1 兩種方法算得剛度系數Tab.1 Stiffness of two kinds methods

4 結束語

本文提出了一種用有限元標定確定液體火箭推進劑貯箱振動簡化模型縱向剛度的方法。通過模擬解析解法建立等效力學模型的過程，設定僅貯液和貯液且施壓兩種工況，根據貯箱形變和液體質心位置變化性質，推導出剛度系數的計算公式，數值計算結果與傳統無矩殼方法基本符合。該方法不僅適于外形規則的貯箱，也能用于難以計算解析解的異形貯箱及貯箱表面加固有筋條等工況，方法的通用性較傳統幾何計算方法等有明顯提高。該方法的分析計算結果可為充液航天器的燃料貯箱結構設計提供理論依據。

［1］ 周思達，劉 莉.運載火箭貯箱流固耦合分析方法綜述［J］.強度與環境，2010，37（3）：52－63.

［2］ WOOD J D.Survey on missile structural dynamics［R］.TRW Space Technology Labs，STL－7102－0041－NU－000，1961.

［3］ PINSON L D.Longitudinal spring constants for liquid－propellant tanks with ellipsoidal ends［M］.Washington DC：National Aeronautics and Space Administration，1964.

［4］ GORMLEY J F，KANA D.Longitudinal vibration of a model space vehicle propellant tank［J］.Journal of Spacecraft and Rockets，1967，4（12）：1585－1591.

［5］ EULITZ W R，GLASER R F，KANA D，et al.Longitudinal vibration of spring－supported cylindrical membrane shells containing liquid［J］.Journal of Spacecraft and Rockets，1968，5（2）：189－196.

［6］ 王其政.結構耦合動力學［M］.北京：中國宇航出版社，1999.

［7］ 付寶連.關于功的互等定理與疊加原理的等價性［J］.應用數學和力學，1985，6（9）：813－818.

［8］ 李菊娥.應用高斯公式應注意的一個問題［J］.高等數學研究，2000（1）：4.

［9］ 王海濤，朱 洪.改進的二分法查找［J］.計算機工程，2006，32（10）：60－62.

［10］ 王志強，王世萍.多邊形裁剪算法［J］.西部電子，1994，5（4）：47－52.