燒蝕熱響應(yīng)計算中的不確定性傳播分析方法研究

章 超,劉 驍,陳江濤,胡星志,國義軍,崔鵬程

(1.空氣動力學(xué)國家重點實驗室,綿陽 621000；2.中國空氣動力研究與發(fā)展中心,綿陽 621000)

0 引 言

物理問題的數(shù)值模擬存在著大量的不確定性，主要來源于建模和模擬兩個過程。建模的不確定性主要是由于人們對現(xiàn)實世界認(rèn)知的局限性，使得模型不能準(zhǔn)確反映真實的物理過程。模擬的不確定性主要體現(xiàn)在數(shù)值計算過程中沒有準(zhǔn)確的描述模型內(nèi)容，存在輸入?yún)?shù)、網(wǎng)格、計算格式、收斂等因素引入的誤差和不確定性[1-2]。以再入式飛行器的燒蝕熱響應(yīng)計算為例。該類飛行器以高超聲速再入地球時，由于受到嚴(yán)酷的氣動加熱，往往需要引入燒蝕熱防護(hù)系統(tǒng)。而提升熱防護(hù)系統(tǒng)設(shè)計的可靠性和精細(xì)化程度依賴于對再入過程中系統(tǒng)結(jié)構(gòu)內(nèi)部溫度分布的準(zhǔn)確預(yù)測[3]。過去五十年來，與燒蝕材料內(nèi)部熱響應(yīng)計算相關(guān)的數(shù)學(xué)建模和計算方法得到了持續(xù)發(fā)展，但由于熱響應(yīng)計算所使用的材料物性參數(shù)存在大量不確定性，且相互關(guān)系不明，綜合影響難以確定。為準(zhǔn)確獲得不確定性影響下的目標(biāo)變量(例如關(guān)鍵部位溫度)的概率分布，需要一套不確定性傳播分析方法。

在常用的不確定度量化方法中，蒙特卡洛(Monte Carlo，MC)方法[4]構(gòu)造簡單、易于實現(xiàn)，被廣泛使用。但是MC方法收斂速度慢，對計算資源需求巨大。多項式混沌(Polynomial chaos，PC)方法[5]是一種更加高效的不確定度量化方法，通過構(gòu)建目標(biāo)變量和不確定性輸入?yún)?shù)多項式形式的響應(yīng)關(guān)系來分析目標(biāo)的不確定性信息，已經(jīng)有了很多成功的應(yīng)用[6-8]。但PC展開的項數(shù)隨著輸入變量的維數(shù)和展開階數(shù)的增加成幾何級數(shù)式急劇增加，易導(dǎo)致“維數(shù)災(zāi)難”現(xiàn)象產(chǎn)生。

在有限的計算資源下，如何通過規(guī)模大小可以接受的樣本數(shù)據(jù)構(gòu)建相對準(zhǔn)確的響應(yīng)模型，分析不確定性參數(shù)的傳播規(guī)律，是各國學(xué)者們一直努力的方向。在機(jī)器學(xué)習(xí)和優(yōu)化設(shè)計等領(lǐng)域廣泛使用的代理模型為不確定度量化分析提供了新的研究思路。代理模型利用有限的樣本，構(gòu)建系統(tǒng)輸出與輸入變量的響應(yīng)關(guān)系，最終目的之一是代替耗時的數(shù)值模擬或者試驗過程，預(yù)測新輸入變量下的輸出。常用的代理模型包括：多項式回歸模型，人工神經(jīng)網(wǎng)絡(luò)，支持向量機(jī)和高斯過程回歸(Gauss process regression，GPR)等。本文使用在工程領(lǐng)域應(yīng)用廣泛的GPR模型。GPR模型，也被稱為Kriging模型,是由南非采礦工程師Krige[9]于1951年提出，最初應(yīng)用在地質(zhì)儲量估計，Sacks等[10]將其進(jìn)一步推廣到確定性計算機(jī)試驗的設(shè)計和分析領(lǐng)域。GPR模型的構(gòu)造與優(yōu)化是基于貝葉斯理論開展的，對非線性函數(shù)具有良好近似能力和獨特的誤差估計功能。目前GPR模型已經(jīng)被廣泛應(yīng)用到不確定性建模和參數(shù)優(yōu)化的研究中[11-13]。本文發(fā)展了基于GPR模型的不確定度量化方法，并以燒蝕熱響應(yīng)模擬中的不確定性輸入?yún)?shù)量化分析為例，驗證了該方法的有效性和高效率。

