游令非，張建國，周霜，杜小松

1. 北京航空航天大學 可靠性與系統(tǒng)工程學院, 北京 100083 2. 北京航空航天大學 可靠性與環(huán)境工程技術(shù)國防重點實驗室, 北京 100083

在航空發(fā)動機適航要求中，將原發(fā)失效能夠引起發(fā)動機危害性影響的部件定義為發(fā)動機限壽件(Engine Life Limited Part,ELLP)，如旋轉(zhuǎn)輪盤和大型旋轉(zhuǎn)封嚴裝置等，其可靠水平直接決定著發(fā)動機整體可靠性和工作性能，設(shè)計中主要通過降低ELLP的失效概率來提高整機的可靠性和安全性[1-2]。而對航空發(fā)動機造成危害性后果的事件有很大部分是由其機械結(jié)構(gòu)的疲勞所導致。對壓氣機盤/渦輪盤進行基于疲勞可靠度計算的概率風險評估成為航空發(fā)動機型號適航取證過程以及保證其結(jié)構(gòu)完整性的關(guān)鍵技術(shù)和重要實施步驟之一。

目前國內(nèi)外針對以輪盤為代表的ELLP疲勞可靠性做了大量研究，美國軍方[3]首先將概率統(tǒng)計理論引入到傳統(tǒng)輪盤疲勞壽命預測方法中，形成了輪盤壽命可靠性設(shè)計方法。Melis等[4]對航空燃氣輪機壓氣機盤的壽命進行了預測和可靠性分析。Gao等[5]應用分布式協(xié)同響應面法對航空發(fā)動機渦輪進行了疲勞損傷的可靠性分析。Zhu 等[6-7]分別建立了渦輪葉片輪盤疲勞可靠性評估的計算試驗框架和概率分析框架。將超速試驗與隨機有限元分析相結(jié)合，量化了試驗數(shù)據(jù)、材料性能和載荷的不確定性以及材料變化和荷載變化對可靠性結(jié)果的影響。Hu 等[8]針對合金渦輪盤的鎳基高溫合金，建立了基于外加機械功密度的概率模型，利用線性異方差函數(shù)對蠕變疲勞壽命中的非常數(shù)偏差進行評價。Wang等[9]探索了渦輪盤低周疲勞壽命預測的區(qū)域可靠性方法以及裂紋的閉合行為。李巖等[10]提出了一種基于Kriging和蒙特卡羅半徑外自適應重要抽樣(Monte Carlo Radius-Outside Adaptive Importance Sampling,MCROAIS)方法混合的結(jié)構(gòu)概率風險評估方法，構(gòu)建了壓氣機盤風險概率模型。陳志英等[11]結(jié)合應力應變場強法和響應面法建立了輪盤疲勞可靠性模型并進行了可靠度計算。除此以外，隨機過程、神經(jīng)網(wǎng)絡(luò)和支持向量機等方法也廣泛地應用到相關(guān)問題中[12-14]。

目前開展ELLP結(jié)構(gòu)概率風險評估的難點在于非線性、小失效概率事件的定量要求驗證問題。但是對于高度非線性極限狀態(tài)函數(shù)以及小失效概率事件，以上方法均有一定的局限:① 階矩法對于非線性較高的問題計算精度不高；② 響應面法與數(shù)值模擬相結(jié)合使計算步驟增加，特別是針對小概率失效事件計算效率低下。因此，一些基于方差抽樣縮減技術(shù)的抽樣方法[15-17]因具有較高的精度和魯棒性成為解決該問題的選擇。

