袁 赫，劉 莉，康 杰

(1 北京理工大學宇航學院， 北京 100081； 2 飛行器動力學與控制教育部重點實驗室， 北京 100081)

0 引言

在運載火箭結構動力學響應分析過程中，存在諸多不確定性，例如外載荷的不確定性、建模時模型的不確定性，這些不確定性因素是客觀存在的，粗略地忽略這些不確定性因素會造成分析誤差。

由于運載火箭結構復雜，目前運載火箭結構動力學響應分析主要采用有限元方法[1]。不確定性分析方法主要分為非采樣方法和采樣方法。非采樣方法包括泰勒展開法[2]、區(qū)間攝動法[3]、隨機伽遼金法[4]等，這些方法雖然高效，但受算法侵入性的局限，對于運載火箭這種大型復雜結構來說適用性差。為考慮復雜結構分析中的不確定性，可將有限元方法和采樣方法結合，經(jīng)典例子是對有限元模型進行蒙特卡羅(MC)仿真[5-6]，后續(xù)出現(xiàn)的擬蒙特卡羅法[7]、馬爾科夫蒙特卡羅法[8]是MC方法的延伸，由于MC方法的非侵入性，使其成為了大型有限元結構不確定性分析簡單有效的方法。但大型有限元結構的結構動力學響應分析代價往往較高，尤其是對帶有多個不確定性變量且計算精度要求高的問題，MC方法的計算代價將無法接受。

多層蒙特卡羅方法(MLMC)是基于方差減小的思想降低計算代價的一種期望估計方法，起先用來求解帶有隨機參數(shù)的偏微分方程[9]，后續(xù)成功的與有限元方法結合應用到工程結構的靜力學[10]、熱力學[11]、可靠性[12]的不確定性分析中，提高了計算效率，但在運載火箭結構動力學不確定性分析領域還鮮有應用。

文中將MLMC方法引入運載火箭結構動力學不確定性分析中，結合有限元方法，以某一具有較大長徑比助推器的捆綁式運載火箭為研究對象，分析了帶有模型不確定性的運載火箭結構在點火起飛工況下的捆綁連桿動力學響應，與MC方法相比，MLMC方法計算效率有顯著的提高。最后，給出了捆綁連桿動響應的統(tǒng)計特性，為運載火箭初步設計階段的強度分析和結構設計提供基礎。

1 不確定性分析方法

1.1 蒙特卡羅仿真

MC方法求解不確定性問題的基本思想是：構造一個概率模型或隨機過程，使問題的解等于其參數(shù)，然后生成N個服從一定分布規(guī)律的采樣點，計算采樣點響應的統(tǒng)計特性，響應期望由式(1)估計：

(1)

MC估計方法的誤差來源于兩種，一種是由有限元網(wǎng)格數(shù)量引起的系統(tǒng)誤差，對于每個采樣點來說隨著有限元模型自由度數(shù)M的增加，其估計值收斂于真值，即

|E[QM-Q]|≤C1M-α

(2)

式中:α為收斂階數(shù)，C1為依賴于M的常數(shù)。

第二種誤差為受有限次采樣數(shù)量影響的采樣誤差，則總誤差可以用均方誤差表示如下[13]:

(3)

由式(3)可知，為了減小估計誤差，最直接的方法是提高采樣數(shù)量或者增加有限元模型網(wǎng)格數(shù)量，但是對于大型結構的結構動力學響應分析來說，這兩種途徑會造成難以接受的計算代價。

1.2 多層蒙特卡羅仿真

中心極限定理證明MC方法提高計算效率的有效方法之一是減小樣本點響應的方差。MLMC方法正是一種基于方差減小的思想來提高計算效率，它通過不同有限元網(wǎng)格數(shù)量的模型組合來估計目標值期望。依據(jù)有限元網(wǎng)格數(shù)量的精細程度，模型可劃分為0到L級，每層模型在每個維度上的網(wǎng)格尺寸比上一層模型精細，精細倍數(shù)為K，文中將K取為2，如圖 1所示的板結構層級有限元模型。

圖1 板結構層級有限元模型

利用期望算子的線性特性，L層模型目標值期望可以表示成0層模型目標值期望加上連續(xù)層級模型目標值之差的期望，即

(4)

等式右側(cè)的無偏估計為：

(5)

假設式(5)中的每層標準蒙特卡羅估計是獨立的，則該估計的均方誤差可以表示為:

(6)

式中：Vl=Var[Yl]。

為了使均方根誤差小于給定誤差ε2，可使式(6)第一項采樣誤差和第二項系統(tǒng)誤差均小于0.5ε2。在采樣誤差的既定要求下，同時使計算代價最小，該問題可以轉(zhuǎn)化為如下數(shù)學問題:

(7)

式中：Cl為計算一次Yl的計算代價。

采用拉格朗日乘子法，將Nl視為連續(xù)變量，從而得到每層模型最優(yōu)采樣數(shù)量。

(8)

隨著模型層級的增加，Ql和Ql-1收斂于真值Q，因此

