失效概率矩獨立全局靈敏度分析的高效算法

蔣獻(xiàn)，王言，孟敏

中國飛行試驗研究院 技術(shù)中心飛機所，西安 710089

靈敏度分析主要衡量輸入變量的不確定性對輸出變量不確定性的影響程度。目前，靈敏度分析主要分為：局部靈敏度分析、區(qū)域靈敏度分析[1-3]以及全局靈敏度分析[4-5]。全局靈敏度分析以其能夠從全局的角度反映輸入變量不確定性對輸出變量不確定性的影響程度而在工程設(shè)計及可靠性評估中廣為應(yīng)用。全局靈敏度分析方法主要有：非參數(shù)法[6]、方差靈敏度分析法[7-9]、基于概率密度函數(shù)的矩獨立靈敏度分析法[10]以及失效概率矩獨立靈敏度分析法[11-12]。目前，方差靈敏度分析法及基于概率密度函數(shù)的矩獨立全局靈敏度分析的計算方法較為成熟，主要包括樣本法[13-15]、矩方法[16-17]以及代理模型法[12,18-20]。

失效概率矩獨立全局靈敏度分析主要衡量的是輸入變量在其整個取值區(qū)域變化時對結(jié)構(gòu)失效概率的平均影響程度。失效概率矩獨立全局靈敏度分析結(jié)果可以起到指導(dǎo)可靠性優(yōu)化設(shè)計的作用。文獻(xiàn)[12]建立了失效概率矩獨立全局靈敏度指標(biāo)與方差全局靈敏度指標(biāo)的關(guān)系?；谖墨I(xiàn)[12]，文獻(xiàn)[21]提出了單層分析法來求解該指標(biāo)，單層分析法是基于樣本的方法且其計算量仍與輸入變量的維數(shù)線性相關(guān)。文獻(xiàn)[22]將極大熵及Nataf變換相結(jié)合提出了一種高效的分析算法，但由于Nataf變換僅適用于正態(tài)分布及小變異系數(shù)下的對數(shù)正態(tài)分布[23]，對非正態(tài)變量及大變異系數(shù)的對數(shù)正態(tài)分布的情況，文獻(xiàn)[22]的方法將產(chǎn)生較大的估計誤差。為避免使用Nataf變換，文獻(xiàn)[17]將極大熵與三點估計結(jié)合，其計算過程仍是雙層嵌套的過程。外層期望的求解采用三點估計，內(nèi)層條件失效概率及無條件失效概率求解過程中的概率密度函數(shù)估計采用極大熵法，該方法的計算量較文獻(xiàn)[22]的高，且計算精度依賴于三點估計對外層期望的估計精度以及極大熵方法中優(yōu)化過程的準(zhǔn)確性。

通過上述對已發(fā)展的失效概率矩獨立全局靈敏度指標(biāo)算法的分析，可以看出極大熵結(jié)合Nataf變換[22]算法的效率較其他算法都高，但Nataf變換的使用限制使得該算法具有了一定的局限性，因此為了繼承該算法的高效性，同時又克服該算法的局限性，本文將乘法降維及Edgeworth級數(shù)展開的思想相結(jié)合，通過Edgeworth級數(shù)展開將輸出的無條件及條件失效概率的求解轉(zhuǎn)化為無條件及條件前四階矩的求解，對于輸出的無條件及條件四階矩的求解，通過乘法降維推導(dǎo)了重復(fù)利用積分網(wǎng)格信息求解的策略，并通過重復(fù)利用積分網(wǎng)格內(nèi)的信息求解外層積分。該算法的計算量(真實物理模型的調(diào)用次數(shù))僅與無條件矩計算中積分網(wǎng)格的建立有關(guān)，條件矩及外層積分的求解都無需額外的計算量，且該算法避免采用Nataf變換，消除了Nataf變換所帶來的限制以及避免了極大熵求解概率密度函數(shù)時的優(yōu)化過程，提高了計算的準(zhǔn)確性，且不失計算的高效性。

本文第1節(jié)簡要回顧了失效概率矩獨立全局靈敏度指標(biāo)的定義。第2節(jié)詳細(xì)介紹了本文所提的基于乘法降維結(jié)合Edgeworth級數(shù)展開的失效概率矩獨立全局靈敏度求解算法。第3節(jié)通過分析航空發(fā)動機渦輪盤以及汽車前軸的失效概率矩獨立全局靈敏度，說明了本文方法的準(zhǔn)確性及高效性。第4節(jié)對本文進(jìn)行了總結(jié)。

