魏大盛，王延榮

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

1 引言

對許多物理和工程問題的描述都涉及接觸現(xiàn)象，以航空發(fā)動機(jī)為例，包括裝配時(shí)的過盈配合、傳動系統(tǒng)中的齒輪傳動、葉盤榫連結(jié)構(gòu)等，接觸的物體在接觸界面上的相互作用較為復(fù)雜，是發(fā)生損傷和破壞的主要原因。

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

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

2 赫茲接觸問題

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

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

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

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

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

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

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

3 樅樹形榫連結(jié)構(gòu)的有限元分析

3.1 接觸應(yīng)力分析

在赫茲接觸問題計(jì)算分析的基礎(chǔ)上，開展了某發(fā)動機(jī)真實(shí)葉盤榫連結(jié)構(gòu)的接觸應(yīng)力分析。其結(jié)構(gòu)如圖5所示，共有三對接觸齒對，接觸面角度約為40°，三對齒的接觸區(qū)長度均接近1.8 mm，接觸區(qū)邊緣圓角半徑約0.675 mm。

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

圖5 樅樹形榫連結(jié)構(gòu)幾何模型Fig.5 Geometrical model of fir tree attachment

以往相關(guān)分析中接觸區(qū)網(wǎng)格稀疏，無法反映細(xì)節(jié)特征，因此本文的有限元模型采用了較密的網(wǎng)格密度。三對齒的接觸區(qū)網(wǎng)格均初步劃分為400個(gè)，邊界層采用規(guī)則網(wǎng)格劃分，內(nèi)部采用自由網(wǎng)格劃分。有限元模型見圖6，單元43 521個(gè)，節(jié)點(diǎn)44 995個(gè)。

圖6 樅樹形榫連結(jié)構(gòu)有限元模型Fig.6 FE model of fir tree attachment

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

3.2 接觸應(yīng)力計(jì)算結(jié)果的收斂性分析

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

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

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

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

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

圖10中給出了網(wǎng)格對收斂的影響趨勢：對于單齒燕尾形榫連結(jié)構(gòu)，沿接觸區(qū)劃分450個(gè)網(wǎng)格可保證計(jì)算收斂；對于三齒樅樹形榫連結(jié)構(gòu)，沿接觸區(qū)劃分1 000個(gè)網(wǎng)格可保證計(jì)算收斂。當(dāng)然，對于不同的計(jì)算模型，網(wǎng)格劃分的密度應(yīng)同幾何參數(shù)(特別是接觸區(qū)長度及接觸區(qū)邊緣圓角半徑)、受力狀態(tài)密切相關(guān)。

3.3 接觸區(qū)計(jì)算網(wǎng)格的改進(jìn)

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

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

表3 不同網(wǎng)格密度模型第三齒應(yīng)力結(jié)果MPaTable 3 Contact stress from different grids on the third tooth

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

圖9 不同模型第三齒接觸區(qū)應(yīng)力σ22分布Fig.9 Distribution of σ22on the third tooth

圖10 計(jì)算收斂性分析Fig.10 Convergence analysis

圖11 改進(jìn)后的計(jì)算網(wǎng)格Fig.11 Improved grids

采用改進(jìn)后的網(wǎng)格進(jìn)行分析，接觸應(yīng)力計(jì)算結(jié)果差別不大，但計(jì)算時(shí)間顯著減少，這對于后續(xù)的三維榫連結(jié)構(gòu)的計(jì)算分析至關(guān)重要。

4 結(jié)論

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

(1)榫連結(jié)構(gòu)接觸區(qū)邊緣存在顯著的應(yīng)力梯度，這是產(chǎn)生微動疲勞磨損的重要因素。

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

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

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

[1]何明鑒.機(jī)械構(gòu)件的微動疲勞[M].北京：國防工業(yè)出版社，1994.

[2]古遠(yuǎn)興.高低周復(fù)合載荷下燕尾榫結(jié)構(gòu)微動疲勞壽命研究[D].南京：南京航空航天大學(xué)，2007.

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

[4]魏大盛，王延榮.榫連結(jié)構(gòu)接觸面幾何構(gòu)形對接觸區(qū)應(yīng)力分布的影響[J].航空動力學(xué)報(bào)，2010，25(2)：407—411.

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

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

[7]Johson K L.接觸力學(xué)[M].徐秉業(yè)，羅學(xué)富，譯.北京：高等教育出版社，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.