作為在可靠度計算中廣泛應用的數(shù)值模擬方法，重要抽樣方法(Important Sampling, IS)[18]將抽樣密度函數(shù)的抽樣中心移到設(shè)計驗算點(以下簡稱驗算點)，可以使更多的樣本點集中在失效域，提高了抽樣的效率。但在未知驗算點(或可靠度指標)的情況下，只能通過蒙特卡羅仿真(Monte Carlo Simulation,MCS)法或解析近似法求解。這極大地限制了它們的適用范圍。為了克服重要抽樣法對驗算點和可靠度指標已知的依賴性，產(chǎn)生了蒙特卡羅自適應重要抽樣(Monte Carlo Adaptive Important Sampling,MCAIS)法的思想。通過自動迭代不斷改進MCS的重要抽樣密度函數(shù)，使得抽樣點越來越趨近于失效域，使抽樣中心逐漸趨近于驗算點。Au和Beck[19]采用基于Metropolis算法的馬爾科夫鏈模擬方法生成失效域樣本，而后運用自適應密度估計方法構(gòu)造重要抽樣密度函數(shù)。Sankaran和Prakasb[20]將自適應重要抽樣法應用在大結(jié)構(gòu)系統(tǒng)的可靠性分析中。黃毅等[21]采用蒙特卡羅方法對中心柱設(shè)計過程中的主要影響因素——低周疲勞進行可靠性評價。陳向前等[22]提出了一種基于樣本概率密度加權(quán)的采樣中心確定方法，該方法在自適應重要抽樣的基礎(chǔ)上增加了有效抽樣中對失效概率貢獻大的樣本出現(xiàn)的概率。馬紀明等[23]提出一種改進的自適應重要抽樣方法，該方法通過失效樣本估計IS密度函數(shù)的初始參數(shù)，提高了算法的收斂時間。戴鴻哲等[24]提出了一種自適應Metropolis算法和快速高斯變換技術(shù)的結(jié)構(gòu)可靠性分析高效自適應重要抽樣方法。但在一些情況下，自適應重要抽樣法并不能得到理想的驗算點近似結(jié)果，比如變量的聯(lián)合概率密度函數(shù)的標準差比較大時，或者概率密度函數(shù)的梯度和極限狀態(tài)函數(shù)的梯度在某一區(qū)域大致相同時等。將會使近似驗算點和真實驗算點間具有較大的誤差。

為了克服以上問題，本文提出一種具有自更新機制的半徑外自適應重要抽樣(Auto Updating Monte Carlo Radius-Outside Adaptive Importance Sampling,AUMCROAIS)疲勞可靠性分析方法，本方法可以兼顧搜索效率和計算精度，在高效地尋找近似驗算點的同時，能較為精確地計算可靠度。本文提出的方法總體思路如下：首先利用MCAIS快速逼近驗算點附近，在MCAIS算法停止后，在以此時的近似驗算點為中心進行極坐標抽樣，并依次構(gòu)造主動學習函數(shù)進行評判，即對極坐標抽樣點先進行極限狀態(tài)函數(shù)的主動學習函數(shù)的構(gòu)造，根據(jù)判據(jù)自動選擇那些離極限狀態(tài)面近的抽樣點，而后針對這些點進行最優(yōu)半徑的主動學習函數(shù)的構(gòu)造并找到最小值點，此時，近似驗算點得到了更新，而后以更新近似驗算點為中心，利用半徑外重要抽樣計算可靠度，繼續(xù)適當縮減方差并重復上述步驟，直到滿足收斂條件，得到最終的可靠度。本方法通過不斷的更新確定出最優(yōu)抽樣半徑，加速失效概率計算的收斂，不受限于自適應抽樣的收斂條件，保證了效率同時提高了精度。AUMCROAIS疲勞可靠性分析方法解決了疲勞可靠性問題中小概率失效和極限狀態(tài)非線性較強的情況，支撐了適航規(guī)章要求(FAR33.75條款)[25]的實施，及對發(fā)動機失效概率事件分析的要求。最后在工程案例上利用本方法分別和蒙特卡羅方法、自適應抽樣法和一階可靠性方法進行對比，說明了本方法的適用性、精確性和高效性。

1 MCROAIS方法原理

圖1 MCROAIS法流程圖
Fig.1 Flow chart of MCROAIS method

(1)

