殷 緣 ，楊 勱 ，信 琦 ，陳云霞

（1.北京航空航天大學可靠性與系統(tǒng)工程學院，北京100191；2.中國船舶工業(yè)系統(tǒng)研究院，北京100094；3.中國航發(fā)沈陽發(fā)動機研究所，沈陽110015）

0 引言

航空發(fā)動機附件機匣是航空發(fā)動機的重要部件，對其力學行為和壽命進行分析以獲得可靠性和壽命指標結果至關重要。附件機匣在工作時，不僅受到軸承載荷、自身和安裝附件重力的影響，同時由于齒輪系統(tǒng)高速轉動產熱以及噴油散熱使得其始終處于較復雜的溫度場中[1]，且由于附件機匣內部齒輪系統(tǒng)嚙合傳動以及在安裝、運行過程中的誤差和故障等因素導致其處于復雜的振動環(huán)境中[2]?？紤]到附件機匣工作環(huán)境的復雜多樣性，以及同時承受機械載荷、振動載荷以及熱載荷，李錦花等結合MASTA軟件和ANSYS軟件計算同時受到軸承載荷和傳動附件安裝彎矩影響下殼體的變形[3]；郭梅等將輪齒嚙合時變剛度和靜傳遞誤差產生的嚙合力作為系統(tǒng)振動激勵，采用有限元方法分析附件機匣系統(tǒng)的振動響應[4]；吳鴻等采用有限元方法建立附件機匣殼體穩(wěn)態(tài)熱分析模型，給出殼體的溫度場分布[5]。然而，上述研究僅包含機械載荷、振動或者溫度場中的某一方面，不能全面反映其在復雜工作環(huán)境和熱固條件下的振動載荷情況。熱、固、振載荷會影響附件機匣應力分布等力學行為，從而影響其疲勞壽命。因此，在熱、固、振條件共同作用下進行附件機匣的力學行為以及壽命建模分析是十分重要的。

本文基于有限元軟件ANSYS Workbench，同時考慮附件機匣殼體自身重力和固定約束條件、軸承載荷、溫度場以及振動載荷作用，在不同工況下進行附件機匣殼體力學行為分析，計算得到附件機匣應力響應PSD譜，并根據(jù)應力響應PSD譜采用雨流循環(huán)計數(shù)方法對附件機匣壽命指標進行估計。

1 力學行為建模及分析

在熱、固、振多載荷耦合作用下，運用ANSYS Workbench對其力學行為進行建模及分析。

1.1 模型定義、邊界以及輸入條件

1.1.1 模型定義

在實際工作環(huán)境中，各種負載懸掛在附件機匣兩側，對其力學行為分析產生影響。因此，在CAD模型上添加質點塊(point mass)來定義負載對附件機匣的慣性影響，如圖1所示。

圖1 模型定義

圖1定義模型中各負載的質量以及在局部坐標系中的位置見表1。

1.1.2 邊界條件

進行ANSYS有限元仿真時，邊界條件指運動邊界上方程組的解應該滿足的條件。針對附件機匣的邊界條件為吊耳處固定約束。

表1 負載附件的質量以及坐標

1.1.3 輸入條件

進行附件機匣熱、固、振耦合作用下的力學行為分析，涉及機械載荷、振動載荷和熱，因此輸入條件有：機匣重力、溫度場、軸承載荷、PSD譜。

根據(jù)ANSYS FLUENT產熱分析，得到附件機匣殼體上各部位溫度，并采用EXCEL擬合得到附件機匣溫度分布函數(shù)。環(huán)境溫度設置為90℃，在附件機匣中部建立局部坐標系，以該坐標系為基準添加溫度，溫度變化圖線如圖2所示。

圖2 附件機匣不同部位的溫度

因此，溫度變化函數(shù)為

式中：x為附件機匣與坐標原點的相對位置；y為該部位的溫度值。

圖3 附件機匣溫度場

加載溫度場到附件機匣上，如圖3所示。附件機匣內部共有10個齒輪軸，每個軸兩端各有1個軸承，僅承受徑向力。采用LMS Virtual Lab.Motion進行整體仿真，選取每個軸承的最大徑向力并加載到附件機匣上，具體數(shù)值見表2。

