航空發(fā)動(dòng)機(jī)限壽件概率失效風(fēng)險(xiǎn)評(píng)估的等效應(yīng)力轉(zhuǎn)化

丁水汀，郭祎瑋，2，李 果，3，周 煜

(1.北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院航空發(fā)動(dòng)機(jī)氣動(dòng)熱力國家級(jí)重點(diǎn)實(shí)驗(yàn)室，北京100191；2.中國航發(fā)商用航空發(fā)動(dòng)機(jī)有限責(zé)任公司，上海200241；3.多倫多大學(xué)航空航天研究所，多倫多M3H5T6；4.北京航空航天大學(xué)交通科學(xué)與工程學(xué)院，北京100191)

1 引言

通常，因失效可能引起航空發(fā)動(dòng)機(jī)危害性事故的部件定義為航空發(fā)動(dòng)機(jī)限壽件[1]，包括發(fā)動(dòng)機(jī)的主要轉(zhuǎn)子結(jié)構(gòu)件和部分靜子結(jié)構(gòu)件，如各類旋轉(zhuǎn)輪盤和承壓機(jī)匣等。航空發(fā)動(dòng)機(jī)適航規(guī)章FAR/CCAR-33.70條款及AC 33.14、AC 33.70-1要求，限壽件必須開展壽命期內(nèi)的概率失效風(fēng)險(xiǎn)評(píng)估工作。開展概率風(fēng)險(xiǎn)評(píng)估可在限壽件設(shè)計(jì)生產(chǎn)階段進(jìn)行，預(yù)測(cè)其是否滿足適航規(guī)章[2]對(duì)失效概率的要求。

有關(guān)限壽件概率失效風(fēng)險(xiǎn)的研究，國外從評(píng)估方法到基礎(chǔ)數(shù)據(jù)庫積累均取得了良好進(jìn)展。美國聯(lián)邦航空局(FAA)與西南研究院等科研機(jī)構(gòu)和GE、普惠等航空發(fā)動(dòng)機(jī)制造商合作，共同開發(fā)了成熟的概率失效風(fēng)險(xiǎn)分析方法。該方法起初由Leverant[3]、McClung[4]等提出，目標(biāo)是解決鈦合金壓氣機(jī)盤材料Hard α夾雜缺陷引起的失效問題。之后，Momin[5]和Wu[6]等進(jìn)一步研究了分析流程，分別從斷裂力學(xué)模型建立、隨機(jī)變量分布及無損探傷檢查等方面對(duì)失效風(fēng)險(xiǎn)分析方法進(jìn)行了最終的完善，形成了達(dá)爾文軟件[3-6]，從而建立了一整套完整的分析理論及流程。

受制于技術(shù)封鎖和經(jīng)驗(yàn)缺乏，國內(nèi)對(duì)限壽件概率失效風(fēng)險(xiǎn)分析的研究目前尚處于起步階段。丁水汀等[7]通過對(duì)適航條款及國內(nèi)外相關(guān)研究的分析，對(duì)概率失效風(fēng)險(xiǎn)評(píng)估方法進(jìn)行了研究；趙合陽等[8]利用GH4133材料低循環(huán)疲勞試驗(yàn)數(shù)據(jù)，對(duì)渦輪盤低循環(huán)疲勞失效進(jìn)行了概率壽命分析；杜少輝等[9]在裂紋擴(kuò)展安全壽命的基礎(chǔ)上，對(duì)輪盤的失效概率進(jìn)行了分析。

目前，概率失效風(fēng)險(xiǎn)評(píng)估主要通過斷裂力學(xué)來分析部件的失效[10]，利用應(yīng)力-強(qiáng)度干涉理論，并考慮材料參數(shù)、缺陷分布及載荷等變量的隨機(jī)性，對(duì)失效概率進(jìn)行定量分析[3，11]。其中，多場(chǎng)耦合的應(yīng)力分布結(jié)果是失效風(fēng)險(xiǎn)分析的關(guān)鍵輸入條件。通常，限壽件在發(fā)動(dòng)機(jī)飛行循環(huán)中承受的應(yīng)力是隨時(shí)間變化的瞬態(tài)結(jié)果，但在失效風(fēng)險(xiǎn)分析中要求應(yīng)力只能以一個(gè)載荷對(duì)的形式來執(zhí)行斷裂力學(xué)計(jì)算[2]。因此，對(duì)于限壽件應(yīng)力分析，必需將飛行循環(huán)內(nèi)的變化應(yīng)力轉(zhuǎn)化為斷裂力學(xué)計(jì)算所需的一個(gè)載荷對(duì)。為此，斷裂力學(xué)計(jì)算時(shí)，通常，選取區(qū)域內(nèi)最大應(yīng)力作為整個(gè)區(qū)域的等效應(yīng)力[12]，即基于局部的等效應(yīng)力轉(zhuǎn)化方法。但這一方法在計(jì)算確定壽命期內(nèi)的概率失效風(fēng)險(xiǎn)時(shí)，得到的結(jié)果通常偏保守。

