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

中國(guó)飛行試驗(yàn)研究院 技術(shù)中心飛機(jī)所，西安 710089

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

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

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

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

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

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

(1)

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

(2)

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

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

(3)

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

(4)

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

(5)

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

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

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

根據(jù)Edgeworth級(jí)數(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固定在某一實(shí)現(xiàn)值時(shí)輸出的前四階條件中心矩，無條件及條件前四階矩的定義式為

(12)

(13)

從式(6)～式(11)可以看出，通過Edgeworth級(jí)數(shù)展開法，輸出的無條件概率密度函數(shù)及條件概率密度函數(shù)的求解轉(zhuǎn)換為求解輸出的無條件及條件前四階中心矩。對(duì)于輸出的無條件及條件前四階中心矩，本文在2.2節(jié)從乘法降維的角度推導(dǎo)了通過重復(fù)利用乘法降維計(jì)算輸出無條件矩中所產(chǎn)生的積分網(wǎng)格中的數(shù)據(jù)同時(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階原點(diǎn)矩可以從一個(gè)n維積分近似轉(zhuǎn)化為n個(gè)一維積分的乘積形式，即

(15)

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

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

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

(16)

(17)

通過輸出的前四階原點(diǎn)矩求得輸出的前四階中心矩，即

(18)

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

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

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

(19)

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

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

(20)

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

(21)

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

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

(22)

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

2.4 計(jì)算量

(23)

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

3 算例分析

3.1 算例1：航空發(fā)動(dòng)機(jī)渦輪盤

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

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

(24)

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

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

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

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

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

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

3.2 算例2：汽車前軸

圖2所示為一汽車前軸示意圖，故危險(xiǎn)截面常發(fā)生在工字梁上，其截面形狀如圖2所示。已知危險(xiǎn)截面的最大正應(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ù)值計(jì)算結(jié)果Table 6 Numerical results of automobile front axle

4 結(jié) 論

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

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

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

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

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

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