劉 健， 張 琨， 李 滕， 辛健強， 姚建堯*,3

(1.重慶大學(xué) 航空航天學(xué)院，重慶 400044；2.中國航天科技集團公司第一研究院 系統(tǒng)工程業(yè)務(wù)部，北京 100076；3.重慶大學(xué) 非均質(zhì)材料力學(xué)重慶市重點實驗室，重慶 400044)

1 引 言

功能梯度材料FGM(functional graded material)是一種由兩種或多種成分組合而成的新型復(fù)合材料。通過改變材料的相對體積分數(shù)、物理狀態(tài)和幾何構(gòu)型，功能梯度材料可以表現(xiàn)出材料性質(zhì)的連續(xù)空間變化[1]。

在功能梯度材料的設(shè)計過程中，最為關(guān)鍵的問題之一是如何描述混合后結(jié)構(gòu)的宏觀屬性。Wakashima等[2]提出了一個自洽的平均場微觀力學(xué)W-T(Wakashima-Tsukamoto)模型，Cho等[3]證明了該模型能得到更準(zhǔn)確的結(jié)果。Kapuria等[4]結(jié)合該模型并通過實驗驗證了該模型的有效性。

已經(jīng)證明FGM具有優(yōu)化溫度場、降低熱應(yīng)力和提高熱阻的能力[5]，但在高溫高壓的工作環(huán)境中，功能梯度材料會出現(xiàn)熱應(yīng)力失配、界面破壞和熱腐蝕退化等失效形式。因此，功能梯度材料的熱傳導(dǎo)和熱力耦合問題是學(xué)者的研究重點，包括數(shù)值算法和解析算法。Chen等[6]用有限元分析了FGM穩(wěn)態(tài)和瞬態(tài)熱傳導(dǎo)問題的靈敏度，給出了響應(yīng)相對于設(shè)計變量的變化率，可以用于FGM的優(yōu)化問題。Sankar等[7]獲得了材料特性隨厚度呈指數(shù)變化的FGM梁中熱應(yīng)力分布的精確解。

然而，功能梯度材料在設(shè)計和加工過程中存在許多不確定性和隨機性[8，9]。概率分析方法是一種考慮熱不確定性的有效方法。屈強等[10]結(jié)合有限元和蒙特卡洛模擬建立了剛性陶瓷瓦熱防護系統(tǒng)概率設(shè)計方法。Zhang等[11]利用蒙特卡洛模擬研究了典型的多層熱防護結(jié)構(gòu)的概率溫度響應(yīng)。辛健強等[12]針對輻射式熱防護系統(tǒng)結(jié)合蒙特卡洛直接抽樣建立了一種概率瞬態(tài)熱分析方法。實踐證明，在熱分析中，蒙特卡洛模擬是功能梯度材料考慮參數(shù)不確定性的一種準(zhǔn)確有效的方法。近年來，代理模型已取代傳統(tǒng)的數(shù)值模擬廣泛應(yīng)用于概率熱傳導(dǎo)和熱應(yīng)力分析。Xin等[13]提出一種結(jié)合拉丁超立方抽樣和響應(yīng)面法的多層熱防護結(jié)構(gòu)概率分析方法。

為確保功能梯度材料在耐熱結(jié)構(gòu)中的可靠性和安全性，對功能梯度材料進行熱力耦合概率分析是很有必要的。Sugano等[14]用解析法研究了初始溫度隨機變化的FGM板的溫度和熱應(yīng)力統(tǒng)計。Chiba等[15]分析了導(dǎo)熱系數(shù)和熱膨脹系數(shù)不確定性對溫度和熱應(yīng)力場的影響。Yang等[16]用幾個獨立的隨機變量對FGM板進行了彈性屈曲分析。Ferrante等[17]研究了組成成分的體積分數(shù)和孔隙率的隨機變化對功能梯度金屬/陶瓷板在恒定熱邊界條件下響應(yīng)的影響。Chiba等[18]分析了用隨機過程表示的熱載荷作用下功能梯度板的熱彈性問題。然而，上述不確定性分析大多都是針對模型某一特定的隨機變量，很少有考慮兩個以上的隨機變量對FGM熱響應(yīng)的影響。而在實際工程問題中，多種不確定性是同時存在的，同時熱響應(yīng)具有時間歷程的性質(zhì)，其輸出是一個多輸出問題。因此，建立一種能夠考慮多種不確定性及多輸出的熱力耦合概率設(shè)計分析方法對功能梯度材料的設(shè)計尤為重要。

