趙洪利 張 猛

(中國民航大學(xué)航空工程學(xué)院 天津 300300)

0 引 言

為了保障飛機(jī)的正常運(yùn)行，航空公司必須為飛機(jī)配備一定量的備件，但這會帶來高額的成本，大約占飛機(jī)運(yùn)營成本的10%左右[1]。作為最為昂貴的備件——備用發(fā)動機(jī)(備發(fā))的數(shù)量一直是航空公司重點關(guān)注的問題，備發(fā)數(shù)量過多會帶來巨額的保障成本，阻礙航空公司的現(xiàn)金流動；備發(fā)數(shù)量過少會造成飛機(jī)停場(AOG)等一系列問題，從而影響航空公司的正常運(yùn)營。因此適當(dāng)?shù)膫浒l(fā)數(shù)量對實現(xiàn)最大化備發(fā)利用率和最小化保障成本至關(guān)重要[2]。

目前，備發(fā)率一般是航空公司用來衡量備發(fā)數(shù)量的一個重要指標(biāo)，是指該機(jī)隊中的備發(fā)數(shù)量與機(jī)隊在翼發(fā)動機(jī)總數(shù)的比值。就備發(fā)數(shù)量的計算問題，大致可分為三種方法：(1) 基于統(tǒng)計學(xué)理論：陳鳳騰等[3]使用改進(jìn)的比率危險模型，結(jié)合高斯公式計算了備發(fā)數(shù)量；白芳等[4]同樣使用高斯公式計算了備發(fā)數(shù)量，通過維修保障成本模型優(yōu)化了備發(fā)數(shù)。(2) 基于算法：趙洪利等[5]基于間隔控制算法，通過固定備發(fā)數(shù)量得出了其對應(yīng)的缺少發(fā)動機(jī)(缺發(fā))時間，并驗證了算法的優(yōu)越性；Sheng等[6]采用遺傳算法，通過仿真尋優(yōu)找到了機(jī)隊備件供應(yīng)的最優(yōu)解決方案以及最小化成本。(3) 基于經(jīng)驗公式：GE公司結(jié)合返廠率(即返廠次數(shù)與飛行時間的比值)、周轉(zhuǎn)周期(即發(fā)動機(jī)維修所用的時間)等一系列因素，給出了計算備發(fā)數(shù)量的經(jīng)驗公式。

但是以上方法均存在考慮發(fā)動機(jī)下發(fā)因素單一，無法給出確切的備發(fā)規(guī)劃等問題，必然會對最終結(jié)果產(chǎn)生影響。針對以上問題，本文使用離散事件仿真的方法，利用MATLAB軟件建立了備發(fā)仿真模型。基于排隊論思想，采用固定增量的時間推進(jìn)機(jī)制，在先進(jìn)先出排隊規(guī)則下，考慮維修廠(Maintenance, Repair and Overhaul, MRO)維修能力的限制，以每三年為一個規(guī)劃期為機(jī)隊制定了15年的備發(fā)規(guī)劃。

1 備發(fā)仿真模型設(shè)計

1.1 模型基本假設(shè)

為準(zhǔn)確地展現(xiàn)發(fā)動機(jī)“在翼—下發(fā)—在翼”的使用流程，在建立仿真模型之前需做出一些基本假設(shè)。

假設(shè)1：發(fā)動機(jī)在使用過程中存在性能衰退、硬件損傷與偶然性損傷多種下發(fā)因素，各因素相互獨立，且故障時間服從一定分布。

假設(shè)2：發(fā)動機(jī)性能恢復(fù)視為非完美維修，性能恢復(fù)后發(fā)動機(jī)排氣溫度裕度(Exhaust Gas Temperature Margin，EGTM)服從一定分布。

假設(shè)3：若有缺發(fā)的情況發(fā)生，采取租賃發(fā)動機(jī)的方式保障航空公司正常運(yùn)營。

假設(shè)4：MRO維修發(fā)動機(jī)時，采取先進(jìn)先出的維修排隊機(jī)制。

1.2 輸入輸出模型參數(shù)

基于以上所有假設(shè)，建立備發(fā)仿真模型，模型輸入輸出參數(shù)包括：

(1) 機(jī)隊規(guī)模：假設(shè)機(jī)隊中共25架全新雙發(fā)飛機(jī)共50臺發(fā)動機(jī)，且分梯次引進(jìn)機(jī)隊，每10天引進(jìn)一架。

(2) 使用情況：根據(jù)某航空公司實際使用情況，設(shè)定飛機(jī)一天運(yùn)行5循環(huán)，一年運(yùn)行1 800循環(huán)，使用壽命為15年。

