李寶玉，張磊剛，裘群海，余雄慶

1. 南京航空航天大學(xué) 航空宇航學(xué)院，南京 210016 2. 中國運(yùn)載火箭技術(shù)研究院，北京 100076

一次二階矩(First Order and Second Moment, FOSM)方法是結(jié)構(gòu)可靠性分析中一種簡單易操作的方法，作為一種實(shí)用方法在工程問題中得到了廣泛應(yīng)用[1-3]。FOSM方法的基本思想就是將非線性的功能函數(shù)進(jìn)行線性化，然后通過輸入變量的一階矩和二階矩來計算線性化后的功能函數(shù)的一階矩和二階矩，使用其求得結(jié)構(gòu)可靠度指標(biāo)，進(jìn)而近似得到功能函數(shù)的失效概率。隨著方法精度的研究與發(fā)展，F(xiàn)OSM方法包括均值一次二階矩(Mean Value FOSM，MVFOSM)方法和改進(jìn)一次二階矩(Advanced FOSM, AFOSM)方法。MVFOSM方法和AFOSM方法的區(qū)別在于二者線性化的點(diǎn)是不同的，前者是在輸入變量的均值點(diǎn)處對功能函數(shù)進(jìn)行線性化，而后者則是在對失效概率貢獻(xiàn)最大的點(diǎn)，即最可能失效點(diǎn)—設(shè)計點(diǎn)處線性化。由于設(shè)計點(diǎn)是對失效概率貢獻(xiàn)最大的點(diǎn)，因此在設(shè)計點(diǎn)處線性展開比在均值點(diǎn)處線性展開對失效概率的近似具有更高的精度[3]。

值得注意的是，不論是AFOSM方法還是MVFOSM方法，對于非線性功能函數(shù)問題，2種方法均需要將非線性功能函數(shù)進(jìn)行線性展開，因此均需要求解功能函數(shù)的導(dǎo)函數(shù)，獲取功能函數(shù)對各輸入變量的梯度信息。對于顯式功能函數(shù)，其導(dǎo)函數(shù)求解較為容易，而對于復(fù)雜工程結(jié)構(gòu)可靠性分析中常見的隱式功能函數(shù)問題，其導(dǎo)函數(shù)則較難求解。有限差分法是求解隱函數(shù)的導(dǎo)函數(shù)的通用方法，但對于高度非線性問題，有限差分法的步長較難確定，而且計算梯度信息所花費(fèi)的計算量相當(dāng)大，特別地，對于以有限元形式表征的功能函數(shù)而言計算代價在工程問題中難以接受。

鑒于AFOSM方法的工程適用性，同時考慮到Kriging代理模型技術(shù)在結(jié)構(gòu)可靠性工程及優(yōu)化領(lǐng)域的應(yīng)用廣泛性，本文研究并提出了一種基于Kriging代理模型的AFOSM可靠性分析方法。在眾多代理模型中，Kriging模型是一種應(yīng)用較為廣泛且精度較高的代理模型[4-9]。例如：Simpson等[6]將Kriging與響應(yīng)面法的計算精度與計算效率進(jìn)行了對比，并將其應(yīng)用于航天飛機(jī)的設(shè)計中；Sakata等[7]改進(jìn)和發(fā)展了抽樣方法，并將Kriging模型應(yīng)用于薄壁管和機(jī)翼等結(jié)構(gòu)的優(yōu)化設(shè)計中；Lucifredi等[8]將動態(tài)Kriging模型與神經(jīng)網(wǎng)絡(luò)進(jìn)行了對比，然后將其應(yīng)用于水電力系統(tǒng)的預(yù)測維護(hù)中。Kriging模型可以表示為張量積基函數(shù)的形式[10-11]，因此可以通過張量積基函數(shù)來解析地推導(dǎo)功能函數(shù)對輸入變量導(dǎo)函數(shù)的表達(dá)式，進(jìn)而利用較少的訓(xùn)練樣本得到較高精度的Kriging代理模型后，直接得到功能函數(shù)對輸入變量梯度信息的近似解，也即以較高的效率為AFOSM提供了較高精度的梯度信息，很好地解決了隱式問題梯度信息難以獲取的難題。

1 改進(jìn)一次二階矩方法

(1)

(2)

設(shè)功能函數(shù)的均值和標(biāo)準(zhǔn)差分別為μg和σg，上述線性極限狀態(tài)方程的可靠度指標(biāo)βR和失效概率Pf可以由式(3)和式(4)精確求解[3]：

(3)

Pf=Φ(-βR)

(4)

式中:Φ(·)為標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)。