本文以某型航空發(fā)動(dòng)機(jī)離心壓氣機(jī)盤這一典型限壽件為例，提出一種將區(qū)域內(nèi)所有節(jié)點(diǎn)在飛行循環(huán)內(nèi)的應(yīng)力轉(zhuǎn)化為等效應(yīng)力的方法，即基于全局的應(yīng)力轉(zhuǎn)化方法。該方法考慮了瞬態(tài)情況，且盡可能反映載荷在過渡過程中的變化，更貼近發(fā)動(dòng)機(jī)真實(shí)運(yùn)轉(zhuǎn)情況，可獲得更為準(zhǔn)確的失效風(fēng)險(xiǎn)評(píng)估結(jié)果。

2 瞬態(tài)載荷計(jì)算

2.1 計(jì)算方法

采用單向流固耦合[13-15]方法對(duì)輪盤類限壽件進(jìn)行溫度及應(yīng)力計(jì)算，以獲得壓氣機(jī)盤隨時(shí)間變化的應(yīng)力載荷分布。即先采用CFX將瞬態(tài)邊界條件代入分析流程計(jì)算瞬態(tài)溫度載荷，隨后將溫度結(jié)果以載荷步的形式導(dǎo)入ANSYS模塊進(jìn)行靜力分析，獲得瞬態(tài)應(yīng)力分布。

2.2 研究模型

以某型航空發(fā)動(dòng)機(jī)離心壓氣機(jī)輪盤為研究對(duì)象(圖1(a))，建立流體域模型(圖1(b))。為減少分析中的計(jì)算量，取輪盤的1/10建立模型(圖1(c))。

圖1 研究模型Fig.1 Research model

2.3 有限元網(wǎng)格

圖2為輪盤子午面網(wǎng)格節(jié)點(diǎn)分布。如圖所示，概率失效風(fēng)險(xiǎn)評(píng)估建立在多場(chǎng)耦合獲取的輪盤載荷分布的基礎(chǔ)之上，因此選用輪盤子午面載荷分布來表征輪盤載荷分布結(jié)果。同時(shí)，分網(wǎng)時(shí)對(duì)有限元網(wǎng)格加以約束，即該截面上網(wǎng)格節(jié)點(diǎn)與流體網(wǎng)格具有完整的節(jié)點(diǎn)對(duì)應(yīng)關(guān)系，以滿足后續(xù)概率失效風(fēng)險(xiǎn)評(píng)估需要。輪盤子午面網(wǎng)格節(jié)點(diǎn)共1 300個(gè)。

圖2 含完整節(jié)點(diǎn)的輪盤子午面Fig.2 Meridian plane of disc with complete nodes

2.4 邊界條件和材料屬性

針對(duì)分析模型，由空氣系統(tǒng)程序獲得典型飛行循環(huán)瞬態(tài)工況的邊界條件。由于本文的主要目的是給出等效應(yīng)力轉(zhuǎn)化方法，因此對(duì)邊界條件進(jìn)行了簡(jiǎn)化以模擬發(fā)動(dòng)機(jī)的工作過程，見圖3。以圖3(a)中轉(zhuǎn)速隨時(shí)間的變化為例，隨著發(fā)動(dòng)機(jī)起動(dòng)，壓氣機(jī)轉(zhuǎn)速逐漸增加，經(jīng)歷一段時(shí)間波動(dòng)后趨于穩(wěn)定，最終工作一段時(shí)間后轉(zhuǎn)速下降。整個(gè)過程歷時(shí)10.0 s，其他邊界條件與轉(zhuǎn)速相匹配，基本模擬了發(fā)動(dòng)機(jī)起動(dòng)-巡航-減速-停機(jī)這4個(gè)主要工作狀態(tài)。需指出，該方法可擴(kuò)展至全部飛行循環(huán)。

壓氣機(jī)盤的材料為鈦合金TC4，其各項(xiàng)屬性見表1[16]。表中，a為熱膨脹系數(shù)，E為彈性模量，ρ為密度，λ為導(dǎo)熱系數(shù)，μ為泊松比。分析中為簡(jiǎn)化起見，不考慮輪盤材料各屬性隨溫度的變化；流體域介質(zhì)為理想可壓縮氣體，且邊界無滑移，壁面光滑[17]。