表2 軸承載荷

隨機振動輸入的加度素功率譜頻率范圍為0～2000 Hz，各頻段振動量值見表3。

表3 該工況下每頻段振動量值

1.2 網格劃分

附件機匣殼體幾何機構較為復雜，建議采用快速自動生成的四面體網格進行劃分。此外，在劃分網格時選擇合適的網格大小對求解的精度以及求解速度十分重要。

結合附件機匣實際尺寸，使用全局8 mm、全局16 mm、附件機匣兩側(平板處)32 mm其余8 mm、兩側16 mm其余8 mm 4種尺寸劃分網格，結果見表4。

表4 網格劃分結果相關參數(shù)

其中，畸變度(skewness)為網格質量的重要參數(shù)，指單元相對其理想形狀的相對扭曲的度量，是1個值在0(極好的)到1(無法接受的)之間的比例因子。根據(jù)表4可知，不同網格尺寸得到的畸變度最大值和平均值基本接近；節(jié)點數(shù)和單元數(shù)為尺寸8 mm的最少，選擇較少節(jié)點和單元的尺寸可以提高運算速度。

對于8 mm所得的單元和節(jié)點數(shù)比16 mm得到的少，其原因在于8 mm時的平滑參數(shù)(transition)設置為fast，而16 mm設置為slow。若8 mm網格尺寸設置為slow則會輸出過多單元，導致計算量太大而不能進行模態(tài)分析；若16 mm設置為fast則網格質量太低不足以比較。而網格尺寸為8 mm時畸變度為0.95的單元所占比例最小，因此8 mm的網格尺寸最理想，其對應的機匣局部網格劃分結果如圖4所示。

圖4 附件機匣網格

1.3 預應力模態(tài)分析

模態(tài)分析用來計算結構的振動特性，包括固有頻率和振型。由于ANSYS隨機振動分析采用模態(tài)疊加法，因此隨機振動分析前必須進行模態(tài)分析。本文中靜力學分析設置的邊界和輸入條件使結構中存在預應力，會導致結構剛度變化進而影響模態(tài)分析及后續(xù)振動分析結果。因此，在進行預應力模態(tài)分析時需考慮機匣載荷和溫度對結構剛度的影響。

設置PSD譜最高頻率為2000 Hz，進行15階模態(tài)分析，固有頻率結果皆小于500 Hz。因此，為保證響應結果的精確性分析，在前45階模態(tài)中，第1、6、15、45階固有頻率結果見表5。

1.4 隨機振動響應輸出

根據(jù)仿真計算的薄弱位置的應力響應PSD譜進行壽命分析。根據(jù)隨機振動理論，響應PSD譜可以由單自由度系統(tǒng)的傳遞函數(shù)H（ω）根據(jù)輸入（激勵）PSD譜通過模態(tài)疊加方法計算得到。第i個自由度的響應PSD譜由以下3部分組成：

動態(tài)部分

擬靜態(tài)部分

協(xié)方差部分

式中：n為振型數(shù)；r1、r2分別為節(jié)點數(shù)和基礎PSD譜數(shù)。

單自由度的傳遞函數(shù)根據(jù)不同的輸入PSD譜和響應類型有不同形式。本文中輸入為加速度PSD譜，輸出為應力PSD譜，其傳遞函數(shù)形式為

式中：ω為外激勵頻率；ωj為第j階模態(tài)的自然圓弧頻率。

采用ANSYS Workbench得到最大等效應力結果見表6。

表6 最大等效應力計算結果

從表6中可知，振動譜沿Z軸方向等效應力最大，如圖5所示。

圖5 等效應力

從圖中可見附件機匣殼體的應力最大點（即薄弱位置），從而能夠輸出薄弱位置的正應力σx、σy、σz和剪應力 τxy、τyz、τzz的響應 PSD 譜。