當(dāng)只對單一輸出或點進行預(yù)測時，代理模型能夠有比較好的預(yù)測效果。而對于一些時間歷程多輸出問題，則需要對每個時刻上的輸出(設(shè)時間維度為N)進行預(yù)測。為提高代理模型的效率和精度，本文采用本征正交分解方法[19]對輸出變量進行減縮。本征正交分解是一種將高維數(shù)據(jù)投影到低維空間以實現(xiàn)數(shù)據(jù)降階的方法。在構(gòu)建代理模型時，只需對包含絕大部分能量的前幾階正交基(模態(tài))的投影系數(shù)進行預(yù)測，從而減少計算量[20]，并能夠提高模型的預(yù)測精度和穩(wěn)定性。

本文基于W-T模型建立了材料屬性與溫度相關(guān)的參數(shù)化有限元模型，給出了功能梯度材料空間分布的方式和整體性能估計的方法。將隨機參數(shù)作為輸入對模型進行熱力耦合分析，以溫度和熱應(yīng)力作輸出，建立了功能梯度材料的概率分析方法。驗證了Kriging代理模型的可行性和準(zhǔn)確性，以及POD降階方法對解決時間歷程多輸出問題的有效性，并分析了不確定性參數(shù)的影響及其靈敏度。最后討論了該方法的有效性。

2 功能梯度材料確定性建模方法

首先給出材料成分在空間的分布和材料屬性隨溫度變化的函數(shù)，并構(gòu)建了基于W-T模型的金屬基陶瓷功能梯度材料模型，用于考慮材料成分的變化和結(jié)構(gòu)整體性能的描述。

2.1 組成成分的空間分布

假設(shè)陶瓷體積分數(shù)和金屬體積分數(shù)由冪關(guān)系給出,

Vc=(z/t)n,Vm=1-Vc

(1,2)

式中下標(biāo)m和c分別表示金屬和陶瓷成分，z為沿厚度方向的坐標(biāo)，t為平板厚度，n為冪指數(shù)，n≥0。這個冪指數(shù)關(guān)系描述了平板是如何通過厚度分級的。圖1給出了相對于不同冪指數(shù)值的陶瓷體積分數(shù)Vc。

圖1 陶瓷體積分數(shù)在板厚方向上的變化

2.2 溫度對材料屬性的影響

ZrO2:

(3)

Ti-6Al-4V:

(4)

材料參數(shù)對溫度的敏感程度和變化規(guī)律不同，為保證結(jié)果的可靠性，準(zhǔn)確定義材料與溫度相關(guān)的函數(shù)在高溫?zé)釓椥匝芯恐酗@得尤為重要。

2.3 基于W-T模型的性能描述

在Mori-Tanaka模型的基礎(chǔ)上，Wakashima等[2]提出的W-T模型廣泛應(yīng)用在功能梯度材料宏觀彈性參數(shù)預(yù)估中。它采用了一種平均場的方法，基于基體中的平均應(yīng)力的概念，近似估計組成成分的相互作用。導(dǎo)出了有效體積模量K和剪切模量G與熱導(dǎo)率k、比熱C和熱膨脹系數(shù)α之間的關(guān)系,

(5)

(6)

式中

(7)

(8)

熱導(dǎo)率k、比熱C和熱膨脹系數(shù)α的表達式為

(9)

C=C1V1+C2V2

(10)

(11)

對于梯度層的區(qū)域來說，某種材料可以作為基體或夾雜物，這就會產(chǎn)生兩個不同的結(jié)果，這兩個結(jié)果與HS模型[21]的上下界相似。采取如式(12)所示的基于體積的加權(quán)平均來估計，這種平均方法能夠在兩個界限之間平滑過渡。

K=KW T 1V1+KW T 2V2

(12)

式中KW T 1是將材料1作為基體代入式(5)得到的結(jié)果，KW T 2是將材料2作為基體代入式(5)得到的結(jié)果。剪切模量G、熱導(dǎo)率k、比熱C和熱膨脹系數(shù)α可通過同樣的方法計算。

3 功能梯度材料熱力耦合概率分析

本文考慮熱導(dǎo)率、比熱及空間分布等13個不確定性參數(shù)對熱應(yīng)力的影響，熱力耦合概率分析流程如圖2所示，步驟如下。

(1) 識別輸入變量及其不確定性，包括溫度相關(guān)的材料屬性和空間分布。

