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

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

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

在重型運載火箭的飛行過程中,球頭-球窩式捆綁機構為了適應芯級與助推器之間的相對彈性變形會發生轉動摩擦,大噸位載荷作用下的轉動摩擦會產生大量熱量導致捆綁機構溫度過高,從而嚴重影響整體強度甚至導致結構失效[1],因此需要精確的熱力耦合有限元分析來輔助結構設計以及結構優化。

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

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

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

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

1 模型修正算法

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

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

1.1 貝葉斯修正框架

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

(1)

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

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

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

(2)

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

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

(3)

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

(4)

式中:σ為寬度系數,一般可取為[10-3,10-1][19];d為統計學距離指標,用來衡量預測樣本與觀測樣本的差異。本文選用歐氏距離作為距離度量指標,其表達式為[20]

(5)

1.2 二次多項式響應面

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

(6)

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

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

(7)

2 虛擬驗證方法

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

2.1 虛擬驗證框架

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

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

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

2.2 通用化理論模型

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

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

(8)

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

(9)

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

(10)

式中,λ為接觸面積修正系數。在此基礎上,假設等效溫升區域如圖2中陰影部分,根據熱力學定律以及集中參數法可以推導出該區域的溫度變化公式,表示為

(11)

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

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

3 數值算例

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

3.1 捆綁機構熱力耦合有限元模型建立

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

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

3.2 第一輪模型修正

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

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

(12)

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

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

圖5 導熱率溫變模型Fig.5 The temperature variation model of thermal conductivity

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

表1 第一輪模型修正的不確定參數

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

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

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

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

代理模型構建完畢后,根據1.1節中介紹的貝葉斯修正框架進行模型修正。設置TMCMC樣本數為2 000,經歷7次迭代后得到各輸入參數的后驗分布如圖7所示。圖7中的9幅圖為每個待修正參數2 000個樣本的后驗分布柱狀圖,區間對應的值越大,說明該區間上的樣本點越多。可以發現,經過修正后的后驗分布非常集中,絕大多數樣本點都聚集在極小的區間內,因此各輸入參數的修正值可以簡單地取為相應區間的中點,即圖7中的“*”號標記點處,具體數值如表3所示。

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

表3 第一輪模型修正輸入參數修正值Tab.3 Updated input parameters after first round calibration

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

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

對修正后的輸入參數進行有限元仿真,得到結構末時刻的溫度場云圖如圖8所示。有限元仿真可以得到結構整體的溫度分布以及最大溫度,這些數據在實際試驗中難以全面測量,但模型修正保證了這些仿真值的正確性,使其能夠在后續用于理論模型標定。

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

3.3 理論模型標定

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

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

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

3.4 第二輪模型修正

第二輪模型修正的對象是在研新型號捆綁機構,其有限元模型基本與前期型號一致,但材料改變為沉淀硬化不銹鋼PH13-8Mo,且其他個別參數更改為:半徑差Δr=0.7 mm,球頭轉動頻率0.8 Hz,外載荷上升至2.4×106N。由于在研型號不具備試驗數據,因此將理論模型計算值作為觀測值對在研型號有限元模型進行驗證及修正。

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

表5 在研型號模型輸入參數Tab.5 Input parameters of new developing model

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

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

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

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

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

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

表7 第二輪修正輸入參數修正值

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

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

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

4 結 論

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

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

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

(3) 基于通用化理論模型,驗證并修正了在研新型號的熱力耦合有限元模型,在不具備試驗數據的情況下,保證了仿真結果的精度,對后續的型號設計、結構優化等環節提供了基礎。