火箭捆綁機構(gòu)熱力耦合模型修正與虛擬驗證方法

賀 揆, 張 薇, 李長龍, 裴長浩

(1. 北京航空航天大學(xué) 可靠性與系統(tǒng)工程學(xué)院,北京 100191;2. 北京宇航系統(tǒng)工程研究所,北京 100076)

在重型運載火箭的飛行過程中,球頭-球窩式捆綁機構(gòu)為了適應(yīng)芯級與助推器之間的相對彈性變形會發(fā)生轉(zhuǎn)動摩擦,大噸位載荷作用下的轉(zhuǎn)動摩擦?xí)a(chǎn)生大量熱量導(dǎo)致捆綁機構(gòu)溫度過高,從而嚴(yán)重影響整體強度甚至導(dǎo)致結(jié)構(gòu)失效[1],因此需要精確的熱力耦合有限元分析來輔助結(jié)構(gòu)設(shè)計以及結(jié)構(gòu)優(yōu)化。

對于熱力耦合有限元建模,袁水林等[2]采用ANSYS軟件建立了某型號球頭-球窩式捆綁機構(gòu)的有限元模型并進(jìn)行了摩擦熱分析,研究了摩擦系數(shù)對捆綁機構(gòu)溫度分布的影響。除進(jìn)行熱力耦合有限元分析以及研究各參數(shù)對結(jié)果的影響外,張薇等[3]還建立了摩擦熱理論模型,并將理論分析結(jié)果與有限元仿真結(jié)果相互驗證。

但是,考慮到熱力耦合這種復(fù)雜的物理過程以及有限元建模中存在的不確定性,有限元分析的結(jié)果與實際情況不可避免地存在著差異。于是為了減小這一差異,提出使用模型修正技術(shù),基于已有的試驗觀測值,通過調(diào)整有限元模型的材料參數(shù)、幾何特征等輸入?yún)?shù),來減小仿真結(jié)果與試驗數(shù)據(jù)的差異,從而提高有限元模型的計算精度[4]。傳統(tǒng)的模型修正主要針對結(jié)構(gòu)動力學(xué)的有限元模型,一般將彈性模量、密度等材料參數(shù)作為輸入?yún)?shù),輸出特征則僅考慮結(jié)構(gòu)的固有頻率或者頻響函數(shù)(frequency response function, FRF)[5]。然而,在考慮熱環(huán)境的有限元模型修正中,熱傳導(dǎo)系數(shù)、比熱容等熱力學(xué)參數(shù)由于存在不確定性也需要作為待修正參數(shù),同時結(jié)構(gòu)溫度則需考慮作為輸出特征。張保強等[6]針對熱傳導(dǎo)模型確認(rèn)挑戰(zhàn)問題,使用貝葉斯方法對熱傳導(dǎo)系數(shù)以及體熱容等參數(shù)進(jìn)行了修正。而對于更復(fù)雜的飛行器結(jié)構(gòu)熱力耦合模型修正問題,現(xiàn)有的研究主要以優(yōu)化的思想來實現(xiàn),通常先構(gòu)造目標(biāo)函數(shù)、確定設(shè)計變量,然后通過遺傳算法[7]、粒子群算法[8-9]以及擬牛頓法[10-11]等優(yōu)化算法進(jìn)行求解。雖然上述優(yōu)化類算法在減小觀測量偏差上確實能起到一定效果,但輸入?yún)?shù)還是難以被修正到精確值[12]。除直接使用優(yōu)化算法以外,何成等[13]還提出了一種分層修正策略,并推導(dǎo)了不確定性參數(shù)均值及其協(xié)方差矩陣的迭代修正方程。目前關(guān)于熱力耦合模型修正的方法還比較單一,所以有必要嘗試更多其他類型的模型修正算法,來解決更復(fù)雜、高維的熱力耦合模型修正問題。

模型修正需要試驗觀測數(shù)據(jù)作為依據(jù),然而在很多實際工程應(yīng)用中,試驗數(shù)據(jù)是比較稀缺的。比如火箭捆綁機構(gòu),其實際工作環(huán)境中載荷巨大、工況嚴(yán)峻,進(jìn)行地面模擬試驗成本高、難度大,這也就造成運載火箭新型號的設(shè)計周期長、設(shè)計效率低。因此,如何充分利用前期型號的已有試驗數(shù)據(jù)來驗證甚至修正在研新型號結(jié)構(gòu)的有限元模型就成為了關(guān)鍵的技術(shù)難題。