對于輸入變量為非正態(tài)分布的情況，Rackwitz和Fiessler[12]提出了一種等價正態(tài)變量算法，簡稱為R-F法，將非正態(tài)變量轉(zhuǎn)化為等價的正態(tài)變量，然后再采用AFOSM法求解可靠度指標(biāo)，進(jìn)而得到失效概率。

2 Kriging代理模型的梯度解析解

Kriging代理模型是一種半?yún)?shù)插值技術(shù)，Sacks等[4]將Kriging模型表達(dá)為線性回歸部分和非參數(shù)部分2部分組成，即

g(x)=fT(x)β+Z(x)

(5)

式中：f(x)=[f1(x)f2(x)…fm(x)]T為輸入向量x的多項式基函數(shù)，可以是0階、一階或二階多項式，提供模型模擬的全局性近似；β=[β1β2…βm]T是需要確定的回歸系數(shù)，m為線性回歸部分基函數(shù)的個數(shù)；Z(x)為隨機(jī)過程實(shí)現(xiàn)的非參數(shù)部分，提供線性回歸部分fT(x)β的偏差，服從N(0,σ2)的正態(tài)分布，協(xié)方差為

Cov[Z(xi),Z(xj)]=σ2R(xi,xj)

i,j=1,2,…,N

(6)

式中：xi與xj為第i和第j個訓(xùn)練樣本點(diǎn)；N表示訓(xùn)練樣本點(diǎn)個數(shù)；σ2是隨機(jī)分布誤差的方差；R(xi,xj)為表征訓(xùn)練樣本點(diǎn)之間的空間相關(guān)性的相關(guān)函數(shù)，表達(dá)為

(7)

式中：n是輸入向量x的維數(shù)；xik和xjk分別為向量xi和xj的第k個分量；τk是相關(guān)性參數(shù)；d決定了函數(shù)在第k個坐標(biāo)方向的光滑性。Kriging模型中的一元相關(guān)函數(shù)有多種形式[4-5]。

未知參數(shù)β和σ2可以通過式(8)和式(9)進(jìn)行估計：

(8)

(9)

(10)

Kriging模型對于結(jié)構(gòu)輸入輸出關(guān)系具有較高的擬合精度，已經(jīng)被廣泛運(yùn)用于結(jié)構(gòu)可靠性分析與優(yōu)化設(shè)計中。關(guān)于模型詳細(xì)的描述，可以參閱文獻(xiàn)[4]。

Kriging模型表達(dá)為張量積基函數(shù)為[5,10]

(11)

式中：λi(i=0,1,…,N)為常系數(shù)。在Kriging模型的多種相關(guān)函數(shù)中，Gaussian相關(guān)函數(shù)被廣泛選擇使用，即hik(xk)=exp(-τk(xk-xik)2)，τk為變量xk對應(yīng)的相關(guān)參數(shù)。式(11)中，Kriging函數(shù)的多項式基函數(shù)取0階，因此，本文所采用的Kriging模型解析表達(dá)式為

(12)

分別對式(12)中等號兩邊輸入變量xq求導(dǎo)數(shù)，推導(dǎo)如下：

(13)

3 基于Kriging解析解的可靠性分析方法

本文提出一種基于Kriging解析解的AFOSM可靠性分析方法，傳統(tǒng)AFOSM方法中的導(dǎo)函數(shù)梯度信息由Kriging模型解析求得。值得注意的是，要采用改進(jìn)的一次二階矩方法求解可靠度指標(biāo)和失效概率，必須先知道設(shè)計點(diǎn)。顯然對于一個給定的非線性功能函數(shù)，其設(shè)計點(diǎn)是不可能預(yù)先知道的，這就必須采用迭代或者直接尋優(yōu)的方法來進(jìn)行問題的求解，一般是通過迭代前后兩次的可靠度指標(biāo)的相對誤差滿足設(shè)定的精度要求為止。對于成熟的AFOSM的迭代求解步驟的詳細(xì)信息可查閱文獻(xiàn)[3]。本文所提方法的簡要步驟如下文所述，流程圖如圖1所示。

1) 鑒于拉丁超立方抽樣(Latin Hypercube Sampling, LHS)[13]樣本具有低偏差特性，采用LHS方法產(chǎn)生N個訓(xùn)練樣本點(diǎn)構(gòu)建Kriging模型。

圖1 可靠性分析流程圖Fig.1 Flow chart of reliability analysis

4) 采用式(4)求解結(jié)構(gòu)失效概率。

4 算例分析

