郭媛媛，孫有朝，李龍彪

南京航空航天大學(xué) 民航學(xué)院，南京 210016

基于蒙特卡羅方法的民用飛機(jī)故障風(fēng)險(xiǎn)評(píng)估方法

郭媛媛，孫有朝*，李龍彪

南京航空航天大學(xué) 民航學(xué)院，南京 210016

為提高航空安全，降低民用飛機(jī)運(yùn)行過程中的故障風(fēng)險(xiǎn)水平，分別建立了單機(jī)風(fēng)險(xiǎn)和機(jī)隊(duì)風(fēng)險(xiǎn)計(jì)算方法，指出單機(jī)風(fēng)險(xiǎn)計(jì)算的關(guān)鍵是確定部件故障分布及故障率，威布爾分布很好地描述了部件故障分布情況。針對(duì)部件故障機(jī)理，給出了部件故障的未檢出概率、部件故障后果的條件概率及相應(yīng)嚴(yán)重性后果，判定風(fēng)險(xiǎn)水平能否接受；機(jī)隊(duì)故障風(fēng)險(xiǎn)在已知單機(jī)風(fēng)險(xiǎn)的基礎(chǔ)上，取決于機(jī)隊(duì)規(guī)模和當(dāng)前飛行循環(huán)水平。基于機(jī)隊(duì)運(yùn)行數(shù)據(jù)預(yù)測(cè)故障飛機(jī)數(shù)量，獲得不同參數(shù)分布、故障飛機(jī)數(shù)量和退役壽命下的機(jī)隊(duì)風(fēng)險(xiǎn)水平。以某型發(fā)動(dòng)機(jī)渦輪盤榫齒裂紋為例，通過蒙特卡羅方法，模擬單機(jī)及機(jī)隊(duì)部件產(chǎn)生故障的時(shí)間，分析部件故障引起的單機(jī)風(fēng)險(xiǎn)水平和機(jī)隊(duì)風(fēng)險(xiǎn)水平，基于運(yùn)行數(shù)據(jù)給出民用飛機(jī)單機(jī)和機(jī)隊(duì)在持續(xù)適航階段的風(fēng)險(xiǎn)水平。

民用飛機(jī)；機(jī)隊(duì)風(fēng)險(xiǎn)；單機(jī)風(fēng)險(xiǎn)；蒙特卡羅方法；累積分布函數(shù)；失效率

民用飛機(jī)作為全球航空運(yùn)輸?shù)闹匾M成部分，保障飛行安全是對(duì)民用飛機(jī)的首要要求。隨著全球航空運(yùn)輸需求增加，其運(yùn)行安全性逐漸受到重視。隨著民用飛機(jī)性能不斷提高，結(jié)構(gòu)愈加復(fù)雜，飛機(jī)部件在運(yùn)營階段發(fā)生故障的概率逐漸增加。因此，民用飛機(jī)風(fēng)險(xiǎn)評(píng)估方法的研究受到越來越多的關(guān)注。

國外最早對(duì)民用飛機(jī)風(fēng)險(xiǎn)進(jìn)行研究[1-5]。美國聯(lián)邦航空管理局于2010年發(fā)布了監(jiān)控安全/分析數(shù)據(jù)(Monitor Safety-Analyze Data，MSAD)條例，對(duì)持續(xù)運(yùn)行安全過程在結(jié)構(gòu)和標(biāo)準(zhǔn)方面給出了指導(dǎo)，要求提出定量的或可發(fā)展為定量的風(fēng)險(xiǎn)分析方法。同時(shí)要求給出可計(jì)算的風(fēng)險(xiǎn)因素，包括機(jī)隊(duì)風(fēng)險(xiǎn)和單機(jī)風(fēng)險(xiǎn)。其中，機(jī)隊(duì)風(fēng)險(xiǎn)為在不采取任何干預(yù)措施時(shí)，機(jī)隊(duì)在壽命期內(nèi)潛在的未檢出故障導(dǎo)致的后果損失。單機(jī)風(fēng)險(xiǎn)是相對(duì)于單架飛機(jī)上，機(jī)組和乘客在每個(gè)飛行單位(飛行小時(shí)或飛行循環(huán))面臨的飛機(jī)故障風(fēng)險(xiǎn)。故障模式為未檢出的結(jié)構(gòu)故障、失控、推進(jìn)系統(tǒng)故障或飛行安全系統(tǒng)故障。分析故障風(fēng)險(xiǎn)，確定災(zāi)難性后果是否需要采取糾正措施并評(píng)估該措施能否有效降低風(fēng)險(xiǎn)。