2.5 瞬態(tài)溫度場(chǎng)計(jì)算結(jié)果

利用CFX軟件對(duì)已繪制好的有限元網(wǎng)格加入相應(yīng)的瞬態(tài)邊界條件進(jìn)行計(jì)算。由于模型取自整盤的1/10，對(duì)流體和固體截面加入周期性對(duì)稱邊界條件。為得到較好的瞬態(tài)收斂結(jié)果，以穩(wěn)定段邊界計(jì)算出的穩(wěn)態(tài)結(jié)果作為初值[18-19]。獲得輪盤工作過程瞬態(tài)溫度場(chǎng)后，將CFX計(jì)算結(jié)果按時(shí)間步導(dǎo)出節(jié)點(diǎn)溫度數(shù)據(jù)。同時(shí)利用2.3節(jié)輪盤子午面的節(jié)點(diǎn)及單元和節(jié)點(diǎn)的對(duì)應(yīng)關(guān)系，通過MATLAB編制程序，實(shí)現(xiàn)每一時(shí)間點(diǎn)對(duì)應(yīng)溫度數(shù)據(jù)的自動(dòng)導(dǎo)入，并繪制各時(shí)刻子午面溫度分布云圖，見圖4。圖中,橫縱坐標(biāo)軸代表子午面上節(jié)點(diǎn)坐標(biāo)，不同顏色表示對(duì)應(yīng)節(jié)點(diǎn)溫度值，下文類似。圖中選取了初始(0 s)、應(yīng)力最大(0.5 s)、應(yīng)力波動(dòng)(1.0 s)、穩(wěn)態(tài)(4.0 s和 8.0s)及終止(10.0 s)6個(gè)最有代表性的時(shí)刻，來表征輪盤工作過程中載荷隨時(shí)間的變化。

圖3 邊界條件隨時(shí)間的變化Fig.3 Boundary conditions change as a function of time

表1 TC4材料屬性Table 1 Material properties of TC4

圖4 輪盤子午面溫度分布云圖Fig.4 Temperature distribution nephogram on meridian surface of compressor disk

2.6 瞬態(tài)應(yīng)力場(chǎng)計(jì)算結(jié)果

應(yīng)力場(chǎng)采用熱固間接耦合分析方法[20]計(jì)算，其中間接法體現(xiàn)在計(jì)算時(shí)按照載荷步循環(huán)重復(fù)進(jìn)行應(yīng)力場(chǎng)計(jì)算。本研究中，每一時(shí)刻的應(yīng)力場(chǎng)計(jì)算均是基于先期獲得的溫度場(chǎng)作為熱邊界條件施加在模型上，并通過MATLAB程序?qū)崿F(xiàn)循環(huán)。圖5給出了各時(shí)刻子午面應(yīng)力分布的有限元云圖。圖6為在發(fā)動(dòng)機(jī)一個(gè)飛行循環(huán)中，一級(jí)壓氣機(jī)盤最大應(yīng)力點(diǎn)處應(yīng)力隨時(shí)間的變化。

圖5 輪盤子午面應(yīng)力分布云圖Fig.5 Stress distribution nephogram on meridian surface of compressor disk

圖6 一級(jí)壓氣機(jī)盤最大應(yīng)力點(diǎn)處應(yīng)力隨時(shí)間的變化Fig.6 Stress changes with time at peak stress

3 全局等效應(yīng)力轉(zhuǎn)化

本文參考傳統(tǒng)安全壽命理論[21-23]，考慮發(fā)動(dòng)機(jī)整個(gè)飛行循環(huán)，在等損傷前提下，提出并采用了一種基于全局的應(yīng)力轉(zhuǎn)化方法獲得概率失效風(fēng)險(xiǎn)評(píng)估的應(yīng)力輸入條件。其具體流程如圖7所示。

圖7 全局應(yīng)力轉(zhuǎn)化流程Fig.7 Global stress transformation process

3.1 雨流計(jì)數(shù)法及線性累計(jì)損傷理論

利用安全壽命理論進(jìn)行低周疲勞應(yīng)力變化分析時(shí)，通常需要利用循環(huán)計(jì)數(shù)法將輪盤不規(guī)則的應(yīng)力隨時(shí)間的變化轉(zhuǎn)化為若干個(gè)應(yīng)力循環(huán)對(duì)。目前工程中使用最廣泛并為美國FAA認(rèn)可的循環(huán)計(jì)數(shù)法是雨流計(jì)數(shù)法[24]。

