陳小磊，郭迎清，閆星輝，姜彩虹

（1.西北工業(yè)大學動力與能源學院，西安710072；2.中航工業(yè)無錫航空動力控制系統(tǒng)研究所，江蘇無錫214100；3.中航工業(yè)沈陽發(fā)動機設計研究所，沈陽110015）

0 引言

隨著新一代航空發(fā)動機對經(jīng)濟性要求的提高，壽命延長控制逐漸成為航空發(fā)動機控制領域的研究熱點。其研究關鍵技術之一是準確計算部件壽命，旨在保證發(fā)動機性能不變或有少量損失的同時，大幅度延長發(fā)動機及其部件壽命[1-2]。航空發(fā)動機部件壽命計算非常復雜，其難點在于如何權(quán)衡發(fā)動機安全運行和部件使用最大化，保證發(fā)動機安全性的同時降低使用成本[3-4]。在發(fā)動機制造過程中，利用先進的工具分析所有部件的失效模式；再利用標準運行任務循環(huán)確定部件最大損傷，并根據(jù)該值計算部件的安全壽命[5]。航空發(fā)動機設計定型后，根據(jù)剩余起降或加減速循環(huán)次數(shù)來確定部件的剩余壽命，通常不再考慮發(fā)動機實際運行環(huán)境[6-8]。使用最大損傷來計算部件安全壽命在一定程度上保證了發(fā)動機運行的安全性，但有很大局限性，即大量部件未完全破損就提前維修或更換，而部分部件在安全壽命結(jié)束前就已經(jīng)完全破損，增加了發(fā)動機運行的危險性[9-10]。發(fā)動機部件的壽命取決于其自身實際負載，如溫度、壓力等，這些因素與發(fā)動機實際運行過程密切相關，外部運行環(huán)境變化、發(fā)動機內(nèi)部噪聲、傳感器或執(zhí)行機構(gòu)等不確定性因素也會影響發(fā)動機實際運行，進而影響部件實際使用環(huán)境，均可導致部件壽命發(fā)生變化。

由于傳統(tǒng)壽命計算方法難以考慮上述不確定性因素，很難分析這些因素對壽命計算的影響。發(fā)動機傳感器測量信息中包含了運行環(huán)境及發(fā)動機內(nèi)部各種不確定性信息，因此可以采用基于傳感器數(shù)據(jù)的壽命計算模型來計算航空發(fā)動機不確定運行環(huán)境下葉片的壽命。

本文以某型軍用渦扇發(fā)動機高壓渦輪導向葉片熱機械損傷壽命為例，研究不確定運行環(huán)境下葉片壽命計算方法。

1 葉片TMF壽命模型及仿真系統(tǒng)

通常利用標準起降或加減速循環(huán)次數(shù)來計算部件的使用壽命，文獻[11]已建立某型軍用渦扇發(fā)動機壽命關鍵部件——高壓渦輪導向葉片在單次加減速循環(huán)過程中的TMF(ThermalMechanicalFatigue)壽命模型，簡單表示為

式中：Nf為葉片安全壽命；TΔmax為葉片前、后緣最大溫度差；Tmetal(TΔmax)為葉片前、后緣最大溫度差下葉片金屬溫度；Pi為發(fā)動機運行時部分截面壓力。

上述參數(shù)均可表示為發(fā)動機運行過程中溫度、壓力的函數(shù)，而溫度和壓力由傳感器測得或由機載模型提供。在無噪聲環(huán)境下，不考慮其他諸如外部運行環(huán)境、加速終值、發(fā)動機退化等因素的變化，葉片每次標稱加減速循環(huán)所經(jīng)歷的損傷相同。但在實際運行過程中，受發(fā)動機內(nèi)部噪聲、執(zhí)行機構(gòu)和傳感器測量噪聲以及運行環(huán)境變化的影響，葉片周圍的溫度和壓力值會發(fā)生變化，從而導致葉片每次循環(huán)所受的損傷不同。因此，在非標稱運行環(huán)境中，將某一固定Nf作為部件的壽命是無效的，甚至會增加發(fā)動機運行的危險性。

