雷世英,孫見忠,劉赫
南京航空航天大學(xué) 民航學(xué)院,南京 211106
民航發(fā)動機熱端部件的故障和失效是造成發(fā)動機性能退化及非計劃換發(fā)的一個重要因素,其很大程度上制約了發(fā)動機整機的使用壽命。按照失效造成的影響嚴(yán)重性,可將熱端部件分為A、B兩大類。A類件又稱為限壽件(Life Limited Part, LLP),包括渦輪軸、渦輪盤等主要旋轉(zhuǎn)件,其失效后會造成機毀人亡的嚴(yán)重后果,工程上一般通過安全性分析給定安全壽命,當(dāng)達(dá)到規(guī)定的使用循環(huán)后進行強制報廢,這種方法由于沒有考慮個體發(fā)動機實際的載荷差異,壽命估計一般較為保守;B類件又稱為關(guān)鍵件,主要是對發(fā)動機性能影響較大的結(jié)構(gòu)件,如壓氣機/渦輪葉片、軸承等,其失效的后果沒有A類件嚴(yán)重,但維修成本較高,發(fā)生故障時還可能造成較大的二次破壞,工程上一般采取視情維修的策略,通過評估剩余壽命來制訂相應(yīng)的維修計劃,因此其壽命預(yù)測的準(zhǔn)確與否對發(fā)動機的安全性和經(jīng)濟性有著重要影響,也是發(fā)動機壽命預(yù)測領(lǐng)域廣泛關(guān)注的問題。如今發(fā)動機制造商和運營商對安全性和經(jīng)濟性要求越來越高,越來越趨向于在確保安全的前提下,盡可能地提高發(fā)動機的可用性,最大限度降低壽命周期運行成本,實現(xiàn)安全性和經(jīng)濟性的最大化。然而實際中要想實現(xiàn)壽命的準(zhǔn)確預(yù)測是相當(dāng)困難的,熱端部件的使用壽命不僅受設(shè)計制造、材料工藝水平的影響,還與發(fā)動機的實際運行環(huán)境、使用方式以及維修模式等有著密切關(guān)聯(lián),這些多方面的因素使得預(yù)測的不確定性大大增加,預(yù)測精度很難滿足實際需要。
目前對于關(guān)鍵件服役壽命的預(yù)測,總體方法上可以分為統(tǒng)計模型法、物理模型法和多模型、多數(shù)據(jù)融合的方法。在不具備零部件設(shè)計信息的情況下,一般使用可靠性分析的方法評估葉片的使用可靠性及剩余壽命,如An等采用貝葉斯方法融合外場可靠性數(shù)據(jù),得到更新后的壽命分布來確定葉片壽命;Zaretsky等根據(jù)發(fā)動機外場使用數(shù)據(jù)分析了葉片的使用壽命,最終推導(dǎo)出一個估算葉片剩余壽命的計算公式;Li等采用模糊理論來表示渦輪葉片壽命預(yù)測中涉及的不確定性,利用線性融合算法對渦輪葉片的可靠性壽命進行區(qū)間評估。上述統(tǒng)計模型法需要提供大量的輸入數(shù)據(jù),但搜集數(shù)據(jù)的時間成本較高,具體的實驗數(shù)據(jù)和現(xiàn)場數(shù)據(jù)也較為缺乏;除此之外,壽命可靠性預(yù)測涉及小樣本問題,得到的結(jié)果反映的是相似使用條件下零部件使用可靠性的平均屬性,難以體現(xiàn)單機狀態(tài)下發(fā)動機自身及運行環(huán)境、運行載荷等方面的差異性。物理模型法則是建立在掌握了詳細(xì)的材料參數(shù)與部件設(shè)計參數(shù)的基礎(chǔ)之上,通過微觀失效機理研究、計算流體力學(xué)與結(jié)構(gòu)有限元法等分析手段,得到關(guān)鍵部位的應(yīng)力、溫度載荷,進而憑借失效機理對應(yīng)的壽命計算模型,對壽命損耗進行評估,如Staroselsky等建立了基于晶體形態(tài)的彈粘塑性模型,考慮了材料滑移層的本構(gòu)行為,并將其耦合到有限元模型中,使用半經(jīng)驗壽命預(yù)測模型估算葉片的剩余壽命;荷蘭NLR中心和德國DLR中心等分別針對發(fā)動機渦輪葉片,通過實驗設(shè)計(Design of Experience, DOE)方法構(gòu)建了葉片溫度和應(yīng)力計算的代理模型,并考慮了蠕變和熱機疲勞兩種損傷模式來預(yù)測葉片剩余壽命;Zhu等通過引入種載荷參數(shù)來考慮渦輪葉片高低周疲勞損傷的耦合效應(yīng),使得預(yù)測誤差的均值和標(biāo)準(zhǔn)差都達(dá)到了較高的精度;Brand?o等等通過逆向工程得到了渦輪葉片的材料成分,借助有限元仿真比較了等溫模型和非等溫模型的蠕變壽命預(yù)測結(jié)果。上述物理模型的預(yù)測方法依托于材料力學(xué)、流體力學(xué)、熱力學(xué)、有限單元法等基本理論,利用發(fā)動機的實際運行參數(shù),獲得葉片關(guān)鍵部位的應(yīng)力和溫度載荷譜,進而借助壽命損耗模型來評估渦輪葉片的剩余壽命。此類方法需要掌握發(fā)動機詳細(xì)的設(shè)計數(shù)據(jù),葉片的應(yīng)力與溫度載荷分析需要進行復(fù)雜的計算流體力學(xué)與結(jié)構(gòu)有限元分析,但在模型計算時做了大量假設(shè),會造成較大的預(yù)測不確定性。而多模型、多數(shù)據(jù)融合的方法則可以充分利用已有的多源異構(gòu)數(shù)據(jù),在葉片失效機理的基礎(chǔ)上綜合考量使用條件、運行環(huán)境等差異,相較于單模型方法能夠提供更高的預(yù)測精度。Pillai等提出了基于有限元代理模型的多源信息融合方法,利用機器學(xué)習(xí)方法融合外場失效等數(shù)據(jù),進一步降低預(yù)測結(jié)果的不確定性;Giesecke等結(jié)合發(fā)動機使用參數(shù)及運行區(qū)域大氣環(huán)境條件建立貝葉斯置信網(wǎng)絡(luò)用來預(yù)測葉片壽命損耗,以輔助決策葉片修理級別;Fu等考慮了微觀結(jié)構(gòu)準(zhǔn)則和蠕變應(yīng)變準(zhǔn)則,并結(jié)合人工神經(jīng)網(wǎng)絡(luò)和改進的投影模型來評估使用中的渦輪葉片的剩余蠕變壽命,開發(fā)了一種基于材料基因工程的渦輪葉片在役蠕變剩余壽命的集成計算方法。
實際上,反映發(fā)動機關(guān)鍵件服役可靠性降級等行為特性的信息存在于多模態(tài)數(shù)據(jù)中,包括產(chǎn)品設(shè)計數(shù)據(jù)與使用維護數(shù)據(jù),如產(chǎn)品幾何參數(shù)、材料數(shù)據(jù)、狀態(tài)監(jiān)測數(shù)據(jù)、離線檢測/檢查報告、故障記錄和大修報告等,貫穿產(chǎn)品壽命周期的不同階段,如何充分利用這些多模態(tài)數(shù)據(jù),實現(xiàn)多領(lǐng)域物理模型與數(shù)據(jù)融合,是提高服役發(fā)動機關(guān)鍵件壽命預(yù)測分析準(zhǔn)確性的重要途徑。本文以民航發(fā)動機高壓渦輪(HPT)葉片蠕變失效為例,在蠕變失效機理分析基礎(chǔ)上,同時考慮服役條件下的使用信息、狀態(tài)參數(shù)及截尾失效時間等數(shù)據(jù),提出的蠕變累積損傷指數(shù)模型實現(xiàn)了多模態(tài)信息的融合,對渦輪葉片服役可靠性進行評估,并預(yù)測其剩余壽命?,F(xiàn)階段外場數(shù)據(jù)獲取有限,因此論文通過仿真數(shù)據(jù)驗證提出方法的可行性和有效性,結(jié)果表明,本文方法比傳統(tǒng)可靠性分析方法預(yù)測精度有顯著提高。
文章的主要內(nèi)容安排如下:第1部分主要介紹渦輪葉片的蠕變失效機理,闡述了蠕變的發(fā)展階段、微觀失效原理以及影響蠕變損傷的常見因素,并引入蠕變損傷評估模型;第2部分主要介紹葉片蠕變壽命計算的整個流程,包括發(fā)動機性能仿真、葉片應(yīng)力計算、葉片溫度計算及蠕變損傷評估等過程;第3部分介紹葉片蠕變累積損傷指數(shù)模型,給出了損傷指數(shù)模型的定義方法以及使用該模型進行剩余壽命預(yù)測的流程;第4部分為模型的仿真驗證部分,使用實際發(fā)動機的狀態(tài)監(jiān)測數(shù)據(jù)進行了全壽命仿真,基于仿真數(shù)據(jù)驗證蠕變累積損傷指數(shù)模型,預(yù)測HPT葉片的可靠性壽命,并與基準(zhǔn)模型進行了結(jié)果對比分析;第5部分是對全文的總結(jié)以及后續(xù)工作的展望。
渦輪葉片作為典型的熱端部件,長期工作在高溫、高壓、高轉(zhuǎn)速的極限環(huán)境下,承受著離心載荷、氣動載荷、溫度載荷、振動載荷及高溫氧化腐蝕的綜合作用,其主要的失效模式包括:熱機疲勞在內(nèi)的低循環(huán)疲勞、振動引起的高循環(huán)疲勞、高溫長時間載荷作用下的蠕變變形和蠕變應(yīng)力斷裂;高溫燃?xì)鉀_刷腐蝕和氧化以及外物損傷等。對于現(xiàn)代高涵道比民航發(fā)動機而言,其飛行過程平穩(wěn)、飛行任務(wù)較為單一、飛行時間較長,這種工況下HPT葉片的蠕變損傷是決定其使用壽命的主要因素之一。同時,對于蠕變壽命的數(shù)值方法研究已經(jīng)有了較為成熟的模型和方法,在工程實踐中得到了廣泛的應(yīng)用。因此以葉片的典型失效模式—蠕變失效為例,來研究累積損傷指數(shù)建模以及外場數(shù)據(jù)融合的方法。
蠕變是指工作在高溫、高應(yīng)力水平條件下的部件,隨著加載時間的延長緩慢地發(fā)生塑性變形的現(xiàn)象。一般情況下,當(dāng)材料的工作溫度大于或等于其0.4~0.6倍熔點溫度,且工作應(yīng)力遠(yuǎn)低于材料的屈服強度時就會產(chǎn)生較為顯著的蠕變變形。如圖1所示,以蠕變時間為橫軸,蠕變導(dǎo)致的材料變形量為縱軸,可以發(fā)現(xiàn)蠕變的發(fā)展一般分為3個階段:① 初始蠕變階段(→),材料已發(fā)生一定形變,初始應(yīng)變?yōu)椋渥兯俾孰S時間逐漸降低;② 二次蠕變階段(→),蠕變速率基本保持恒定;③ 三次蠕變階段(→),蠕變速率隨時間逐漸增加,直到蠕變時間發(fā)生蠕變斷裂,此時變形量為。
圖1 蠕變發(fā)展階段Fig.1 Development stages of creep
渦輪葉片的蠕變壽命對溫度的升高十分敏感,研究表明,葉片的溫度每升高10~15 ℃,其蠕變使用壽命將減少一半。同時,蠕變是一個依賴時間的熱變形過程,葉片暴露在高溫下的時間越長,其累積變形量越大。對于發(fā)動機渦輪葉片而言,工程上一般將蠕變損傷作為葉片應(yīng)力、葉片溫度和持續(xù)時間的函數(shù)。
葉片設(shè)計時所要求的持久強度一般在數(shù)萬小時,通過長時間的試驗獲得相關(guān)的蠕變數(shù)據(jù)是不現(xiàn)實的,而想要精確地建立起預(yù)測蠕變損傷和剩余壽命的模型就越發(fā)困難。工程上一般使用由實驗驗證的蠕變持久方程來預(yù)測葉片的蠕變持久強度。目前,國內(nèi)外較為成熟的葉片設(shè)計規(guī)范和手冊普遍使用拉森-米勒(L-M)參數(shù)方程、葛-唐吾(G-D)方程、曼森-索柯普(M-S)方程、曼森-哈弗特(M-H)方程等方法進行蠕變壽命的計算和實驗驗證。選用工程上常用的L-M參數(shù)方程預(yù)測渦輪葉片的蠕變持久強度,該方法基于材料的蠕變試驗數(shù)據(jù),使用拉森米勒參數(shù)來表征材料在不同溫度和應(yīng)力組合下的蠕變斷裂行為。針對葉片特定的某點,拉森米勒參數(shù)滿足:
(1)
式中:為該點的溫度載荷,K;為該點的應(yīng)力水平,MPa;為蠕變斷裂時間,h;為介于17~23之間的常數(shù),一般情況取=20。同時,應(yīng)力水平與拉森米勒參數(shù)之間的關(guān)系還可以用一個玻爾茲曼函數(shù)來表示:
(2)
式中:、、、是玻爾茲曼函數(shù)的參數(shù)。式(2) 可以看作是拉森米勒方程在全應(yīng)力和溫度范圍內(nèi)的擴展,模型的參數(shù)可以由材料在不同溫度下的蠕變斷裂實驗數(shù)據(jù)和抗拉強度數(shù)據(jù)通過插值得到。由式(1)可以推導(dǎo)出在一定的應(yīng)力、溫度載荷作用下,持續(xù)時間內(nèi)的蠕變損傷為
(3)
由于目前航空公司實際數(shù)據(jù)樣本比較難以獲取,基于某機隊歷史飛行數(shù)據(jù),借助發(fā)動機性能仿真模型仿真得到整機關(guān)鍵站位的溫度、壓力、流量等性能參數(shù),進一步基于HPT葉片的幾何參數(shù)、材料參數(shù),構(gòu)建起葉片截面平均溫度和參考點應(yīng)力的數(shù)學(xué)計算經(jīng)驗?zāi)P?,得到每個航班飛行任務(wù)中葉片關(guān)鍵部位的應(yīng)力譜和溫度譜,并基于L-M參數(shù)方程計算每個航班的蠕變損傷量。
蠕變損傷計算涉及到溫度和應(yīng)力,需要進行相應(yīng)的熱分析和應(yīng)力分析。按照分析維度的不同可劃分為:考慮面平均或體平均的0D分析、考慮直線分布的1D分析、考慮平面分布的2D分析、考慮體分布的3D分析。隨著分析維度的提高,模型中包含的單元類型和數(shù)量越來越多,建模的復(fù)雜程度逐漸提高,求解結(jié)果也會越來越趨近于真實情況??梢姺治龅木S度從根本上決定了分析結(jié)果的準(zhǔn)確性,但單純增加分析維度會迅速延長模型的求解時間,模型的輸入也越來越趨向于精細(xì)化,便利性大大降低。參考發(fā)動機部分設(shè)計公式,葉片關(guān)鍵部位溫度和應(yīng)力載荷歷程計算用0D和1D的計算方法進行簡化,在提高計算速度的同時,保證了一定的計算精度。所使用的蠕變損傷集成計算流程如圖2所示,其中包含了用于中間環(huán)節(jié)計算和參數(shù)傳遞的子模型,同時也顯示了不同子模型所需要的輸入數(shù)據(jù)。
從圖2中可以看出,蠕變損傷綜合評估模型由4個子模型組成:
1) 性能仿真模型,仿真發(fā)動機HPT葉片進出口相關(guān)的氣動參數(shù)和熱參數(shù),為葉片的溫度和應(yīng)力計算模型提供輸入條件。
2) 葉片溫度計算模型,對葉片進行溫度分析,并為葉片應(yīng)力計算模型提供輸入條件。
3) 葉片應(yīng)力計算模型,計算葉片參考點的應(yīng)力。
4) 蠕變損傷計算模型,根據(jù)給定的飛參數(shù)據(jù)計算葉片蠕變累積損傷。
圖2 蠕變損傷集成計算流程Fig.2 Integrated calculation process of creep damage
進行葉片的溫度和應(yīng)力的計算時需要提供葉片進出口的速度、溫度、壓力、流量等參數(shù)作為仿真計算輸入的邊界條件。這些參數(shù)往往不能由傳感器直接測量,工程中通常使用一個能夠反映發(fā)動機整機性能的模型來近似替代。基于發(fā)動機性能軟件Turbofan Simulator,通過使用動態(tài)鏈接庫的方法編程調(diào)用,建立能夠匹配實際發(fā)動機的性能參數(shù)的部件級性能仿真模型。性能仿真實際上是對發(fā)動機的氣路進行分析,通過調(diào)整相關(guān)的健康參數(shù),使得發(fā)動機性能模型的輸出與實測的氣路參數(shù)相匹配。其中非線性模型由于精度較高的優(yōu)勢,最先是由Stamatis等引進到工業(yè)燃?xì)廨啓C的氣路分析當(dāng)中,之后這種方法被推廣到了航空發(fā)動機上,隨著應(yīng)用技術(shù)的成熟,最終開發(fā)出了民航發(fā)動機性能模型商業(yè)軟件Turbofan Simulator。Turbofan Simulator能夠?qū)崿F(xiàn)海拔高度、環(huán)境溫度、飛行馬赫數(shù)、控制變量等運行參數(shù)同發(fā)動機整機各部位的溫度、壓力、流量等性能參數(shù)的快速匹配,同時可以調(diào)節(jié)表征發(fā)動機氣路部件性能退化的健康指數(shù),模擬發(fā)動機服役條件下的性能退化行為。性能模型包含了風(fēng)扇(Fan)、低壓壓氣機(Low Pressure Compressor, LPC)、高壓壓氣機(High Pressure Compressor, HPC)、燃燒室、高壓渦輪(High Pressure Turbine, HPT)和低壓渦輪(Low Pressure Turbine, LPT)5個單元體,結(jié)構(gòu)如圖3所示。
性能模型輸入?yún)?shù)包括飛行高度、標(biāo)準(zhǔn)大氣溫差、飛行馬赫數(shù)等;蠕變損傷集成計算模型需要的性能模型輸出參數(shù)包括高壓轉(zhuǎn)子轉(zhuǎn)速、高壓壓氣機進出口壓力、HPT進出口溫度等,詳細(xì)數(shù)據(jù)見表1。
圖3 發(fā)動機性能模型結(jié)構(gòu)示意圖Fig.3 Structural diagram of engine performance model
表1 性能仿真模型輸入和部分輸出參數(shù)
應(yīng)力計算模型主要計算葉片參考點的應(yīng)力,根據(jù)對應(yīng)力成分的分析,將兩種主要的應(yīng)力來源納入應(yīng)力計算模型當(dāng)中:由發(fā)動機轉(zhuǎn)動的離心載荷引起的離心應(yīng)力和由氣體彎矩引起的氣動應(yīng)力。
為了較為全面地考核葉片的應(yīng)力分布情況,針對某型民航發(fā)動機HPT渦輪葉片,首先通過逆向工程獲取其計算機輔助設(shè)計模型,如圖4所示。
圖4 逆向工程得到的HPT葉片三維數(shù)模Fig.4 Three-dimensional model of HPT blade obtained by reverse engineering
進一步采用了一種偽二維解析的方法,按葉片高度取渦輪葉片的4個截面:葉根、25%葉高、50%葉高、75%葉高;根據(jù)葉片可能出現(xiàn)蠕變失效的危險點位置,每個葉片截面上分別取前緣最遠(yuǎn)點(Leading Edge, LE)、后緣最遠(yuǎn)點(Trailing Edge, TE)、葉背最遠(yuǎn)點(Blade Back, BB)作為計算參考點,并使用三維建模軟件對葉片截面的幾何尺寸進行測繪,截面選取的位置和計算參考點位置示意如圖5 和圖6所示。
圖5 葉片截面示意圖Fig.5 Schematic diagram of blade section
圖6 葉片截面尺寸示意圖Fig.6 Schematic diagram of blade section size
取葉根截面、50%葉高截面和葉尖截面,將葉片分為“葉根-葉中”和“葉中-葉尖”兩部分,由于本文中的葉片型號沒有葉冠,因此可近似認(rèn)為葉尖截面的離心應(yīng)力為0、25%葉高截面和75%葉高截面的應(yīng)力可由葉根和葉中截面的計算結(jié)果進行插值得到。
葉根和50%葉高截面的平均離心力和離心應(yīng)力分別為
=
(4)
(5)
式中:為葉片材料密度;為計算截面和葉尖截面的平均面積;為截面到葉尖的高度;為旋轉(zhuǎn)角速度;為“截面—葉尖”部分的重心到旋轉(zhuǎn)中心的距離;為對應(yīng)的截面面積。葉根截面和葉中截面的壓力、彎矩分別為
(6)
=∑()
(7)
式中:為截面截取的葉片段繞旋轉(zhuǎn)軸一周在垂直于旋轉(zhuǎn)軸方向投影形成的環(huán)形面積;為葉片段的平均靜壓差;為葉片的個數(shù);為葉片段的重心到葉根的距離。
流經(jīng)葉片進口和出口的氣流速度變化會導(dǎo)致氣體沿軸向和切向的動量變化,葉根和葉中截面軸向和切向的動量力以及動量彎矩分別為
(8)
(9)
=∑()
(10)
=∑(Tansec)
(11)
式中:為投影的環(huán)狀面積單位面積的質(zhì)量流量;和分別為軸向和切向速度差。葉片截面的方向彎矩和截面參考點的彎矩應(yīng)力分別為
=(+)sin+cos
(12)
=(+)cos-sin
(13)
(14)
式中:為葉片安裝角;和為葉片段產(chǎn)生的彎矩沿主軸方向的投影分量;和分別為截面參考點到截面慣性主軸的距離;和分別為截面的最大和最小慣性矩。葉片截面尺寸的定義如圖6所示,、、分別代表葉片截面前緣、葉背最遠(yuǎn)點、尾緣到軸的距離;、、分別代表葉片截面前緣、葉背最遠(yuǎn)點、尾緣到軸的距離;為葉片安裝角。因此,截面參考點的總應(yīng)力可表示為
=+
(15)
溫度計算模型主要用來計算葉片的基體溫度,本文中涉及的葉片型號主要冷卻方式為氣膜冷卻,采用1D的溫度計算模型不僅能反映蠕變失效的點位,結(jié)合上文中的應(yīng)力計算結(jié)果和后續(xù)的蠕變損傷評估結(jié)果,還能反映不同位置處溫度和應(yīng)力的不同組合對蠕變損傷的影響。類比于應(yīng)力計算模型,計算溫度時將葉片按照葉片長度的0~25%、25%~50%、50%~75%、75%~100%依次劃分為4段,從葉根到葉尖,把前一段的冷卻氣體出口溫度作為后一段的冷卻氣體入口溫度,計算得到葉片葉根、25%葉高、50%葉高、75%葉高共4個截面的平均溫度。
計算時首先需要指定徑向溫度分布因子(Radial Temperature Distribution Factor, RTDF),RTDF的值一般不超過0.2。同時在進行溫度計算前,需要做如下假設(shè):
1) 由于葉片的旋轉(zhuǎn)效應(yīng),葉片最高溫度的位置會從葉片中部向葉尖區(qū)域轉(zhuǎn)移(圖7),因此假設(shè)最大燃?xì)鉁囟任挥诰嚯x葉根75%葉高的位置。
2) 葉根和葉尖的燃?xì)鉁囟葹樽畹腿細(xì)鉁囟取?/p>
3) 從75%葉高到葉根的溫度降低是線性的。
4) 從75%葉高到葉尖的溫度降低是線性的。
最大和最小燃?xì)鉁囟葹?/p>
=+
(16)
(17)
式中:為渦輪進口溫度;為燃燒室溫升;為徑向溫度分布因子。葉片段的平均金屬溫度為
=-(-)
(18)
(19)
(20)
(21)
圖7 徑向溫度分布Fig.7 Radial temperature distribution
使用上述的溫度、應(yīng)力計算模型,每次計算可以得到葉片4個截面、每個截面3個參考點共12個點的溫度和應(yīng)力情況,由于葉片的葉型尺寸和工作條件決定了葉片的蠕變失效位置,因此蠕變失效最先出現(xiàn)在葉片的哪個位置并不固定,表2給出了某航班爬升工況下的溫度和應(yīng)力計算結(jié)果。
從表2中可以看出,截面溫度隨著葉片跨度的增加而逐漸升高,截面應(yīng)力隨著葉片跨度的增加而逐漸減?。辉谕唤孛鎯?nèi),前緣參考點應(yīng)力最大,尾緣參考點應(yīng)力次之,葉背參考點應(yīng)力最小,但應(yīng)力水平差距較小。綜合式(1)和式(2)可計算出表2中溫度和應(yīng)力對應(yīng)的蠕變壽命,見表3。
表2 溫度和應(yīng)力計算結(jié)果Table 2 Temperature and stress calculation results
表3 蠕變壽命計算結(jié)果Table 3 Creep life calculation results
從蠕變壽命計算結(jié)果中可以看出,50%葉高截面的前緣蠕變壽命最短,可以作為整個葉片的危險點。由式(3)可求出所有時刻對應(yīng)的單位時間危險點蠕變損傷,之后遵循Minner線性法則,對一個飛行循環(huán)持續(xù)時間內(nèi)的蠕變損傷進行累加,即可得到某個飛行循環(huán)對應(yīng)的累積蠕變損傷:
(22)
式中:為在第個飛行循環(huán)、第個載荷狀態(tài)下葉片的蠕變持久壽命;為在第個飛行循環(huán)、第個載荷狀態(tài)下的持續(xù)時間;即為第個飛行循環(huán)對應(yīng)的累積蠕變損傷。
在飛機服役過程中,即使有兩臺發(fā)動機的飛行任務(wù)剖面完全一致,甚至于渦輪葉片表面的載荷歷程完全一致,但由于實際的運行環(huán)境和操作條件等因素不可能完全相同,最終導(dǎo)致兩臺發(fā)動機的葉片使用壽命存在較大差異。傳統(tǒng)的基于統(tǒng)計分布的方法雖然可以擬合出葉片失效樣本整體的壽命分布,但很難基于歷史的環(huán)境載荷對單機狀態(tài)進行有效的壽命評估;而對于非設(shè)計單位而言,葉片材料數(shù)據(jù)、三維模型等資料難以獲取,使得精確的有限元仿真計算難以進行,其計算結(jié)果很難反映出實際情況下的不確定性因素對壽命的影響。
對于發(fā)動機來說,飛參是最常見的發(fā)動機動態(tài)時序數(shù)據(jù),可以從機載設(shè)備上直接獲取,其次它作為傳感器數(shù)據(jù),本身就集成了運行環(huán)境、運行載荷等因素的綜合作用,直接反映了發(fā)動機的使用歷程,體現(xiàn)了單機運行條件。將這種動態(tài)數(shù)據(jù)納入分析的范圍也可以更好地解釋發(fā)動機個體剩余壽命的差異性。由第1部分蠕變損傷評估模型可知,葉片的溫度和應(yīng)力直接決定了葉片的蠕變損傷,而高壓轉(zhuǎn)子轉(zhuǎn)速直接決定了葉片的離心應(yīng)力;排氣溫度通常與葉片溫度有著較好的相關(guān)性,這兩種時序數(shù)據(jù)作為常見的狀態(tài)監(jiān)測數(shù)據(jù)記錄在機載快速存儲設(shè)備中,非常容易獲取。因此基于動態(tài)時序數(shù)據(jù)和,定義了一種HPT葉片的蠕變累積損傷指數(shù)模型,用于表征使用環(huán)境載荷等參數(shù)對蠕變壽命分布的影響。通過對Ⅰ型截尾數(shù)據(jù)累積損傷指數(shù)失效分布相關(guān)參數(shù)的最大似然估計,得到基于累積損傷指數(shù)的失效分布函數(shù),并在第4節(jié)對葉片壽命進行評估。
累積損傷指數(shù)(Cumulative damage index, CDI)描述了動態(tài)環(huán)境載荷因素對觀測樣本失效時間分布的影響,環(huán)境載荷歷程稱為樣本的動態(tài)協(xié)變量過程。假設(shè)隨機變量表示樣本的失效時間,隨機變量表示截尾數(shù)據(jù)的失效標(biāo)志位(=1表示觀察結(jié)束前樣本已失效;=0表示在觀察結(jié)束時樣本未失效),()表示樣本空間內(nèi)的個體在一段時間間隔下的協(xié)變量集合。=(,,())為可觀測到的組合隨機變量,對于樣本內(nèi)的第個樣本個體,具體數(shù)據(jù)可表示為(,,()),=1,2,…,,其中為樣本總量;表示第個樣本個體的失效時間;表示第個樣本個體的失效標(biāo)志位;()表示第個樣本個體失效時間內(nèi)的協(xié)變量過程,則時間內(nèi)的累積損傷指數(shù)定義為
(23)
圖8 葉片溫度與排氣溫度散點圖Fig.8 Scatter diagram of blade temperature and exhaust temperature
圖9 葉片應(yīng)力與高壓轉(zhuǎn)子轉(zhuǎn)速平方散點圖Fig.9 Scatter diagram of blade stress and high pressure rotor speed
(24)
式中:、為線性回歸的回歸系數(shù);、為常數(shù)項。忽略、并將式(24)代入式(2)和式(3)中可得:
(25)
忽略式(25)中的次要常數(shù)項,并作進一步化簡
(26)
式中:、、為相關(guān)的系數(shù)。綜合上述分析,使用協(xié)變量和的函數(shù)形式來定義蠕變損傷指數(shù):
(27)
式中:和是損傷指數(shù)的相關(guān)未知參數(shù);=8為給定的常數(shù)項。進而一段時間內(nèi)的蠕變累積損傷指數(shù)()可表示為
(28)
在失效時間內(nèi),每個樣本依據(jù)隨機的協(xié)變量過程(∞)產(chǎn)生對應(yīng)的累積損傷指數(shù),當(dāng)其達(dá)到一定的閾值時樣本失效,假定的累積分布函數(shù)(,)服從威布爾分布,轉(zhuǎn)化為極值分布后滿足:
(29)
(30)
式中:為極值分布函數(shù);為極值分布概率密度函數(shù);=(,)′為極值分布的特征參數(shù)。與失效時間之間的關(guān)系為
(31)
故的分布函數(shù)和概率密度函數(shù)分別為
(;,)=((;,());)
(32)
(;,)=′(();)((;,());)
(33)
渦輪葉片的服役壽命普遍在上萬小時,實際當(dāng)中一般統(tǒng)計截止時間內(nèi)的飛行循環(huán)數(shù)作為葉片的壽命信息,如果截止時間內(nèi)葉片失效,則失效標(biāo)志位記為1;否則失效標(biāo)志位記為0。失效標(biāo)志位和截止時間內(nèi)的葉片壽命信息一同構(gòu)成失效樣本,這類數(shù)據(jù)也稱為截尾數(shù)據(jù)。截尾數(shù)據(jù)中的失效時間是一個含有連續(xù)部分和離散部分的混合隨機變量,綜合式(31)和式(32)可以推導(dǎo)出其對應(yīng)的似然函數(shù)為
(34)
式中:=(,)′和=(,)′是似然函數(shù)中的待估計參數(shù)。結(jié)合式(28)、式(31)、式(32)和式(34),通過迭代求解可以估計出參數(shù)、、和的值。
葉片的剩余壽命是基于當(dāng)前的使用壽命定義的,本身是一個隨機變量。若定義葉片剩余壽命分布(Distribution of Remaining Life, DRL)為
(;)=(<≤+|>,())>0
(35)
則第臺發(fā)動機基于內(nèi)的歷史協(xié)變量過程()=()的條件下,經(jīng)過飛行循環(huán)的剩余壽命分布可表示為
(;)=((<≤+|>,
(),(,+)))
(36)
式中:(,)=(():<<)是~飛行循環(huán)內(nèi)葉片的協(xié)變量歷史;是關(guān)于(,+)|()=()的條件期望。而基于當(dāng)前的歷史協(xié)變信息預(yù)測未來的協(xié)變過程十分復(fù)雜,且直接得到剩余壽命分布的數(shù)學(xué)表達(dá)式較為困難。對于一臺服役狀態(tài)下的某機隊的民用航空發(fā)動機而言,一旦投入到航線開始運行,通常需要重復(fù)執(zhí)行相同的航班任務(wù),其歷史飛行航線結(jié)構(gòu)就相對較為固定,因此可以從歷史航班和對應(yīng)的協(xié)變量信息中抽取樣本,作為未來飛行航班信息的預(yù)測,從而使用蒙特卡洛數(shù)值仿真的方法,對剩余壽命分布進行評估。蒙特卡洛仿真流程如下:
2) 得到更新的協(xié)變量歷史數(shù)據(jù):
7)重復(fù)過程6)次,得到:
要驗證提出方法的有效性,不僅需要葉片的壽命統(tǒng)計數(shù)據(jù),還需要對應(yīng)的發(fā)動機工況監(jiān)測數(shù)據(jù),目前收集的外場樣本數(shù)據(jù)有限,在某機隊航班歷史數(shù)據(jù)基礎(chǔ)上,采用第2節(jié)所述渦輪葉片蠕變壽命仿真方法,通過仿真獲取渦輪葉片的壽命數(shù)據(jù)樣本及其歷史航班的動態(tài)協(xié)變量數(shù)據(jù),驗證提出的模型與方法。
收集的某機隊的歷史飛行數(shù)據(jù)涵蓋了1~12月的航班數(shù)據(jù),設(shè)計轉(zhuǎn)速為14 600 r/min,其中一個航班的和歷程,如圖10所示。
圖11為渦輪葉片蠕變壽命周期仿真流程,其中FN=1,2,…,100為發(fā)動機編號。在滿足飛機起飛安全性的前提下,以相對于發(fā)動機正常起飛推力較小的推力起飛的方式叫做減推力起飛,也稱靈活推力起飛?,F(xiàn)代高涵道比民航發(fā)動機一般具有較高的推力儲備,因此采用減推力起飛不僅完全可行,這種起飛方式不僅能夠延長發(fā)動機的使用壽命,還能夠節(jié)省燃油,降低運行成本。因
圖10 排氣溫度和高壓轉(zhuǎn)子轉(zhuǎn)速歷程Fig.10 Course of exhaust temperature and high pressure rotor speed history
圖11 蠕變壽命周期仿真Fig.11 Creep life cycle simulation
此考慮到飛機的減推力起飛,在進行仿真時將發(fā)動機的額定推力等級考慮在內(nèi),對歷史航班按照4個推力等級劃分,其燃油流量分別為滿推力的100%、97%、94%和91%,分別對應(yīng)T0,T1,T2,T3起飛階段減推力的4個等級;同時考慮到壽命評估的不確定因素,可將不確定性來源總體上劃分為三類:一是葉片設(shè)計階段的材料參數(shù)分布、計算機輔助設(shè)計誤差等引起的不確定性;二是葉片制造階段的材料晶體缺陷、驗收測量精度等引起的不確定性;三是葉片使用階段的自然環(huán)境載荷、發(fā)動機性能退化等引起的不確定性。這些不確定性相互作用,最終導(dǎo)致了服役葉片的實際可用時間與設(shè)計時間不符??紤]到以上不確定性來源的影響,使用一個隨機變量~(0,)代表這些因素的綜合作用,其中為隨機變量方差,由于隨機變量的均值影響的是全壽命樣本整體飛行循環(huán)的偏移,這種數(shù)據(jù)偏移的影響理論上可以通過分布的標(biāo)準(zhǔn)化過程進行消除,因此為了方便后續(xù)處理和分析,隨機變量的均值取為0。壽命周期仿真時,每次從相應(yīng)的樣本空間隨機抽取航班,記錄下飛參數(shù)據(jù)和計算得到的蠕變損傷,則不確定性因素影響下的航班損傷為′=(1+),將′作為此次抽取的實際損傷進行累積,當(dāng)損傷累積到1時停止抽樣,每進行1次這樣的抽樣循環(huán)便得到了一臺發(fā)動機的全壽命樣本,按照上述方法,每個推力等級下各仿真25臺全壽命樣本,4個推力等級共包含100個全壽命樣本,對應(yīng)的壽命分布如圖12所示。
其中1臺發(fā)動機全壽命周期的蠕變損傷量分布,如圖13所示,可以看到發(fā)動機的單個航班蠕變損傷量大約分布于10~10之間,同一臺發(fā)動機不同航班的蠕變損傷呈現(xiàn)較強的分散性,損傷量差距高達(dá)10倍,這主要是不同航班的運行載荷不同導(dǎo)致的結(jié)果。
圖12 HPT葉片壽命分布直方圖Fig.12 Distribution histogram of HPT blade lifetime
取其中一個高損傷與一個低損傷航班,并分別統(tǒng)計起飛爬升階段的數(shù)據(jù),從圖14和圖15中可以看出,蠕變損傷較大的航班對應(yīng)的相對排氣溫度和相對高壓轉(zhuǎn)子轉(zhuǎn)速水平普遍高于蠕變損傷較小的航班,說明服役條件下的載荷水平差異是引起不同航班蠕變損傷差異性的重要因素。
分別統(tǒng)計1臺發(fā)動機主要飛行階段的累積損傷以及相應(yīng)的持續(xù)時間占比,統(tǒng)計結(jié)果如圖16所示。5階段為飛機爬升階段,葉片處于高負(fù)載狀態(tài),此階段持續(xù)時間所占比重較小,約為總時間的12.5%左右,但累積蠕變損傷卻占總損傷的80%以上;6階段為飛機巡航階段,葉片處于低負(fù)載狀態(tài),此階段的時間占比達(dá)到了50%以上,但累積損傷占總損傷不到10%,說明葉片在高溫和高應(yīng)力水平下的單位時間蠕變損傷效應(yīng)已經(jīng)超過了飛行歷程內(nèi)蠕變損傷的時間累積效應(yīng)。
圖13 航班蠕變損傷統(tǒng)計情況Fig.13 Statistics of flight creep damage
圖14 起飛爬升階段相對排氣溫度歷程Fig.14 Relative exhaust temperature course during takeoff and climb phase
圖15 起飛爬升階段相對高壓轉(zhuǎn)子轉(zhuǎn)速歷程Fig.15 Relative high pressure rotor speed course during takeoff and climb phase
圖16 飛行階段損傷統(tǒng)計Fig.16 Damage statistics during flight phase
仿真得到的100臺發(fā)動機樣本中,隨機抽取50臺作為失效樣本,這部分發(fā)動機保留壽命周期的飛行航班數(shù)據(jù)和對應(yīng)的動態(tài)協(xié)變量信息,失效標(biāo)志位記為=1(=1,2,…,50);剩余50臺作為未失效樣本,對飛行循環(huán)進行隨機截尾處理,記錄截尾時刻之前的航班及對應(yīng)的協(xié)變量信息,失效標(biāo)志位記為=0(=51,52,…,100),以此模擬機隊發(fā)動機葉片實際使用壽命的統(tǒng)計情況,其中前20臺發(fā)動機截尾數(shù)據(jù)分布如圖17所示。
圖17 發(fā)動機截尾數(shù)據(jù)Fig.17 Censored engine data
使用R語言編寫似然函數(shù),并調(diào)用optim優(yōu)化函數(shù)進行估計參數(shù)的迭代求解,其中,優(yōu)化方法設(shè)置為BFGS,BFGS方法屬于Quasi-Newton方法,其在最優(yōu)值附近的計算收斂速度快,迭代次數(shù)少。經(jīng)過546次迭代計算,得到式(28)中的參數(shù)估計值:
將、反代可以求出每個航班的累積損傷指數(shù)。將每臺發(fā)動機所有航班累積損傷指數(shù)累加,即可得到每臺發(fā)動機截尾時刻對應(yīng)的總損傷指數(shù),計算結(jié)果的統(tǒng)計情況如圖18和圖19所示。
由圖18和圖19可見:① 失效樣本的累積損傷指數(shù)普遍大于未失效樣本的累積損傷指數(shù),且兩種損傷指數(shù)分界較為明顯,可以很好地區(qū)分失效樣本和未失效樣本,說明本文定義的蠕變累積損傷指數(shù)可以很好的表征HPT葉片的實際損傷狀態(tài)。② 蠕變累積損傷指數(shù)與飛行循環(huán)數(shù)沒有明顯的強相關(guān)性,僅僅借助飛行循環(huán)數(shù)或飛行時
圖18 蠕變累積損傷指數(shù)散點圖Fig.18 Scatter plot of creep cumulative damage index
圖19 蠕變累積損傷指數(shù)直方圖Fig.19 Histogram of creep cumulative damage index
間來表征HPT葉片的累積損傷并不合適,低飛行循環(huán)數(shù)累積損傷指數(shù)可能較大,葉片已經(jīng)發(fā)生蠕變失效,而高飛行循環(huán)數(shù)累積損傷指數(shù)可能較小,葉片仍可以正常服役。這是因為HPT葉片蠕變累積損傷的不僅僅取決于飛行時間,還取決于經(jīng)歷的溫度、應(yīng)力等環(huán)境載荷參數(shù)。③ 失效樣本的累積損傷指數(shù)呈一定的分散性,失效閾值并不是一個確定值,這是因為即使同一型號的HPT葉片,不同個體在“設(shè)計—制造—使用”中不可避免地引進不確定性,這些不確定性的相互作用正是導(dǎo)致服役條件下失效葉片累積損傷指數(shù)發(fā)散的重要原因。
式(30)中的估計參數(shù)值和95%置信區(qū)間邊界值見表4,其中為極值分布的均值,ln()為極值分布標(biāo)準(zhǔn)差以e為底數(shù)的對數(shù)值。從而得到基于累積損傷指數(shù)的失效概率圖(圖20),可以看到,累積損傷指數(shù)基本符合威布爾分布。
對50個未失效樣本的未來協(xié)變過程進行蒙特卡洛仿真,仿真設(shè)置為步長=80 000,單輪仿真次數(shù)=1 000,仿真輪次=100。仿真計算得到的其中四臺發(fā)動機葉片剩余壽命可靠性分布如圖21所示。
表4 蠕變壽命計算結(jié)果Table 4 Creep life calculation results
圖20 隨累積損傷指數(shù)變化的失效概率圖Fig.20 Diagram of failure probability changing with cumulative damage index
圖21 剩余壽命可靠性分布Fig.21 Distribution of remaining life reliability
上述分布表示了從統(tǒng)計截止時間起,未來葉片累積失效率隨飛行循環(huán)數(shù)的變化過程。由于協(xié)變量歷史信息的不同,圖中四臺發(fā)動機剩余壽命可靠性分布具有明顯差異,這也體現(xiàn)了實際服役條件下葉片的個體差異性。針對剩余壽命分布,工程上常用式(37)所示的平均剩余壽命(Mean Residual Life, MRL)作為剩余壽命預(yù)測值。
(37)
式中:()為可靠度函數(shù);為當(dāng)前截止時間。
在得到了葉片的壽命分布函數(shù)的估計參數(shù)后,可以根據(jù)分布類型推導(dǎo)出對應(yīng)的可靠度分布函數(shù),之后代入式(37)中即可計算出與分布對應(yīng)的平均剩余壽命。一般情況下通過直接統(tǒng)計截尾數(shù)據(jù)的方式也可以得到使用壽命的分布函數(shù)。把這種傳統(tǒng)的處理方法記為基準(zhǔn)模型,分別使用本文定義的累積損傷指數(shù)模型與基準(zhǔn)模型按照式(37) 進行平均剩余壽命計算,部分結(jié)果如表5所示,其中使用壽命指的是截止時間的實際使用壽命,剩余壽命是在已知全壽命的基礎(chǔ)上去除使用壽命后剩余的飛行循環(huán)數(shù),這兩部分?jǐn)?shù)據(jù)對應(yīng)著全壽命仿真中模擬的未失效樣本。其中30臺發(fā)動機的預(yù)測偏差如圖22所示。
結(jié)合表5與圖22中的預(yù)測結(jié)果可以看出,基于失效時間分布的剩余壽命預(yù)測偏差大約在-30%~1 200%的范圍內(nèi),說明預(yù)測誤差具有較大的分散性,基本不具備參考價值;使用累積損傷指數(shù)進行的預(yù)測偏差大約在-20%~20%之間分
表5 平均剩余壽命預(yù)測結(jié)果對比Table 5 Comparison of MRL prediction results
圖22 預(yù)測結(jié)果偏差Fig.22 Forecast result deviation
布,說明預(yù)測誤差分散性小,結(jié)果總體在可接受范圍內(nèi),可作為維修計劃定制的參考剩余壽命。
1) 針對民航發(fā)動機渦輪葉片的剩余壽命預(yù)測,提出了一種數(shù)據(jù)驅(qū)動與失效物理融合的可靠性剩余壽命預(yù)測方法,基于服役條件下的歷史工況監(jiān)測數(shù)據(jù),結(jié)合發(fā)動機性能仿真模型、葉片應(yīng)力和溫度計算模型對葉片表面參考點的溫度和應(yīng)力歷程進行計算;在此基礎(chǔ)上使用拉森米勒參數(shù)法對葉片蠕變損傷進行評估;之后針對蠕變損傷定義了一種累積損傷指數(shù)模型,融合歷史協(xié)變量信息排氣溫度和高壓轉(zhuǎn)子轉(zhuǎn)速對葉片剩余壽命進行了可靠性評估,并通過數(shù)據(jù)仿真分析驗證了本文建立模型的可行性,為渦輪葉片可靠性評估提供了一種新的思路。
2) 針對服役條件下的HPT葉片蠕變壽命預(yù)測,使用本文的模型方法得到的預(yù)測結(jié)果,預(yù)測偏差大約在±20%范圍內(nèi)分布,而傳統(tǒng)的基于失效時間分布的可靠性分析方法,得到的預(yù)測偏差約在-30%~1 200%之間,說明提出的模型具有更好的適用性,在壽命預(yù)測精度上也更具有優(yōu)勢。
3) 反映發(fā)動機關(guān)鍵件服役可靠性降級等行為特性的信息存在于多模態(tài)數(shù)據(jù)中,包括產(chǎn)品設(shè)計數(shù)據(jù)與使用維護數(shù)據(jù),如產(chǎn)品幾何參數(shù)、材料數(shù)據(jù)、狀態(tài)監(jiān)測數(shù)據(jù)、離線檢測/檢查報告、故障記錄和大修報告等,貫穿產(chǎn)品壽命周期的不同階段,如何充分利用這些多模態(tài)數(shù)據(jù),實現(xiàn)多領(lǐng)域物理模型與數(shù)據(jù)融合,是提高服役發(fā)動機關(guān)鍵件壽命預(yù)測分析準(zhǔn)確性的重要途徑。以民航發(fā)動機HPT葉片蠕變失效為例,在蠕變失效機理分析基礎(chǔ)上,同時考慮服役條件下的使用信息、狀態(tài)參數(shù)及截尾失效時間等數(shù)據(jù),提出的蠕變累積損傷指數(shù)模型實現(xiàn)了多模態(tài)信息的融合,對渦輪葉片服役可靠性進行評估,并預(yù)測其剩余壽命。
4) 本文中涉及的葉片損傷正向計算中,每個計算環(huán)節(jié)都會不可避免地引進計算的不確定性,這些不確定性既有認(rèn)知的不確定性,也有材料、設(shè)計等固有的不確定性,這是正向計算所不可避免的。因此“數(shù)據(jù)驅(qū)動與失效物理融合”是對可靠性壽命評估方面創(chuàng)新性的思考,它主要體現(xiàn)在累積損傷指數(shù)模型的建立過程以及外場使用數(shù)據(jù)、截尾壽命數(shù)據(jù)的融合過程中。損傷指數(shù)定義的方式基于一系列的正向的計算基礎(chǔ)上的推導(dǎo),而發(fā)動機性能模型、溫度計算模型、應(yīng)力計算模型和失效模式對應(yīng)的損傷計算物理模型都與失效物理相對應(yīng);之后通過分析蠕變失效機理以及物理模型,推導(dǎo)了一種新的蠕變損傷累積指數(shù)模型,在此初始模型的基礎(chǔ)上進一步建立似然函數(shù),融合外場使用數(shù)據(jù)、截尾壽命數(shù)據(jù)來估計模型參數(shù),從而實現(xiàn)了在失效物理機理基礎(chǔ)上對外場數(shù)據(jù)的融合。模型直接用到的是飛行監(jiān)測數(shù)據(jù)與外場失效時間數(shù)據(jù),而由于累積損傷模型基于失效模式定義,因此預(yù)測結(jié)果在物理層面的解釋性更好。
現(xiàn)階段外場數(shù)據(jù)獲取有限,因此通過仿真數(shù)據(jù)驗證了提出方法的可行性和有效性。下一步的工作將基于目前的研究繼續(xù)深入,從以下幾個方面開展:一是從熱機械疲勞、高溫氧化等其他失效機理入手建立對應(yīng)的累積損傷指數(shù)模型,驗證本文方法的通用性和魯棒性;二是研究有效的動態(tài)協(xié)變量預(yù)測方法,來代替目前使用的蒙特卡洛仿真過程,以取得更好的參數(shù)估計結(jié)果和預(yù)測精度;三是跟蹤預(yù)測部件的實時狀態(tài),對本文的模型方法進行實際數(shù)據(jù)的驗證。