MCROAIS雖然相比MCS提高了抽樣效率，但是MCROAIS依然存在一些缺陷，表現(xiàn)在：① 小概率情形下且變量的標準差較大時，將會導致β在近驗算點附近變化緩慢，將提前結(jié)束迭代過程；② 概率密度函數(shù)的梯度和極限狀態(tài)函數(shù)的梯度在某一區(qū)域大致相同，導致β變化緩慢，將提前結(jié)束迭代過程；③ 若非線性較高時同時收斂偏離驗算點，每次迭代時標準差的減小將不會使近似驗算點的收斂路線得到修正，反而將加速迭代結(jié)束，得到不精確的結(jié)果；④ 若收斂準則嚴格(即提高精度要求)，則計算效率將大大增加。為此，在第2節(jié)中將討論如何基于MCROAIS法加以改進并構(gòu)造新的抽樣方法的以同時兼顧抽樣效率和計算精度。

2 基于MCROAIS法的極坐標抽樣和主動學習函數(shù)

通過第1節(jié)分析可知，MCROAIS法并不能在任何極限狀態(tài)函數(shù)和變量分布，特別是小概率事件下得到理想的結(jié)果，本節(jié)將通過極坐標抽樣和主動學習函數(shù)的構(gòu)造，使其在滿足自適應抽樣收斂條件下，繼續(xù)構(gòu)造和識別離驗算點近的樣本點，從而達到自更新的目的。極坐標抽樣即以每次自適應迭代的結(jié)果得出的近似驗算點為中心繼續(xù)構(gòu)造抽樣，進而用主動學習函數(shù)依次對貼近極限狀態(tài)和最優(yōu)半徑進行篩選，更新出最優(yōu)點，用以構(gòu)造最優(yōu)半徑或者下一次迭代的起始點。極坐標抽樣和主動學習函數(shù)相結(jié)合的方法提高了抽樣效率，同時也考慮了計算精度。

2.1 極坐標抽樣

以單次自適應抽樣的近似驗算點為中心進行極坐標抽樣，極坐標抽樣基于布朗運動的軌跡概率分布，其特點是分別基于極角和極徑分布進行抽樣，且樣本點數(shù)控制在很小的范圍內(nèi)，大大地提高了抽樣效率同時也使樣本點更多地接近近似驗算點一側(cè)。下面以二維布朗運動加以簡要說明，高維同理。設(shè)x方向和y方向的位移X和Y都是獨立的隨機變量，符合相同的正態(tài)分布N(0,σ2)。X和Y的聯(lián)合分布是fX和fY的乘積：

(2)

式中：σ2為方差。

位移R2=X2+Y2是隨機變量的函數(shù)，其分布函數(shù)為

(3)

式中：r為極徑；θ為極角；D是x2+y2≤r2的區(qū)域，通過變成極坐標求得結(jié)果。對式(3)求導即得到概率密度分布:

(4)

(5)

(6)

則極徑r和極角θ的聯(lián)合分布為

pRΘ(r,θ)=

(7)

由于無規(guī)行走的極徑和極角是相互獨立的隨機量，故2個邊緣分布可以直接求積分，則在近似驗算點x*對應的近似最優(yōu)半徑β*下極徑和極角的分布為

(8)

(9)

圖2 極坐標抽樣
Fig.2 Polar coordinate sampling

2.2 主動學習函數(shù)的構(gòu)造

(10)

3 AUMCROAIS可靠度計算方法

AUMCROAIS法通過對最優(yōu)樣本的不斷更新，達到對小概率事件和強非線性極限狀態(tài)函數(shù)的真實驗算點的近似，從而得到最優(yōu)半徑，繼而用MCROAIS法進行可靠度求解。其包括內(nèi)循環(huán)和外循環(huán)2部分。

3.1 自更新機制

外循環(huán)包括以下幾個部分：以每次改進的自適應抽樣的近似驗算點為中心進行的極坐標抽樣(pR(β*),pΘ(θ))，而后依次針對“近極限狀態(tài)面”和“最優(yōu)半徑”進行主動學習函數(shù)的構(gòu)造和判別，首先建立近極限狀態(tài)面主動學習函數(shù)LG(xi)并判別找出那些近極限狀態(tài)面的點xi(i=1,2,…,N)，再建立最優(yōu)半徑主動學習函數(shù)Lβ*(xi)并根據(jù)判別，找出近極限狀態(tài)面的點中半徑最小的點，即最優(yōu)驗算點x*，并基于更新的驗算點進行半徑外重要抽樣可靠度計算R。最終進行可靠度收斂判別，若不滿足要求，則進行下一次外循環(huán)直至滿足外循環(huán)的可靠度收斂判別?；谝陨喜襟E可同時兼顧效率和精度，同時構(gòu)建雙主動學習函數(shù)，給極限狀態(tài)建立了統(tǒng)一的評判標準，有助于快速高效地識別出最優(yōu)驗算點。