與國外相比，國內(nèi)關(guān)于民用飛機(jī)風(fēng)險(xiǎn)方面的研究起步較晚，主要在飛行安全與航空器技術(shù)的風(fēng)險(xiǎn)評(píng)估方面進(jìn)行了一定的相關(guān)性研究[6-13]。采用蒙特卡羅仿真方法模擬發(fā)動(dòng)機(jī)多故障發(fā)生情況，計(jì)算各故障模式的風(fēng)險(xiǎn)因子以及每次飛行風(fēng)險(xiǎn)；建立航空發(fā)動(dòng)機(jī)部件故障風(fēng)險(xiǎn)仿真模型、發(fā)動(dòng)機(jī)故障危險(xiǎn)等級(jí)劃分、危險(xiǎn)系數(shù)確定和風(fēng)險(xiǎn)因子計(jì)算的方法，預(yù)測(cè)發(fā)動(dòng)機(jī)故障的發(fā)生情況，評(píng)估發(fā)動(dòng)機(jī)運(yùn)營階段的故障風(fēng)險(xiǎn)，制定合理的可降低風(fēng)險(xiǎn)方案，分析維修周期對(duì)風(fēng)險(xiǎn)因子的影響，給出合理的降低風(fēng)險(xiǎn)方案；借助概率統(tǒng)計(jì)和安全性分析方法，判斷事件對(duì)單架飛機(jī)和整個(gè)機(jī)隊(duì)的影響程度；提出結(jié)合極值理論與Copula模型量化評(píng)估平尾結(jié)冰條件下飛行風(fēng)險(xiǎn)概率的方法；提出一種基于Kriging和蒙特卡羅半徑外重要抽樣混合的結(jié)構(gòu)概率風(fēng)險(xiǎn)評(píng)估方法；基于蒙特卡羅飛行仿真實(shí)驗(yàn)對(duì)結(jié)冰條件下人-機(jī)-環(huán)系統(tǒng)，構(gòu)建了飛行風(fēng)險(xiǎn)發(fā)生的判定條件。隨著中國民機(jī)項(xiàng)目的蓬勃發(fā)展，中國自主設(shè)計(jì)與生產(chǎn)的飛機(jī)將越來越多地獲得適航認(rèn)證并投入運(yùn)營，但是缺乏基于MSAD提出的對(duì)持續(xù)運(yùn)行安全上單機(jī)及機(jī)隊(duì)風(fēng)險(xiǎn)評(píng)估的相關(guān)研究，本文對(duì)單機(jī)風(fēng)險(xiǎn)及機(jī)隊(duì)風(fēng)險(xiǎn)定量計(jì)算給出了詳細(xì)分析過程，并以算例進(jìn)行延伸驗(yàn)證，動(dòng)態(tài)監(jiān)測(cè)機(jī)隊(duì)運(yùn)行狀態(tài)，確保民用飛機(jī)機(jī)隊(duì)運(yùn)行安全。

本文基于機(jī)隊(duì)運(yùn)行數(shù)據(jù)，建立了單機(jī)與機(jī)隊(duì)機(jī)械部件的風(fēng)險(xiǎn)評(píng)估模型，雙參數(shù)威布爾分布可擬合民用飛機(jī)上各種類型的機(jī)械部件壽命數(shù)據(jù)[14]。統(tǒng)計(jì)表明，雙參數(shù)威布爾分布對(duì)發(fā)動(dòng)機(jī)故障模式分布具有很好的擬合性，在確定部件故障數(shù)據(jù)的情況下，能夠很好地反映部件使用壽命與故障率之間的關(guān)系[6]。同時(shí)針對(duì)航空發(fā)動(dòng)機(jī)零部件等高可靠性的產(chǎn)品，其故障統(tǒng)計(jì)通常為小樣本，利用威布爾分布在處理小樣本數(shù)據(jù)時(shí)優(yōu)勢(shì)明顯[15]。本文采用威布爾分布模型預(yù)測(cè)航空發(fā)動(dòng)機(jī)故障風(fēng)險(xiǎn)，指出影響風(fēng)險(xiǎn)變化趨勢(shì)的關(guān)鍵因素，通過威布爾分布分別描述單機(jī)和機(jī)隊(duì)的部件失效率及累積分布函數(shù)，以蒙特卡羅方法仿真部件故障時(shí)間。預(yù)測(cè)壽命期內(nèi)飛機(jī)運(yùn)行過程中飛機(jī)部件故障的風(fēng)險(xiǎn)水平及機(jī)隊(duì)中故障飛機(jī)數(shù)量，確定了機(jī)隊(duì)風(fēng)險(xiǎn)評(píng)估方法。以渦輪盤榫齒故障為例，建立故障未檢出概率流程圖，給出故障不安全后果的因果鏈、條件概率及相應(yīng)死亡率，計(jì)算不同運(yùn)行狀態(tài)的單機(jī)及機(jī)隊(duì)故障風(fēng)險(xiǎn)，評(píng)估在壽命期內(nèi)單機(jī)風(fēng)險(xiǎn)及機(jī)隊(duì)風(fēng)險(xiǎn)是否處于可接受水平，分析維修間隔及退役壽命的制定是否合理，在綜合比較安全要求和經(jīng)濟(jì)效益的情況下，確定維修間隔、退役壽命及單機(jī)和機(jī)隊(duì)的運(yùn)行風(fēng)險(xiǎn)水平。

1 單機(jī)與機(jī)隊(duì)風(fēng)險(xiǎn)評(píng)估模型