為解決上述難題,本文提出一種熱力耦合模型修正與虛擬驗證方法,并將以上方法用于修正某在研型號火箭捆綁機構(gòu)的熱力耦合有限元模型。首先,基于某前期型號的試驗數(shù)據(jù),使用基于距離判別的貝葉斯修正方法對前期型號有限元模型進(jìn)行模型修正;然后,通過修正后的前期型號有限元模型標(biāo)定通用化的摩擦熱理論模型,建立溫度與結(jié)構(gòu)幾何尺寸、材料參數(shù)以及工況參數(shù)的關(guān)系公式;最后,將在研型號的參數(shù)代入通用理論模型得到理論溫度值,并將理論值作為依據(jù)進(jìn)行新一輪模型修正,從而得到精確的在研型號有限元模型。

1 模型修正算法

這部分介紹本文所采用的模型修正方法,將用于修正后續(xù)熱力耦合有限元模型。目前常用的模型修正方法包括優(yōu)化類修正、協(xié)方差修正、區(qū)間類修正以及貝葉斯修正。優(yōu)化類修正主要基于模型輸入輸出參數(shù)間的梯度信息,為傳統(tǒng)的確定性修正算法;協(xié)方差與區(qū)間修正除參數(shù)均值外還考慮了參數(shù)的方差以及范圍信息,涉及了模型修正中的不確定性分析;貝葉斯修正是一類典型的隨機修正方法,利用待修正參數(shù)的先驗分布與后驗分布描述修正前后的參數(shù)隨機特征。上述各類方法的詳細(xì)比較與特征可見Broggi等[14-17]研究。本文將使用貝葉斯修正算法,一方面因為該方法適用于解決高維問題,往往能夠極大縮減預(yù)測值與實際值的誤差,同時還具有很高的計算效率;另一方面,熱力耦合模型修正問題面臨輸入?yún)?shù)不確定性差異巨大且相互耦合等難點,貝葉斯修正方法得到的輸入?yún)?shù)后驗分布有助于在統(tǒng)計意義下估計最優(yōu)解。

同時,熱力耦合的單次有限元分析十分耗時,而模型修正算法又需要大量地調(diào)用有限元模型,因此本文將通過二次多項式響應(yīng)面的方法建立代理模型,從而盡可能減少計算成本。

1.1 貝葉斯修正框架

本文所采用的貝葉斯修正框架將結(jié)合近似貝葉斯計算方法(approximate Bayesian computation, ABC)以及過渡馬爾科夫鏈蒙特卡洛算法(transitional Markov chain Monte Carlo, TMCMC),該框架的核心為概率統(tǒng)計理論中的貝葉斯原理,表示為

(1)

式中:x為輸入?yún)?shù),即待修正參數(shù);Yexp為試驗觀測樣本;P(x)為待修正參數(shù)的先驗分布;P(x|Yexp)為待修正參數(shù)后驗分布,也就是依據(jù)觀測值修正得到的結(jié)果;P(Yexp)為歸一化參數(shù),本質(zhì)是觀測數(shù)據(jù)的概率積分,其確保后驗分布概率累計積分為1;PL(Yexp|x)為似然函數(shù),本質(zhì)是現(xiàn)有試驗觀測數(shù)據(jù)基于輸入?yún)?shù)樣本的條件概率。

貝葉斯修正框架中的一大難題就是歸一化參數(shù)P(Yexp)難以計算,因此提出使用TMCMC算法,通過迭代過程,從一系列中間概率密度函數(shù)采樣去逐漸逼近后驗分布。第j代中間概率密度函數(shù)可表示為

Pj∝PL(Yexp|x)βjP(x)

(2)

式中,βj為似然函數(shù)的權(quán)重指數(shù),其中,j=1,2,…,m,m為總迭代次數(shù),在迭代開始時β0=0,式(2)右邊即為先驗分布,后隨著迭代逐漸逼近βm=1,等式右邊也就變?yōu)楹篁灧植?。每次迭代步中βj都根據(jù)上一個迭代步中的采樣樣本計算得到。TMCMC算法中,馬爾科夫鏈以及Metropolis-Hastings采樣方法不斷地從具有高似然值的樣本中傳播出新的樣本,從而實現(xiàn)了從復(fù)雜的后驗分布采樣。該方法的更多細(xì)節(jié)可見Ching等[18]研究。

