航空發(fā)動機渦輪部件蠕變失效風險分析

李昊燃，陳 健

（中國航發(fā)商用航空發(fā)動機有限責任公司 設計研發(fā)中心，上海 200241）

1 引言

典型航空發(fā)動機部件處于復雜的高溫、高載荷耦合環(huán)境中，蠕變是引起其原發(fā)失效的機制之一。雖然傳統(tǒng)部件的設計按照確定性設計原則進行，然而對某個特定部件而言，無論載荷、幾何還是材料性能和缺陷，均存在一定的分散性，由此導致該部件在預期壽命周期內存在一定的失效風險。CCAR33.75安全分析條款[1]中，要求申請人必須對發(fā)動機進行安全性分析，以評估預期可能發(fā)生的所有失效的后果。因此，有必要開展航空發(fā)動機的蠕變風險研究。雖然確定性的設計方法和準則可以通過確保發(fā)動機關鍵部件在壽命周期內的結構完整性，實現蠕變風險控制，但近年來，隨著概率計算方法的進步，概率風險計算成為發(fā)展的一種趨勢[2-3]，概率設計以及基于風險進行管理逐漸獲得應用[4-6]，概率風險計算可以為航空發(fā)動機產品的設計制造以及使用管理提供指導性意見。

對于用概率方法進行蠕變風險計算，文獻[7]在1997 年提出曲面擬合結合積分降階法解決渦輪盤蠕變可靠性問題。文獻[8]在此基礎上進行改進，采用有限元計算方法結合Monte Carlo 法以及人工神經網絡方法進行渦輪盤蠕變可靠性分析。文獻[9]主要考慮渦輪葉片壽命消耗主要是蠕變和低周循環(huán)疲勞，對短壽命發(fā)動機高壓渦輪葉片使用壽命進行預測。文獻[10]等建立了考慮多軸應力影響的含缺陷高溫構件蠕變壽命預測方法。文獻[11-12]等對鈦合金、粉末高溫合金的保載疲勞開展了分析和實驗研究，揭示了蠕變與低周的交互作用。文獻[13]提出了考慮材料蠕變應變模型參數分散性及應力松弛效應的渦輪盤蠕變-疲勞壽命可靠性分析方法。

文獻[14]用有限元分析以及Monte Carlo法對鎳基單晶葉片可靠性進行研究。文獻[15]基于中國航空材料手冊試驗數據以及蠕變持久壽命預測理論建立蠕變壽命預測模型，并運用異常數據處理、參數估計和擬合優(yōu)度檢驗等方法確定數據分布類型，計算任意可靠度下的持久蠕變壽命。文獻[16]將應力和應變能密度作為隨機變量，利用試驗結果構建了鋼材的蠕變-疲勞壽命可靠性模型。文獻[17]采用響應面法結合Monte Carlo法對考慮蠕變損傷的渦輪盤低周疲勞壽命進行可靠性評估。

綜上所述，國內外學者在蠕變壽命模型、蠕變-疲勞壽命交互模型以及壽命可靠性方面開展了系列研究，但考慮應力松弛的蠕變可靠性分析，特別是針對航空發(fā)動機零部件的實際應用分析較為鮮見。結合航空發(fā)動機典型零部件服役載荷特點，對適用于航空發(fā)動機的蠕變風險分析方法進行研究。在有限元分析的基礎上，采用響應面法擬合蠕變壽命函數，通過Monte Carlo法進行蠕變壽命可靠度計算及蠕變風險分析。并以某型民用航空發(fā)動機渦輪葉片為例，開發(fā)蠕變風險計算程序，并開展算例分析。

2 蠕變風險分析方法

主要介紹蠕變風險計算的流程，以轉速ω作為輸入的隨機變量，對于危險點的蠕變壽命可靠度進行計算。進行彈塑性蠕變軸對稱問題分析可以采用有限差分法，也可以采用有限元法。兩種方法計算結果相差不大，實際應用中通常使用ANSYS進行有限元計算。

首先，用有限差分法或ANSYS有限元計算出隨機變量ω對應的一組應力值σ，接著將計算得到的應力值σ和給定溫度T帶入熱強參數綜合方程計算出一組蠕變壽命tc，進而求得一組對應的蠕變損傷dc，計算壽命周期內累積損傷Dc，給定蠕變臨界損傷DCR，構建功能函數：

用隨機法進行模擬抽樣，計算給定壽命的可靠度。

2.1 蠕變持久壽命概率模型

