張 忠，韓 麗，任 方，李海波，秦朝紅

(北京強(qiáng)度環(huán)境研究所可靠性與環(huán)境工程技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京100076)

1 引言

運(yùn)載火箭全箭動(dòng)力學(xué)特性是姿態(tài)控制系統(tǒng)設(shè)計(jì)、載荷分析、POGO分析等的基礎(chǔ)和依據(jù)，全箭動(dòng)特性建模和分析是必不可少的［1］。運(yùn)載火箭的結(jié)構(gòu)變得越來越復(fù)雜，呈現(xiàn)出空間模態(tài)特征，僅僅依靠傳統(tǒng)的梁模型無法很好的反映其動(dòng)力學(xué)特性，需建立復(fù)雜的三維模型，導(dǎo)致相應(yīng)的計(jì)算成本增加，有時(shí)所建立的有限元模型甚至?xí)猩习偃f個(gè)自由度，如何提高計(jì)算效率成為解決該問題的關(guān)鍵。子結(jié)構(gòu)模態(tài)綜合法是解決大型復(fù)雜結(jié)構(gòu)特征值計(jì)算問題的有效方法，模態(tài)綜合法將整個(gè)系統(tǒng)劃分為多個(gè)子結(jié)構(gòu)，每個(gè)子結(jié)構(gòu)獨(dú)立建模，然后再組配成整體結(jié)構(gòu)，其優(yōu)點(diǎn)是可以縮減系統(tǒng)自由度，提高計(jì)算效率。根據(jù)子結(jié)構(gòu)界面約束條件的不同，可以把子結(jié)構(gòu)位移表達(dá)式分為三類:①自由界面模態(tài)綜合法［2-6］;②約束界面模態(tài)綜合法［7-10］;③混合界面模態(tài)綜合法［3，11］，邱吉寶對(duì)模態(tài)綜合法進(jìn)行了系統(tǒng)的研究，并在文獻(xiàn)［11］中對(duì)動(dòng)態(tài)子結(jié)構(gòu)法進(jìn)行了總結(jié)。

運(yùn)載火箭貯箱內(nèi)裝有的液體燃料，約占運(yùn)載火箭質(zhì)量的80%，液體建模的正確與否直接影響全箭動(dòng)力學(xué)建模的準(zhǔn)確性。目前工程上液體建模多采用集中質(zhì)量方法和虛質(zhì)量方法［12］。集中質(zhì)量方法可應(yīng)用于運(yùn)載火箭梁模型或三維模型，將液體當(dāng)作單純的集中質(zhì)量處理，其缺點(diǎn)是無法反映液體運(yùn)動(dòng)與貯箱的柔性變形之間的耦合關(guān)系;為了解決集中質(zhì)量建模的問題，大型有限元軟件NASTRAN提供了液體虛質(zhì)量建模方法，但工程上應(yīng)用不多，其主要原因是虛質(zhì)量建模會(huì)給出一個(gè)與貯箱自由度規(guī)模相當(dāng)?shù)臐M質(zhì)量陣，導(dǎo)致計(jì)算效率大大降低，甚至無法進(jìn)行大規(guī)模模型的計(jì)算。

為了解決上述問題，本文提出了采用動(dòng)態(tài)子結(jié)構(gòu)和虛質(zhì)量相結(jié)合的方法，在保證計(jì)算精度的同時(shí)降低了求解規(guī)模，提高了計(jì)算效率。

2 基本理論

2.1 虛質(zhì)量方法

虛質(zhì)量建模是建立在液體無旋、無粘、不可壓假定之下的，液體的作用最終轉(zhuǎn)化為附加到貯箱壁上的附加質(zhì)量，這些附加質(zhì)量的計(jì)算并不是按照液體的實(shí)際質(zhì)量，而是按照液體運(yùn)動(dòng)方程，將液體內(nèi)部的運(yùn)動(dòng)方程轉(zhuǎn)化到邊界上，最終形成一個(gè)邊界自由度上的質(zhì)量陣，虛質(zhì)量建模方法能夠克服集中質(zhì)量方法的弊端［12］。

假設(shè)流體是無粘、不可壓縮的理想流體，流場是擬靜態(tài)場，結(jié)構(gòu)是線彈性小變形結(jié)構(gòu)，根據(jù)流體力學(xué)模型，利用Helmholtz方法求解Laplace方程，借助邊界元方法離散技術(shù)進(jìn)行數(shù)值求解，可獲得液體對(duì)結(jié)構(gòu)的集中質(zhì)量。

液體速度勢(shì)及壓力場如公式(1)(2)所示。

式中，vi——任意結(jié)點(diǎn)ri處速度向量;

pi——任意面Aj上的壓力;

Aj——結(jié)構(gòu)體表面上微元面積;

σj——j結(jié)點(diǎn)處的流速向量;