(2) 采用拉丁超立方抽樣方法生成獨立的輸入變量。利用參數(shù)化有限元模型和熱載荷條件，計算功能梯度平板的瞬態(tài)熱特性，得到其溫度分布。

(3) 將每個時刻溫度場以載荷的形式施加到結(jié)構(gòu)上進行瞬態(tài)結(jié)構(gòu)分析，輸出每個時刻的最大熱應(yīng)力。

(4) 通過輸入變量和輸出的熱應(yīng)力構(gòu)建Kriging代理模型。

(5) 利用Kriging代理模型進行蒙特卡洛模擬，得到功能梯度平板的熱應(yīng)力概率特性，并獲得每個輸入變量對輸出結(jié)果的靈敏度。

圖2 FGM熱力耦合概率分析方法流程

analysis method flow

3.1 不確定性有限元模型

本算例中，平板的尺寸大小為30 mm×16 mm，保持底面恒定300 K，上表面施加熱流載荷q=300 kW/m2，左右邊界保持絕熱。首先進行一次確定性熱分析，結(jié)果表明，功能梯度平板在300 s之前就已達到穩(wěn)定，上表面最高溫度為1281.84 K，如圖3所示。在瞬態(tài)熱分析的基礎(chǔ)上進行瞬態(tài)結(jié)構(gòu)分析，設(shè)置下表面簡支的邊界條件，輸出每個時刻最大von Mises stress。

根據(jù)式(3,4)，將材料屬性隨溫度變化的函數(shù)中的常數(shù)項作為隨機參數(shù)輸入變量，隨機參數(shù)的離散程度通過變異系數(shù)(COV)來控制，其定義為變量標(biāo)準(zhǔn)差與均值之比。根據(jù)文獻[13]的研究，材料屬性的變異系數(shù)設(shè)置為0.02，另一隨機輸入變量為材料的空間分布冪指數(shù)n，變異系數(shù)為0.01，不確定性輸入?yún)?shù)列入表1。

3.2 POD降階的Kriging模型及其精度分析

由于熱力耦合分析的輸出是時間歷程上的熱應(yīng)力，即多輸出問題。因此每一個樣本的輸出為一向量，總輸出為一矩陣。樣本矩陣可表示為

(13)

式中Pi j為來自第i個訓(xùn)練樣本的材料屬性和空間分布的第j個隨機輸入?yún)?shù)，子矩陣Q為訓(xùn)練樣本的輸出，Qi j為第i個訓(xùn)練樣本的第j時刻的最大熱應(yīng)力。m為隨機參數(shù)的個數(shù)，b為訓(xùn)練樣本的個數(shù)，nt為時間采集步數(shù)。因此，矩陣S的每一行代表一個訓(xùn)練樣本。

使用傳統(tǒng)的Kriging模型預(yù)測時間歷程上的多輸出問題時，無法保證熱應(yīng)力在時域上的連續(xù)性。針對這一問題，本文首先通過POD方法提取時間歷程上熱應(yīng)力的特征信息，從而降低待預(yù)測數(shù)據(jù)的維度，進而構(gòu)建Kriging代理模型。

3.2.1 本征正交分解

通過POD將大量相互依賴的變量減少到少量的不相關(guān)變量，同時盡可能多地保留原始變量的變化。POD對輸出矩陣Q降階過程如下。

令矩陣V=QT,V的第i列向量vi為第i組參數(shù)下熱應(yīng)力在時間歷程上的變化。熱應(yīng)力隨參數(shù)和時間變化的關(guān)系可以由式(14)逼近，

(14)

(15)

式中(·,·)為內(nèi)積運算。Φ為向量集合V張成空間的一組規(guī)范正交基，那么這組基中任何一個基都可以用原始向量之間的線性疊加來表示，即

(16)

特征值的大小表征了該特征值對應(yīng)的POD基所包含矩陣V特征的多少。一般情況下，往往前幾個POD基就能包含90%甚至99%以上的廣義能。因此，對于本文多輸出問題，只需要從中抽取一定數(shù)量快照進行本征正交分解，分析少數(shù)廣義能比重較大的POD基就能得到系統(tǒng)的主導(dǎo)特征。

取Φ的前J基向量按列歸一化得到Ψ，降階系數(shù)可通過式(17)求得

A=ΨTV

(17)

因此通過POD降階后，式(13)的形式可轉(zhuǎn)換為

(18)