(3) 周轉(zhuǎn)周期與MRO維修能力：發(fā)動機(jī)周轉(zhuǎn)周期為75天且MRO可同時維修10臺發(fā)動機(jī)。

(4) 備發(fā)數(shù)量：備發(fā)數(shù)量起始設(shè)為1臺，若不能滿足保障率要求時增加備發(fā)數(shù)量。

模型輸出參數(shù)包括：

(1) 總?cè)卑l(fā)時間Tlack。

(2) 各規(guī)劃期的保障率ω：

(1)

式中：M為機(jī)隊發(fā)動機(jī)數(shù)量;P為規(guī)劃期長度;Tlack1是各規(guī)劃期內(nèi)的缺發(fā)時間。

1.3 仿真方法與流程

離散事件仿真系統(tǒng)是狀態(tài)僅在離散時間點上發(fā)生變化的仿真方法。事件被初始化后，通過處理在仿真時間軸上不連續(xù)點處發(fā)生的帶有時間戳的事件，進(jìn)行所有特定計算。所使用的MATLAB軟件是一款強(qiáng)大的數(shù)學(xué)計算軟件，已被廣泛地運(yùn)用在算法開發(fā)、程序編寫等方面。在MATLAB仿真平臺上，利用.m文件進(jìn)行程序編寫，建立離散事件備發(fā)仿真模型，流程如圖1所示。

圖1 備發(fā)仿真流程

(1) 初始化各項輸入?yún)?shù)、Tlack計數(shù)器和Tlack1計數(shù)器。

(2) 編寫性能衰退、硬件損傷和偶然性損傷三種下發(fā)因素隨機(jī)數(shù)產(chǎn)生函數(shù)并調(diào)用，使用sort函數(shù)將產(chǎn)生的隨機(jī)數(shù)排序，給定機(jī)隊所有發(fā)動機(jī)故障下發(fā)時間戳。

(3) 時鐘推進(jìn)1個循環(huán)，遍歷所有發(fā)動機(jī)判斷是否到達(dá)下發(fā)時間戳。若是，再次推進(jìn)1個循環(huán)；若否，發(fā)動機(jī)進(jìn)廠維修并判斷MRO工位情況。同時使用switch函數(shù)判斷故障類型并記錄。

(4) 若MRO有工位，則工位數(shù)-1，判斷備發(fā)情況；若MRO無工位，發(fā)動機(jī)排隊等待工位出現(xiàn)。

(5) 若有備發(fā)，則備發(fā)數(shù)-1，重新產(chǎn)生故障下發(fā)時間，并判斷發(fā)動機(jī)壽命情況；若無備發(fā)，記錄缺發(fā)開始時刻。等待維修完成后，記錄缺發(fā)結(jié)束時刻，并統(tǒng)計缺發(fā)時間Tlack。

(6) 若壽命不終止，時鐘再次推進(jìn)1個循環(huán)；若終止，輸出各規(guī)劃期保障率和總?cè)卑l(fā)時間Tlack，單次仿真結(jié)束，輸出故障記錄。

(7) 多次重復(fù)仿真，使用plot函數(shù)作圖查看各規(guī)劃期內(nèi)的保障率平均收斂情況。經(jīng)測試，仿真次數(shù)為500次時收斂情況較好。

2 下發(fā)策略與性能恢復(fù)

導(dǎo)致發(fā)動機(jī)下發(fā)的因素是很復(fù)雜的，本文考慮發(fā)動機(jī)的性能衰退、硬件損傷和偶然性損傷三種主要因素競爭下發(fā)，即單次發(fā)動機(jī)下發(fā)時間戳t為：

t=min{t1,t2,t3}

(2)

式中：t1為性能衰退下發(fā)時間隨機(jī)數(shù)；t2為硬件損傷下發(fā)時間隨機(jī)數(shù)；t3為偶然性損傷下發(fā)時間隨機(jī)數(shù)。

2.1 性能衰退偽隨機(jī)數(shù)

2.1.1性能衰退模型建立

發(fā)動機(jī)的性能衰退，一般會在發(fā)動機(jī)EGTM的衰退情況上有所反映。EGTM定義為發(fā)動機(jī)排氣溫度(Exhaust Gas Temperature，EGT)紅線值與在海平面拐點溫度下全功率起飛時的EGT之差[7]。隨著發(fā)動機(jī)的使用， EGTM值會逐漸減小直至到達(dá)警戒值。本文將采用隨機(jī)wiener過程，對發(fā)動機(jī)EGTM衰退過程建模。

發(fā)動機(jī)EGTM在t時刻的衰退量為X(t),假設(shè)該衰退量服從wiener過程{X(t),t≥0}有：