1 基于失效概率的矩獨立全局靈敏度指標(biāo)的定義

對于一個分析模型Y=g(X),Y為結(jié)構(gòu)或系統(tǒng)的輸出,X=[X1,X2,…,Xn]為結(jié)構(gòu)或系統(tǒng)模型的隨機輸入，n為模型的輸入變量個數(shù),g(·)為模型的極限狀態(tài)函數(shù)。失效概率的計算式為

(1)

式中：F為失效域，定義失效域為F={X:g(X)≤ 0}；fX(x)為輸入變量的聯(lián)合概率密度函數(shù)，x為變量X的實現(xiàn)。因此，式(1)亦可表達(dá)為極限狀態(tài)函數(shù)的概率密度函數(shù)fY(y)在失效域的積分，即

(2)

式中：y為變量Y的實現(xiàn)。

為衡量輸入變量對結(jié)構(gòu)失效概率的影響，文獻(xiàn)[11]建立了失效概率矩獨立全局靈敏度指標(biāo)，其衡量的是無條件失效概率與條件失效概率的絕對差異的平均值，即

(3)

式中：i為第i個輸入變量；Pf|Xi為Xi固定在其分布內(nèi)某一實現(xiàn)值時輸出的條件失效概率；條件失效概率Pf|Xi可以表示為極限狀態(tài)函數(shù)的條件概率密度函數(shù)fY|Xi(y)在失效域內(nèi)的積分，即

(4)

將式(2)及式(4)代入到式(3)中,可得

(5)

從式(5)可以看出，若能高效準(zhǔn)確地求解輸出的無條件及條件概率密度函數(shù)，以及外層期望，即可高效準(zhǔn)確地分析得到失效概率矩獨立全局靈敏度。在第2節(jié)中，本文詳細(xì)介紹了所提的基于Edgeworth級數(shù)展開結(jié)合乘法降維的輸出無條件概率密度及條件概率密度函數(shù)的求解算法。

2 基于乘法降維及Edgeworth級數(shù)展開的失效概率矩獨立全局靈敏度高效求解算法

2.1 基于前四階矩的Edgeworth級數(shù)展開法求解無條件及條件輸出概率密度函數(shù)

根據(jù)Edgeworth級數(shù)展開法[24]，輸出的無條件概率密度函數(shù)可以表示為

fY(y)≈

(6)

式中:

(7)

(8)

輸出的條件概率密度函數(shù)可以表示為

fY|Xi(y)≈

(9)

式中:

(10)

(11)

其中：α1g、α2g、α3g、α4g為輸出的無條件前四階中心矩；α1g|Xi、α2g|Xi、α3g|Xi、α4g|Xi為Xi固定在某一實現(xiàn)值時輸出的前四階條件中心矩，無條件及條件前四階矩的定義式為

(12)

(13)

從式(6)～式(11)可以看出，通過Edgeworth級數(shù)展開法，輸出的無條件概率密度函數(shù)及條件概率密度函數(shù)的求解轉(zhuǎn)換為求解輸出的無條件及條件前四階中心矩。對于輸出的無條件及條件前四階中心矩，本文在2.2節(jié)從乘法降維的角度推導(dǎo)了通過重復(fù)利用乘法降維計算輸出無條件矩中所產(chǎn)生的積分網(wǎng)格中的數(shù)據(jù)同時求解輸出條件前四階中心矩及外層期望的計算式。

2.2 乘法降維求解輸出的無條件前四階中心矩

根據(jù)文獻(xiàn)[16]，模型Y=g(x)可以近似表達(dá)為如下的一維變量函數(shù)的乘積形式，即

g(x)≈[g(μ)]1-n·

(14)

式中:μ=[μ1,μ2,…,μn]是輸入變量X的均值向量，k(k=1，2，…，n)為輸入變量的次序。

因此，通過式(14)輸出的αp階原點矩可以從一個n維積分近似轉(zhuǎn)化為n個一維積分的乘積形式，即

(15)

表1 高斯積分網(wǎng)格Table 1 Gaussian integration grid

根據(jù)高斯積分，式(15)可以由式(16)計算得到，即