1 燒蝕熱響應(yīng)數(shù)值模擬方法

炭化燒蝕材料廣泛應(yīng)用于再入式高超聲速飛行器的熱防護(hù)系統(tǒng)。該類材料在受到外部加熱時，其中的樹脂會發(fā)生熱解反應(yīng)并產(chǎn)生熱解氣體，熱解通常會在一定溫度范圍內(nèi)進(jìn)行，形成熱解區(qū)。材料熱解后留下多孔的炭化層并在表面發(fā)生燒蝕，熱解氣體流經(jīng)炭化層引射到表面上，與表面燒蝕產(chǎn)物一起對氣動加熱起阻塞作用，材料的這種熱解反應(yīng)不僅影響表面燒蝕速度，而且還影響材料內(nèi)部的溫度分布。

假設(shè)熱解型防熱材料由三個組分構(gòu)成,密度分別為ρA,ρB,ρC，體積分?jǐn)?shù)為XA,XB,XC，則固體材料的密度為：

ρs=XAρA+XBρB+XCρC

(1)

每一個組分的熱解過程各由一個化學(xué)動力學(xué)方程式描述：

(2)

假設(shè)完全未熱解的原始材料的密度為ρv，完全炭化后的材料密度為ρc，引入炭化率：

(3)

防熱材料為多孔介質(zhì)，其孔隙率為：

φ=φv(1-α)+φcα

(4)

熱解產(chǎn)生的熱解氣體在多孔介質(zhì)中的流動符合達(dá)西定理：

(5)

其中，Γ為材料的滲透率，υg為熱解氣體的黏度，P為氣體壓力，Vg為熱解氣體流動的速度。熱解氣體的流動滿足動量守恒方程：

(6)

材料內(nèi)部的能量輸運(yùn)過程滿足能量守恒方程：

(7)

本文中燒蝕熱響應(yīng)計算所采用的求解器是基于開源有限元計算框架MOOSE開發(fā)的，以有限元方法作為控制方程的數(shù)值求解方法。接下來以能量方程為例對控制方程的離散做簡要介紹。有限元方法以變分或加權(quán)余量法為基礎(chǔ)，通過分塊逼近思想，對偏微分方程進(jìn)行求解。對于能量守恒方程，先將待求變量T寫為插值函數(shù)與待求系數(shù)乘積的形式

(8)

將式(8)代入到能量守恒方程中，并將式中所有項移動到等號同一端，再乘以檢驗函數(shù)后在控制單元內(nèi)積分，令檢驗函數(shù)等于插值函數(shù)，即Wi=Ni，可得到：

(9)

對式(9)中的擴(kuò)散項運(yùn)用高斯積分公式：

(10)

并將溫度和溫度梯度的表達(dá)式寫成矩陣乘積的形式：

(11)

可得到能量守恒方程的Galerkin弱格式：

(12)

接下來再基于等參變換，構(gòu)造三維插值函數(shù)。通過坐標(biāo)變換，將任意全局坐標(biāo)系下的非標(biāo)準(zhǔn)單元變換為局部坐標(biāo)下的標(biāo)準(zhǔn)單元，為此需要建立坐標(biāo)轉(zhuǎn)換關(guān)系式：

(13)

根據(jù)等參變換，待求變量與坐標(biāo)變換使用同樣的插值函數(shù)，以變量T為例，在局部坐標(biāo)系下：

(14)

本文采用拉格朗日插值多項式來構(gòu)造插值函數(shù)。插值函數(shù)的形式確定之后，可將全局坐標(biāo)下的單元有限元方程變換到局部坐標(biāo)系下。先寫出插值函數(shù)的變換式：

(15)

同時，全局坐標(biāo)下的積分項dxdydz變換為：

dxdydz=|J|dξdηdζ

(16)

其中，|J|為雅可比矩陣的行列式。最后可得到單元有限元方程在局部坐標(biāo)下的表達(dá)式：

(17)

對于動量方程和質(zhì)量守恒方程，可根據(jù)同樣的方式進(jìn)行離散。獲得單元有限元方程后，將其中的矩陣按照一定的規(guī)則放入到總體矩陣，可得到總體方程。對上述控制方程組，最終獲得的總體方程為瞬態(tài)非線性方程組。本文對其在時間方向采用二階Euler后差(BDF2)進(jìn)行時間積分，并結(jié)合Newton-GMRES迭代對所得到的大型稀疏矩陣進(jìn)行求解。上述計算方法及求解器的驗證與考核詳見文獻(xiàn)[3]。