應力響應PSD譜如圖6（a）所示。已知功率譜密度值-頻率值關系曲線下面積的開方即為均方根值（RMS）[6]。從圖 6（a）中可見，σx響應 PSD譜的 RMS皆遠大于其余應力，并且各應力頻率主要集中于90～150 Hz，將該頻率區(qū)間圖像放大，如如圖 6（b）所示，發(fā)現(xiàn)σz和σyz的RMS遠小于其余應力。

圖6 應力響應PSD譜

2 壽命建模及分析

2.1 雨流循環(huán)計數(shù)方法

根據(jù)第1.4節(jié)輸出的響應PSD譜，采用雨流循環(huán)計數(shù)法（RFC）計算壽命。RFC是循環(huán)計數(shù)法的1種。循環(huán)技術法指在已知峰值概率密度函數(shù)的情況下，對應力峰值進行循環(huán)計數(shù)得到幅值信息，然后進行壽命計算。常用的循環(huán)計數(shù)法還有變程計數(shù)法（RC）和水平穿越計數(shù)法（NB）[7]。

雨流計數(shù)法是公認最好的疲勞損傷估計方法，但其計算過程比較繁瑣，準確解析式很難給出[8-9]。然而，在任意平穩(wěn)高斯過程中，線性損傷準則有效的情況下，雨流計數(shù)損傷總是處于2種損傷值之間[10-11]

式中為雨流計數(shù)損傷；DRC為變程計數(shù)損傷；DNB為水平穿越計數(shù)損傷。

從式（6）中可知，估計雨流計數(shù)損傷可以通過在其上下界中找到1個合適的中間點。因此，D.Benasciutti和R.Tovo提出雨流計數(shù)損傷可以由和加權線性組合得到[12]

變程計數(shù)法和水平穿越計數(shù)法[13]計算疲勞損傷的近似公式為

式中：υ0為譜密度函數(shù)Sx（ω）的平均上跨率

式中：λm為單邊譜密度函數(shù)Sx（ω）的譜距

方差由 λ0決定

C和k為S-N曲線中的參數(shù)

權重b的選擇取決于響應的譜密度函數(shù)。D.Benasciutti和R.Tovo給出b的表達式為

其中

對于譜密度函數(shù) Sx（ω），α1/α2為其帶寬參數(shù)，取值為[0,1][14]。

2.2 線性疲勞累計損傷理論

根據(jù)線性疲勞累積損傷理論，計算附件機匣的總損傷量從而給出其疲勞壽命[15]。

式中：D為總損傷量；Di為單個循環(huán)造成的損傷。

線性疲勞累計損傷理論假定損傷量D=1時試件將發(fā)生疲勞破壞。因此，疲勞壽命為

2.3 疲勞壽命計算

根據(jù)得到的響應PSD譜（圖6），運用雨流循環(huán)計數(shù)法并采用MATLAB編程計算附件機匣的疲勞壽命。由于σx響應PSD譜RMS遠大于其余應力，根據(jù)σx響應PSD譜計算附件機匣的損傷以及疲勞壽命，結果見表7。

根據(jù)6個應力的響應PSD譜計算其單位時間內的損傷，再疊加計算所有功率譜造成的單位時間損傷，從而得到附件機匣總損傷以及壽命結果，見表8。

表7 損傷及疲勞壽命計算結果

從表7、8中可知，2種計算方式的結果非常接近，表明 RMS非常小的 σy、σz，τxy、τyz、τzz對附件機匣的損傷以及壽命計算結果影響很小。

3 結論

本文綜合考慮附件機匣復雜的工作環(huán)境，包括自身重力和固定約束條件、軸承載荷、溫度場以及振動載荷，分析計算附件機匣殼體在機械載荷、溫度并施加隨機振動譜得到的應力響應PSD譜，結合雨流循環(huán)計數(shù)方法計算附件機匣殼體的疲勞壽命。在壽命計算時，對比單獨采用RMS最大的應力PSD譜和采用所有應力PSD譜計算得到的壽命結果分別為3.25828×106min和3.25823×106min，相差僅為50 min。因此，可以選用RMS最大的PSD譜計算附件機匣的疲勞壽命。