eij——從j點(diǎn)到i點(diǎn)的單位向量

ρ——流體密度

ri——結(jié)點(diǎn)位置

將式(1)及式(2)積分得到公式(3)(4)。

附加質(zhì)量陣如式(5)所示。

獲得附加質(zhì)量后，質(zhì)量矩陣會(huì)變成一個(gè)相對(duì)滿的陣，這給求解特征值問題帶來一定的困難，需采用動(dòng)態(tài)子結(jié)構(gòu)方法進(jìn)行模型縮減。

2.2 動(dòng)態(tài)子結(jié)構(gòu)模態(tài)綜合方法

考慮液體附加質(zhì)量后的結(jié)構(gòu)動(dòng)力學(xué)方程如式(6)所示。

其中，M ＝M0+MA，M0為結(jié)構(gòu)質(zhì)量陣，MA為液體附加質(zhì)量陣，K為剛度陣，F(xiàn)為載荷向量。

將動(dòng)力學(xué)方程按保留自由度和非保留自由度分塊后得式(7)。

式中，下標(biāo)m、s分別代表交界面和非交界面自由度，并且對(duì)于求解系統(tǒng)的特征值問題，除了界面約束力Fm外，沒有其他任何力，即Fs=0。

約束模態(tài)綜合法的子結(jié)構(gòu)坐標(biāo)變換矩陣由邊界約束的子結(jié)構(gòu)主模態(tài)以及約束模態(tài)組成。由模態(tài)展開定理可以得到式(8)。

將式(7)代入式(6)，并前乘TT可以得到以模態(tài)坐標(biāo)描述的剛度矩陣、質(zhì)量矩陣、模態(tài)力向量如式(10)～(12)所示。

式(13)即為包含所需低頻模態(tài)信息的動(dòng)力學(xué)方程。通過以上一系列的求解，有限元模型的自由度大大減少。

可以將廣義模型當(dāng)作物理模型一樣看待完成耦合分析。耦合分析過程是根據(jù)各子結(jié)構(gòu)交界面的位移和力的協(xié)調(diào)條件完成整體耦合結(jié)構(gòu)運(yùn)動(dòng)方程的組裝，然后求解整體運(yùn)動(dòng)方程的動(dòng)力學(xué)問題。

兩個(gè)獨(dú)立子結(jié)構(gòu)a、b的動(dòng)力學(xué)方程根據(jù)式(13)，可以寫成式(14)～(15)。

采用Excel 2013軟件建立數(shù)據(jù)庫，通過SPSS19.0統(tǒng)計(jì)軟件進(jìn)行數(shù)據(jù)處理，檢驗(yàn)水準(zhǔn)取p=0.05，所有數(shù)據(jù)均以±sd表示。采用單因素方差分析(one-way ANOVA)，方差齊時(shí)組間比較采用Duncan法，方差不齊時(shí)組間比較采用DunnettsT3法。p＜0.05為差異有統(tǒng)計(jì)學(xué)意義。

可以將廣義模型當(dāng)作物理模型一樣看待，應(yīng)用式(16)(17)所示的協(xié)調(diào)條件和平衡條件

得到組合體的整體動(dòng)力學(xué)方程為如式(18)所示。

式(18)即為耦合系統(tǒng)的整體動(dòng)力學(xué)方程，求解(18)時(shí)即可得到整體耦合系統(tǒng)的模態(tài)信息。

對(duì)于多級(jí)子結(jié)構(gòu)耦合分析，只需按照以上步驟依次進(jìn)行計(jì)算即可。

3 算例分析

3.1 全箭結(jié)構(gòu)有限元建模與分析

有限元模型主體結(jié)構(gòu)采用梁—?dú)とS模型，推進(jìn)劑采用集中質(zhì)量進(jìn)行模擬，只計(jì)質(zhì)量，不計(jì)轉(zhuǎn)動(dòng)慣量，并通過RBE3與貯箱殼單元連接;助推與芯級(jí)連接采用梁(桿)單元進(jìn)行模擬，并釋放旋轉(zhuǎn)自由度模擬鉸接;發(fā)動(dòng)機(jī)采用等效梁單元模擬，所建立的有限元模型如圖1所示。

圖1 有限元模型Fig.1 FEM model

首先，開展不包含液體推進(jìn)劑的動(dòng)態(tài)子結(jié)構(gòu)模態(tài)綜合計(jì)算分析。對(duì)于運(yùn)載火箭和航天器組成的系統(tǒng)級(jí)模型，通過在星箭界面施加固定界面約束將該系統(tǒng)分成兩部分，分別作為運(yùn)載火箭子結(jié)構(gòu)和航天器子結(jié)構(gòu)，保留航天器和運(yùn)載火箭子結(jié)構(gòu)的低階模態(tài)，運(yùn)用上述動(dòng)態(tài)子結(jié)構(gòu)法進(jìn)行綜合，得到整體結(jié)構(gòu)的計(jì)算結(jié)果，并與整體有限元解進(jìn)行比較，驗(yàn)證其計(jì)算精度。