單機(jī)風(fēng)險(xiǎn)RI是機(jī)組和乘客在每個(gè)飛行小時(shí)/循環(huán)面臨的飛機(jī)故障風(fēng)險(xiǎn)。單機(jī)風(fēng)險(xiǎn)針對(duì)單架飛機(jī)，對(duì)每個(gè)飛行循環(huán)產(chǎn)生故障的概率進(jìn)行模擬。機(jī)隊(duì)風(fēng)險(xiǎn)RF是衡量機(jī)隊(duì)風(fēng)險(xiǎn)的重要指標(biāo)，是機(jī)隊(duì)在壽命期內(nèi)由于某部件發(fā)生故障導(dǎo)致的運(yùn)行風(fēng)險(xiǎn)。單機(jī)及機(jī)隊(duì)故障風(fēng)險(xiǎn)評(píng)估流程分為3步，如圖1所示。

圖1 機(jī)隊(duì)與單機(jī)故障風(fēng)險(xiǎn)評(píng)估過程Fig.1 Processes of fleet and individual failure risk assessment

步驟1通過安全分析監(jiān)控?cái)?shù)據(jù)確定機(jī)隊(duì)運(yùn)行過程中存在潛在的不安全狀態(tài)時(shí)，收集機(jī)隊(duì)相關(guān)數(shù)據(jù)，包括當(dāng)前運(yùn)行數(shù)據(jù)、維修/運(yùn)行數(shù)據(jù)、模型數(shù)據(jù)、設(shè)計(jì)數(shù)據(jù)和歷史數(shù)據(jù)，如飛機(jī)數(shù)量、各飛機(jī)當(dāng)前飛行循環(huán)、飛機(jī)退役壽命、當(dāng)前已檢出故障飛機(jī)數(shù)量和災(zāi)難性后果的致死率等。確定機(jī)隊(duì)情況，包括研究狀況、不安全后果及故障模型。確定飛機(jī)產(chǎn)生故障的參數(shù)，包括機(jī)隊(duì)飛行循環(huán)、飛機(jī)數(shù)量和不安全狀態(tài)或后果，為機(jī)隊(duì)故障風(fēng)險(xiǎn)評(píng)估提供數(shù)據(jù)支持。

步驟2結(jié)合部件故障概率、故障數(shù)量、維修間隔和未檢出概率，確定分布模型的形狀參數(shù)、壽命參數(shù)、故障率函數(shù)和風(fēng)險(xiǎn)分布函數(shù)。通過故障因果鏈確定不安全后果的條件概率和致死率。計(jì)算單架飛機(jī)某部件在壽命期內(nèi)的故障率和單機(jī)風(fēng)險(xiǎn)；給出預(yù)計(jì)故障飛機(jī)數(shù)，計(jì)算機(jī)隊(duì)風(fēng)險(xiǎn)。

步驟3確定風(fēng)險(xiǎn)評(píng)估方法有效性，判定故障引起的機(jī)隊(duì)風(fēng)險(xiǎn)是否在可接受水平。若是，無需采取措施，繼續(xù)監(jiān)控飛行數(shù)據(jù)[3]；若否，針對(duì)風(fēng)險(xiǎn)評(píng)估結(jié)果，分析維修間隔和退役壽命的制定是否合理，采取最優(yōu)糾正措施，重新評(píng)估機(jī)隊(duì)故障風(fēng)險(xiǎn)，執(zhí)行步驟2。

1.1 單機(jī)風(fēng)險(xiǎn)模型

通過確定飛機(jī)部件在運(yùn)行期間是否發(fā)生故障，預(yù)測(cè)在壽命期內(nèi)由于該故障導(dǎo)致機(jī)上個(gè)體遭受致命傷害的概率[16-17]。單機(jī)風(fēng)險(xiǎn)為

RI=h(tI)μIψIfI

(1)

式中：h(tI)為部件工作到時(shí)刻tI時(shí)尚未故障，tI時(shí)刻后的下一單位時(shí)間發(fā)生故障的概率；μI為部件故障導(dǎo)致飛機(jī)發(fā)生不安全后果的未檢出概率；ψI為飛機(jī)災(zāi)難性后果發(fā)生概率；fI為死亡率。發(fā)動(dòng)機(jī)部件故障概率h(tI)為

(2)

式中：β為形狀參數(shù)；η為特征壽命。查閱文獻(xiàn)確定單架飛機(jī)部件雙參數(shù)取值；收集機(jī)隊(duì)運(yùn)行數(shù)據(jù)，通過極大似然估計(jì)確定機(jī)隊(duì)部件雙參數(shù)取值。

文獻(xiàn)[18]利用專家經(jīng)驗(yàn)確定運(yùn)輸類飛機(jī)部件故障檢查效果及檢出概率，據(jù)此判斷損傷、故障是否明顯，檢查手段是否有效，是否受到疲勞損傷影響等確定未檢出概率μ值。如圖2所示，存在5種原因使故障未檢出，未檢出概率越低表明部件在導(dǎo)致不安全后果前被檢出的概率越大。

步驟1缺陷或者故障是否明顯。若部件由潛在故障發(fā)展為功能故障的過程中被機(jī)組或維修檢查人員容易檢出，未檢出概率應(yīng)設(shè)為較小值μ=0.01；若否，執(zhí)行步驟2。

