樅樹形榫連結構接觸應力的有限元分析及建模研究

魏大盛，王延榮

(北京航空航天大學 能源與動力工程學院，北京100083)

1 引言

對許多物理和工程問題的描述都涉及接觸現象，以航空發動機為例，包括裝配時的過盈配合、傳動系統中的齒輪傳動、葉盤榫連結構等，接觸的物體在接觸界面上的相互作用較為復雜，是發生損傷和破壞的主要原因。

航空發動機榫連結構，是典型的接觸結構，在熱、機械載荷的作用下，接觸部位通常存在高應力梯度，隨著結構使用壽命的日益增長，接觸疲勞、微動疲勞問題在工程中逐漸得到重視[1～3]。文獻[4]對燕尾形榫連結構的接觸應力進行了深入分析，準確計算了高應力梯度位置的接觸應力分布，并探討了網格密度對計算結果的影響，改善了國內以往針對榫連結構分析時計算結果精度較低的情況[5，6]。

相比于壓氣機的燕尾形榫連結構，渦輪的樅樹形榫連結構除承受離心載荷之外，還承受較高的熱載荷作用，同時結構復雜，接觸對數目增加，這都增加了接觸分析的難度。正因如此，本文首先對典型赫茲接觸問題進行計算分析，并同解析解進行對照，以明確采用有限元法處理接觸問題的流程及注意事項；其次，針對某發動機三齒樅樹形榫連結構的接觸問題展開分析，重點考查建模方式對接觸應力求解精度的影響。

2 赫茲接觸問題

接觸問題的高度非線性使其求解難度大大提高，目前有限元法是處理接觸問題的首選方法，但其主要有以下幾個難點：解不易收斂；解收斂，但合理性及精度不好評價；求解時計算量過大，尤其是三維問題。因此，在開展榫連結構接觸應力分析之前，首先采用有限元法對具有解析解的赫茲接觸問題進行分析[7]，以此明確有限元法求解接觸問題的一般流程。

以兩平行軸彈性圓柱體的接觸問題作為算例，其幾何形式及接觸體內部的應力分布見圖1。采用有限元程序MSC.MARC建立的相應的分析模型見圖2，有限元解及解析解同時在圖3、圖4中給出(圖中，formula表示解析解，fem表示有限元解)，表1則給出了幾個關鍵數值。由接觸面(Y=0)以及對稱面(X=0)的應力計算結果來看，接觸區內部應力的解析解同有限元解吻合很好，接觸表面應力的解析解同有限元解則略有差異，由數值計算誤差所致。

以上計算分析表明，只要模型建立合適，求解方法得當，有限元法在處理接觸問題時可以給出足夠高的精度。

圖1 赫茲接觸形式及應力分布Fig.1 Hertz contact and stress distribution

圖2 赫茲接觸的有限元模型Fig.2 Finite element model of hertz contact

圖3 接觸面(Y=0)上的應力分布Fig.3 Stress distribution on the contact surface

圖4 接觸體內部沿對稱面(X=0)的應力分布Fig.4 Inner stress on the symmetrical plane

3 樅樹形榫連結構的有限元分析

3.1 接觸應力分析

在赫茲接觸問題計算分析的基礎上，開展了某發動機真實葉盤榫連結構的接觸應力分析。其結構如圖5所示，共有三對接觸齒對，接觸面角度約為40°，三對齒的接觸區長度均接近1.8 mm，接觸區邊緣圓角半徑約0.675 mm。

表1 赫茲接觸的有限元解及解析解的對比Table 1 Numerical solution by FEM and analytic one

圖5 樅樹形榫連結構幾何模型Fig.5 Geometrical model of fir tree attachment

以往相關分析中接觸區網格稀疏，無法反映細節特征，因此本文的有限元模型采用了較密的網格密度。三對齒的接觸區網格均初步劃分為400個，邊界層采用規則網格劃分，內部采用自由網格劃分。有限元模型見圖6，單元43 521個，節點44 995個。

圖6 樅樹形榫連結構有限元模型Fig.6 FE model of fir tree attachment

計算時邊界條件主要考慮了離心力、溫度場以及未參與有限元建模的葉片部分的拉力，其作用于葉片上端面。樅樹形榫槽的接觸對數量較多，增加了接觸分析時收斂的難度。因此計算中采用了固定步長時間增量，增量步為100；同時載荷也采用了逐步增加的方式，其中轉速隨時間成正比增加，而拉力則與時間成非線性變化，具體對應關系可由離心力公式推導得出。計算結果表明，第三對接觸齒接觸區的下邊緣點接觸應力最大，見圖7和圖8。當然，計算時采用的幾何模型較為精確，這同實際加工時具有尺寸公差的情況不同，不同齒距偏差組合對不同榫槽/榫頭齒對接觸次序、接觸應力大小都有較大影響，因而接觸應力最大點的位置在實際情況中可能有所變化。

3.2 接觸應力計算結果的收斂性分析

由于接觸區邊緣存在高應力梯度，同時網格密度對接觸區邊緣的應力值也有較大影響[8]，因此，為精確計算接觸區邊緣應力，將接觸區網格分別細化2倍及3倍，以此考查有限元結果的收斂性。設接觸區四邊形單元邊長為LE，接觸區邊緣圓角半徑為RF，則參數LE/RF可用于度量接觸區網格密度。不同網格密度有限元模型的參數見表2。

