韋 笑，楊瓊梁，張美艷，唐國(guó)安

（1.復(fù)旦大學(xué) 力學(xué)與工程科學(xué)系，上海 200433；2.上海宇航系統(tǒng)工程研究所，上海 201109）

0 引言

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

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

1 貯箱和推進(jìn)劑縱向振動(dòng)的彈簧振子模型

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

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

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

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

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

由式（3）、（4），可確定振子模型中的二個(gè)參數(shù)

振子模型中3個(gè)彈簧串－并聯(lián)后的剛度就是貯箱貯液段的縱向剛度K，則第三個(gè)剛度參數(shù)

在微小動(dòng)態(tài)擾動(dòng)前提下，無論結(jié)構(gòu)是否均質(zhì)或各向異性，位移都可用柔度系數(shù)表示為外載荷的線性函數(shù)

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

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

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

與式（6）類似，第三個(gè)彈簧的剛度系數(shù)

對(duì)柔度系數(shù)，可用通用程序用有限元方法計(jì)算。設(shè)計(jì)兩組計(jì)算工況：

則，可得式（7）中的柔度系數(shù)為

2 液體質(zhì)心在貯箱變形后位置計(jì)算

建立貯箱貯液部分的有限元模型如圖2（a）所示。此處僅示意了無加強(qiáng)筋的箱體結(jié)構(gòu)，但由有限元方法的特點(diǎn)可知，對(duì)更復(fù)雜的箱體結(jié)構(gòu)變形也能進(jìn)行建模計(jì)算。有限元計(jì)算能給出模型中每個(gè)節(jié)點(diǎn)的位移分量，且可認(rèn)為模型中每個(gè)單元在變形后的位置是確定的。

定義液體的質(zhì)心坐標(biāo)為

其中積分區(qū)域?yàn)樽冃魏笠后w所占據(jù)部分的三維空間。

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

根據(jù)高斯公式，有

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

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

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

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

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

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

3 運(yùn)載火箭縱向振動(dòng)的貯箱建模應(yīng)用

為驗(yàn)證本文方法，以橢球底柱形貯箱為例進(jìn)行數(shù)值計(jì)算。該貯箱貯液部分表面分布節(jié)點(diǎn)6 041個(gè)，共被劃分為6 000個(gè)有限元網(wǎng)格。取貯箱數(shù)據(jù)為：筒段內(nèi)貯液高度h＝8m；筒段半徑R＝1.675m；后底短軸r＝1.047m；短長(zhǎng)軸之比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。算例中貯箱內(nèi)加注的液體為水。計(jì)算柔度時(shí)，取Mp＝76 575.6kg，x··p＝9.8m／s2；F0＝750 441N。有限元計(jì)算程序選用MSC.Nastran。用本文有限元計(jì)算方法計(jì)算各項(xiàng)數(shù)值結(jié)果為：c11＝2.435 6×10－9N／m，c21＝7.346 9×10－10N／m，c22＝4.759 6×10－9N／m。代入式（9）、（10）所得振子模型的3個(gè)縱向等效剛度系數(shù)見表1，與用文獻(xiàn)［2］無矩殼方法所得相應(yīng)的剛度系數(shù)比較，兩種方法的差異見表1。

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

4 結(jié)束語

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

［1］ 周思達(dá)，劉 莉.運(yùn)載火箭貯箱流固耦合分析方法綜述［J］.強(qiáng)度與環(huán)境，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］ 王其政.結(jié)構(gòu)耦合動(dòng)力學(xué)［M］.北京：中國(guó)宇航出版社，1999.

［7］ 付寶連.關(guān)于功的互等定理與疊加原理的等價(jià)性［J］.應(yīng)用數(shù)學(xué)和力學(xué)，1985，6（9）：813－818.

［8］ 李菊娥.應(yīng)用高斯公式應(yīng)注意的一個(gè)問題［J］.高等數(shù)學(xué)研究，2000（1）：4.

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

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