步驟2當(dāng)前執(zhí)行的檢查措施是否有效。若是，進(jìn)一步判斷潛在故障是否在檢查后發(fā)展為功能故障，若是，應(yīng)適當(dāng)縮短該部件的維修計(jì)劃間隔，并設(shè)μ=0.06；若否，執(zhí)行步驟3。

圖2 確定未檢出故障概率流程Fig.2 Flowchart for determining non-detection failure probability

步驟3檢查是否起作用。若結(jié)構(gòu)設(shè)計(jì)符合安全壽命理論，如檢查脆性材料起落架結(jié)構(gòu)完整度不能很好地起到檢查監(jiān)控作用，在發(fā)生之前很難檢出，需按時(shí)移除替換，設(shè)μ=0.99；若結(jié)構(gòu)設(shè)計(jì)符合損傷容限理論，如渦輪盤，對(duì)其檢查起作用，執(zhí)行步驟4。

步驟4是否受疲勞損傷故障影響。相鄰結(jié)構(gòu)部件在相似應(yīng)力水平下運(yùn)行并發(fā)生裂紋，若是，判斷故障是由多個(gè)部件(如相鄰部件，相似結(jié)構(gòu)的縱梁)或是由多個(gè)點(diǎn)(如一個(gè)結(jié)合點(diǎn)的相鄰緊固件)增長,這些裂紋造成災(zāi)難性后果前難以被檢出[3]；若否，執(zhí)行步驟5。

步驟5加載路徑設(shè)計(jì)如果是單條，故障極不可能在災(zāi)難性后果發(fā)生前被檢出，μ=0.9；如果是多條，未檢出概率低于單條加載路徑，μ=0.3。

Pai為部件發(fā)生第a種故障，發(fā)展至i狀態(tài)的條件概率。給出相應(yīng)部件故障因果鏈，如圖3所示。

圖3 部件故障因果鏈Fig.3 Causal chain of parts failure

該部件故障導(dǎo)致2種不安全后果。由當(dāng)前故障狀態(tài)順次發(fā)展最終至不安全后果，相應(yīng)死亡率為fI1、fI2。該部件發(fā)生不安全后果的風(fēng)險(xiǎn)為

1.2 機(jī)隊(duì)風(fēng)險(xiǎn)模型

在機(jī)隊(duì)運(yùn)行過程中，存在通過無損檢測(cè)未檢出的故障概率，預(yù)測(cè)機(jī)隊(duì)從當(dāng)前到退役過程中的故障飛機(jī)數(shù)量，由此發(fā)生系統(tǒng)故障甚至危及飛行安全的災(zāi)難性事故的可能性，給出嚴(yán)重性等級(jí)fF,機(jī)隊(duì)風(fēng)險(xiǎn)RF為

RF=εFμFψFfF

(3)

式中：εF為故障飛機(jī)預(yù)測(cè)數(shù)量；μF為未檢出故障概率；ψF為發(fā)生災(zāi)難性后果的條件概率，因果鏈條件概率乘積之和，與故障類型和飛行循環(huán)相關(guān)；fF為嚴(yán)重性等級(jí)。

ε為機(jī)隊(duì)中每架飛機(jī)從當(dāng)前飛行循環(huán)到退役飛行循環(huán)發(fā)生故障的條件概率之和[18]，以累積分布函數(shù)表示：

(4)

式中：N為機(jī)隊(duì)飛機(jī)數(shù)量；Ai為第i架飛機(jī)預(yù)測(cè)故障概率；F(t)為累積分布函數(shù)，t為飛行循環(huán)，tR為退役飛行循環(huán)，ti為當(dāng)前飛行循環(huán)，tR≥ti。

機(jī)隊(duì)中飛機(jī)設(shè)計(jì)型號(hào)存在差異，同一部件可能不完全相同，據(jù)極大似然估計(jì)函數(shù)，由形狀參數(shù)獲得壽命參數(shù)，即η=f(ti,β),如式(6)所示。

(5)

(6)

式中：m為機(jī)隊(duì)故障飛機(jī)數(shù)。當(dāng)β確定時(shí)，故障飛機(jī)數(shù)越少，機(jī)隊(duì)特征壽命越長。由式(4)～式(6)可得Ai為

取β=2,4，tR/η=2/3,1,2，Ai隨ti/tR變化趨勢(shì)如圖4所示。

圖4 β與tR/η對(duì)Ai的影響Fig.4 Effect of β and tR/η on contribution to Ai