考慮到外部運行環(huán)境、噪聲等各種因素均會對Nf產(chǎn)生影響，因此將Nf看作隨機變量，利用概率統(tǒng)計理論來確定Nf。這里定義Nfn為葉片在標稱運行環(huán)境下的安全壽命，Nfa為非標稱運行環(huán)境下的安全壽命。那么標稱運行環(huán)境下使用n次循環(huán)后，葉片失效概率為

式中：Ffn為Nfn的累積分布函數(shù)；Pfn為葉片在非標稱運行環(huán)境下的失效概率。

研究目的是尋找Pfn（nc）和Pfa（nc）之間的關系，以確定在非標稱運行環(huán)境下，葉片失效概率及等效使用壽命。

文中構(gòu)建的仿真系統(tǒng)結(jié)構(gòu)如圖1所示。

圖1 仿真系統(tǒng)結(jié)構(gòu)

圖中，engine為某型軍用渦扇發(fā)動機模型；ECU為其對應的變增益PID控制器；LifeModel為葉片壽命計算模塊。仿真系統(tǒng)輸入包括外部環(huán)境（高度和大氣溫度）、油門桿角度；輸出為葉片在加速循環(huán)過程中的總應變差及在此應變差下葉片TMF可使用壽命。

2 葉片失效概率及等效使用壽命估算

因噪聲呈隨機分布特征，葉片壽命難以使用某一確定值來衡量，本文使用概率統(tǒng)計的方法來描述葉片使用壽命情況。工業(yè)中常用的產(chǎn)品壽命分布類型有：指數(shù)分布、正態(tài)分布、對數(shù)正態(tài)分布和威布爾分布等[12]。其中威布爾分布因其通用性強得到廣泛應用，指數(shù)分布、正態(tài)分布及瑞利分布等可看作是威布爾分布的特例。航空發(fā)動機部件壽命分布也主要使用威布爾分布[13-14]。

2.1 威布爾分布

威布爾分布是隨機變量分布之一，又稱韋伯分布、韋氏分布，由瑞典物理學家WallodiWeibull于1939年引進，是可靠性分析及壽命檢驗的理論基礎。威布爾分布在可靠性工程中被廣泛應用，尤其適用于機電類產(chǎn)品的磨損累計失效的分布。因易利用概率紙推斷出其分布參數(shù)，被廣泛應用于各種壽命試驗的數(shù)據(jù)處理。

目前，2參數(shù)的威布爾分布主要用于滾動軸承的壽命試驗以及高應力水平下的材料疲勞試驗，3參數(shù)的威布爾分布用于低應力水平的材料及某些零件的壽命試驗，具有比對數(shù)正態(tài)分布更大的適用性。其概率密度函數(shù)為

故障累積概率為

式中：β 為形狀參數(shù)，影響威布爾分布曲線的形狀，在“威布爾概率紙”上稱為威布爾斜率；η 為尺度參數(shù)，影響分布的離散程度；t0為位置參數(shù)，影響分布曲線起點的位置。

這3個參數(shù)通常由試驗確定。若t0=0，則為2參數(shù)的威布爾分布。

2.2 威布爾參數(shù)的估計

當前影響威布爾分布在壽命分布領域廣泛應用的主要因素是如何精確估計上述3個參數(shù)。資料顯示，威布爾分布參數(shù)通常由試驗確定，即通過一定次數(shù)的試驗，記錄產(chǎn)品失效的時間或者循環(huán)次數(shù)，然后利用威布爾圖來確定參數(shù)值。

令式（4)中t0=0，則式中只剩下2個參數(shù)β 和η，通過移相、取對數(shù)，最終得到

可簡寫成