本節(jié)將所提的基于Kriging解析解的AFOSM法應(yīng)用于數(shù)值及工程算例中以驗證所提方法的正確性。4.1節(jié)利用簡單的數(shù)值算例來驗證所推導(dǎo)的Kriging梯度解析函數(shù)的正確性。4.2節(jié)將所提的基于Kriging解析解的AFOSM法應(yīng)用于某航空發(fā)動機(jī)渦輪盤模型的可靠性分析中，以驗證所提可靠性分析方法的正確性與高效性。在此基礎(chǔ)上，4.3節(jié)將其應(yīng)用于以有限元隱式形式表征的某導(dǎo)彈舵面結(jié)構(gòu)的可靠性分析中，進(jìn)一步驗證本文所提方法的正確性與工程適用性。

4.1 數(shù)值算例

考慮2個非線性功能函數(shù)g1(x)和g2(x)，均包含4個基本隨機(jī)變量xi(i=1,2,3,4)，且均服從正態(tài)分布xi～N(1,0.12)。功能函數(shù)表示為

(14)

易得2個功能函數(shù)對各輸入變量偏導(dǎo)數(shù)在其均值點(diǎn)處函數(shù)值的精確解分別為

(15)

運(yùn)用式(13)所推導(dǎo)的Kriging梯度解析函數(shù)進(jìn)行上述導(dǎo)數(shù)值求解，分別抽取一定數(shù)量的訓(xùn)練樣本建立2個功能函數(shù)的Kriging代理模型，代理模型建立后，將輸入變量均值點(diǎn)代入可直接得到相應(yīng)的梯度信息如表1所示。

由表1所示的結(jié)果可以看出，本文所推導(dǎo)的Kriging梯度解析解具有較高的精度，通過對比驗證了所推導(dǎo)的基于Kriging代理模型的梯度解析公式的正確性。同時可看出，借助于Kriging代理模型的優(yōu)勢，僅需要少量的訓(xùn)練樣本即可以獲得高精度的結(jié)果，計算代價小。

表1 數(shù)值算例的Kriging梯度解析結(jié)果

Table 1 Results of Kriging based gradient analysis of numerical example

項目導(dǎo)函數(shù)值誤差/%訓(xùn)練樣本Ng1(x)/x11.0000g1(x)/x23.0000g1(x)/x35.0000g1(x)/x47.000020g2(x)/x12.0030.15g2(x)/x25.9990.017g2(x)/x39.9970.03g2(x)/x414.0020.01425

4.2 航空發(fā)動機(jī)渦輪盤算例

本算例分別抽取50、100、300個訓(xùn)練樣本點(diǎn)建立Kriging代理模型，獲取相應(yīng)的梯度解析結(jié)果，進(jìn)而采用AFOSM法進(jìn)行渦輪盤可靠性分析。從表3可以看出，3種情況下，采用本文所提的基于Kriging解析解的AFOSM方法的可靠性分析結(jié)果均與Monte Carlo(MC)方法比較吻合，驗證了該方法的正確性。隨著訓(xùn)練樣本的增加，Kriging模型的精度也隨之增加，相應(yīng)的梯度解析解精度也隨之增加。然而，隨著訓(xùn)練樣本數(shù)量的增加，可靠性分析結(jié)果的精度提高程度很小，可見Kriging模型本身已經(jīng)具有較高的精度，然而其可靠性分析結(jié)果的結(jié)果受限于AFOSM方法本身的精度。

圖2 某型航空發(fā)動機(jī)渦輪盤示意圖[14]Fig.2 Schematic diagram of an aero-engine turbine disc[14]

表2 某型航空發(fā)動機(jī)渦輪盤模型輸入變量分布信息

Table 2 Input variable distribution information of an aero-engine turbine disc

輸入變量分布類型均值變異系數(shù)ρ/(kg·m-3)正態(tài)8 2400.1C/(kg·m)正態(tài)5.670.1A/m2正態(tài)6.2×10-30.1J/m4正態(tài)1.22×10-40.1n/(r·s-1)正態(tài)2000.05

表3 某型航空發(fā)動機(jī)渦輪盤模型可靠性分析結(jié)果

表3中所示方法的計算量也顯示出基于Kriging解析解的AFOSM法的計算效率，針對一般工程問題，只需少量的計算成本即可得到較高精度的可靠性分析結(jié)果，很大程度降低了計算代價。

4.3 飛行器舵面算例