內(nèi)循環(huán)即根據(jù)本算法的特點和需要改進的自適應抽樣，旨在快速使近似驗算點x*收斂于真實驗算點x*附近。在迭代過程中，為使近似驗算點快速收斂，內(nèi)循環(huán)中自適應抽樣的收斂準則可適當放大，同時，由于第1次外循環(huán)已經(jīng)到達真實驗算點附近，所以第2次以后的循環(huán)中(如果需要)，近似驗算點的收斂步長將適當減小，則：在第1次循環(huán)中σhV(v)i=σhV(v)i-1×1，hV(v)為重要抽樣密度函數(shù)，第2次及以后的循環(huán)中σhV(v)i=σhV(v)i-1×0.9。同時，為加速內(nèi)循環(huán)的收斂，將γ設(shè)為1.02，以此來提高內(nèi)循環(huán)的計算精度。

3.2 算法流程

步驟4若沒有樣本點落在失效域內(nèi)，令βi=β**并把hV(v)移到Y(jié)**，否則比較β*>γβi-1，若不等式成立，令βi=β**，σhV(v)i=σhV(v)i-1×0.9(k≠1時)。i=i+1，當|βi-βi-1|<0.01時，轉(zhuǎn)到步驟2，直到滿足要求，記錄此時的Y=x*，σhV(v)k=σhV(v)i×0.9。

步驟5以x*為中心，進行極坐標抽樣(pR(r),pΘ(θ))，樣本點記為xi(i=1,2,…,N)。

步驟10根據(jù)下式確定第k次外循環(huán)下ELLP的失效概率:

4 應用實例

將本文提出的方法應用于某型發(fā)動機低壓壓氣機輪盤，通過與MCS、MCROAIS法和FORM法進行對比，驗證本文方法的高效率、魯棒性和仿真精度。

針對影響發(fā)動機安全性的關(guān)鍵件壓氣機輪盤的主要失效模式—低循環(huán)疲勞斷裂，本文選取某型發(fā)動機低壓三級輪盤為典型關(guān)鍵件，開展發(fā)動機輪盤疲勞可靠性計算。

4.1 壓氣機輪盤疲勞試驗及仿真驗證

采用立式旋轉(zhuǎn)試驗系統(tǒng)開展疲勞壽命試驗，試驗在3 mmHg真空度下進行[10]，如圖4所示。立式旋轉(zhuǎn)系統(tǒng)控制輪盤的轉(zhuǎn)速由低速500 r/min至高速9 550 r/min之間呈三角波式周期性變化，周期為30 s，如圖5所示。

通過該試驗可確定其疲勞危險位置，試驗轉(zhuǎn)速控制在500-9 550-500 r/min，共完成55 300次試驗器循環(huán)，約461 h，完成114 265次試驗器循環(huán)后，對該盤進行無損探傷檢查，在銷釘孔6點鐘方向處出現(xiàn)了穿透性裂紋并延伸至盤體，分析可知低壓三級輪盤銷釘孔裂紋屬于疲勞斷裂，起源于銷釘孔的內(nèi)端面。

同時，利用有限元方法對轉(zhuǎn)速為9 550 r/min的試驗結(jié)果進行仿真驗證。渦輪盤的可靠性主要取決于結(jié)構(gòu)危險點的應力應變分布。由于結(jié)構(gòu)上的形狀和負載完全對稱，因此本例中考慮輪盤的1/37，整體模型見圖6，1/37模型見圖7。為還原真實約束，仿真時添加空心銷軸，軸材料為3Cr13，輪盤材料為TC11。結(jié)合結(jié)構(gòu)循環(huán)對稱的幾何特性，對有限元模型施加以下位移約束：