在燒蝕熱響應(yīng)數(shù)值模擬中包含了模型構(gòu)造、不確定性輸入?yún)?shù)、數(shù)值求解過程等諸多不確定性因素。本文主要分析與材料物性相關(guān)的輸入?yún)?shù)的不確定性傳播。

2 不確定性傳播分析方法

針對不確定輸入?yún)?shù)的不確定性傳播問題，本文發(fā)展了基于代理模型的分析流程，包括以下步驟：

(1)不確定性因素的量化；

(2)通過確定性計算生成樣本；

(3)構(gòu)建代理模型；

隨著數(shù)碼技術(shù)和互聯(lián)網(wǎng)應(yīng)用的進(jìn)步，攝影已深深的融入社會生活，數(shù)碼音樂播放器，移動電話等數(shù)碼產(chǎn)品開始配備攝影功能，拍攝的圖片可以無線傳播，攝影開始多元化發(fā)展。

(4)代理模型校驗及可能的模型參數(shù)優(yōu)化調(diào)整；

(5)敏感性分析及輸出響應(yīng)統(tǒng)計特性。

流程如圖1所示，對模型的精度有要求時可以迭代步驟(2)～(4)。

第2.1節(jié)和2.2節(jié)將簡單介紹使用的Kriging代理模型以及模型校驗與敏感性分析方法，結(jié)合概率統(tǒng)計方法，得到輸出響應(yīng)的統(tǒng)計特性。

2.1 Kriging代理模型及校驗

Kriging代理模型是一種插值模型，插值結(jié)果是已知樣本的線性加權(quán)，

(18)

其中，y(i)為樣本點，數(shù)量為n，給出加權(quán)系數(shù)w=[w(1),w(2),…,w(n)]T的表達(dá)式，即可得到其他位置的預(yù)估值。引入統(tǒng)計學(xué)假設(shè)：將未知函數(shù)看成是某個高斯靜態(tài)隨機(jī)過程的具體實現(xiàn)。該靜態(tài)隨機(jī)過程定義為

Y(x)=β0+Z(x)

(19)

其中，β0為未知常數(shù)，也稱為全局趨勢模型，Y(x)表示數(shù)學(xué)期望值；Z(x)表示的是均值為零、方差為σ2(σ2(x)≡σ2,?x)的靜態(tài)隨機(jī)過程。在參數(shù)空間的不同位置，這些隨機(jī)變量存在一定的相關(guān)性，通過協(xié)方差表示為：

Cov[Z(x),Z′(x)]=σ2R(x,x′)

(20)

式中：R(x,x′)為只與空間距離相關(guān)的“相關(guān)函數(shù)”，距離為零時取值為1，距離無窮大時取值為0，相關(guān)性的大小隨距離的增大而減小?；谝陨霞僭O(shè)，Kriging模型尋找最優(yōu)加權(quán)系數(shù)w，使得均方差

(21)

最小，并滿足無偏差條件

(22)

(23)

式(23)中，R是由相關(guān)函數(shù)組成的“相關(guān)矩陣”，r為未知點與已知樣本點之間相關(guān)函數(shù)組成的“相關(guān)矩陣”，F(xiàn)是全1的列矩陣，全局趨勢模型的表達(dá)式為β0=(FTR-1F)-1FTR-1ys。相關(guān)函數(shù)的選擇及模型參數(shù)的訓(xùn)練可以參考文獻(xiàn)[14]。

對代理模型的校驗可以采用留出法或交叉驗證法，模型精度可以采用R2指標(biāo)來評判，

(24)

如果存在額外的數(shù)據(jù)量足夠的校驗數(shù)據(jù)，可以采用分布比較的方法。假設(shè)(x,yv)是校驗數(shù)據(jù)，相應(yīng)的代理模型響應(yīng)是(x,ySM)，定義均值指標(biāo)dE測量它們之間的差異，

dE(yv,ySM)=E(|yv-ySM|)

(25)

其中，E為期望算子。還可以通過將校驗數(shù)據(jù)和代理模型響應(yīng)值表示成累計概率分布的形式，直觀地顯示出兩者差異。

2.2 輸出響應(yīng)的統(tǒng)計特性和敏感性分析

(26)