圖3為某型飛行器舵面結(jié)構(gòu)示意圖[15-16]，其由5條桁條、6條翼肋、蒙皮以及用于和機(jī)體相連接的固定端所組成。各個部件的材料均為某鈦合金材料，材料的楊氏模量為E，泊松比為ν。翼肋與桁條的厚度均為trs，蒙皮的厚度為tsk。在飛行器的飛行過程中，舵面用于調(diào)節(jié)以及穩(wěn)定飛行姿態(tài)，其所承受的載荷主要為氣動載荷，大小為P，垂直作用于上蒙皮。舵面結(jié)構(gòu)材料的彈性模量、泊松比，桁條、翼肋、蒙皮的厚度以及所承受的氣動載荷均為服從正態(tài)分布的隨機(jī)變量，詳細(xì)分布信息見表4。以舵面結(jié)構(gòu)在承受氣動載荷情況下的最大變形不超過δc=8 mm為門限值，建立結(jié)構(gòu)系統(tǒng)的功能函數(shù)為

g(E,ν,P,tsk,trs)=δc-max(δ)

(16)

由于此算例為以有限元形式表征的隱式問題，舵面結(jié)構(gòu)在氣動載荷作用下的位移無法解析求解，本節(jié)采用有限元分析軟件MSC.Patran/Nastran對結(jié)構(gòu)的變形進(jìn)行求解，有限元模型如圖4[15-16]所示。一次響應(yīng)計算時間較長，若采用Monte Carlo法進(jìn)行可靠性分析，其計算代價是無法接受的。為了結(jié)果對比，在建立Kriging模型后，采用MC方法(Kriging_MC)針對建立的代理模型進(jìn)行結(jié)構(gòu)可靠性分析，以此結(jié)果作為參照解進(jìn)行對比。

如圖5所示，抽取200個訓(xùn)練樣本建立舵面結(jié)構(gòu)功能函數(shù)的Kriging代理模型，進(jìn)而抽取106個樣本代入建立的代理模型中進(jìn)行舵面結(jié)構(gòu)的可靠性分析，以此結(jié)果作為參照解。在抽取80個訓(xùn)練樣本的情況下，采用本文所提方法的可靠性分析結(jié)果與參照解較為接近，再次驗證了所提結(jié)構(gòu)可靠性分析方法的正確性，結(jié)果的差異性主要來源于AFOSM方法本身對于高度非線性問題的適用性。然而值得指出的是，若直接采取數(shù)值模擬法進(jìn)行有限元隱式問題的可靠性分析的計算代價是無法承受的，進(jìn)行200次有限元模型分析同樣需要花費(fèi)較長時間，計算成本相對較高，本文方法在80次模型調(diào)用的情況下，所得結(jié)果精度較高，進(jìn)而進(jìn)一步驗證了本文所提方法的工程適用性。

圖3 某型飛行器舵面結(jié)構(gòu)示意圖[15-16]Fig.3 Schematic diagram of an aircraft rudder structure[15-16]

表4 某型飛行器舵面結(jié)構(gòu)輸入變量分布信息

Table 4 Input variable distribution information of an aircraft rudder structure

輸入變量分布類型均值變異系數(shù)E/MPa正態(tài)1176000.1ν正態(tài)0.30.1P/MPa正態(tài)0.160.1tsk/mm正態(tài)10.05trs/mm正態(tài)30.05

圖4 某型飛行器舵面結(jié)構(gòu)有限元模型[15-16]Fig.4 Finite element model for an aircraft rudder structure[15-16]

表5 某型飛行器舵面結(jié)構(gòu)可靠性分析結(jié)果

Table 5 Reliability analysis results of an aircraft rudder structure

方法失效概率計算量Kriging_ANA_AFOSM0.008980Kriging_MC0.0092200

5 結(jié) 論

1) 針對工程廣泛應(yīng)用的AFOSM法對于隱式功能函數(shù)問題存在梯度信息較難求解的現(xiàn)狀，尤其是對于以有限元形式表征的工程復(fù)雜問題，梯度計算難度大、計算代價大，本文提出了基于Kriging模型梯度解析解的AFOSM可靠性分析方法，為AFOSM法提供了高精確的梯度信息，有效地解決了AFOSM方法的工程應(yīng)用難題。

2) 鑒于Kriging模型在工程應(yīng)用中的優(yōu)勢以及考慮到Kriging模型可以很好地表征為張量積基函數(shù)的形式，本文研究推導(dǎo)出了基于Kriging代理模型的功能函數(shù)對輸入變量導(dǎo)函數(shù)的解析表達(dá)式，Kriging模型建立好后可直接得到梯度解析結(jié)果。

3) 本文提出的基于Kriging代理模型梯度解析解的AFOSM法拓寬了AFOSM方法的應(yīng)用范疇，成為AFOSM方法的補(bǔ)充，為工程復(fù)雜問題的結(jié)構(gòu)可靠性分析提供了可供選擇的方法。算例驗證了基于Kriging模型的梯度解析結(jié)果具有較高精度，因此，基于Kriging模型解析解的AFOSM方法的精度與適用性與AFOSM本身密切相關(guān)。