這是1條不通過零點，斜率為β 的直線。通常采用雙對數(shù)坐標紙（如圖2所示）來確定威布爾參數(shù)。圖中縱坐標表示失效發(fā)生累積率，橫坐標表示失效發(fā)生的數(shù)量，這里為葉片的使用壽命。將試驗數(shù)據(jù)在圖中標出，并擬合成1條直線，則直線的斜率就是威布爾分布的參數(shù)β，而另一參數(shù)η 為直線上縱坐標為0.632所對應的橫坐標值。

圖2 威布爾圖

由于未進行實物試驗，因此設定η 為標稱運行環(huán)境下發(fā)動機標稱加速循環(huán)計算得到的葉片可使用壽命，也即η=34900，而β=4，根據(jù)“10%法則”[15]，設定安全壽命為3490次標稱循環(huán)時，葉片失效概率為0.01%。在工程應用過程中，這些參數(shù)可以使用上文方法進行估計。

在標稱運行環(huán)境下發(fā)動機運行 次加減速循環(huán)后，葉片失效概率可表示為

2.3 非標稱運行環(huán)境下葉片失效概率計算方法

為了分析非標稱運行環(huán)境下葉片的失效概率，需要統(tǒng)計出在額定運行循環(huán)nc下葉片失效次數(shù)。為了與標稱運行環(huán)境進行對比，這里選取nc=3490。計算Pfa（3490）的前提是必須知道Nfa的分布情況。發(fā)動機運行環(huán)境與葉片壽命沒有直接聯(lián)系，很難確定其分布情況。為了避免該問題，可以利用MonteCarlo仿真來解決非標稱運行環(huán)境下葉片失效概率計算問題[16-17]。

在Matlab/Simulink下，設定發(fā)動機運行環(huán)境，進行仿真并記錄葉片失效次數(shù)。當仿真次數(shù)很多（比如104）后，葉片失效比例就接近葉片真實失效概率。用表示第m 次仿真部件壽命，經(jīng)過M 次仿真后，葉片失效概率Pfa（3490）可表示為

其中ψ 定義為

很難直接利用MonteCarlo仿真來實現(xiàn)上述方法。因為Pfa（3490）的值很小，在標稱運行環(huán)境中，Pfa（3490）約為0.01%，在非標稱運行環(huán)境中，Pfa（3490）也在這個數(shù)量級上。需要仿真近10000次，葉片才會出現(xiàn)1次失效；同時也只有這1次仿真可以用于計算失效概率，大幅度降低了仿真效率。只有仿真次數(shù)足夠多時，MonteCarlo仿真的失效比例才接近于真實失效概率，仿真次數(shù)將遠多于10000次。

針對上述問題，考慮到在標稱運行環(huán)境下葉片失效概率容易求出，本文提出利用標稱運行環(huán)境作為基準值。假設Nfa與Nfn有如下線性關系

Pfa（3490）可通過下式求得

令neq=3490λ，表示非標稱運行環(huán)境中的等效使用壽命。λ值已知，即可求得Pfa（3490）。由于葉片壽命Nfa為隨機變量，因此λ 也是隨機變量，需要利用統(tǒng)計法來取得λ 的估計值。將式中的Nfn用3490代替，則

根據(jù)Nfa求出λ，利用統(tǒng)計知識可得

式中：p（nfa）為Nfa的分布情況，式中由于出現(xiàn)Nfa，因此無法直接求得，需要再次使用MonteCarlo仿真，得到

式中：λm為第m 次MonteCarlo仿真中比例系數(shù)。式（14）與式（8）的區(qū)別在于：所有的仿真均有用，可以大幅減少仿真次數(shù)，同時保證MonteCarlo仿真的準確性。具體實現(xiàn)步驟如下：

（1）確定仿真外部環(huán)境，如溫度的變化；

（2）定義仿真次數(shù)M，另m=1:M，進行M 次仿真，收集并利用式（12）計算λm；

（3）根據(jù)式（14）計算出估計λˉ；

（4）計算出等效使用壽命neq；

（5）根據(jù)式（11）計算出非標稱情況下葉片失效概率。

3 不同運行環(huán)境下葉片壽命計算

3.1 標準運行環(huán)境