以損傷作為蠕變持久壽命考核量進行蠕變持久壽命預測及可靠性分析。目前材料手冊中[18]常用的蠕變壽命預測方程主要有蠕變持久方程和熱強參數綜合方程兩種，美國渦噴渦扇發(fā)動機設計通用規(guī)范、我國國軍標和發(fā)動機設計規(guī)范都推薦采用拉森-米勒（L-M）方程、葛-唐吾（G-D）方程、曼森-索柯普（M-S）方程和曼森-哈弗特（M-H）方程這四種持久方程，通常選擇其中短期蠕變試驗數據擬合最好的方程最為最終的蠕變持久方程。選用熱強參數綜合方程進行蠕變壽命概率分析，考慮到材料蠕變性能存在一定的分散性，在熱強參數綜合方程中添加隨機參數進行可靠性分析。構建如下概率模型：

式中：σ—持久強度；Ci—待定系數；P—熱強參數；η—材料蠕變應力隨機參數。

上述四種持久方程相對四個的熱強參數方程分別為：

式中：C5、T0—材料常數。

將隨機變量ω對應的隨機應力值σ及溫度T帶入熱強參數綜合方程可計算出蠕變持久壽命tc。

2.2 蠕變損傷計算

蠕變損傷計算采用Robinson法則，將應力和溫度按時間步長ΔT進行細分，劃分為n個子區(qū)間。以第i個子區(qū)間為例，在該時間段內產生的蠕變損傷定義為：

式中：ΔT—時間步長；Di—該子區(qū)間內蠕變損傷增量；tci—熱強參數綜合方程計算出的每個區(qū)間的單點蠕變壽命。

則整個循環(huán)內蠕變累積損傷為：

式中：n—總的時間步數，子區(qū)間內取應力和溫度的最大值計算蠕變壽命的應力和溫度。通過上述方法可計算全壽命周期內蠕變累積損傷。

2.3 蠕變持久壽命可靠度計算

根據Miner累積損傷理論，對航空發(fā)動機典型部件最危險點進行蠕變持久壽命可靠性分析。進行輪盤蠕變持久壽命可靠性分析時，采用瞬時蠕變累積損傷近似計算方法，可以求出應力σ以及蠕變持久壽命預測模型隨機參數η相對應的蠕變累積損傷dc，進而計算出蠕變壽命周期N次飛行循環(huán)內的累積損傷Dc。因此理論上可對各隨機變量進行大量模擬抽樣，計算各隨機變量組下的蠕變累積損傷，對蠕變累積損傷進行統(tǒng)計分析可以得出航空發(fā)動機典型部件關鍵點的可靠度。

由于上述蠕變累積損傷計算過程涉及不同時刻點的大量有限元計算，工程中難以實現，故在進行給定壽命的航空發(fā)動機典型部件蠕變持久壽命可靠性評估時，通過對各關鍵點進行有限元分析，得到應力及溫度隨時間變化的歷程，然后進行蠕變損傷計算，進而采用Monte-Carlo法獲得關鍵點蠕變持久壽命可靠度，重復該方法完成整個航空發(fā)動機典型部件的蠕變風險計算。具體主要計算步驟如下：

（1）在有限元計算得到的應力數組的基礎上，對參數η進行隨機化，計算單個位置點蠕變累積損傷dc。

（2）進行有限次抽樣，計算相應樣本點對應的蠕變累積損傷值Dc。

（3）給定蠕變臨界損傷DCR，根據強度干涉理論，可以得出蠕變損傷臨界失效函數，如式（1）所示。

（4）采用Monte-Carlo法對對應力和隨機參數η的隨機抽樣模擬抽樣，可得蠕變損傷臨界失效函數Z的分布，Z≥1時認為結構失效，進而可確定給定壽命的可靠度。

3 蠕變風險計算流程

基于蠕變風險計算的理論和流程，開發(fā)蠕變風險計算工具，利用隨機抽樣法，對關鍵位置點應力和熱強參數綜合方程中隨機參數η進行Ns次抽樣。

每一次抽樣均以計算位置點應力、溫度隨時間變化的歷程為輸入，在經歷服役要求的循環(huán)數后，計算該點的累積蠕變損傷。

在完成Ns次抽樣模擬后，統(tǒng)計其中的失效樣本，計算得到該點在蠕變條件下的失效概率，即完成給定循環(huán)數下計算位置點可靠度計算。重復該方法對不同壽命周期（循環(huán)數）內關鍵位置點的可靠度進行計算，完成該關鍵點的風險計算。若進行部件定量風險分析，對關鍵位置點風險計算流程及算例分析進行詳細說明。若要進行部件蠕變風險分析，則需通過部件各關鍵位置點的風險計算為部件的定量風險分析提供依據。