當(dāng)ti/tR=0時(shí)，Ai最大，新飛機(jī)有最高的損傷率。當(dāng)ti/tR=1時(shí)，Ai為零，飛機(jī)退役不再運(yùn)行，故障發(fā)生率為零。當(dāng)η/tR≥1時(shí)，退役壽命低于特征壽命，機(jī)隊(duì)初始運(yùn)行時(shí)損傷率Ai較低，隨飛機(jī)飛行循環(huán)增加不斷降低，較少在退役之前發(fā)生故障。η/tR=1.5，β=4，對(duì)應(yīng)的曲線，在飛機(jī)交付運(yùn)行時(shí)，Ai-max僅為0.06，表明即使是剛交付的新飛機(jī)，在整個(gè)壽命期僅有6%的概率發(fā)生故障；當(dāng)η/tR≤1時(shí)，退役壽命高于特征壽命，機(jī)隊(duì)至少63.2%的飛機(jī)在退役前發(fā)生故障。η/tR相同時(shí)，Ai隨β增大而升高。

2 發(fā)動(dòng)機(jī)渦輪盤榫齒故障風(fēng)險(xiǎn)評(píng)估

2.1 單機(jī)風(fēng)險(xiǎn)評(píng)估

某型渦輪發(fā)動(dòng)機(jī)維修間隔為1 000 h，在使用中由于Ⅱ級(jí)渦輪盤榫齒裂紋斷裂導(dǎo)致葉片甩出，曾發(fā)生多次事故[19-20]，文獻(xiàn)[21]指出渦輪盤榫齒在工作時(shí)承受葉片離心力引起的拉伸應(yīng)力、彎曲應(yīng)力、剪切應(yīng)力、齒面上的擠壓應(yīng)力及溫度分布不均勻引起的熱應(yīng)力，容易出現(xiàn)疲勞斷裂。通過故障數(shù)據(jù)分析，渦輪盤榫齒裂紋故障概率服從威布爾分布。由文獻(xiàn)[6]確定雙參數(shù)β=2，η=1 418/h。由式(2)獲得渦輪榫齒故障率曲線h(tI)，如圖5所示， 發(fā)動(dòng)機(jī)運(yùn)行時(shí)間越長， 故障率上升速率越快，在到達(dá)維修間隔1 000 h時(shí)，渦輪盤榫齒故障率為7.0×10-4/h。

圖5 渦輪盤榫齒故障率Fig.5 Failure rate of dovetail gear in turbine disk

由圖2可知，通過現(xiàn)有的電位法探傷等檢查方法檢查渦輪盤榫齒裂紋故障時(shí)存在漏檢的情況，不容易被檢出，在造成災(zāi)難性后果前難以被檢出，結(jié)構(gòu)設(shè)計(jì)符合損傷容限理論，渦輪榫齒在相似應(yīng)力水平上無相似結(jié)構(gòu)運(yùn)行，不受廣泛的疲勞故障影響，具有多條加載路徑設(shè)計(jì)，μ值取0.3。

榫齒裂紋因果鏈如圖6所示，渦輪盤榫齒產(chǎn)生裂紋概率P1與當(dāng)前運(yùn)行循環(huán)有關(guān)。當(dāng)ti/tR分別接近于20%、40%、60%、80%時(shí)，P1分別為0.8、0.5、0.1、0.01。通過蒙特卡羅方法，對(duì)渦輪盤榫齒開始出現(xiàn)裂紋故障時(shí)間，模擬100 000次，tR=1 000 h，β=2，η=1 418/h， 預(yù)測(cè)渦輪盤榫齒在運(yùn)行661 h時(shí)生成可探測(cè)到的裂紋，ti/tR=661/1 000= 66%，則P1=0.1。由榫齒裂紋導(dǎo)致發(fā)生不安全后果的各狀態(tài)條件概率由歷史統(tǒng)計(jì)數(shù)據(jù)獲得。渦輪盤榫齒故障導(dǎo)致4種不安全后果中的一種或多種：空中解體，墜毀，機(jī)身基本完整、存在人員傷亡，未正常起飛滑行、偏離跑道。4種不安全后果的嚴(yán)重性等級(jí)為死亡數(shù)與總?cè)藬?shù)之比，如圖6所示。

圖6 渦輪盤榫齒裂紋因果鏈Fig.6 Causal chain of crack on dovetail gear in turbine disk

ψIfI=1.6×10-4,1.37×10-4,4×10-9,3×10-7,2.45×10-7,1×10-11,5×10-9,4.9×10-8,5×10-11,1.5×10-8。

∑(ψIfI)=2.98×10-4，為單架飛機(jī)當(dāng)前由于榫齒裂紋故障產(chǎn)生的所有不安全后果概率。

單機(jī)風(fēng)險(xiǎn)RI=h(tI)μIψIfI=h(tI)×0.3×2.98×10-4=8.94×10-5h(tI),如圖7所示。該發(fā)動(dòng)機(jī)在下次維修之前，渦輪盤榫齒故障導(dǎo)致的后果風(fēng)險(xiǎn)小于部件故障可接受風(fēng)險(xiǎn)水平10-7/飛行小時(shí)[2]。

圖7 單機(jī)故障風(fēng)險(xiǎn)水平Fig.7 Failure level of individual risk

2.2 機(jī)隊(duì)風(fēng)險(xiǎn)評(píng)估

某機(jī)隊(duì)在役飛機(jī)有507架，當(dāng)前飛行循環(huán)如圖8所示。檢出故障飛機(jī)數(shù)m取1、3、5、7、9，tR=50 000飛行循環(huán),70 000飛行循環(huán)，分析退役時(shí)間不同對(duì)機(jī)隊(duì)運(yùn)行安全的影響。