Mαpg≈gαp-nαp(μ)·

(16)

(17)

通過輸出的前四階原點矩求得輸出的前四階中心矩，即

(18)

2.3 乘法降維求解輸出的條件矩

根據(jù)式(14)可以推導(dǎo)出輸出的條件矩的近似計算式為

gαp-nαp(μ)gαp(xi,μ-i)·

(19)

根據(jù)高斯積分準(zhǔn)則，式(19)可計算得

Mαpg|Xi(xi)=gαp-nαp(μ)gαp(xi,μ-i)·

(20)

fY|Xi(y)=MXi(g(xi,μ-i))

(21)

式中：MXi(·)是通過Edgeworth級數(shù)展開求解概率密度函數(shù)的過程，該過程無需額外的模型調(diào)用次數(shù)。

對于式(5)的外層一維積分，本文仍采用高斯積分求得。因此，式(5)可通過式(22)進(jìn)行估計:

(22)

乘法降維求解矩的過程中對隨機變量的分布形式?jīng)]有限制，Edgeworth級數(shù)展開法近似概率密度函數(shù)的過程中也對隨機變量的分布形式?jīng)]有限制，因此本文所提算法在求解失效概率矩獨立全局靈敏度中對隨機輸入變量的分布形式?jīng)]有限制。

2.4 計算量

(23)

式中：p為輸入變量中對稱分布的個數(shù)。

3 算例分析

3.1 算例1：航空發(fā)動機渦輪盤

渦輪盤是航空發(fā)動機關(guān)鍵轉(zhuǎn)動部件之一，在起動和加速過程中承受著巨大的離心力和熱應(yīng)力。加上機構(gòu)形狀復(fù)雜，在工作過程中容易出現(xiàn)應(yīng)力集中部位，如榫槽槽底、銷釘孔等。在工作一段時間后，可能在這些部位出現(xiàn)裂紋故障，如圖1所示。某型航空發(fā)動機渦輪盤在工作時所受載荷為F=Cω2/2π+2ρω2J，ρ、C、ω和J分別是質(zhì)量密度、系數(shù)、轉(zhuǎn)動角速度和截面慣性矩，ω=2πn′,n′為轉(zhuǎn)動頻率。渦輪盤的強度極限為σs，截面積為A。建立渦輪盤破裂失效的極限狀態(tài)函數(shù)為

g(σs,ρ,C,A,J,n′)=σsA/F

(24)

6個輸入變量的分布參數(shù)如表2所示。根據(jù)乘法降維中求解一維變量積分的高斯數(shù)值積分法則建立求解整數(shù)矩的高斯積分網(wǎng)格如表3所示。根據(jù)式(23)計算的構(gòu)建這高斯積分網(wǎng)格所需的模型調(diào)用次數(shù)為27次。通過重復(fù)利用表3中的信息，可以將輸出的無條件前四階矩、輸出的條件前四階矩及指標(biāo)求解過程中的一維外層積分同時求得，所得失效概率矩獨立全局靈敏度指標(biāo)如表4所示。表4給出了雙層MCS(Monte Carlo Simulation)法、文獻(xiàn)[22]提出的極大熵結(jié)合Nataf變換法、文獻(xiàn)[17]提出的三點估計結(jié)合極大熵法以及本文所提算法的計算結(jié)果。從表4的分析結(jié)果中可以看出，極大熵結(jié)合Nataf變換法在求解質(zhì)量密度ρ及系數(shù)C的失效概率矩獨立靈敏度指標(biāo)

圖1 航空發(fā)動機渦輪盤模型裂紋示意圖Fig.1 Diagram of crack of aero-engine turbine disk model

表2 基本變量分布參數(shù)Table 2 Distribution parameters of input variables