(9)

式(8)表明采樣數(shù)量Nl是l的減函數(shù)，也就意味著多層蒙特卡羅估計的采樣數(shù)量集中在低層模型也就是計算代價小的模型上。

式(6)第二項系統(tǒng)誤差決定了模型的最高層級L，假設|E[Yl]|的收斂速度的階為O(K-α)，則

(10)

1)|E[Ql-Q]|≤C12-αl；

2)Vl≤C22-βl；

3)Cl≤C32γl。

因此，存在一個正數(shù)C4，對于任意ε≤e-1，存在最高層級L和相應的采樣數(shù)量使得多層蒙特卡羅估計的均方誤差小于給定誤差，并且計算代價有界，如式(11)所示。

(11)

2 運載火箭模型建立

運載火箭建模方法主要有3種：集中質(zhì)量-梁建模方法、局部三維建模方法、三維建模方法。其中集中質(zhì)量-梁模型和局部三維模型針對實際結構做了大量的簡化，而三維模型能更好的模擬運載火箭結構的剛度特性，尤其是對于超靜定結構，能更好的模擬捆綁結構的傳力情況[14]。文中基于Nastran建立運載火箭三維有限元模型，模型數(shù)據(jù)引用文獻[15]，芯級全長52 m，直徑3.35 m，4個助推器全長26 m，直徑2.25 m，芯級與助推器之間采用前、中、后3個捆綁結構，以下分別敘述各部分結構建模方法。

2.1 主體結構建模

這里的主體結構指運載火箭整流罩、級間段、助推器等部段，大部分為薄壁加筋結構。為了簡便起見，這里采用當量厚度的薄殼來模擬，建模原則為拉壓、彎曲、扭轉(zhuǎn)剛度與原來相同。在建模時首先在站點處建立結構質(zhì)量點，四周按照箭體半徑布置相應的殼單元，質(zhì)量點與殼之間采用RBE3單元連接，將集中質(zhì)量平均分配到殼上。

2.2 液體推進劑建模

由于液體推進劑沒有剛度，建模時只需模擬其質(zhì)量特性。貯箱橫向彎曲變形時，因液體推進劑無黏性，推進劑不會隨著貯箱壁面一起轉(zhuǎn)動；在貯箱縱向變形時，液體只跟隨箱底一起運動。因此，建模時將貯箱節(jié)點附近的推進劑按集中質(zhì)量法等效到該節(jié)點上，縱向質(zhì)量效應只作用在貯箱底部節(jié)點上，橫向質(zhì)量效應作用在貯箱的各個站點，貯箱柱段和貯箱下底節(jié)點的耦合質(zhì)量矩陣分別為式(12)和式(13)[14]。

(12)

(13)

其中∑me為貯箱內(nèi)推進劑總質(zhì)量；me為節(jié)點附近區(qū)域推進劑質(zhì)量。

2.3 捆綁連接結構建模

捆綁連接結構是捆綁式火箭非常重要的連接方式，文中模型采用三捆綁連接方式，圖 2(a)給出了由3根連桿(兩根直桿，一根斜拉桿)組成的前、中輔助捆綁結構。建模時首先確定桿在助推器上的連接中心，依據(jù)桿的空間角度確定芯級的連接中心，在兩個連接中心之間建立桿單元。為避免捆綁連接處的剛度過于薄弱，在芯級和助推器捆綁點處建立加強框，桿的連接中心用RBE2單元固接到加強框上。主捆綁是典型的球頭-球窩結構，中間沒有轉(zhuǎn)動約束，建模時采用兩根梁單元模擬球頭和球窩，兩根梁采用多點約束連接，并釋放轉(zhuǎn)動自由度來模擬轉(zhuǎn)動效應。加強框的剛度一般由地面點火實驗確定，在運載火箭設計初期，該結構的剛度不能準確的給出，只能憑經(jīng)驗給出一定的范圍，具有較大的不確定性。

圖2 捆綁連接結構

2.4 發(fā)動機結構建模

發(fā)動機結構包括發(fā)動機支架及發(fā)動機本身，發(fā)動機支架的作用是將發(fā)動機推力傳遞到火箭上。建模時首先確定發(fā)動機的質(zhì)心位置，在該位置和發(fā)動機對接面站點之間建立梁單元模擬支架，支架質(zhì)量等效到發(fā)動機上。對于發(fā)動機本身來說，由于文中分析的是關于全箭整體的工況，不需要發(fā)動機復雜的結構形式，建模時采用集中質(zhì)量單元模擬其質(zhì)量特性即可，建立原則為其質(zhì)量大小和位置與真實結構相等。依據(jù)上述建模方法，運載火箭的三維模型如圖 3所示。

3 結構動力學不確定性分析

3.1 計算工況