不同網格密度有限元模型的接觸應力計算結果見表3，表4則給出了相鄰網格密度計算結果的相對誤差。據此，可根據下式對收斂性做出判斷[8]：

圖7 三對齒接觸壓力σ22比較Fig.7 Contact stress component σ22on each tooth

圖8 接觸壓力σ22分布Fig.8 Distribution of contact stress component σ22

式中n表示第n個計算模型。根據此收斂判據不斷對接觸區網格進行細化，當兩次計算結果相差不超過5%時可認為計算收斂。由表3中的計算結果可知，隨著網格密度的增加，峰值應力計算結果提高，計算逐漸收斂。網格2同網格1的計算結果相差最大9%，說明接觸區網格密度為400時的模型計算結果并未收斂，而網格3同網格2的計算結果相差在3%以下，網格2的計算結果已經達到收斂。圖9中給出了樅樹形榫槽第三齒接觸區應力計算結果的分布及收斂趨勢。

圖10中給出了網格對收斂的影響趨勢：對于單齒燕尾形榫連結構，沿接觸區劃分450個網格可保證計算收斂；對于三齒樅樹形榫連結構，沿接觸區劃分1 000個網格可保證計算收斂。當然，對于不同的計算模型，網格劃分的密度應同幾何參數(特別是接觸區長度及接觸區邊緣圓角半徑)、受力狀態密切相關。

3.3 接觸區計算網格的改進

接觸是典型的非線性問題，求解不易收斂，計算時間較長。而前面的分析表明，傳統榫連結構接觸區邊緣存在較高的應力梯度，其應力計算結果依賴于網格密度，較為稠密的網格劃分可以滿足計算精度，但同時計算時間會顯著增加。因此，有必要對前面采用的計算網格進行改進，以提高計算效率。對于本次研究中的榫連結構，其高應力梯度區域，也即需要較密網格劃分的部分，僅存在于葉/盤接觸區的邊緣部分，因此僅對該部分的網格進行細化即可。改進后的計算網格如圖11所示。

表2 不同網格密度的有限元模型Table 2 FE models of different grids

表3 不同網格密度模型第三齒應力結果MPaTable 3 Contact stress from different grids on the third tooth

表4 不同模型計算誤差分析Table 4 Error analysis of different FE model

圖9 不同模型第三齒接觸區應力σ22分布Fig.9 Distribution of σ22on the third tooth

圖10 計算收斂性分析Fig.10 Convergence analysis

圖11 改進后的計算網格Fig.11 Improved grids

采用改進后的網格進行分析，接觸應力計算結果差別不大，但計算時間顯著減少，這對于后續的三維榫連結構的計算分析至關重要。

4 結論

本文采用有限元法對經典赫茲接觸問題以及三齒樅樹形榫連結構的接觸應力展開了分析，并重點研究了網格疏密程度對數值解的影響及由此帶來的收斂性問題，得出了一些具有工程意義的結論：

(1)榫連結構接觸區邊緣存在顯著的應力梯度，這是產生微動疲勞磨損的重要因素。

(2)以往對于榫連結構接觸應力的分析，其精度明顯不夠，這無疑將使微動疲勞壽命的計算精度降低；增加高應力梯度位置的網格密度，有助于計算精度的提高。

(3)樅樹形榫槽的危險位置出現在第三齒接觸區下邊緣點位置，但此結論將受到齒距公差等現實條件的影響。

本文以接觸應力分析作為重點，研究表明有限元建模方式對接觸應力結果有較大影響，密集的網格劃分可以準確刻畫接觸區應力的細節特征，但同時也使得計算時間急劇增加，如何同時保證計算精度及效率，相應的建模方式還需深入探討。另外，榫連結構幾何參數對接觸應力的影響，微動疲勞壽命的評估，也將作為下一步的研究工作。

[1]何明鑒.機械構件的微動疲勞[M].北京：國防工業出版社，1994.

[2]古遠興.高低周復合載荷下燕尾榫結構微動疲勞壽命研究[D].南京：南京航空航天大學，2007.

[3]劉道新，劉 軍，劉元鏞.微動疲勞裂紋萌生位置及形成方式研究[J].工程力學，2007，24(3)：42—47.

[4]魏大盛，王延榮.榫連結構接觸面幾何構形對接觸區應力分布的影響[J].航空動力學報，2010，25(2)：407—411.

[5]溫衛東，高德平.榫頭/榫槽接觸問題邊界元分析[J].航空動力學報，1992，7(2)：117—120.

[6]鄭旭東，蔚奪魁，王兆豐，等.某型航空發動機渦輪葉片和輪盤榫齒裂紋故障力學分析[J].航空發動機，2005，31(3)：35—38.

[7]Johson K L.接觸力學[M].徐秉業，羅學富，譯.北京：高等教育出版社，1991.

[8]Sinclair G B，Cormier N G，Griffin J H，et al.Contact Stresses in Dovetail Attachments:Finite Element Modeling[J].Journal of Engineering for Gas Turbines and Power，2002，124：182—189.