圖1 蠕變風險計算流程Fig.1 Process of Creep Risk Calculation

在部件的各個區(qū)域指定不同數量的樣本，Monte-Carlo法可以估計樣本的失效概率。每個循環(huán)中的一組應力為隨機變量σ，在符合蠕變規(guī)律的基礎上隨機生成。考慮材料內部缺陷分布不同，為了不失一般性，在傳統(tǒng)蠕變壽命預測模型的基礎上，添加隨機參數η，進行可靠性分析。具體計算步驟如下：

（1）根據有限元分析結果，確定單個循環(huán)內時間數組t、應力數組σ、溫度數組T。

（2）查閱材料性手冊，確定熱強參數方程的系數C1～C5、熱強參數限制值PMax、PMin，以及溫度限制值TMax、TMin。

（3）定義需要計算不同的壽命周期數組N，其中任意一個壽命用N(h)表示。

（4）定義Monte-Carlo模擬次數Ns，其中任意一次抽樣用k表示。初始失效樣本數Nf=0。

（5）變量隨機化。

隨機抽取初始應力歷程中最大應力點σrand，基于蠕變規(guī)律可按比例得到本次循環(huán)中其他時間點對應的應力值。生成的這組隨機應力記為σrand。

熱強參數綜合方程中添加材料蠕變應力隨機參數η，通常可假設其服從正態(tài)分布[19]。

（6）確定計算點在給定的壽命周期N(h)內是否失效。

調用蠕變損傷計算子程序。輸入參數為：蠕變時間點數組t，應力數組σrand，溫度數組T。

返回參數為：單次飛行循環(huán)的蠕變損傷dc。

計算壽命周期內N(h)蠕變累積損傷Dc。

如果Dc<1，則不發(fā)生失效；否則發(fā)生失效，Nf=Nf+1，直到k=Ns。

（7）重復（2）～（6）步，直到完成Ns次抽樣。

計算得到該關鍵位置點經歷N(h)后的失效概率Pf_h=Nf/Ns。

（8）計算不同壽命周期N對應的失效概率，得到關鍵點在不同壽命周期內的失效概率。

4 蠕變風險算例分析

渦輪葉片作為發(fā)動機的重要零部件，處在高溫、高轉速、高振環(huán)境中，蠕變是其主要失效形式之一。在進行蠕變壽命計算時，包括尺寸、裝配容差、機械及溫度載荷、材料力學性能、材料蠕變性能等輸入數據的不確定性，往往會會對蠕變壽命造成影響，單用一條中值曲線來描述蠕變的整個過程顯然是不科學的，因此需要開展風險分析。

以DD6材料的渦輪葉片為例進行蠕變風險計算。DD6是我國第二代鎳基單晶高溫合金，具有高溫強度高、綜合性能好、組織穩(wěn)定及鑄造工藝性能好的優(yōu)點。該合金適合于制作1100℃以下工作的具有復雜內腔的燃氣渦輪工作葉片等高溫零件。以DD6渦輪葉片為例，通過自主編寫的程序進行蠕變風險計算，為葉片熱機耦合環(huán)境下相對失效風險定量分析工作提供基礎。

利用Monte-Carlo 法實現蠕變風險計算，通過大量抽樣，計算給定的壽命周期下（循環(huán)數）關鍵點的失效概率，對渦輪葉片的關鍵點進行可靠性分析。在此基礎上，計算不同壽命周期（循環(huán)數）下關鍵點的失效概率，為渦輪葉片的風險定量分析提供依據。通過自主開發(fā)的程序進行蠕變風險計算，具體計算過程如下：

4.1 有限元模型

渦輪葉片有限元模型圖，如圖2所示。葉片施加離心力和溫度載荷，約束葉片榫頭一側節(jié)點軸向位移，以模擬擋環(huán)的軸向約束，在榫頭接觸面上施加等效接觸壓力。提取葉片關鍵位置點，如圖2所示，考慮關鍵位置點應力松弛后應力歷程進行蠕變壽命分析。

圖2 渦輪葉片有限元模型Fig.2 Finite Element Model of the Turbine Blade

4.2 載荷譜簡化