式中:z為標(biāo)準(zhǔn)正態(tài)分布分位數(shù)，s為樣本標(biāo)準(zhǔn)偏差，置信水平表示為100(1-α)。α通常取值0.1或者0.05。

敏感性分析是不確定度量化傳播的一個重要步驟，其中全局敏感性分析能夠衡量模型輸入變量的不確定性在其整個分布域內(nèi)對輸出響應(yīng)的不確定性的影響[15]。本文基于代理模型對輸入不確定性參數(shù)進(jìn)行敏感性分析?；舅悸啡缦拢?/p>

1)重采樣生成大量的樣本;

2)基于代理模型快速計算輸出響應(yīng)值；

3)基于敏感性分析方法分析每個輸入?yún)?shù)對輸出響應(yīng)的貢獻(xiàn)。

本文使用了Sobol指標(biāo)[16]來量化每個輸入變量對系統(tǒng)輸出方差的貢獻(xiàn)程度。Sobol指標(biāo)通過Sobol分解得到，是一種基于方差分解的全局敏感性分析方法，其優(yōu)點是考慮了隨機(jī)輸入?yún)?shù)在整個取值范圍內(nèi)對輸出響應(yīng)的貢獻(xiàn)以及隨機(jī)輸入?yún)?shù)的交互作用。

許多時候，工程上對于輸出響應(yīng)符合某種特定的分布存在期待，可以引入假設(shè)檢驗中的擬合優(yōu)度檢驗來解決這個問題。常用的擬合優(yōu)度檢驗方法有χ2檢驗、柯式檢驗法等。

至此給出了基于代理模型的不確定度量化的理論基礎(chǔ)和主要過程。接下來，將以一維燒蝕熱響應(yīng)問題為例，分析其不確定性的傳播過程。

3 燒蝕熱響應(yīng)數(shù)值模擬中的不確定性傳播分析

3.1 問題描述及樣本數(shù)據(jù)處理

本節(jié)選取了一維燒蝕熱響應(yīng)算例，算例來自4th Abaltion Workshop[17]官方網(wǎng)站提供的標(biāo)準(zhǔn)算例。本算例計算所使用的材料，TACOT，是Abaltion Workshop主辦方為方便數(shù)據(jù)對比而提供的一種假想材料，其物性參數(shù)的不確定度無法通過試驗手段獲得，計算條件如圖2所示。

圖2 計算條件

選擇25 mm厚度位置的溫度值作為目標(biāo)輸出，引入9個材料物性參數(shù)作為不確定性輸入，如表1所示。

表1 輸入?yún)?shù)及其不確定度范圍

圖3 確定性數(shù)值仿真結(jié)果直方圖

從圖3可以看出，30個樣本點的數(shù)據(jù)呈現(xiàn)較大的散布，分布無明顯規(guī)律，500個樣本點的數(shù)據(jù)呈現(xiàn)出較強(qiáng)的正態(tài)分布特征，對500個樣本數(shù)據(jù)進(jìn)行擬合優(yōu)度檢驗。

由表2得到檢驗統(tǒng)計量

表2 擬合優(yōu)度檢驗數(shù)據(jù)表

(27)

從結(jié)果對比可以發(fā)現(xiàn)，30個樣本數(shù)據(jù)和500個樣本數(shù)據(jù)的統(tǒng)計特性相差較大。需要指出的是這里并不說明500個樣本數(shù)據(jù)的結(jié)果一定是更準(zhǔn)確的，但充分說明MC方法統(tǒng)計結(jié)果非常依賴于樣本數(shù)據(jù)的規(guī)模。在工程問題中，計算資源有限的情況下需要更加高效的不確定度量化方法。

3.2 代理模型構(gòu)建及校驗

本文基于GPR模型利用30個樣本數(shù)據(jù)構(gòu)建代理模型，模型基函數(shù)為常數(shù)，核函數(shù)為“Matern 5/2”，在使用代理模型進(jìn)行不確定性分析之前，首先驗證代理模型的準(zhǔn)確度，由第2.1節(jié)所述，采用兩種方式。

1)30樣本點自身的交叉驗證；采用留一交叉驗證，R2=0.99。顯示泛化誤差很小，模型的準(zhǔn)確度很高。

