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

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

1 引言

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

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

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

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

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

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

2.1 計算方法

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

2.2 研究模型

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

圖1 研究模型Fig.1 Research model

2.3 有限元網(wǎng)格

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

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

2.4 邊界條件和材料屬性

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

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

2.5 瞬態(tài)溫度場計算結(jié)果

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

圖3 邊界條件隨時間的變化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)力場計算結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 轉(zhuǎn)化過程

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

表2 每個點雨流計數(shù)程序運行結(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)下壓氣機輪盤子午面的全局和局部等效應(yīng)力的分布云圖?？煽闯?，與最大應(yīng)力時刻的局部等效應(yīng)力相比，全局等效應(yīng)力的應(yīng)力水平顯著下降，但分布規(guī)律相似。

4 概率失效風(fēng)險分析計算

4.1 計算方法

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

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

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

4.2 計算結(jié)果

將全局等效應(yīng)力分布與基于局部的應(yīng)力分布分別代入概率失效風(fēng)險分析流程，得到兩種不同應(yīng)力轉(zhuǎn)化方法下的壓氣機盤失效概率，如圖12所示。結(jié)果表明，使用循環(huán)數(shù)達到20 000時，基于局部應(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)化方法具有更小的失效概率。盡管這一過程需要后期的驗證性試驗數(shù)據(jù)支持，但上述典型示例表明轉(zhuǎn)換方法不同對失效概率有著顯著且重要的影響。

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

5 結(jié)論

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

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