動態子結構方法在全箭動特性分析中的應用研究

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

(北京強度環境研究所可靠性與環境工程技術重點實驗室，北京100076)

1 引言

運載火箭全箭動力學特性是姿態控制系統設計、載荷分析、POGO分析等的基礎和依據，全箭動特性建模和分析是必不可少的［1］。運載火箭的結構變得越來越復雜，呈現出空間模態特征，僅僅依靠傳統的梁模型無法很好的反映其動力學特性，需建立復雜的三維模型，導致相應的計算成本增加，有時所建立的有限元模型甚至會有上百萬個自由度，如何提高計算效率成為解決該問題的關鍵。子結構模態綜合法是解決大型復雜結構特征值計算問題的有效方法，模態綜合法將整個系統劃分為多個子結構，每個子結構獨立建模，然后再組配成整體結構，其優點是可以縮減系統自由度，提高計算效率。根據子結構界面約束條件的不同，可以把子結構位移表達式分為三類:①自由界面模態綜合法［2-6］;②約束界面模態綜合法［7-10］;③混合界面模態綜合法［3，11］，邱吉寶對模態綜合法進行了系統的研究，并在文獻［11］中對動態子結構法進行了總結。

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

為了解決上述問題，本文提出了采用動態子結構和虛質量相結合的方法，在保證計算精度的同時降低了求解規模，提高了計算效率。

2 基本理論

2.1 虛質量方法

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

假設流體是無粘、不可壓縮的理想流體，流場是擬靜態場，結構是線彈性小變形結構，根據流體力學模型，利用Helmholtz方法求解Laplace方程，借助邊界元方法離散技術進行數值求解，可獲得液體對結構的集中質量。

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

式中，vi——任意結點ri處速度向量;

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

Aj——結構體表面上微元面積;

σj——j結點處的流速向量;

eij——從j點到i點的單位向量

ρ——流體密度

ri——結點位置

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

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

獲得附加質量后，質量矩陣會變成一個相對滿的陣，這給求解特征值問題帶來一定的困難，需采用動態子結構方法進行模型縮減。

2.2 動態子結構模態綜合方法

考慮液體附加質量后的結構動力學方程如式(6)所示。

其中，M ＝M0+MA，M0為結構質量陣，MA為液體附加質量陣，K為剛度陣，F為載荷向量。

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

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

約束模態綜合法的子結構坐標變換矩陣由邊界約束的子結構主模態以及約束模態組成。由模態展開定理可以得到式(8)。

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

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

可以將廣義模型當作物理模型一樣看待完成耦合分析。耦合分析過程是根據各子結構交界面的位移和力的協調條件完成整體耦合結構運動方程的組裝，然后求解整體運動方程的動力學問題。

兩個獨立子結構a、b的動力學方程根據式(13)，可以寫成式(14)～(15)。

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

可以將廣義模型當作物理模型一樣看待，應用式(16)(17)所示的協調條件和平衡條件

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

式(18)即為耦合系統的整體動力學方程，求解(18)時即可得到整體耦合系統的模態信息。

對于多級子結構耦合分析，只需按照以上步驟依次進行計算即可。

3 算例分析

3.1 全箭結構有限元建模與分析

有限元模型主體結構采用梁—殼三維模型，推進劑采用集中質量進行模擬，只計質量，不計轉動慣量，并通過RBE3與貯箱殼單元連接;助推與芯級連接采用梁(桿)單元進行模擬，并釋放旋轉自由度模擬鉸接;發動機采用等效梁單元模擬，所建立的有限元模型如圖1所示。

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

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

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

表1 不同方法模態頻率計算結果對比Table 1 Comparison of modal frequencies obtained by different methods

3.2 考慮液體推進劑的全箭模態計算分析

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

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

集中質量方法將液體當作單純的集中質量處理，虛質量方法對附加質量的計算并不是按照液體的實際質量，而是按照液體運動方程，將液體內部的運動方程轉化到邊界上，最終形成一個邊界自由度上的質量陣。虛質量法對液體推進劑的模擬更加真實，克服了集中質量方法的弊端。

圖2 助推器一階彎曲模態Fig.2 1st bending mode of booster

圖3 助推器二階彎曲模態Fig.3 2nd bending mode of booster

圖4 助推器一階縱向模態Fig.4 1st Longitudinal mode of booster

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

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

表2 全箭模態計算結果對比Table 2 Modal frequencies of Launch Vehicle

圖5 全箭模態振形Fig.5 Modal shapes of launch vehicle

4 結論

本文對動態子結構方法在全箭動特性分析中的應用進行了研究，提出了虛質量法和動態子結構方法相結合的液體推進劑建模方法，并采用不同液體推進劑建模方法對全箭動特性進行了計算分析，得到的主要結論如下:

(1)動態子結構方法可以在保證計算精度的同時，極大的提高了計算效率，比較適合于全箭此類復雜結構的動力學問題研究;

(2)動態子結構與虛質量相結合的液體推進劑建模方法，一方面解決了計算效率問題，另一方面也解決了傳統集中質量方法的弊端，可以實現縱、橫、扭一體化建模;

(3)對于彎曲模態，虛質量方法和集中質量法均可以很好的描述液體質量的運動狀態，因此計算結果一致;對于縱向模態，計算結果與模態形狀、液體表現出的質量特性有關。

［1］林宏，羅恒，潘忠文，等.運載火箭動特性有限元模型修正技術研究［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］王文亮，杜作潤，陳康元.模態綜合技術短評和一種新的改進［J］.航空學報，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］邱吉寶，向樹紅，張正平.計算結構動力學［M］.合肥:中國科學技術大學出版社，2009:346-446.

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