X(t)=X(0)+λt+σB(t)

(3)

T=inf{t1|X(t1)≥Df}

(4)

式中：Df=S-L,Df是EGTM衰退閾值，S是發(fā)動機(jī)EGTM的起始數(shù)值，L是EGTM的警戒值;t1是故障下發(fā)時間；inf表示集合的下限值。此時該隨機(jī)wiener過程的失效函數(shù)為[9]：

(5)

式中：Φ表示標(biāo)準(zhǔn)正態(tài)分布的累積概率密度函數(shù)。

2.1.2模型參數(shù)計算與偽隨機(jī)數(shù)生成

表1是某航空公司機(jī)隊10臺某系列發(fā)動機(jī)實際使用過程中的EGTM衰退數(shù)據(jù)。

表1 某機(jī)隊發(fā)動機(jī)EGTM樣本數(shù)據(jù)

由于該隨機(jī)wiener過程的失效函數(shù)較為復(fù)雜，難以直接生成符合該函數(shù)的隨機(jī)數(shù)，本文采用牛頓迭代法求解性能衰退下發(fā)時間。牛頓法是一種用來求非線性方程解的迭代算法，在給定初值合理的情況下，可以很快求解出方程解[11]，本文中所用迭代公式如下：

(6)

2.2 硬件損傷偽隨機(jī)數(shù)

2.2.1硬件損傷模型建立

航空發(fā)動機(jī)在使用過程中，主要的故障包括低周疲勞、熱機(jī)械疲勞、蠕變、腐蝕等[12]，這些損傷通常會損耗部件壽命，造成發(fā)動機(jī)下發(fā)。本文采用威布爾分布來表征該種故障的發(fā)生，其失效函數(shù)為：

(7)

式中：t2為故障時間；β為形狀參數(shù)；η為尺度參數(shù)。

2.2.2模型參數(shù)計算與偽隨機(jī)數(shù)生成

根據(jù)某機(jī)隊發(fā)動機(jī)實際下發(fā)數(shù)據(jù)顯示，發(fā)動機(jī)硬件損傷常見的故障分為三種：高壓壓氣機(jī)(HPC)葉片損傷、燃燒室(Combustor)損傷和高壓渦輪(HPT)一級葉片損傷，分別提取故障下發(fā)時間如表2所示。

表2 實際發(fā)動機(jī)下發(fā)時間

根據(jù)表2中數(shù)據(jù)采用中位秩回歸參數(shù)估計法估算模型參數(shù)，結(jié)果如表3所示。

表3 模型參數(shù)

使用反變換法，計算威布爾分布的失效時間t2為：

(8)

式中：F(t2)是區(qū)間[0,1]內(nèi)的均勻分布隨機(jī)數(shù)，同樣由rand函數(shù)產(chǎn)生，計算生成一系列硬件損傷下發(fā)時間偽隨機(jī)數(shù)t2。

2.3 偶然性損傷偽隨機(jī)數(shù)

2.3.1偶然性損傷模型建立

發(fā)動機(jī)的偶然性損傷，是指發(fā)動機(jī)在使用過程中由于鳥擊、雷擊等突發(fā)事件，造成發(fā)動機(jī)故障下發(fā)，且這類故障具有無記憶性的特點。所以本文采用指數(shù)分布模型，對發(fā)動機(jī)偶然性損傷下發(fā)時間建模，其失效函數(shù)為：

F(t3)=1-exp(-λ3t3)

(9)

式中：t3為故障時間點；λ3為故障率。

2.3.2模型參數(shù)計算與偽隨機(jī)數(shù)生成

根據(jù)最新的《鳥擊航空器信息分析報告》，能夠造成發(fā)動機(jī)下發(fā)的偶然性損傷故障率λ3=4.335×10-6/cycle。

同樣使用反變換法，生成偶然性損傷下發(fā)時間偽隨機(jī)數(shù)t3，具體過程如上文，這里不再贅述。

2.4 性能恢復(fù)EGTM衰退閾值偽隨機(jī)數(shù)

(10)

(11)

式中：c為修正系數(shù)，保證Δ在區(qū)間[0,Df]內(nèi)累積概率為1；F(Δ)為該截斷正態(tài)分布的累積概率。

采用正態(tài)分布擬合該型號發(fā)動機(jī)性能恢復(fù)后EGTM衰退閾值Δ具體如表4所示，根據(jù)式(10)計算c=1.015。同樣使用牛頓法，迭代計算性能恢復(fù)后EGTM衰退閾值偽隨機(jī)數(shù)。

表4 性能恢復(fù)后EGTM衰退閾值