[1]呂亞國，劉振俠，路彬，等.航空發(fā)動機附件機匣熱分析研究[J].潤滑與密封,2011,36（10）：62-80.LU Yaguo,LIU Zhenxia,LU Bin,et al.Thermal analysis of aeroengine accessory gearbox[J].Lubrication Engineering,2011，36（10）：62-80.(in Chinese)

[2]王桂華，劉海年.航空發(fā)動機成附件振動環(huán)境試驗剖面確定方法研究[J].推進技術，2013，34（8）：1101-1107.WANG Guihua,LIU Hainian.Study on formulating method for vibration environment test profiles of aeroengine accessories[J].Journal of Propulsion Technology，2013，34（8）：1101-1107.(in Chinese)

[3]李錦花，史妍妍，張茂強，等.航空發(fā)動機附件機匣殼體變形分析[J].航空發(fā)動機，2013，39（3）：59-72.LI Jinhua,SHI Yanyan,ZHANG Maoqiang,et al.Analysis of accessory gearbox housing distortion for aeroengine[J].Aeroengine，2013，39（3）：59-72.(in Chinese)

[4]郭梅，陳聰慧，王建軍，等.發(fā)動機附件機匣結構系統(tǒng)振動特性[J].航空動力學報，2013，28（7）：1607-1612.GUO Mei,CHEN Chonghui,WANG Jianjun.Vibration characteristics of accessory gearbox structure system of engine [J].Journal of Aerospace Power，2013，28（7）：1607-1612.(in Chinese)

[5]吳鴻，李國權，齊樂華.基于ANSYS的附件機匣抽油口潤滑油溫度計算[J].航空發(fā)動機，2006，32（3）：31-35.WU Hong,LI Guoquan,QI Lehua.Lubrication temperature calculation of the oil outlet for accessory gearbox based on ANSYS[J].Aeroengine，2006，32（3）：31-35.(in Chinese)

[6]鄭志國，王宇峰.隨機振動中的參數(shù)介紹及計算方法[J].電子產品可靠性與環(huán)境試驗，2009，26（7）：45-48.ZHENG Zhiguo,WANG Yufeng.Introduction of parameters in random vibration and their calculation[J].Electronic Product Reliability and Environmental Testing，2009，26（7）：45-48.(in Chinese)

[7]王明珠.結構疲勞振動壽命分析方法研究[D].南京:南京航空航天大學,2009.WANG Mingzhu.Research on life analysis method for structure vibration fatigue[D].Nanjing：Nanjing University of Aeronautics and Astronautics,2009.

[8]Petrucci G，Zuccarello B.On the estimation of the fatigue cycle distribution from spectral density data [J].Proceedings of the Institution of Mechanical Engineers,Part C:Journal of Mechanical Engineering Science，1999，213（8）：819-831.

[9]Petrucci G,Di Paola M,Zuccarello B.On the characterization of dynamic properties of random processes by spectral parameters[J].Journal of Applied Mechanics，2000，67（3）：519-526.

[10]Rychlik I.Note on cycle counts in irregular loads[J].Fatigue and Fracture and Engineering Materials and Structures，1993，16（4）：377-390.

[11]Frendahl M,Rychlik I.Rainflow analysis:Markov method[J].International Journal of Fatigue，1993，15（4）：265-272.

[12]Benasciutti D,Tovo R.Comparison of spectral methods for fatigue analysis of broad-band Gaussian random processes[J].Probabilistic Engineering Mechanics，2006，21（4）：287-299.

[13]Rychlik I.On the narrow-band approximation for expected fatigue damage[J].Probabilistic Engineering Mechanics，1993，8（1）：1-4.

[14]Lutes L D,Sarkani S.Stochastic analysis of structural and mechanical vibrations[M].Oxford：Butterworth-Heinemann,2003：261-272.

[15]姚起杭，姚軍.工程結構的振動疲勞問題[J].應用力學學報，2006，23（1）：12-15.YAO Qihang,YAO Jun.Vibration fatigue in engineering structures[J].Chinese Journal of Applied Mechanics，2006，23（1）：12-15.(in Chinese)