式中J?nt，即輸出矩陣經(jīng)過POD降階后維度大大降低。不僅減少了后續(xù)Kriging模型的構(gòu)建過程，降階后的輸出矩陣還包含了各時間歷程輸出之間的關(guān)系，這對與時間相關(guān)的結(jié)果預(yù)測很重要。

3.2.1 代理模型構(gòu)建及精度分析

通過計算得到500個樣本構(gòu)建降階的Kriging模型，流程如圖4所示。

對比POD降階的Kriging模型和傳統(tǒng)的Kriging模型的結(jié)果，如圖5所示?？梢钥闯?，POD降階的Kriging模型預(yù)測時間歷程上的多輸出問題效果相較傳統(tǒng)的要好。傳統(tǒng)的Kriging導(dǎo)致了預(yù)測結(jié)果存在明顯的非正常波動。而通過POD降階后的模型很好地解決了這一問題，降階后的矩陣提取了輸出矩陣的絕大部分信息，考慮了相鄰向量之間的關(guān)系，從而能夠與原始數(shù)據(jù)吻合得比較好。

圖4 降階的Kriging模型構(gòu)建流程

圖5 降階的Kriging模型與Kriging及真實值的結(jié)果對比

with Kriging and original values

為了評價代理模型的精度，定義歸一化均方根誤差NRMSE為

(19)

式中s為樣本數(shù)，PV為預(yù)測值，RV為真實值。

如圖6所示，當(dāng)訓(xùn)練樣本達到360時，模型開始收斂。通過構(gòu)建的Kriging模型對整個加載時間歷程中的熱應(yīng)力進行預(yù)測，結(jié)果如圖7所示。

3.3 熱力耦合概率分析

由于每個樣本在時間歷程上均有一個多輸出結(jié)果，在應(yīng)力-時間圖中表現(xiàn)為每一條曲線代表一個樣本結(jié)果。所有樣本各個時刻最大熱應(yīng)力值如圖8所示，圖中粗曲線為該樣本的最大和最小熱應(yīng)力值，星條線為確定性結(jié)構(gòu)分析的結(jié)果，其他樣本結(jié)果均分布在這兩條極值曲線中間?？梢钥闯?，局部最大熱應(yīng)力的最大值和最小值相差約60 MPa，而穩(wěn)定時達到了近100 MPa的偏差。由此可知，輸入?yún)?shù)的微小隨機性能會引起較大的熱應(yīng)力偏差，在工程應(yīng)用中尤其需要重點關(guān)注不確定性對熱應(yīng)力的影響。

圖6 NRMSE隨著訓(xùn)練樣本的變化

圖7 預(yù)測熱應(yīng)力和原始熱應(yīng)力的對比

圖8 500個樣本各時刻最大熱應(yīng)力

同樣地，為了研究不確定性輸入?yún)?shù)對熱應(yīng)力的影響，分別計算每個隨機參數(shù)與輸出的相關(guān)系數(shù)，如圖9所示。可以看出，在整個時間歷程中，對最大熱應(yīng)力響應(yīng)影響最大的是K1和α1，即ZrO2的熱導(dǎo)率和熱膨脹系數(shù)，其中ZrO2的熱膨脹系數(shù)對熱應(yīng)力的影響是正相關(guān)的，而熱導(dǎo)率的影響則是負相關(guān)的。其次則為E1,α2和n，其他變量都集中在0附近，對結(jié)果影響較小。

圖9 隨機輸入變量的靈敏度

4 結(jié) 論

根據(jù)本文研究，可以得到以下結(jié)論。

(1) 本文提出的概率分析方法能夠同時考慮多種隨機因素對功能梯度材料熱響應(yīng)的影響，并驗證了該方法的準(zhǔn)確性和有效性，給功能梯度材料的結(jié)構(gòu)設(shè)計和靈敏度分析提供一定指導(dǎo)。

(2) 通過靈敏度分析可知，對于各時刻的最大熱應(yīng)力，影響最大的是陶瓷材料的熱導(dǎo)率和熱膨脹系數(shù)，其次為金屬材料的熱膨脹系數(shù)和陶瓷材料的彈性模量以及空間分布冪系數(shù)。因此，在功能梯度材料的工程應(yīng)用中需要重點關(guān)注材料參數(shù)以及材料空間分布的不確定性對結(jié)果的影響。

(3) 基于POD的Kriging模型對多輸出的結(jié)構(gòu)響應(yīng)的預(yù)測具有較好的結(jié)果。POD降階法對時間歷程上的多輸出問題能夠大大減少計算量，優(yōu)化預(yù)測結(jié)果。