貝葉斯修正框架中的第二個難點就是似然函數(shù)PL(Yexp|x)的計算,似然函數(shù)依照定義可寫為

(3)

由式(3)可以看到,似然函數(shù)的計算需要對每一個試驗樣本估計出基于待修正參數(shù)的條件概率,而每一次估計都需要大量地調(diào)用數(shù)值模型,因此計算成本極高。為解決這個難題,提出使用ABC方法來簡化似然函數(shù),表示為

(4)

式中:σ為寬度系數(shù),一般可取為[10-3,10-1][19];d為統(tǒng)計學(xué)距離指標(biāo),用來衡量預(yù)測樣本與觀測樣本的差異。本文選用歐氏距離作為距離度量指標(biāo),其表達(dá)式為[20]

(5)

1.2 二次多項式響應(yīng)面

多項式響應(yīng)面是最常用的代理模型之一,它在參數(shù)空間內(nèi)構(gòu)造輸入輸出間的多項式函數(shù)關(guān)系,從而逼近真實有限元模型,原理簡單、直接。對于熱力耦合有限元模型,輸出為溫度、應(yīng)力等物理量,輸入則為力學(xué)性能參數(shù)、熱性能參數(shù)等等,相應(yīng)的輸入輸出關(guān)系也就往往呈現(xiàn)單峰值、低階非線性、輸入?yún)?shù)變化范圍小等特點,因此本文使用二次多項式響應(yīng)面來構(gòu)建熱力耦合有限元分析的代理模型。對于第個輸出特征,基于二次多項式的代理模型可構(gòu)造為

(6)

式中,β為模型的待定系數(shù),通過最小二乘法估計得到。

在代理模型構(gòu)建完成后,還需對其精度進(jìn)行驗證。于是本文選擇均方根誤差(root mean square error, RMSE)作為代理模型的檢驗指標(biāo),表達(dá)式為

(7)

2 虛擬驗證方法

虛擬驗證方法重點解決新型號捆綁機構(gòu)設(shè)計中沒有任何試驗數(shù)據(jù)的情況,利用已有的捆綁機構(gòu)摩擦熱試驗歷史數(shù)據(jù)以及熱力耦合分析的通用理論公式,建立新舊型號之間的關(guān)聯(lián),從而驗證并修正在研新型號的有限元模型。

2.1 虛擬驗證框架

對于一個系列的火箭捆綁機構(gòu),其結(jié)構(gòu)外形、接觸方式以及運動形式基本一致,不同的型號僅在幾何尺寸、工況參數(shù)以及材料參數(shù)上存在一定范圍內(nèi)的差異。因此考慮建立通用化的摩擦熱理論模型,使得理論計算結(jié)果與有限元仿真結(jié)果在不同輸入?yún)?shù)下仍能保持一致。結(jié)合有限元模型、模型修正技術(shù)以及通用化理論模型,可以構(gòu)造出虛擬驗證的基本原理框架,如圖1所示。

圖1 虛擬驗證框架Fig.1 Virtual verification framework

在該虛擬驗證框架中:①需要對前期、在研兩型號的捆綁機構(gòu)進(jìn)行初步的有限元建模,其中部分輸入?yún)?shù)因認(rèn)知不完全存在不確定性;②以前期型號的摩擦熱試驗數(shù)據(jù)為觀測值,使用1.1節(jié)介紹的貝葉斯修正算法對前期型號有限元模型開展第一輪模型修正;③建立通用化理論模型,利用修正后的前期有限元模型的仿真結(jié)果去標(biāo)定理論模型的待定系數(shù);④將在研型號的輸入?yún)?shù)代入理論模型,將理論計算結(jié)果作為觀測值,在初步驗證有限元仿真結(jié)果合理后,對在研型號有限元模型進(jìn)行第二輪模型修正;⑤得到高精度的在研型號有限元模型,從而輔助在研型號后續(xù)的設(shè)計研制工作。

2.2 通用化理論模型