文中的研究對象為帶有四助推器的超靜定捆綁運載火箭三維模型，計算工況為豎立在發(fā)射臺上時發(fā)動機點火時刻，發(fā)動機開機曲線采用50噸級氫氧發(fā)動機開機試驗數(shù)據(jù)[16]，對全箭做瞬態(tài)響應分析，計算捆綁連桿中受載情況最嚴苛的中捆綁直桿的動響應(軸力)極值。不確定性變量為輔助捆綁點的加強框剛度，服從均值為7.4×1011Pa，標準差為1.51×1011Pa的正態(tài)分布。

3.2 多層蒙特卡羅應用思路

1)模型層級劃分

按1.2節(jié)給出的層級模型網(wǎng)格細化方法，將運載火箭三維模型分為4個等級(L=3)，每一層級模型兩個站點間的四邊形網(wǎng)格在徑向和軸向較上一層級細化2倍。表 1列出了各層級模型每兩個站點間網(wǎng)格數(shù)量，圖 3給出了各層級火箭有限元模型。

表1 層級模型網(wǎng)格數(shù)量

圖3 運載火箭層級有限元模型

首先，加強框剛度取均值對各層模型做開機工況計算。圖 4給出了Cl的變化情況，其對數(shù)值正比于2γl，參數(shù)γ取決于計算機性能，這里γ≈2.4。橙色圖線為各層模型的中捆綁直桿的動響應極值，隨著模型層級的增加，網(wǎng)格的不斷細化，該值趨于收斂，從工程應用的角度來說，最高層級L=3已經(jīng)足夠精確。

2)具體實施步驟

步驟一：給出每層模型初始采樣量N0，隨機生成N0個服從正態(tài)分布的加強框剛度，按生成的剛度修改模型；

步驟二：從l=0開始計算Ql、Ql-1，當l=0時Ql-1=0；

步驟三：計算方差Vl=Var[Yl]；

步驟四：根據(jù)式(8)計算每層模型最優(yōu)采樣數(shù)量Nl；

步驟五：定義ΔNl=Nl-N0，若ΔNl>0，重復步驟一，此時令樣本數(shù)量為ΔNl重新計算Ql，Ql-1；

步驟六：依據(jù)式(10)進行收斂性檢查，確定最高層級L，若不滿足式(10)條件，令L=L+1重回步驟二。

圖4 Cl和Ql隨層級的變化

3.3 結果分析

1)計算效率對比

本節(jié)從計算時長和采樣數(shù)量兩方面說明MLMC方法在運載火箭結構動力學不確定性分析的高效性。圖 5比較了MC方法和MLMC方法在達到指定相對均方誤差θ時的實際仿真時間，結果顯示，在不同θ下，MLMC方法計算的時長始終小于MC方法，計算效率提高了70～80倍。

圖5 計算時長對比

表2列出了在不同θ下MLMC方法和MC方法的采樣數(shù)量，以及MLMC方法在L層模型上的當量采樣數(shù)量，分析結果可知隨著模型層級的增加，采樣數(shù)量Nl大幅度下降；當精度要求提高時，每層的采樣數(shù)量接近于等幅度增加，這是由式(8)決定的；雖然MLMC方法在各層級模型上的總采樣數(shù)量總體高于MC方法，但是其絕大部分的采樣點集中在計算代價非常小的0層模型上，這解釋了圖 5顯示的MLMC方法計算時長較MC方法大幅度降低現(xiàn)象，對于MC方法，上述結果在θ=2.5×10-5時的值省去了，這是由于其計算代價巨大，難以接受。

表2 MLMC和MC方法采樣數(shù)量對比

2)計算結果分析

圖 6給出了MC方法Ql的方差和θ=6.25×10-5時MLMC方法Yl的方差情況，隨著層級的增加，Yl的方差大幅度下降，說明MLMC方法能明顯減小方差，參數(shù)β可由Yl方差圖線斜率估計得到，Ql的方差大致保持常數(shù)。因此，MLMC方法的L層模型響應的方差特性可以由低層模型表示，而其均值特性可以通過MLMC方法給出的式(4)進行逼近，應用這個思想，該算例中的中捆綁直桿動響應極值的均值為4.108×104N，方差為1.74×104N2。

圖6 方差與層數(shù)之間的關系

3)參數(shù)估計

表 3列出了MLMC采樣的參數(shù)估計結果，在這四種情況下均β>γ，按定理1的結論，MLMC在本算例中的計算代價應正比于θ-2，這符合圖 5所給出的結果。

表3 MLMC方法參數(shù)估計

4 結論

文中將MLMC方法引入到運載火箭結構動力學不確定性分析中。考慮模型的不確定性，分析了發(fā)動機點火工況下運載火箭捆綁連桿響應。結果表明，相較于MC方法，MLMC方法的計算效率有很大的提高，在文中給出的計算精度下，計算效率提高了70～80倍；同時，MLMC方法通過低層模型模擬不確定性傳播的方差特性，通過低高層模型的結合逼近不確定性傳播的期望特性，能高效地獲取包含均值和方差的捆綁連桿響應極值統(tǒng)計特性，為運載火箭初步設計階段考慮不確定性的強度分析和結構設計提供基礎。