1) 以發(fā)動機輪盤的軸線為z軸建立圓柱坐標系，在該圓柱坐標系內(nèi)，如圖8所示，約束A、B兩面上全部節(jié)點沿旋轉(zhuǎn)方向θ的位移。

2) 在總體笛卡爾坐標系下，約束圖8中銷釘兩端面——C、D面，以及輪盤E、F面上全部節(jié)點沿x方向——即沿銷釘軸向的位移。

葉片載荷施加在銷軸上，載荷合力為29 189.6 N(這一載荷為航空發(fā)動機研究所相關(guān)部門給出的葉片產(chǎn)生的等效載荷)，方向為笛卡爾坐標系y方向。參數(shù)詳情見表1。

圖3 AUMCROAIS方法流程
Fig.3 Flowchart of AUMCROAIS method

圖4 低壓壓氣機輪盤試驗組裝
Fig.4 Test assembly of low pressure compressor disk

圖5 500～9 550 r/min轉(zhuǎn)速的變化示意圖
Fig.5 Schematic diagram of change of rotation speed between 500 and 9 550 r/min

圖6 低壓壓氣機輪盤模型
Fig.6 Low pressure compressor disk model

圖7 1/37輪盤和空心銷模型
Fig.7 Model of 1/37 disk and hollow pin

圖8 各面約束
Fig.8 All-sided constraints

選取表1中的數(shù)值作為仿真參數(shù)，利用ANSYS18.1進行仿真。模擬結(jié)果如圖9所示。由結(jié)果可知應力應變水平最高處位于輪盤與銷軸連接處(位置在銷釘孔下部內(nèi)端面，節(jié)點號1535)，此位置即結(jié)構(gòu)的破壞點，最大應力為791.64 MPa。從而驗證了試驗的正確性和準確性，下面將用試驗結(jié)果進行發(fā)動機輪盤第一級疲勞應力作用下(9 550 r/min)的疲勞可靠性計算。

表1 渦輪盤參數(shù)Table 1 Turbine disk parameters

圖9 9 550 r/min轉(zhuǎn)速下輪盤的應力云圖
Fig.9 Stress cloud chart of disk at 9 550 r/min rotate speed

4.2 壓氣機輪盤疲勞可靠性分析

lgNlife=lgK-mlg(Smax-703.84)

(11)

式中：Nlife為疲勞載荷循環(huán)次數(shù)；σa為循環(huán)疲勞應力幅值；Smax為應力場強；m和K為材料參數(shù)；壓氣機輪盤材料為TC11，由《飛機結(jié)構(gòu)金屬材料力學性能手冊》[26]得到其材料參數(shù)K=6.278 4×1015，m=4.736，利用應力-強度干涉理論構(gòu)建模型為

K(Smax-703.84)-m-nmax

(12)

K(Smaxi-703.84)-m-nmaxi

(13)

同時，由TC11進行了材料級試驗[8]，統(tǒng)計得到TC11材料參數(shù)m和K的均值、變異系數(shù)分別為μm=4.628、COVm=0.01，μK=6.135 26×1015，COVK=0.015，服從正態(tài)分布。

第1級應力場強的均值和變異系數(shù)為μSmax1=791.64，COVSmax1=0.1，服從極值分布。利用本文中提出的AUMCROAIS疲勞可靠性分析方法，對第一級循環(huán)應力作用下輪盤疲勞失效概率與應力循環(huán)次數(shù)的關(guān)系進行求解，同時用MCS, FORM法與MCROAIS法進行比較。其中由于MCS抽樣規(guī)模較大，認為其結(jié)果是最接近真實結(jié)果的并以MCS的結(jié)果作為精確度對比。

在發(fā)動機輪盤第1級疲勞應力作用下得到疲勞失效概率Pr隨載荷循環(huán)次數(shù)的關(guān)系如圖10所示，載荷循環(huán)從10 000～55 000，4種方法的可靠度指標計算結(jié)果見表2，可以看出隨著疲勞載荷循環(huán)次數(shù)的增加，失效概率逐漸升高。