模態(tài)綜合計(jì)算結(jié)果與整體有限元模型計(jì)算結(jié)果對(duì)比列入表1中。從表1中可以看出，在低頻范圍內(nèi)，具有較高的計(jì)算精度。除了計(jì)算精度外，計(jì)算效率也是工程廣泛關(guān)注的一個(gè)問題。對(duì)于本算例，整體有限元模型具有30多萬個(gè)自由度，而經(jīng)過一系列縮聚后其模態(tài)綜合模型的自由度不足1000個(gè)，計(jì)算效率的提升是不可言喻的。

表1 不同方法模態(tài)頻率計(jì)算結(jié)果對(duì)比Table 1 Comparison of modal frequencies obtained by different methods

3.2 考慮液體推進(jìn)劑的全箭模態(tài)計(jì)算分析

本文采用動(dòng)態(tài)子結(jié)構(gòu)法+虛質(zhì)量縮聚方法來降低自由度數(shù)目，進(jìn)行了全箭的液體推進(jìn)劑虛質(zhì)量模擬方法研究。首先，將芯二級(jí)、芯一級(jí)、助推器的氧箱、燃箱提交計(jì)算生成廣義模型，然后采用超單元的方法將廣義模型與級(jí)間段、箱間段等其余部分有限元模型進(jìn)行耦合分析，得到全箭動(dòng)特性，最后與傳統(tǒng)集中質(zhì)量方法進(jìn)行對(duì)比。

對(duì)單獨(dú)助推器的零秒狀態(tài)進(jìn)行了模態(tài)計(jì)算。如圖2、圖3所示，通過彎曲模態(tài)的對(duì)比發(fā)現(xiàn)，兩種計(jì)算方法的計(jì)算結(jié)果雖有一定的差異，但是均在工程接受范圍之內(nèi)。但是通過圖4的縱向模態(tài)對(duì)比發(fā)現(xiàn)，兩種方法在計(jì)算結(jié)果差異較大，傳統(tǒng)方法計(jì)算得到的縱向模態(tài)要明顯高于虛質(zhì)量法。這是因?yàn)楦鶕?jù)液體推進(jìn)劑無粘性的特性，貯箱橫向彎曲變形時(shí)，除自由液面附近的推進(jìn)劑外，其余大部分推進(jìn)劑跟隨箭體結(jié)構(gòu)一起平動(dòng)，以上兩種方法均滿足此假設(shè)，因此計(jì)算結(jié)果一致。而在縱向變形時(shí)時(shí)，推進(jìn)劑僅跟隨箱底運(yùn)動(dòng)，虛質(zhì)量法對(duì)液體推進(jìn)劑運(yùn)動(dòng)的模擬與實(shí)際情況相符。而傳統(tǒng)RBE3液體推進(jìn)劑質(zhì)量建模方法相當(dāng)于液體質(zhì)量沿貯箱縱向均勻分布并隨箱壁一起運(yùn)動(dòng)，與實(shí)際情況不符，因此計(jì)算結(jié)果誤差較大。

集中質(zhì)量方法將液體當(dāng)作單純的集中質(zhì)量處理，虛質(zhì)量方法對(duì)附加質(zhì)量的計(jì)算并不是按照液體的實(shí)際質(zhì)量，而是按照液體運(yùn)動(dòng)方程，將液體內(nèi)部的運(yùn)動(dòng)方程轉(zhuǎn)化到邊界上，最終形成一個(gè)邊界自由度上的質(zhì)量陣。虛質(zhì)量法對(duì)液體推進(jìn)劑的模擬更加真實(shí)，克服了集中質(zhì)量方法的弊端。

圖2 助推器一階彎曲模態(tài)Fig.2 1st bending mode of booster

圖3 助推器二階彎曲模態(tài)Fig.3 2nd bending mode of booster

圖4 助推器一階縱向模態(tài)Fig.4 1st Longitudinal mode of booster

采用虛質(zhì)量方法對(duì)全箭的動(dòng)特性也進(jìn)行了計(jì)算，計(jì)算結(jié)果列入表2中，從表中可以看出，與單獨(dú)助推器計(jì)算得到的結(jié)論有所不同，兩種計(jì)算方法得到的全箭模態(tài)結(jié)果基本一致，甚至縱向的計(jì)算結(jié)果也基本一致。這與全箭的質(zhì)量特性和模態(tài)形狀有關(guān)，縱向模態(tài)如圖5所示，其模態(tài)形狀主要表現(xiàn)為助推器的上下運(yùn)動(dòng)，并且其縱向變形很小，推進(jìn)劑耦合效應(yīng)不明顯，此時(shí)助推器縱向表現(xiàn)為整體質(zhì)量效應(yīng)，兩種方法的對(duì)液體質(zhì)量的描述均合理，因此得到的計(jì)算結(jié)果一致。