由式(6)，β=2，m=1，tR=50 000飛行循環(huán)時(shí),有

由式(4)可得

通過蒙特卡羅方法，對(duì)渦輪盤榫齒產(chǎn)生裂紋時(shí)間，模擬100 000次。

當(dāng)tR=50 000飛行循環(huán)，β=2，η=635 821時(shí)，模擬結(jié)果如圖9所示。

預(yù)測(cè)渦輪盤榫齒在運(yùn)行33 311飛行循環(huán)時(shí)生成可探測(cè)的裂紋，ti/tR=66.7%，P1=0.1。

圖8 機(jī)隊(duì)飛行循環(huán)統(tǒng)計(jì)圖Fig.8 Fleet data of flight cycle

圖9 預(yù)測(cè)渦輪盤榫齒故障的飛行循環(huán)Fig.9 Predicted flight cycle of cracking on dovetail gear in turbine disk

(ψFfF)max=0.000 16

∑(ψFfF)=0.000 298

RF=εμFψFfF=

2.1×0.3×0.000 298=0.000 188

當(dāng)tR=70 000飛行循環(huán)時(shí)，ti=33 636飛行循環(huán)，ti/tR=48.05%，P1=0.5。

(ψFfF)max=0.000 6

∑(ψFfF)=0.001 488

RF=εμFψFfF=

5.1×0.3×0.001 488=0.002 28

相應(yīng)η、ε和RF,如表1所示。

分別取tR=50 000,70 000，生成圖10所示曲線RF=g(m),給出風(fēng)險(xiǎn)隨故障飛機(jī)數(shù)變化趨勢(shì)，民用飛機(jī)機(jī)隊(duì)運(yùn)行風(fēng)險(xiǎn)可接受水平為0.02[2]。風(fēng)險(xiǎn)水平與退役壽命密切相關(guān)，退役壽命越長，機(jī)隊(duì)故障風(fēng)險(xiǎn)越高。

表1 機(jī)隊(duì)故障風(fēng)險(xiǎn)參數(shù)(β=2)Table 1 Failure risk parameters of fleet(β=2)

圖10 機(jī)隊(duì)故障風(fēng)險(xiǎn)水平Fig.10 Failure level of fleet risk

取(β,m)=(2,9),(2,1)。由表1或圖11獲得預(yù)測(cè)飛機(jī)出現(xiàn)渦輪盤榫齒故障數(shù)。表1中β=2，m=9，tR=50 000飛行循環(huán)時(shí)，預(yù)測(cè)故障飛機(jī)數(shù)量為18.8架，當(dāng)ti=tR時(shí)，存在9+18.8≈28架飛機(jī)可能出現(xiàn)渦輪盤榫齒故障；由式(3)，F(xiàn)(t)隨飛行循環(huán)變化趨勢(shì)如圖11所示，tR=50 000飛行循環(huán)，β=2，m=9，ti=tR時(shí)，機(jī)隊(duì)累積故障率為5.41%，因此有507×5.41%≈28架飛機(jī)可能出現(xiàn)渦輪盤榫齒故障。

圖11 威布爾累積分布函數(shù)Fig.11 Weibull cumulative distribution functions

3 結(jié) 論

威布爾分布很好地建立了單機(jī)運(yùn)行故障率曲線和機(jī)隊(duì)隨飛行循環(huán)變化的風(fēng)險(xiǎn)水平。對(duì)發(fā)動(dòng)機(jī)渦輪盤榫齒故障進(jìn)行單機(jī)和機(jī)隊(duì)風(fēng)險(xiǎn)計(jì)算，單機(jī)運(yùn)行風(fēng)險(xiǎn)取決于部件的故障率和未檢出概率；機(jī)隊(duì)運(yùn)行風(fēng)險(xiǎn)由機(jī)隊(duì)當(dāng)前運(yùn)行規(guī)模及飛行循環(huán)確定。

1) 單機(jī)渦輪盤榫齒故障風(fēng)險(xiǎn)水平低于可接受水平10-7/飛行小時(shí)，但隨著運(yùn)行時(shí)間接近退役壽命，故障風(fēng)險(xiǎn)水平接近10-7/飛行小時(shí)。

2) 機(jī)隊(duì)渦輪盤榫齒故障風(fēng)險(xiǎn)在故障飛機(jī)增至9架，服役壽命為70 000飛行循環(huán)時(shí)，風(fēng)險(xiǎn)水平超過機(jī)隊(duì)風(fēng)險(xiǎn)可接受水平0.02。

3) 若當(dāng)前退役壽命下對(duì)應(yīng)風(fēng)險(xiǎn)值已超出上限值，則需根據(jù)風(fēng)險(xiǎn)計(jì)算結(jié)果采取相應(yīng)糾正措施。

[1] CABLER J M.8110.107A Monitor safety/analyze data[S]. Washington, D.C.: Federal Aviation Administration, 2012.