本文所使用的通用化理論模型為球頭-球窩式捆綁機構(gòu)摩擦熱計算模型,下面將列出部分重要的計算公式及推導(dǎo)過程。首先,建立球頭-球窩結(jié)構(gòu)的Hertz接觸模型,具體外形如圖2所示,計算得到球頭-球窩接觸形式的嚙合深度δ,表示為

圖2 球頭-球窩Hertz接觸模型Fig.2 Hertz contact model of ball-and-socket mechanism

(8)

式中:Fz為等效集中載荷; Δr為球頭與球窩的半徑差;E為材料彈性模量。同時,計算得到等效接觸面的橫向長度a,即圖2中AD段的長度,可寫為

(9)

于是,接觸面上的平均接觸壓力可以寫為

(10)

式中,λ為接觸面積修正系數(shù)。在此基礎(chǔ)上,假設(shè)等效溫升區(qū)域如圖2中陰影部分,根據(jù)熱力學(xué)定律以及集中參數(shù)法可以推導(dǎo)出該區(qū)域的溫度變化公式,表示為

(11)

式中:α為等效接觸區(qū)域質(zhì)量修正系數(shù);γ為熱散失項修正系數(shù);b為AB段長度;θ為球頭轉(zhuǎn)動角度;f為轉(zhuǎn)動頻率;T0為初始溫度;ρ為材料密度;h0為BD段長度;CV為材料比熱容。

通用化理論模型中包含3個待定系數(shù),分別為λ、α以及γ。因此在構(gòu)建模型時,需要以第一輪修正后的有限元模型作為依據(jù)去標(biāo)定這三個參數(shù)。本文采用分部策略,先通過有限元仿真的接觸壓力計算結(jié)果去標(biāo)定式(10)中的λ,然后再通過溫度仿真結(jié)果去標(biāo)定式(11)中的α與γ。

3 數(shù)值算例

本章將使用前文中提到的方法建立某型號運載火箭捆綁機構(gòu)的精確熱力耦合有限元模型。

3.1 捆綁機構(gòu)熱力耦合有限元模型建立

球頭球窩式捆綁機構(gòu)的外形如圖2所示,球接觸面的標(biāo)稱半徑為75 mm,存在半徑差Δr,即球頭球面半徑為(75-Δr)mm。為模擬捆綁機構(gòu)在實際工作情況中的轉(zhuǎn)動摩擦情況,在有限元仿真中令球窩固定,令球頭繞球心在XZ平面上進(jìn)行周期性的往復(fù)擺動,同時均布載荷持續(xù)作用在球頭上表面。仿真使用熱力耦合分析,單元類型設(shè)置為C3D8RT。網(wǎng)格尺寸整體設(shè)置為5 mm,接觸面區(qū)域加密至2 mm。通過演示算例可以看到該結(jié)構(gòu)的運動及傳熱效果,截取三個時刻的仿真結(jié)果如圖3所示。

圖3 演示算例仿真結(jié)果Fig.3 Simulation results of demonstration case

3.2 第一輪模型修正

第一輪模型修正的對象為某前期型號的熱力耦合有限元模型,其外形與3.1節(jié)中的建模一致。對于該型號,半徑差設(shè)置為0.5 mm,接觸面摩擦系數(shù)0.15,球頭轉(zhuǎn)動頻率2 Hz,最大轉(zhuǎn)角2°,外載荷等效為集中力,大小2×106N,初始溫度20 ℃。該型號捆綁機構(gòu)的材料為合金結(jié)構(gòu)鋼30CrMnSiNi2A。

首先,對待修正參數(shù)進(jìn)行選擇。有限元模型的輸入?yún)?shù)包括:材料力學(xué)性能參數(shù)、材料熱學(xué)性能參數(shù)以及熱環(huán)境參數(shù),這些參數(shù)在一定程度上都具有不確定性,在建立有限元模型時難以確定其數(shù)值。參考文獻(xiàn)[13],彈性模量以及導(dǎo)熱率隨溫度的改變會發(fā)生明顯變化,因此有必要先單獨對這兩個參數(shù)進(jìn)行處理。分別取溫度為20 ℃以及500 ℃時對應(yīng)的彈性模量以及導(dǎo)熱率,分別記為E20,E500,k20和k500, 并將這兩個參數(shù)的溫變特性用一次函數(shù)來描述,得到表達(dá)式為

(12)