這里的標準運行環(huán)境是指不考慮氣候差異所導致的發(fā)動機進口溫度變化，僅考慮由發(fā)動機內(nèi)部噪聲、材料差異所引起的葉片壽命變化。根據(jù)文獻[18]，添加溫度、轉(zhuǎn)速、壓力傳感器和執(zhí)行機構(gòu)噪聲，見表1。

表1 主要傳感器及執(zhí)行機構(gòu)噪聲

根據(jù)第2.3節(jié)中給出的方法，確定標準運行環(huán)境下與標稱運行環(huán)境之間的比例系數(shù)λˉ的值。設定運行環(huán)境ALT=0km，Tin=288.15K，各參數(shù)噪聲分布見表1，材料屬性不確定性包含在壽命模型中。仿真5000次，結(jié)果如圖3（a）、（b）所示。分別記錄第m 次仿真下渦輪導葉壽命根據(jù)式（12）計算λm，再利用式（14）計算出估計。

每次計算得到的葉片總應變差如圖3（a）所示，由于噪聲等不確定因素，壽命模型每次計算得到的葉片總應變差均不相同，但基本分布在1個固定的范圍內(nèi)（0.0063~0.00655mm）。計算得到的葉片可使用壽命分布情況如圖3（b）所示，計算得到λˉ=1.00636，在標稱運行環(huán)境下運行3490次后，葉片等效使用壽命次循環(huán)。即在標準運行環(huán)境下運行3490次循環(huán)后，相當于在標稱運行環(huán)境下運行3512次循環(huán)。

圖3 不同環(huán)境下葉片總應變差及可使用壽命分析結(jié)果

則根據(jù)式（11）計算葉片失效概率

即在標準運行環(huán)境下運行3490次循環(huán)后，葉片失效概率為0.0103%。

3.2 典型運行環(huán)境

在實際使用中，很難保證飛機每次均在標準運行環(huán)境下起飛或加速，運行高度和大氣溫度均會變化，需要研究這些差異對壽命計算的影響。設定典型運行環(huán)境為ALT=0~1km，Tin=288.15±10K，即高度和溫度在上述范圍內(nèi)隨機變化，噪聲幅值與標準運行環(huán)境下相同，其他設定同上節(jié)，其仿真結(jié)果如圖3（c）、（d）所示。

在典型運行環(huán)境下，受運行高度及大氣溫度的變化影響， 葉片上的總應變在更大的范圍（0.0062~0.0066mm）內(nèi)變化，計算得到=1.0326。在典型運行環(huán)境下運行3490次后，葉片等效使用壽命為3490×=3604次循環(huán)。即在典型運行環(huán)境下運行3490次循環(huán)后，相當于在標稱運行環(huán)境下運行3604次循環(huán)。

則根據(jù)式（11）得

即在典型運行環(huán)境下運行3490次循環(huán)后，葉片失效概率為0.0114%。

3.3 高溫運行環(huán)境

考慮到部分實際使用的飛機主要運行航線分布在近赤道附近，大部分起飛或加速時外界大氣溫度比典型運行環(huán)境下的高，為研究這種高溫運行環(huán)境對葉片壽命的影響，設定高溫運行環(huán)境ALT=0~1km，Tin=268.15±10K，噪聲幅值同標準運行環(huán)境下的，其仿真結(jié)果如圖3（e）、（f）所示。

在高溫運行環(huán)境下，由于大氣溫度升高，葉片上的平均總應變比典型運行環(huán)境下的大，主要分布于0.0064~0.0070mm，此時計算得到=1.5694，這樣當在高溫運行環(huán)境下運行3490次后，葉片等效使用壽命為3490×=5477次循環(huán)。即在高溫運行環(huán)境下運行3490次循環(huán)后，相當于在標稱運行環(huán)境下運行3604次循環(huán)。

則根據(jù)式（11）

即在高溫運行環(huán)境下運行3490次循環(huán)后，葉片失效概率為0.0607%。

3.4 低溫運行環(huán)境