代理模型響應(yīng)與確定性計算結(jié)果的差值分布如圖5所示，90%的校驗數(shù)據(jù)與代理模型輸出響應(yīng)的差值在1.5 K以內(nèi)。圖6給出了代理模型響應(yīng)輸出與確定性計算結(jié)果的累計概率分布比較，可以看出差別較大的輸出結(jié)果主要集中在輸出溫度小于350 K的區(qū)域。綜合圖4～6的校驗結(jié)果，30個樣本生成的GPR代理模型具有很高的精度，能夠很好地代替數(shù)值模擬過程。

圖4 代理模型響應(yīng)與確定性計算結(jié)果比較

圖5 代理模型響應(yīng)與確定性計算結(jié)果差值

圖6 代理模型響應(yīng)與確定性計算累計分布比較

3.3 敏感性分析及輸入?yún)?shù)優(yōu)化

采用拉丁超立方[18](Latin hypercube sampling，LHS)采樣10000次，采用代理模型進(jìn)行快速計算，再利用Sobol指標(biāo)進(jìn)行敏感性分析。從表3可以看出,對輸出響應(yīng)影響最大的是輸入?yún)?shù)1和5，參數(shù)2,3影響明顯，參數(shù)4和7有一定的影響，參數(shù)6,8,9對輸出響應(yīng)完全沒有影響。

表3 不同物性參數(shù)Sobol指標(biāo)

在本文的一維燒蝕熱響應(yīng)問題中，若將9個參數(shù)優(yōu)化為6個參數(shù)(參數(shù)1，2，3，4，5，7)，對計算結(jié)果沒有影響。若進(jìn)一步優(yōu)化為4個(參數(shù)1，2，3，5)，采用同樣的30個樣本生成新的代理模型，R2=0.95，在500個確定性計算中的校驗中，均值誤差為0.03%，期望指標(biāo)dE=6.00 K,為均值的1.6%。減少不確定性參量的個數(shù)，可以減少構(gòu)建代理模型所需要的樣本數(shù)量，從而減少確定性計算的耗費(fèi)。通過對一維燒蝕熱響應(yīng)輸入?yún)?shù)的敏感性分析可以得到對此類問題影響最大的輸入?yún)?shù)為原始材料導(dǎo)熱率、炭化材料導(dǎo)熱率、原始材料比熱和原始材料密度。需要說明的是，實際上孔隙率對于材料熱導(dǎo)率和密度具有關(guān)聯(lián)關(guān)系，在本例中不顯著可能和響應(yīng)溫度位置、參數(shù)不確定度的范圍等因素有關(guān)。這也說明，代理模型的快速計算分析并不考慮物理本質(zhì)，是一種高效的替代方法。

3.4 輸出響應(yīng)統(tǒng)計特性

圖7給出了輸出響應(yīng)的直方圖，與500個樣本數(shù)據(jù)結(jié)果也很接近。因此可以認(rèn)為由30個樣本數(shù)據(jù)構(gòu)建的代理模型給出的目標(biāo)變量不確定性統(tǒng)計信息與500個樣本數(shù)據(jù)的結(jié)果基本一致，但計算量大大減少。由式(26)可以給出響應(yīng)均值置信水平為99%的置信區(qū)間為(355.64 K，356.08 K)。圖8給出了輸出響應(yīng)的概率圖，從圖中可以看出置信水平為90%時的置信區(qū)間為(341.78 K，370.62 K)，區(qū)間寬度占比為8%。

圖7 輸出響應(yīng)分布直方圖

圖8 輸出響應(yīng)概率圖

4 結(jié) 論

針對數(shù)值模擬中存在較多不確定性因素的問題，本文提出了研究不確定性傳播的流程。以一維熱燒蝕響應(yīng)數(shù)值模擬中物性參數(shù)的不確定性研究為例，給出了完整的針對多個輸入?yún)?shù)的不確定傳播分析方法，識別出了對輸出響應(yīng)影響較大的不確定性參數(shù)，給出了輸出響應(yīng)的統(tǒng)計特性。

這一不確定性分析流程的關(guān)鍵在于高效地獲得高準(zhǔn)確度的代理模型，這對試驗設(shè)計、代理模型的選擇和優(yōu)化提出了較高要求，是后續(xù)研究的重點。另外，本文忽略了除輸入?yún)?shù)不確定性以外的其他不確定性，也沒有考慮輸入?yún)?shù)間的關(guān)聯(lián)關(guān)系，如何統(tǒng)籌考慮各種不確定性因素以及相互關(guān)系對目標(biāo)輸出的影響，也是后續(xù)研究的重點之一。