為證明用一次函數(shù)描述溫變特性的合理性,針對高強度鋼30CrMnSiNi2A,通過JMat Pro數(shù)據(jù)庫軟件,計算一組不同溫度下的彈性模量以及導(dǎo)熱率作為參考值,繪制如圖4與圖5所示。取E20,E500,k20和k500構(gòu)造一次函數(shù)也繪制在圖4與圖5中,現(xiàn)考察該擬合函數(shù)對參考數(shù)據(jù)點的擬合精度。參考式(7),計算一次函數(shù)上數(shù)值相對于對應(yīng)參考數(shù)值的百分比形式的RMSE,彈性模量的擬合RMSE為0.59%,導(dǎo)熱率的為0.12%,因此認(rèn)為使用一次函數(shù)來擬合溫變特性是具有足夠高精度的。

圖4 彈性模量溫變模型Fig.4 The temperature variation model of elasticity modulus

圖5 導(dǎo)熱率溫變模型Fig.5 The temperature variation model of thermal conductivity

除此以外,材料參數(shù)中的密度、比熱容與線膨脹系數(shù)以及描述結(jié)構(gòu)散熱情況的表面換熱系數(shù)與輻射率都存在不確定性,因此也被選為待修正參數(shù)。最終,9個待修正參數(shù)及其變化范圍列在表1所示。

表1 第一輪模型修正的不確定參數(shù)

下一步,對輸出特征進(jìn)行選擇。對于熱力耦合分析,關(guān)注的輸出主要為應(yīng)力分布以及結(jié)構(gòu)升溫情況。在前期型號的摩擦熱試驗中,觀測量由一系列布在結(jié)構(gòu)外部的傳感器測得,包括實時的應(yīng)力與溫度。在模型修正的輸出特征選擇上,應(yīng)保證與試驗觀測量一致。因此,選擇球頭、球窩上4個測點的溫度及應(yīng)力值共8個物理量作為輸出特征,記為Si與Ti,i=1,2,3,4。4個測點的具體位置如圖6所示。

圖6 測點位置參考Fig.6 Location reference of measurement points

在模型修正中需大量調(diào)用有限元分析,因此對上述輸入輸出構(gòu)建代理模型。進(jìn)行試驗設(shè)計獲取多組有限元仿真結(jié)果,對上述9個待修正參數(shù)設(shè)計9因素3水平正交表,共27組試驗,記為L27(39)。低、中、高三個水平分別對應(yīng)著輸入?yún)?shù)變化范圍的下界、中點以及上界。于是用這27組數(shù)據(jù)計算式中的待定系數(shù)β,從而得到前期型號有限元模型的二次多項式響應(yīng)面代理模型。使用拉丁超立方采樣從取值范圍內(nèi)隨機得到5組驗證樣本,分別代入有限元模型與代理模型得到輸出響應(yīng)。對5組試驗樣本的有限元仿真值與代理模型擬合值進(jìn)行比較,并依照式計算對于每個輸出特征的RMSE,計算如表2所示?？梢钥吹?8個特征中均方根誤差最高達(dá)到3.00%,平均誤差為1.59%,精度能滿足要求,代理模型可用于后續(xù)模型修正。

表2 代理模型關(guān)于8個輸出特征的均方根誤差

代理模型構(gòu)建完畢后,根據(jù)1.1節(jié)中介紹的貝葉斯修正框架進(jìn)行模型修正。設(shè)置TMCMC樣本數(shù)為2 000,經(jīng)歷7次迭代后得到各輸入?yún)?shù)的后驗分布如圖7所示。圖7中的9幅圖為每個待修正參數(shù)2 000個樣本的后驗分布柱狀圖,區(qū)間對應(yīng)的值越大,說明該區(qū)間上的樣本點越多?？梢园l(fā)現(xiàn),經(jīng)過修正后的后驗分布非常集中,絕大多數(shù)樣本點都聚集在極小的區(qū)間內(nèi),因此各輸入?yún)?shù)的修正值可以簡單地取為相應(yīng)區(qū)間的中點,即圖7中的“*”號標(biāo)記點處,具體數(shù)值如表3所示。

圖7 第一輪模型修正后驗分布柱狀圖Fig.7 Posterior distributions after first round calibration