通過以上討論可以看出，采用不同方法計(jì)算結(jié)果的差異與所關(guān)心的模態(tài)形狀、液體表現(xiàn)出的質(zhì)量特性有關(guān)。如果液體質(zhì)量的耦合效應(yīng)不明顯，集中質(zhì)量方法可以作為虛質(zhì)量法的一個(gè)特例，此時(shí)兩種方法計(jì)算得到的結(jié)果一致(如全箭動(dòng)特性的計(jì)算)，為了簡化可以采用集中質(zhì)量方法進(jìn)行計(jì)算。但是，如果耦合效應(yīng)不可忽略，集中質(zhì)量方法可能無法正確模擬液體推進(jìn)劑的運(yùn)動(dòng)(如單獨(dú)助推器的計(jì)算)，此時(shí)建議采用虛質(zhì)量法進(jìn)行計(jì)算。

表2 全箭模態(tài)計(jì)算結(jié)果對(duì)比Table 2 Modal frequencies of Launch Vehicle

圖5 全箭模態(tài)振形Fig.5 Modal shapes of launch vehicle

4 結(jié)論

本文對(duì)動(dòng)態(tài)子結(jié)構(gòu)方法在全箭動(dòng)特性分析中的應(yīng)用進(jìn)行了研究，提出了虛質(zhì)量法和動(dòng)態(tài)子結(jié)構(gòu)方法相結(jié)合的液體推進(jìn)劑建模方法，并采用不同液體推進(jìn)劑建模方法對(duì)全箭動(dòng)特性進(jìn)行了計(jì)算分析，得到的主要結(jié)論如下:

(1)動(dòng)態(tài)子結(jié)構(gòu)方法可以在保證計(jì)算精度的同時(shí)，極大的提高了計(jì)算效率，比較適合于全箭此類復(fù)雜結(jié)構(gòu)的動(dòng)力學(xué)問題研究;

(2)動(dòng)態(tài)子結(jié)構(gòu)與虛質(zhì)量相結(jié)合的液體推進(jìn)劑建模方法，一方面解決了計(jì)算效率問題，另一方面也解決了傳統(tǒng)集中質(zhì)量方法的弊端，可以實(shí)現(xiàn)縱、橫、扭一體化建模;

(3)對(duì)于彎曲模態(tài)，虛質(zhì)量方法和集中質(zhì)量法均可以很好的描述液體質(zhì)量的運(yùn)動(dòng)狀態(tài)，因此計(jì)算結(jié)果一致;對(duì)于縱向模態(tài)，計(jì)算結(jié)果與模態(tài)形狀、液體表現(xiàn)出的質(zhì)量特性有關(guān)。

［1］林宏，羅恒，潘忠文，等.運(yùn)載火箭動(dòng)特性有限元模型修正技術(shù)研究［J］.載人航天，2011，(6):30-34.

［2］Hou S N.Review of modal synthesis techniques and a new approach［J］.Shock and Vibration Bulletin，1969，40(4):25-39.

［3］MacNeal R H.A hybrid method of component mode synthesis［J］.Computers＆ Structures，1971，1(4):581-601.

［4］Rubin S.Improved component-mode representation for structural dynamic analysis［J］.AIAA Journal，1975，13(8):995-1006.

［5］Craig R，Chang C J.Free-interface methods of substructure coupling for dynamic analysis［J］.AIAA Journal，1976，14(11):1633-1635.

［6］王文亮，杜作潤，陳康元.模態(tài)綜合技術(shù)短評(píng)和一種新的改進(jìn)［J］.航空學(xué)報(bào)，1979，3:32-51.

［7］Hurty W C.Dynamic analysis of structural systems using component modes［J］.AIAA Journal，1965，3(4):678-685.

［8］Suarez L E，Singh M P.Improved fixed interface method for modal synthesis［J］.AIAA Journal，1992，30(12):2952-2958.

［9］Kubomura K.A theory of substructure modal synthesis［J］.Journal of Applied Mechanics，1982，49(4):903-909.

［10］Jezequel L，Seito H D.Component modal synthesis methods based on hybrid models，part i:Theory of hybrid models and modal truncation methods［J］.Journal of Applied Mechanics，1994，61(1):100-108.

［11］邱吉寶，向樹紅，張正平.計(jì)算結(jié)構(gòu)動(dòng)力學(xué)［M］.合肥:中國科學(xué)技術(shù)大學(xué)出版社，2009:346-446.

［12］楊劍，張璞，陳火紅.新編MD Nastran有限元實(shí)例教程［M］.北京:機(jī)械工業(yè)出版社，2007:315-335.