[2] U.S. Department of Transportation. Advisory circular 39-8 Continued airworthiness assessments of powerplant and auxiliary power unit installations of transport category airplanes[S]. Washington, D.C.: Federal Aviation Administration, 2003.

[3] Transport Airplane Directorate. PS-ANM-25-05 Transport airplane risk assessment methodology handbook[S]. Washington, D.C.: Federal Aviation Administration, 2011.

[4] BRISTOL T. Safety management system manual[S]. Washington, D.C.: Federal Aviation Administration, 2014.

[5] GALILEO T, MICAELA D. Risk assessment techniques for civil aviation security[J]. Reliability Engineering and System Safety, 2011, 96(8): 892-899.

[6] 李龍彪, 畢蘇藝, 孫有朝. 一種基于故障統(tǒng)計(jì)數(shù)據(jù)的發(fā)動(dòng)機(jī)風(fēng)險(xiǎn)預(yù)測(cè)[J]. 南京航空航天大學(xué)學(xué)報(bào), 2015, 47(4): 559-565.

LI L B, BI S Y, SUN Y C. Risk prediction of aero-engine based on failure statistics data[J]. Journal of Nanjing University of Aeronautics and Astronautics, 2015, 47(4): 559-565 (in Chinese).

[7] HEDO J M, Val M R. Assessment of narrow-body transport airplane evacuation by numerical simulation[J]. Journal of Aircraft, 2011, 48(5): 1785-1794.

[8] 徐慶宏, 孫有朝, 李龍彪. 航空發(fā)動(dòng)機(jī)部件故障風(fēng)險(xiǎn)評(píng)估方法研究[J]. 機(jī)械科學(xué)與技術(shù), 2016, 35(8): 1291-1296.

XU Q H, SUN Y C, LI L B. Study on risk assessment method for aero-engine component failure[J]. Mechanical Science and Technology for Aerospace Engineering, 2016, 35(8): 1291-1296 (in Chinese).

[9] 郭媛媛, 李龍彪, 胡宇群, 等. 民用飛機(jī)運(yùn)行風(fēng)險(xiǎn)評(píng)估方法[J]. 航空工程進(jìn)展, 2016, 7(2): 253-258.

GUO Y Y, LI L B, HU Y Q, et al. Operation risk assessment method of civil aircraft[J]. Advances in Aeronautical Science and Engineering, 2016, 7(2): 253-258 (in Chinese).

[10] 王健名, 徐浩軍, 薛源, 等. 基于極值理論的平尾結(jié)冰飛行風(fēng)險(xiǎn)評(píng)估[J]. 航空學(xué)報(bào), 2016, 37(10): 3011-3022.