表3 第一輪模型修正輸入?yún)?shù)修正值Tab.3 Updated input parameters after first round calibration

將上述輸入?yún)?shù)的修正值再代入有限元模型中,驗證各輸出特征的仿真值與試驗值是否足夠接近。計算得到輸出特征的試驗值、仿真值以及誤差列在表4中,其中應(yīng)力單位為MPa,溫度單位為℃。由表4可得,輸出特征中最大的預(yù)測誤差僅為2.45%,各輸出特征的平均誤差更是只有1.13%,表明修正后的有限元模型已經(jīng)十分接近于實際情況。

表4 第一輪修正輸出特征及誤差

對修正后的輸入?yún)?shù)進(jìn)行有限元仿真,得到結(jié)構(gòu)末時刻的溫度場云圖如圖8所示。有限元仿真可以得到結(jié)構(gòu)整體的溫度分布以及最大溫度,這些數(shù)據(jù)在實際試驗中難以全面測量,但模型修正保證了這些仿真值的正確性,使其能夠在后續(xù)用于理論模型標(biāo)定。

圖8 第一輪修正后末時刻溫度云圖Fig.8 Temperature nephogram at final moment after first round calibration

3.3 理論模型標(biāo)定

基于第一輪修正后的有限元模型,對通用化理論模型進(jìn)行標(biāo)定。首先,通過接觸壓力標(biāo)定式(10)中的接觸面積修正系數(shù)λ。式(10)中計算得到的是理論接觸面上的平均接觸壓力,而實際上,在存在半徑差且不考慮彈性形變時,球頭-球窩接觸為點接觸,因此不妨就直接使用有限元仿真得到的球頭接觸面中心最大接觸壓力來進(jìn)行標(biāo)定。利用3.2節(jié)中的結(jié)果,最終標(biāo)定接觸面積修正系數(shù)λ=0.426。

同理,通過球頭頂點處的溫度仿真值,標(biāo)定式(11)中的等效接觸區(qū)域質(zhì)量修正系數(shù)α以及熱散失項修正系數(shù)γ。由于待標(biāo)定系數(shù)有兩個,因此需要構(gòu)造大于等于2組的仿真結(jié)果來作為標(biāo)定的依據(jù)。在不改變結(jié)構(gòu)外形以及材料的基礎(chǔ)上,考慮只改變轉(zhuǎn)動頻率來構(gòu)造多組數(shù)據(jù),并通過Levenberg-Marquardt方法減小仿真值與理論值的差異,將得到的最優(yōu)解作為α以及γ的標(biāo)定值。最終標(biāo)定等效接觸區(qū)域質(zhì)量修正系數(shù)α=0.104,熱散失項修正系數(shù)γ=0.043。將上述標(biāo)定數(shù)值代入理論模型,得到不同轉(zhuǎn)動頻率下溫度變化的理論計算值與仿真值,如圖9所示。可以看到,理論計算值與仿真值差異較小,因此可以認(rèn)為該理論模型能代表修正后有限元模型的計算結(jié)果。

圖9 不同轉(zhuǎn)動頻率下仿真與理論溫度值對比Fig.9 Comparison of temperatures calculated by FEM and theoretical model under different rotational frequencies

3.4 第二輪模型修正

第二輪模型修正的對象是在研新型號捆綁機構(gòu),其有限元模型基本與前期型號一致,但材料改變?yōu)槌恋碛不讳P鋼PH13-8Mo,且其他個別參數(shù)更改為:半徑差Δr=0.7 mm,球頭轉(zhuǎn)動頻率0.8 Hz,外載荷上升至2.4×106N。由于在研型號不具備試驗數(shù)據(jù),因此將理論模型計算值作為觀測值對在研型號有限元模型進(jìn)行驗證及修正。

類似3.2節(jié)中的步驟,先確定模型的輸入?yún)?shù),根據(jù)在研型號的材料PH13-8Mo,總結(jié)各輸入?yún)?shù)如表5所示。其中,由于彈性模量、密度以及比熱容是理論模型的輸入?yún)?shù),因此不將其考慮為待修正參數(shù),而直接認(rèn)為是確定值。所以,第二輪修正的待修正參數(shù)分別為20 ℃及500 ℃時的導(dǎo)熱率、線膨脹系數(shù)、表面換熱系數(shù)以及輻射率5個物理量。