部分實際使用的飛機主要運行航線分布在近極地附近，其大部分起飛或加速時外界大氣溫度將比典型運行環(huán)境的低，本節(jié)中設定運行環(huán)境ALT=0~1 km，Tin=268.15±10K，噪聲幅值同標準運行環(huán)境下的，研究這種低溫運行環(huán)境對葉片壽命影響的仿真結(jié)果如圖3（g）、（h）所示。

在低溫運行環(huán)境下，由于大氣溫度下降，葉片上的平均總應變要比典型運行環(huán)境下的小，主要分布于0.0059~0.0065mm，此時計算得到λˉ=2599，當在低溫運行環(huán)境下運行3490次后，葉片等效使用壽命為次循環(huán)。即在低溫運行環(huán)境下運行3490次循環(huán)后，相當于在標稱運行環(huán)境下運行2599次循環(huán)。

則根據(jù)式（11）得

即在低溫運行環(huán)境下運行3490次循環(huán)后，葉片失效概率為0.0031%。

不同運行環(huán)境下葉片失效概率及3490次循環(huán)等效使用壽命見表2。

表2 不同環(huán)境下葉片壽命及失效概率

從表2中可見，在標稱運行環(huán)境下，不考慮噪聲影響，當運行循環(huán)數(shù)為3490次時，葉片失效概率為0.01%；而在考慮噪聲等不確定因素影響的標準運行環(huán)境中，運行同樣循環(huán)數(shù)，葉片失效概率增大至0.0103%；在典型運行環(huán)境中，由于運行高度及大氣溫度的變化，葉片等效使用壽命增加至3604次循環(huán)，此時運行3490次循環(huán)后葉片失效概率為0.0114%；而在高溫運行環(huán)境下，等效使用壽命高達5477次循環(huán)，葉片失效概率為0.0607%，但在冷天環(huán)境下，運行3490次循環(huán)僅相當于在標稱運行環(huán)境下運行2559次循環(huán)，因此其3490次循環(huán)葉片失效概率僅為0.0031%。

4 結(jié)論

本文利用MonteCarlo仿真和統(tǒng)計學理論中的威布爾分布分析了部件在各種運行環(huán)境中的壽命及失效概率，為航空發(fā)動機壽命延長控制研究提供了包含運行環(huán)境信息的部件壽命計算模型，從而在全壽命期實現(xiàn)壽命延長控制。主要結(jié)論如下：

（1）仿真結(jié)果顯示，運行環(huán)境不確定使葉片計算壽命有大幅度變化，單純將某一固定值作為發(fā)動機部件的安全壽命是不合理的。

（2）計算結(jié)果表明，不考慮和考慮噪聲等不確定因素下，渦輪葉片失效概率相差3%；高、低溫環(huán)境下，對渦輪葉片壽命及失效概率都有很大影響。

該方法可擴展到實際葉片負載和加速循環(huán)不確定性影響研究中，為發(fā)動機壽命監(jiān)視提供更加準確的部件使用壽命估計值。同時也可以根據(jù)該估計值來調(diào)整壽命延長控制策略，盡可能延長部件壽命。

[1]陳小磊，郭迎清，杜憲.航空發(fā)動機全壽命期自適應壽命延長控制[J].推進技術，2013，31（1）：107-114.CHEN Xiaolei，GUO Yingqing，DU Xian.Adaptive life extending control of aircraft engine in the whole life[J].Journal of Propulsion Technology，2013，31（1）：107-114（in Chinese）.

[2]陳小磊，郭迎清，張書剛.航空發(fā)動機壽命延長控制綜述[J].航空發(fā)動機，2013，39（1）：17-22.CHEN Xiaolei，GUO Yingqing，ZHANG Shugang.Summary of life extending control for an aeroengine[J].Aeroengine，2013，39（1）：17-22.（in Chinese）

[3]Vittal S,Hajela P,Joshi A.Review of approaches to gas turbine life management[R].AIAA-2004-4372.