WANG J M, XU H J, XUE Y, et al. Flight risk evaluation of tailplane icing based on extreme value theory[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(10): 3011-3022 (in Chinese).

[11] 李巖, 張曙光, 宮綦. 一種改進(jìn)的航空發(fā)動(dòng)機(jī)結(jié)構(gòu)概率風(fēng)險(xiǎn)評(píng)估方法[J]. 航空學(xué)報(bào), 2016, 37(2): 597-608.

LI Y, ZHANG S G, GONG Q. An improved probabilistic risk assessment method of structural parts for aeroengine[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(2): 597-608 (in Chinese).

[12] 薛源, 徐浩軍, 胡孟權(quán). 結(jié)冰條件下人-機(jī)-環(huán)系統(tǒng)的飛行風(fēng)險(xiǎn)概率[J]. 航空學(xué)報(bào), 2016, 37(11): 3328-3339.

XUE Y, XU H J, HU M Q. Flight risk probability of pilot-aircraft-environment system under icing conditions[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(11): 3328-3339 (in Chinese).

[13] 孫楊慧, 楊坤, 侯乃先, 等. 渦輪盤概率風(fēng)險(xiǎn)的評(píng)估方法[J]. 航空動(dòng)力學(xué)報(bào), 2016, 31(9): 2195-2202.

SUN Y H, YANG K, HOU N X, et al. Probabilistic risk assessment method for turbine disk[J]. Journal of Aerospace Power, 2016, 31(9): 2195-2202 (in Chinese).

[14] 凌丹. 威布爾分布模型及其在機(jī)械可靠性中的應(yīng)用研究[D]. 成都: 電子科技大學(xué), 2010: 3-24.

LING D. Research on Weibull distribution and its applications in mechanical reliability engineering[D]. Chengdu: University of Electronic Science and Technology of China, 2010: 3-24 (in Chinese).

[15] 趙洪利, 劉宇文. 基于蒙特卡羅模擬的航空發(fā)動(dòng)機(jī)故障風(fēng)險(xiǎn)預(yù)測(cè)[J]. 北京航空航天大學(xué)學(xué)報(bào), 2015, 41(3): 545-550.

ZHAO H L, LIU Y W. Forecasting for aero-engine failure risk based on Monte Carlo simulation[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(3): 545-550 (in Chinese).

[16] Society of Automotive Engineers. ARP5150 safety assessment of transport airplanes in commercial service[S]. SAE, 2003.

[17] ABERNETHY R B. The new weibull handbook[M]. 5th ed. Florida: Society of Automotive Engineers International, 2010: 2-11.

[18] VIOLETTE M G, SAFARIAN P, HAN N, et al. Transport airplane risk analysis[J]. Journal of Aircraft, 2015, 52(2): 395-402.

[19] 何衛(wèi)鋒, 李應(yīng)紅, 聶祥樊, 等. 激光沖擊葉片榫頭變形控制與疲勞試驗(yàn)[J]. 航空學(xué)報(bào), 2014, 35(7): 2041-2048.

HE W F, LI Y H, NIE X F, et al. Deformation control and fatigue test of blade tenon by laser shock peening[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(7): 2041-2048 (in Chinese).

[20] 孫瑞杰, 閆曉軍. 渦輪葉片榫齒部位疲勞/蠕變?cè)囼?yàn)的新特點(diǎn)[J]. 航空動(dòng)力學(xué)報(bào), 2007, 22(3): 419-424.

SUN R J, YAN X J. New characteristics of fatigue-creep tests on serration of turbine blades[J]. Journal of Aerospace Power, 2007, 22(3): 419-424 (in Chinese).

[21] 劉濱春, 王慧. 某型發(fā)動(dòng)機(jī)渦輪盤榫齒裂紋分析[J]. 工程與試驗(yàn), 2012, 52(2): 26-27.

LIU B C, WANG H. Analysis of serration crack in turbine disc for engineer[J]. Engineering and Test, 2012, 52(2): 26-27 (in Chinese).

FailureriskassessmentmethodofcivilaircraftbasedonMonteCarlomethod

GUOYuanyuan,SUNYouchao*,LILongbiao

CollegeofCivilAviation,NanjingUniversityofAeronauticsandAstronautics,Nanjing210016,China

Inordertoimproveaviationsafetyandreducetherisklevelofpartsfailureintheprocessofcivilaircraftoperation,theindividualandfleetriskshavebeencalculated.Thekeypointsofindividualriskarehowtoconfirmthefailuredistributionandrate,whichturnsoutthatWei-bulldistributionisapplicable.Accordingtothefailuremechanismofpart,thenon-detectionprobabilityandtheconditionalprobabilitythedetectwillleadstoanunsafeoutcomeaswellasrelatedinjuryratehavebeenlisted,inordertodeterminetherisklevel.Basedontheindividualrisk,thekeypointsoffleetriskarefleetscaleandcurrentflightcircles.Thenumberoffailureflightshasbeenpredictedbasedonoperationaldata.Theriskvaluehasbeenacquiredunderdifferentoperatedsituationsuchasparameterdistribution,numberoffailureflightsaswellasretirementcycles.ThefailuretimeofindividualandfleethasbeensimulatedbyMonteCarlomethodinthecaseofserrationcrackinturbinediscofaero-engine.Therisklevelsofindividualandfleetaboutthisparthavebeenanalyzed,andthetwokindsrisklevelshavebeengivenatcontinuedairworthinessstagebasedonoperationaldata.

civilaircraft;fleetrisk;individualrisk;MonteCarlomethod;cumulativedistributionfunction;failurerate

2017-01-12；Revised2017-03-02；Accepted2017-04-18;Publishedonline2017-05-031644

URL：http://hkxb.buaa.edu.cn/CN/html/20171015.html

s：JointlyfundofNationalNaturalScienceFoundationofChinaandCivilAviationAdministrationofChina(U1333119);DefenseIndustrialTechnologyDevelopmentProgram(JCKY2013605B002);CivilAircraftSpecialFoundationofMinistryofIndustryandInformationTechnology(MJ-F-2011-33);FundofShanghaiEngineeringResearchforCivilAircraftHealthMonitoring(GCZX-2015-05);TechnicalProjectFoundationofStateAdministrationofScience,TechnologyandIndustryforNationalDefense(Z052013B003)

.E-mailsunyc@nuaa.edu.cn

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2017.221126

V239

A

1000-6893(2017)10-221126-09

2017-01-12；退修日期2017-03-02；錄用日期2017-04-18；< class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2017-05-031644

http://hkxb.buaa.edu.cn/CN/html/20171015.html

國家自然科學(xué)基金與中國民用航空局聯(lián)合資助(U1333119)； 國防基礎(chǔ)科研計(jì)劃(JCKY2013605B002)； 工信部民機(jī)專項(xiàng)(MJ-F-2011-33)； 上海民用飛機(jī)健康監(jiān)控工程技術(shù)研究中心基金(GCZX-2015-05)； 國防科工局技術(shù)基礎(chǔ)科研項(xiàng)目(Z052013B003)

.E-mailsunyc@nuaa.edu.cn

郭媛媛，孫有朝，李龍彪．基于蒙特卡羅方法的民用飛機(jī)故障風(fēng)險(xiǎn)評(píng)估方法J．航空學(xué)報(bào),2017,38(10):221126.GUOYY,SUNYC,LILB.FailureriskassessmentmethodofcivilaircraftbasedonMonteCarlomethodJ.ActaAeronauticaetAstronauticaSinica,2017,38(10):221126.

(責(zé)任編輯：張晗)