通過ANSYS對渦輪葉片模型進行有限元計算，得到一個飛行循環(huán)內渦輪葉片關鍵位置點的等效應力和溫度譜。典型民用航空發(fā)動機載荷譜包括地面慢車、起飛、爬升、巡航、空中慢車、進近、反推等若干飛行段，參考當前文獻通常將其簡化[20]，考慮各飛行段對于蠕變的貢獻，簡化一個飛行循環(huán)內應力及對應的溫度隨時間變化。選擇、一些關鍵的時間點描述應力、溫度隨時間變化的歷程，本算例應力取Mises應力。

考慮到應力松弛對蠕變壽命影響較大，因此，參考文獻[16]，考慮有限元循環(huán)加載前51個循環(huán)的應力松弛，對應力進行修正，并對應力和溫度數據進行歸一化處理。修正后的應力、溫度隨時間變化歷程，如圖3所示。蠕變風險計算輸入數據，如表1所示。

圖3 應力/溫度譜Fig.3 Stress/Temperature Spectrum

表1 渦輪葉片關鍵位置點蠕變風險計算數據Tab.1 Creep Risk Calculation Data for the Key Position Point of Turbine Blade

為了確保計算結果的準確性，這里蠕變風險計算程序在求解時設置了參數限制，如表2所示。計算應力、溫度對應的壽命時，如果應力計算得到的熱強參數P的最大值超過PMax或者溫度超過TMax，則不容許外插，壽命設為0；如果應力計算得到的熱強參數P小于PMin或者溫度小于TMin，直接用PMin、TMin代替當前應力或者溫度值。

表2 參數限制Tab.2 Parameter Limits

4.3 材料蠕變性能

查閱《航空發(fā)動機設計用材料數據手冊（第四冊）》[21]，通過試驗數據擬合獲得熱強參數方程中的系數，如表3所示。

表3 熱強參數方程相關參數Tab.3 Related Parameters of the Comprehensive Equation of Thermal Strength Parameters

4.4 給定壽命周期內可靠度計算

確定關鍵點的輸入數據、相關參數、限制值后，程序進行不同壽命周期（循環(huán)數）內可靠度的計算。根據強度干涉理論，蠕變臨界損傷DCR=1，蠕變損傷臨界失效函數如式（1）所示，當Z≥1時失效。

由于不同轉速下關鍵點處對應不同的應力值，因此應力σ為隨機變量，計算時應對σ進行隨機抽樣。由于本算例缺乏足夠的真實應力數據，因此，假設應力σ服從正態(tài)分布，以單次應力σ為中值，以5MPa為標準差，采用Monte-Carlo法對應力σ進行抽樣。參考文獻[19]通常假設熱強參數綜合方程中隨機參數η服從正態(tài)分布N(0,0.014)，用Monte-Carlo法對隨機參數η進行抽樣。

進行10000次抽樣后，統(tǒng)計得出Z≥1的概率，從而獲得給定壽命周期（循環(huán)數）內蠕變損傷臨界失效函數Z的概率密度函數，并計算得到給定壽命周期（循環(huán)數）關鍵位置點的可靠度。

4.5 可靠度風險計算結果

程序對關鍵點在不同的壽命周期內進行了蠕變風險計算，得到不同壽命周期對應的蠕變失效概率，如表4所示。蠕變失效概率與壽命周期關系，如圖4所示。

圖4 壽命周期與失效概率關系圖Fig.4 Relationship Between Life Cycle and Creep Failure Probability

表4 關鍵位置點蠕變風險計算結果Tab.4 Results of Creep Risk Calculation at the Key Position Point

由于本算例僅考慮有限元循環(huán)加載前51 個循環(huán)的應力松弛，計算結果偏保守，如果考慮更多循環(huán)的應力松弛結果將更接近實際情況，后續(xù)將進一步開展該方面的研究。

5 結論

（1）傳統(tǒng)部件的設計按照確定性設計原則進行，然而實際情形中，無論載荷、幾何還是材料性能和缺陷，均存在一定的分散性，導致該部件在預期壽命周期內存在一定的失效風險，因此在設計中引入概率方法開展蠕變風險分析是十分必要的。

（2）發(fā)動機典型零部件如渦輪葉片，蠕變在總損傷中占比較大，而在不同壽命周期內由于蠕變的影響發(fā)生失效的概率不同，基于風險的研發(fā)和壽命管理是一種發(fā)展趨勢，量化的蠕變風險分析技術可以有效指導航空發(fā)動機零部件的設計、制造和使用管理等工作。