3 模型驗證

圖2是做500仿真的平均總?cè)卑l(fā)時間點，該值收斂于45 100循環(huán)。由于總?cè)卑l(fā)時間點直接影響保障率大小，所以這里對假設(shè)H0:Tlack=45 100 cycle做t檢驗，以驗證模型的正確性。

圖2 平均總?cè)卑l(fā)時間點仿真結(jié)果

重新仿真產(chǎn)生10組Tlack數(shù)據(jù)如表5所示，取置信度為95%，雙邊檢驗臨界值t0.025,9=2.26。經(jīng)計算t0=2.008

表5 10次Tlack抽樣數(shù)據(jù)

4 仿真結(jié)果與分析

4.1 仿真結(jié)果

一般航空公司對于保障率的要求在95%～97%之間[15]。本文取97%的保障率標(biāo)準(zhǔn)，1臺備發(fā)時各規(guī)劃期內(nèi)的保障率仿真結(jié)果如圖3所示。

圖3 1臺備發(fā)保障率仿真結(jié)果

仿真結(jié)果顯示，1臺備發(fā)只能滿足1～3年內(nèi)的機(jī)隊運(yùn)營，其余規(guī)劃期內(nèi)保障率均在97%以下。想要滿足整個壽命階段的保障率要求必須要增加備發(fā)數(shù)量，因此依次增加備發(fā)數(shù)量重新仿真，結(jié)果如圖4-圖6所示。

圖4 2臺備發(fā)保障率仿真結(jié)果

圖5 3臺備發(fā)保障率仿真結(jié)果

圖6 4臺備發(fā)保障率仿真結(jié)果

依據(jù)以上所有仿真結(jié)果，當(dāng)4臺備發(fā)時滿足各階段的保障率要求，現(xiàn)給出各個規(guī)劃期內(nèi)的最佳備發(fā)數(shù)如表6所示。

表6 各規(guī)劃期內(nèi)最佳備發(fā)數(shù)量

4.2 結(jié)果分析

由于1～3年屬于發(fā)動機(jī)的首次壽命階段，發(fā)動機(jī)在翼時間相對較長，因此少有下發(fā)，僅需要1臺備發(fā)；4～6年內(nèi)，由于發(fā)動機(jī)初始壽命不夠分散，發(fā)動機(jī)的首次下發(fā)多發(fā)生在該階段，出現(xiàn)集中下發(fā)的現(xiàn)象，因此所需備發(fā)數(shù)量最多為4臺；7～15年內(nèi)，發(fā)動機(jī)已進(jìn)入成熟期，各階段下發(fā)次數(shù)以及發(fā)動機(jī)在翼時間已接近穩(wěn)定，因此保障率幾乎維持不變，所需備發(fā)數(shù)也減少至3臺。

目前航空公司的備發(fā)率一般維持在5%～8%。對比本文最多備發(fā)數(shù)為4臺即備發(fā)率為8%的仿真結(jié)果，具有合理性，且減少了1～3年備發(fā)數(shù)量過多造成的成本浪費，4～6年也可以在購買3臺備發(fā)的基礎(chǔ)上額外租發(fā)1臺以滿足備發(fā)需求，使備發(fā)規(guī)劃更為靈活。

4.3 維修能力敏感度分析

設(shè)定備發(fā)數(shù)量為4臺，將MRO維修能力改為5臺和15臺，對比維修能力為10臺時各規(guī)劃期的保障率，具體結(jié)果如圖7所示。

圖7 維修能力敏感度分析圖

由圖7可得出以下結(jié)論：

(1) 當(dāng)維修能力為5臺時，第2個規(guī)劃期即4～6年的保障率很低，小于維修能力為10臺和15臺同階段保障率，說明5臺的維修能力是不能滿足該機(jī)隊的維修需求的。

(2) 維修能力為10臺和15臺在各個規(guī)劃期內(nèi)的保障率相差都不大，因此10臺的維修能力即可滿足該機(jī)隊的維修需求，工位不夠而出現(xiàn)維修等待的現(xiàn)象很少發(fā)生。

5 結(jié) 語

本文采用離散系統(tǒng)仿真的方法，依據(jù)排隊理論建立了備發(fā)仿真模型。分規(guī)劃期給出了機(jī)隊15年內(nèi)不同時間段的最佳備發(fā)數(shù)量，使備發(fā)規(guī)劃更加細(xì)致準(zhǔn)確且靈活，減少了備發(fā)數(shù)量過多造成的成本浪費，同時對MRO的維修能力進(jìn)行評估。對比航空公司現(xiàn)行的備發(fā)率，本文結(jié)論較為準(zhǔn)確。