表5 在研型號模型輸入?yún)?shù)Tab.5 Input parameters of new developing model

對于輸出特征,由于理論模型只能計算出球頭頂點的溫度值,因此取該位置上t=2 s,4 s, 6 s, 8 s, 10 s共5個時刻的溫度作為輸出特征。根據(jù)已確定的輸入輸出,進(jìn)行試驗設(shè)計,并構(gòu)建在研型號的代理模型。對5個待修正參數(shù)建立5因素3水平的正交表,共18組試驗,記為L18(35)。依照正交表進(jìn)行18組有限元仿真,并根據(jù)仿真結(jié)果建立二次多項式響應(yīng)面代理模型。同樣地,再獲得5組驗證樣本,并計算代理模型對于輸出特征的RMSE如表6所示,平均RMSE僅為0.68%,代理模型誤差較小,可用于后續(xù)模型修正。

表6 代理模型關(guān)于5個時刻溫度的均方根誤差

將18組試驗得到的末時刻球頭頂點溫度值的分布繪制如圖10所示,可以看到溫度的變化范圍為[271.50, 308.58]。將相應(yīng)輸入?yún)?shù)代入理論模型后,得到理論計算值為277.12 ℃,落在試驗樣本的輸出變化區(qū)間內(nèi)。由此,可以初步說明在研型號有限元仿真結(jié)果具有一定合理性,但仍需模型修正過程進(jìn)一步提高模型精度。

圖10 仿真樣本分布圖Fig.10 Distribution of simulation samples

進(jìn)一步地,依照理論計算值對在研型號有限元模型進(jìn)行修正,在TMCMC算法中設(shè)置樣本數(shù)為2 000,經(jīng)過6次迭代,得到后驗分布如圖11所示。由圖11可知,第二次模型修正得到的參數(shù)取值區(qū)間得到了明顯的縮減,但與第一輪修正的結(jié)果不同,修正后樣本在區(qū)間內(nèi)分布得相對分散。于是,在確定修正值時,使用核密度估計方法擬合出各參數(shù)后驗分布的概率密度函數(shù),并將概率密度函數(shù)最大值對應(yīng)的橫坐標(biāo)作為修正值, 圖11中的星形標(biāo)記點處,具體數(shù)值如表7所示。

圖11 第二輪模型修正后驗分布柱狀圖Fig.11 Posterior distributions after second round calibration

表7 第二輪修正輸入?yún)?shù)修正值

將該組修正值代入有限元模型,計算得到結(jié)構(gòu)溫度云圖如圖12所示。各輸出特征的理論計算值以及有限元仿真值分別列在表8,可以看到平均誤差為4.94%。由圖12、表8可知,第二輪模型修正依然能顯著減小在研模型有限元模型與理論模型的差距,而此時的有限元仿真結(jié)果經(jīng)過虛擬驗證也更具可靠性。

圖12 第二輪修正后末時刻溫度云圖Fig.12 Temperature nephogram at final moment after first round calibration

表8 第二輪修正輸出特征及誤差

4 結(jié) 論

針對缺少試驗數(shù)據(jù)的在研新型號火箭捆綁機構(gòu)有限元模型修正,本文提出考慮建模不確定性的熱力耦合模型修正方法以及基于通用化理論模型的虛擬驗證方法,充分利用前期型號的已有試驗數(shù)據(jù),驗證并修正了在研新型號熱力耦合有限元模型。得到結(jié)論如下:

(1) 本文通過貝葉斯修正框架對前期型號的熱力耦合有限元模型進(jìn)行了模型修正,最終修正后的仿真結(jié)果與實際試驗值相差小,證明了前期型號有限元模型與實際情況貼合較好,仿真值具有較高精度。

(2) 標(biāo)定了通用化熱力耦合理論模型,理論計算值能反映修正后有限元模型的仿真結(jié)果。對于同一系列不同型號的捆綁機構(gòu),該理論模型計算值也具有較高的通用性參考價值。

(3) 基于通用化理論模型,驗證并修正了在研新型號的熱力耦合有限元模型,在不具備試驗數(shù)據(jù)的情況下,保證了仿真結(jié)果的精度,對后續(xù)的型號設(shè)計、結(jié)構(gòu)優(yōu)化等環(huán)節(jié)提供了基礎(chǔ)。