基于多層蒙特卡羅的火箭結(jié)構(gòu)動力學(xué)不確定性分析*

袁 赫，劉 莉，康 杰

(1 北京理工大學(xué)宇航學(xué)院， 北京 100081； 2 飛行器動力學(xué)與控制教育部重點(diǎn)實(shí)驗(yàn)室， 北京 100081)

0 引言

在運(yùn)載火箭結(jié)構(gòu)動力學(xué)響應(yīng)分析過程中，存在諸多不確定性，例如外載荷的不確定性、建模時(shí)模型的不確定性，這些不確定性因素是客觀存在的，粗略地忽略這些不確定性因素會造成分析誤差。

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

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

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

1 不確定性分析方法

1.1 蒙特卡羅仿真

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

(1)

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

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

(2)

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

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

(3)

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

1.2 多層蒙特卡羅仿真

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

圖1 板結(jié)構(gòu)層級有限元模型

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

(4)

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

(5)

假設(shè)式(5)中的每層標(biāo)準(zhǔn)蒙特卡羅估計(jì)是獨(dú)立的，則該估計(jì)的均方誤差可以表示為:

(6)

式中：Vl=Var[Yl]。

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

(7)

式中：Cl為計(jì)算一次Yl的計(jì)算代價(jià)。

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

(8)

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

(9)

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

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

(10)

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

2)Vl≤C22-βl；

3)Cl≤C32γl。

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

(11)

2 運(yùn)載火箭模型建立

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

2.1 主體結(jié)構(gòu)建模

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

2.2 液體推進(jìn)劑建模

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

(12)

(13)

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

2.3 捆綁連接結(jié)構(gòu)建模

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

圖2 捆綁連接結(jié)構(gòu)

2.4 發(fā)動機(jī)結(jié)構(gòu)建模

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

3 結(jié)構(gòu)動力學(xué)不確定性分析

3.1 計(jì)算工況

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

3.2 多層蒙特卡羅應(yīng)用思路

1)模型層級劃分

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

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

圖3 運(yùn)載火箭層級有限元模型

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

2)具體實(shí)施步驟

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

步驟二：從l=0開始計(jì)算Ql、Ql-1，當(dāng)l=0時(shí)Ql-1=0；

步驟三：計(jì)算方差Vl=Var[Yl]；

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

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

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

圖4 Cl和Ql隨層級的變化

3.3 結(jié)果分析

1)計(jì)算效率對比

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

圖5 計(jì)算時(shí)長對比

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

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

2)計(jì)算結(jié)果分析

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

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

3)參數(shù)估計(jì)

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

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

4 結(jié)論

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