中出現(xiàn)了嚴(yán)重的錯誤，這是源于Nataf變換對大變異系數(shù)的對數(shù)正態(tài)分布的聯(lián)合概率密度函數(shù)的近似不夠準(zhǔn)確。三點估計結(jié)合極大熵變換的方法得到的結(jié)果與雙層MCS得到的結(jié)果較為接近，但其模型調(diào)用次數(shù)為390次。本文算法僅通過27次模型調(diào)用，就得到了與雙層MCS較為接近的結(jié)果，其計算量僅為三點估計結(jié)合極大熵變換算法的6.92%。表4中也給出了雙層MCS法、文獻(xiàn)[22]所提的極大熵結(jié)合Nataf變換法、文獻(xiàn)[17]提出的三點估計結(jié)合極大熵法和基于本文所提方法在分析失效概率矩獨立全局靈敏度指標(biāo)的計算耗時。從表4中可以看出本文所提算法的計算耗時較已有對比算法都少，主要原因是本文所提算法在概率密度函數(shù)估計中避免了優(yōu)化過程且能重復(fù)利用無條件矩求解過程中的積分網(wǎng)格內(nèi)的信息，同時求得輸出的條件矩及指標(biāo)求解過程中的外層積分。雙層MCS方法的計算時間最長，極大熵結(jié)合Nataf變換算法較三點估計結(jié)合極大熵算法的計算時間長，其原因在于極大熵結(jié)合Nataf變換的算法中不僅存在優(yōu)化過程估計概率密度函數(shù)的耗時，同時還存在Nataf變換中求解非線性方程的耗時。因此，從計算耗時的角度再次說明了本文所提算法的高效性。6個隨機輸入變量的基于失效概率的矩獨立全局靈敏度的排序為n′>ρ>A>σs>J>C，該排序可用于指導(dǎo)該結(jié)構(gòu)的可靠性優(yōu)化設(shè)計。

表3 航空發(fā)動機渦輪盤的高斯積分網(wǎng)格Table 3 Gaussian integration grid of aero-engine turbine disk

表4 航空發(fā)動機渦輪盤的數(shù)值計算結(jié)果Table 4 Numerical results of aero-engine turbine disk

3.2 算例2：汽車前軸

圖2所示為一汽車前軸示意圖，故危險截面常發(fā)生在工字梁上，其截面形狀如圖2所示。已知危險截面的最大正應(yīng)力為σ=M/Wx和τ=T/Wρ，其中M和T分別為前軸所受的彎矩和轉(zhuǎn)矩，Wx和Wρ分別為結(jié)構(gòu)的截面系數(shù)和極截面系數(shù)，且有

(25)

Wρ=0.8bt2+0.4[a3(h-2t)/t]

(26)

式中：a、b、h、t為工字梁幾何參數(shù)。

圖2 汽車前軸示意圖Fig.2 Schematic diagram of automobile front axle

表5 輸入變量的分布參數(shù)Table 5 Distribution parameters of input variables

表6 汽車前軸的數(shù)值計算結(jié)果Table 6 Numerical results of automobile front axle

4 結(jié) 論

本文針對可靠性優(yōu)化設(shè)計中較為關(guān)鍵的可靠性全局靈敏度分析指標(biāo)，即失效概率矩獨立全局靈敏度指標(biāo)，提出了高效準(zhǔn)確的計算方法。本文推導(dǎo)了如何重復(fù)利用乘法降維估計輸出無條件前四階矩的輸入-輸出樣本來估計輸出的條件前四階矩以及外層的輸入變量的一維積分。通過理論推導(dǎo)及算例分析，得到如下結(jié)論：

1) 本文所提算法的模型調(diào)用次數(shù)僅為乘法降維估計輸出無條件前四階矩的模型調(diào)用次數(shù)。

2) 本文所提算法的計算量與極大熵結(jié)合Nataf變換的計算量一致，但避免使用Nataf變換從而克服了Nataf變換使用過程中所帶來的計算誤差以及避免了Nataf變換的局限性。對于極大熵結(jié)合Nataf變換法無法使用的大變異系數(shù)的對數(shù)正態(tài)分布情況，本文所提算法得到了較為精確的結(jié)果。

3) 較極大熵結(jié)合Nataf變換，本文所提算法也無需優(yōu)化過程，從而避免了優(yōu)化結(jié)果的不穩(wěn)定性對靈敏度分析結(jié)果的影響。

4) 在計算耗時上，本文算法的計算耗時對比于已有的對比算法都要少，主要原因是本文所提算法無需優(yōu)化擬合輸出的概率密度函數(shù)以及Nataf變換，從而避免了優(yōu)化過程的耗時以及Nataf變換中求解非線性方程過程的耗時。

5) 由于本文算法僅需明確輸入輸出關(guān)系，因此，不論顯式輸入輸出關(guān)系還是隱式輸入輸出關(guān)系，本文所提算法都適用。