利用雨流計(jì)數(shù)法獲得與瞬態(tài)載荷等效的一組載荷對(duì)后，通過查尋對(duì)應(yīng)材料的S-N曲線可得到每個(gè)載荷對(duì)對(duì)應(yīng)的循環(huán)壽命Ni。由于計(jì)數(shù)后每個(gè)載荷對(duì)的應(yīng)力比不同，需要利用Goodman方程進(jìn)行統(tǒng)一：

式中：σa為任意實(shí)際應(yīng)力循環(huán)對(duì)的應(yīng)力幅值，σm為平均應(yīng)力，σb為材料靜強(qiáng)度。

根據(jù)線性累計(jì)損傷理論，若構(gòu)件在恒幅交變載荷Si作用下對(duì)應(yīng)循環(huán)次數(shù)ni，則疲勞損傷Di為：

對(duì)于每個(gè)轉(zhuǎn)化后的等幅載荷，若破壞準(zhǔn)則為：

可以認(rèn)為，若每一個(gè)等幅載荷單獨(dú)加載的疲勞壽命之和D1與一個(gè)載荷D經(jīng)過相同多次加載后的疲勞壽命D2相同，則載荷可以等效。

最后，反查S-N曲線即可得到等效的應(yīng)力載荷對(duì)[21-23]。至此，通過線性累計(jì)損傷理論即可把雨流計(jì)數(shù)后的一組載荷對(duì)轉(zhuǎn)化為可用于斷裂力學(xué)計(jì)算的一個(gè)載荷。

3.2 轉(zhuǎn)化過程

對(duì)壓氣機(jī)盤進(jìn)行瞬態(tài)分析時(shí)，需要對(duì)盤-盤腔模型每個(gè)節(jié)點(diǎn)在整個(gè)發(fā)動(dòng)機(jī)飛行循環(huán)中的應(yīng)力變化進(jìn)行雨流計(jì)數(shù)處理。為提高效率，通過C語言編制“四點(diǎn)法雨流計(jì)數(shù)”程序[25-26]完成轉(zhuǎn)換工作。該程序主要由兩部分組成：①簡(jiǎn)化數(shù)據(jù)，剔除不會(huì)對(duì)材料造成疲勞損傷的、相鄰的、應(yīng)力值相同的點(diǎn)，提取瞬態(tài)應(yīng)力的波峰和波谷；②依照四點(diǎn)法編寫雨流計(jì)數(shù)程序，其計(jì)算結(jié)果為每個(gè)點(diǎn)所受的與瞬態(tài)應(yīng)力等效的恒幅載荷，如表2所示。表中，恒幅載荷最大應(yīng)力幅SI和恒幅載荷最小應(yīng)力幅SI′構(gòu)成了Si。根據(jù)材料手冊(cè)確定應(yīng)力轉(zhuǎn)化中所使用的S-N曲線，且材料屈服極限為973 MPa，應(yīng)力比為-1。S-N曲線可轉(zhuǎn)化為如圖8中便于程序計(jì)算的標(biāo)準(zhǔn)冪函數(shù)形式[27]。圖9為全局應(yīng)力轉(zhuǎn)化程序框圖。首先將雨流計(jì)數(shù)法得到的一組載荷對(duì)利用Goodman方程統(tǒng)一為應(yīng)力比為-1后，再導(dǎo)入S-N曲線得到對(duì)應(yīng)的壽命，進(jìn)而利用線性累計(jì)損傷理論計(jì)算全局等效應(yīng)力，最后通過程序得到整個(gè)飛行循環(huán)輪盤的全局等效應(yīng)力分布。

表2 每個(gè)點(diǎn)雨流計(jì)數(shù)程序運(yùn)行結(jié)果Table 2 Rain-flow program operation results

圖8 冪函數(shù)擬合S-N曲線Fig.8 Power function fittingS-Ncurve

圖9 全局應(yīng)力轉(zhuǎn)化程序框圖Fig.9 Global stress transformation program diagram

3.3 轉(zhuǎn)化結(jié)果

圖10給出了經(jīng)轉(zhuǎn)化的典型飛行循環(huán)下壓氣機(jī)輪盤子午面的全局和局部等效應(yīng)力的分布云圖。可看出，與最大應(yīng)力時(shí)刻的局部等效應(yīng)力相比，全局等效應(yīng)力的應(yīng)力水平顯著下降，但分布規(guī)律相似。