[4]劉廷毅.航空發(fā)動機研制全壽命管理研究及建議[J].航空發(fā)動機，2012，38（1）：1-6.LIU Tingyi.Research and suggestion of life cycle management for aeroengine development[J].Aeroengine，2012，38（1）：1-6.（in Chinese）

[5]蘇清友，孔瑞蓮，陳筱雄.航空渦噴、渦扇發(fā)動機主要零部件定壽指南[M].北京：航空工業(yè)出版社，2004：89-97.SU Qingyou,KONG Ruilian,CHEN Xiaoxiong.The life guide of the main parts of aviation turbojet and turbofan engine[M].Beijing：Aviation Industry Press，2004：89-97.（in Chinese）

[6]Suarez E L，Duffy M J，Seto D，et al.Advanced life prediction systems for gas turbine engines[R].AIAA-2003-4985.

[7]Zhang J.Fatigue life prediction under random loading using distributional stress-life relationship[R].AIAA-99-1600.

[8]Hsiung H C，Dunn A J，Loh D L.Stress analysis and life prediction of gas turbine blade[R].AIAA-88-3218.

[9]蔡向暉，東巖，李偉.航空發(fā)動機單機技術狀態(tài)壽命控制及應用[J].航空發(fā)動機，2009，35（4）：16-19.CAI Xianghui，DONG Yan，LI Wei.Life control and application of aeroengine for single engine technical condition [J].Aeroengine，2009，35（4）：16-19.（in Chinese）

[10]Grelotti R A，Glanovsky J L.Usage-based life prediction and fleet anagement for gas turbine engines[R].AIAA-2010-2972.

[11]陳小磊，郭迎清，陸軍.基于修改加速控制規(guī)律的航空發(fā)動機壽命延長控制[J].航空動力學報，2011，26（9）：2116-2121.CHEN Xiaolei，GUO Yingqing，LU Jun.Life extending control of aircraft engine based on adjusting acceleration schedule[J].Journal of Aerospace Power，2011,26（9）：2116-2121.（in Chinese）

[12]凌丹.威布爾分布模型及其在機械可靠性中的應用研究[D].成都：電子科技大學,2011.LING Dan.The application of Weibull distribution model in mechanical reliability[D].Chengdu：University of Electronic Science and Technology,2011.（in Chinese）

[13]薛向珍，李育錫，王三民.某直升機主減速器傳動系統(tǒng)的壽命與可靠性計算方法[J].航空動力學報，2011，26（3）：635-641.XUE Xiangzhen，LI Yuxi，WANG Sanmin.Method of lifetime and reliability of some helicopter's main reducer[J].Journal of Aerospace Power，2011，26（3）：635-641.（in Chinese）

[14]余國林.小子樣飛機系統(tǒng)使用可靠性評估方法研究與應用[D].南京：南京航空航天大學,2005.YU Guolin.Study and application of the evaluation technology on the small sample aircraft system operational reliability[D].Nanjing：Nanjing University of Aeronautics and Astronautics，2005.（in Chinese）

[15]航空發(fā)動機設計手冊總編委會.航空發(fā)動機設計手冊（第18冊）：葉片輪盤及主軸強度分析[M].北京:航空工業(yè)出版社,2001：205-235.Editorial Board of Aeroengine Design Manual.Aeroengine design manual（section 18）：strength analysis of vane wheel and spindle[M].Beijing：Aviation Industry Press，2001：205-235.（in Chinese）

[16]Liu J S.Monte Carlo strategies in scientific computing[M].New York:Springer,2008：53-57.

[17]肖剛，李天柁，余梅.動態(tài)系統(tǒng)可靠性仿真的五種蒙特卡羅方法[J].計算物理，2001,18（2）:173-176.XIAO Gang，LI Tianduo，YU Mei.Five Monte Carlo methods for simulation of dynamic reliability system [J].Chinese Journal of Computation Physics，2001,18（2）:173-176.（in Chinese）

[18]Behbahani A，Semega K.Sensing challenges for controls and PHM in the hostile operating conditions of modern turbine Engine[R].AIAA-2008-5280.