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

(1.空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室,綿陽 621000；2.中國空氣動(dòng)力研究與發(fā)展中心,綿陽 621000)

0 引 言

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

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

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

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

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

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

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

(1)

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

(2)

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

(3)

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

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

(4)

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

(5)

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

(6)

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

(7)

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

(8)

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

(9)

對(duì)式(9)中的擴(kuò)散項(xiàng)運(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)

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

(15)

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

dxdydz=|J|dξdηdζ

(16)

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

(17)

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

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

2 不確定性傳播分析方法

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

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

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

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

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

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

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

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

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

2.1 Kriging代理模型及校驗(yàn)

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

(18)

其中，y(i)為樣本點(diǎn)，數(shù)量為n，給出加權(quán)系數(shù)w=[w(1),w(2),…,w(n)]T的表達(dá)式，即可得到其他位置的預(yù)估值。引入統(tǒng)計(jì)學(xué)假設(shè)：將未知函數(shù)看成是某個(gè)高斯靜態(tài)隨機(jī)過程的具體實(shí)現(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ù)”，距離為零時(shí)取值為1，距離無窮大時(shí)取值為0，相關(guān)性的大小隨距離的增大而減小。基于以上假設(shè)，Kriging模型尋找最優(yōu)加權(quán)系數(shù)w，使得均方差

(21)

最小，并滿足無偏差條件

(22)

(23)

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

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

(24)

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

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

(25)

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

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

(26)

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

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

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

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

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

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

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

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

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

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

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

圖2 計(jì)算條件

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

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

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

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

由表2得到檢驗(yàn)統(tǒng)計(jì)量

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

(27)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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