4 概率失效風(fēng)險(xiǎn)分析計(jì)算

4.1 計(jì)算方法

圖10 兩種等效應(yīng)力分布云圖Fig.10 Two types of equivalent stress distribution nephograms

根據(jù)文獻(xiàn)[28]中的概率失效風(fēng)險(xiǎn)計(jì)算方法，采用分區(qū)評(píng)估方式，并結(jié)合應(yīng)力-強(qiáng)度干涉原理、線彈性斷裂力學(xué)及概率統(tǒng)計(jì)建立分析模型；應(yīng)用Monte Carlo法、Runge-Kutta法和Newman法，確定模型的數(shù)值解法[29-31]。在此基礎(chǔ)上，可以將概率失效風(fēng)險(xiǎn)評(píng)估的流程，編制為失效風(fēng)險(xiǎn)的計(jì)算程序。其中，限壽件的應(yīng)力分布將影響斷裂力學(xué)計(jì)算中裂紋擴(kuò)展驅(qū)動(dòng)力的大小從而影響輪盤失效風(fēng)險(xiǎn)，限壽件的溫度分布變化會(huì)影響材料的裂紋擴(kuò)展速率和斷裂韌性，如圖11所示。因此，當(dāng)失效概率在某一應(yīng)力水平(圖中為500 MPa)之下時(shí)，可認(rèn)為不發(fā)生失效；當(dāng)應(yīng)力超過這一水平時(shí)，失效概率增長率隨應(yīng)力的增加而增加，直到完全失效為止。

圖11 失效概率隨應(yīng)力和溫度分布的變化Fig.11 The change of failure probability with the stress and temperature distribution

4.2 計(jì)算結(jié)果

將全局等效應(yīng)力分布與基于局部的應(yīng)力分布分別代入概率失效風(fēng)險(xiǎn)分析流程，得到兩種不同應(yīng)力轉(zhuǎn)化方法下的壓氣機(jī)盤失效概率，如圖12所示。結(jié)果表明，使用循環(huán)數(shù)達(dá)到20 000時(shí)，基于局部應(yīng)力轉(zhuǎn)化方法獲得的失效概率為2.2 10-11次/飛行循環(huán)，而基于全局應(yīng)力轉(zhuǎn)化的為4.8 10-12次/飛行循環(huán)。這表明，與基于局部應(yīng)力轉(zhuǎn)化的方法相比，采用FAA所推薦的基于全局應(yīng)力轉(zhuǎn)化方法具有更小的失效概率。盡管這一過程需要后期的驗(yàn)證性試驗(yàn)數(shù)據(jù)支持，但上述典型示例表明轉(zhuǎn)換方法不同對(duì)失效概率有著顯著且重要的影響。

圖12 不同應(yīng)力轉(zhuǎn)化方法下的壓氣機(jī)盤失效概率Fig.12 Failure probability of compressor disk under two different stress transformation methods

5 結(jié)論

通過對(duì)概率失效風(fēng)險(xiǎn)評(píng)估輸入條件的兩類轉(zhuǎn)化方法進(jìn)行對(duì)比，在相同飛行循環(huán)(20 000)內(nèi)，采用全局應(yīng)力轉(zhuǎn)化方法獲得的等效應(yīng)力分布作為輸入條件代入概率失效風(fēng)險(xiǎn)評(píng)估流程時(shí)獲得的失效概率更小(全局 4.8 10-12次/飛行循環(huán)，局部 2.2 10-11次/飛行循環(huán))，表明不同等效應(yīng)力轉(zhuǎn)化方法對(duì)概率失效風(fēng)險(xiǎn)有顯著影響。同時(shí)，也揭示出基于局部應(yīng)力轉(zhuǎn)化方法獲得的失效概率是否過于保守，是我國民機(jī)型號(hào)未來適航取證過程中需仔細(xì)考慮的問題。

本文的研究成果對(duì)概率失效風(fēng)險(xiǎn)評(píng)估中關(guān)鍵輸入?yún)?shù)應(yīng)力分布數(shù)據(jù)的確定給出了定量參考，后續(xù)可以通過協(xié)調(diào)行業(yè)內(nèi)相關(guān)發(fā)動(dòng)機(jī)研究所，以真實(shí)壓氣機(jī)盤為分析對(duì)象，獲取工作狀態(tài)下輪盤各項(xiàng)邊界條件的變化數(shù)據(jù)，以失效風(fēng)險(xiǎn)分析的數(shù)值計(jì)算結(jié)果對(duì)試驗(yàn)結(jié)果進(jìn)行評(píng)估。