通過將本文方法的仿真效率、魯棒性以及仿真精度與MCS方法進行比較，從而驗證本文提出方法的適用性。

1) 仿真效率：通過MCROAIS收斂速度對比，可以驗證本方法的仿真效率。如圖11所示，抽取12 000次、24 000次2種情況作為對比其可靠度指標β的收斂過程，可以看出本方法在抽樣效率上的優(yōu)越性，由于結(jié)合了極坐標抽樣和主動學習函數(shù)，使近似驗算點能迅速接近真實驗算點x*，不斷更新最優(yōu)抽樣半徑；另外，在抽樣次數(shù)上與MCS和MCROAIS進行對比，在每一級循環(huán)下，MCS需108次抽樣；MCROAIS分別需32次和25次迭代，每次1 000個樣本點，則MCROAIS共需32 000和25 000個樣本點；而本方法則分別需要2次外循環(huán)，11次內(nèi)循環(huán)； 2次外循環(huán)，8次內(nèi)循環(huán)；每次內(nèi)循環(huán)需抽樣1 000個樣本點，外循環(huán)36個樣本點，則本方法分別需要11 072和8 072個樣本點。同時, 在與FORM法收斂速度對比上可以看出，AUMCROAIS的收斂速度與FORM基本相同，且在一定程度上要先于FORM到達真實MPP點附近。綜上所述，本方法在仿真效率上得到了很大的提升。

圖10 第一級疲勞應力作用下輪盤失效概率與 疲勞應力循環(huán)次數(shù)的關(guān)系
Fig.10 Relationship between failure probability of disk and number of fatigue stress cycles under the first-order fatigue stress

表2 可靠度指標結(jié)果對比Table 2 Comparisons of reliability index results

圖11 可靠度指標收斂速度對比
Fig.11 Comparison of convergence rate of reliability index

2) 魯棒性：針對高非線性的風險概率模型的失效概率仿真，驗證本方法的魯棒性。如圖10所示，圍繞精確解(MCS結(jié)果)，AUMCROAIS和MCROAIS均有所波動，但總體來說本文方法波動較小，更具參考性。同時，F(xiàn)ORM法的結(jié)果隨著應力循環(huán)次數(shù)的增加越來越偏離真實值，而本方法依然在小范圍內(nèi)波動，說明本方法將在因高非線性而導致FORM偏離真實值或者不適用的情況下依然能給出近似解。

3) 仿真精度：通過仿真結(jié)果誤差分析，來驗證本文方法的仿真精度。仿真精度見表2，可以看出可靠度指標總體誤差不大于0.052 84，總體結(jié)果和FORM法基本相近，相比MCROAIS精度得到了很大提升。

從計算結(jié)果對比中可以看出，在不降低仿真精度和仿真效率得情況下，AUMCROAIS得到了較滿意的結(jié)果，解決了在小概率失效和非線性情況下MCROAIS并不能得到穩(wěn)定的精確解的問題，同時由于極坐標抽樣和主動學習函數(shù)的高效的抽樣和判別能力，使本方法的仿真效率大大提高，因此通過上述分析驗證了本文所提方法的較高的計算效率、魯棒性和計算精度。

5 結(jié) 論

本文針對小概率失效事件和強非線性可靠性分析問題，綜合考慮了抽樣效率和抽樣精度，針對航空發(fā)動機結(jié)構(gòu)疲勞可靠性分析問題，提出了AUMCROAIS方法，經(jīng)具體案例驗證表明：

1) 本文提出的AUMCROAIS方法實施方便，計算誤差較小，對類似的小概率失效事件和強非線性可靠性分析問題具有一定的指導意義，例如本文中航空發(fā)動機結(jié)構(gòu)疲勞可靠性分析案例中本方法的計算誤差最大不超過0.052 84，和MCS方法計算結(jié)果貼合度較高。

2) 本文提出的AUMCROAIS方法可以結(jié)合仿真，在完善隨機參數(shù)信息的情況下，可得到較精確的疲勞可靠度計算結(jié)果，同時計算效率大大提高，總算次數(shù)約為MCS方法的1/3左右。