弧齒錐齒輪有限元建模與接觸分析

侯祥穎，方宗德，鄧效忠，蔡香偉，寧程豐

（1.西北工業大學機電學院，陜西西安710072；2.河南科技大學機電工程學院，河南洛陽471023）

弧齒錐齒輪有限元建模與接觸分析

侯祥穎1，方宗德1，鄧效忠2，蔡香偉1，寧程豐1

（1.西北工業大學機電學院，陜西西安710072；2.河南科技大學機電工程學院，河南洛陽471023）

為了使有限元軟件能夠對齒面微米量級變化的齒輪進行分析，提出了基于高精度建模的有限元分析方法。將齒面的理論計算公式轉化成數值計算算法，按照有限元理論和有限元軟件的編程語法將齒輪實體進行節點坐標求解和網格單元劃分，可以求得一對裝配好的有限元齒輪模型。高精度建模的高重合度弧齒錐齒輪算例中，將有限元計算結果與實驗結果對比，發現兩者的一致性：彎曲應力曲線形狀一致，變化規律相近；接觸印痕區域一致，變化趨勢相同。高精度建模有限元分析能夠清晰的反映齒面微米量級改變導致的印痕和應力的變化。

弧齒錐齒輪；有限元；接觸分析；彎曲應力；接觸印痕

弧齒錐齒輪是一種不完全共軛局部點接觸的齒輪副，齒面形狀和加工機床的調整都比較復雜。目前國外已經較多地將弧齒錐齒輪應用于船艦傳動系統，國內在船艦上的應用相對較少，這方面的研究正在進行。弧齒錐齒輪傳動平穩，是動力傳輸的關鍵部件，在強度、動態性能與可靠性方面等方面都有很高的要求。鄧效忠等［1］曾提出高重合度弧齒錐齒輪的設計方法，采用LTCA理論，通過有限元應力影響矩陣法對其相關性能進行了詳細討論［2］。高重合度弧齒錐齒輪是通過齒面修形改善其嚙合性能的典型，不同的弧齒錐齒輪齒面精度會使弧齒錐齒輪有不同的傳動性能。本文以一對高重合度弧齒錐齒輪為例，通過高精度齒面建模和裝配，導入有限元分析軟件分析，并與實驗結果進行了對比驗證分析。

1 高重合度弧齒錐齒輪

文獻［3?5］等利用NURBS曲面函數對齒面進行重構來建立弧齒錐齒輪模型。但是此方法計算求得模型的齒面與理論齒面仍有一定誤差；在對一對齒輪進行有限元接觸分析時，采用手動裝配，可能存在干涉和裝配誤差等現象，產生計算誤差或錯誤。對于齒面修形處于微米量級的齒輪來說，上述方法建模精度遠遠達不到要求。采用曲面擬合方法的誤差分布如圖1所示，誤差可能達到十幾甚至幾十微米，且有很大隨機性，分布不均勻；而且此方法需要進行三維建模，然后再將模型導入有限元軟件進行實體分割、網格劃分等操作，工作量大，可重復性差。

在保證齒輪基本參數（如表1所示）不變的條件下，采用局部綜合法使接觸路徑傾斜和改變齒面的失配量，設計了3種重合度不同的弧齒錐齒輪。（詳見文獻［2］）高重合度設計方法不改變齒輪輪坯基本參數，只對小輪齒面加工參數進行調整，而保持大輪的加工參數不變，通過小輪齒面的改變使接觸路徑傾斜，接觸路徑延長，以達到增大重合度的目的。小輪齒面的修形量往往只有十幾微米甚至幾微米，故要求有限元模型要有極高的精度才能反映出齒面變化的影響。

圖1 NURBS曲面重構誤差分析Fig.1 Error analysis of NURBS surface reconstruction

表1 弧齒錐齒輪參數Table 1 Basic parameters of spiral bevel gear

本文算例采取的傾斜角分別為20°、50°、80°接觸路徑，通過TCA分析得到傳動誤差如圖2所示，圖2（a）理論設計重合度只有1.25左右，圖2（b）理論設計重合度增加至1.75左右，圖2（c）理論設計重合度增加至2.0以上，傳動誤差曲線明顯可以看出存在三齒嚙合區。

圖2 3種不同齒輪的傳動誤差Fig.2 Transmission error of three different designs

2 弧齒錐齒輪有限元模型建立

2.1 通過數值計算求解有限元模型

弧齒錐齒輪齒面較為為復雜，其加工原理和方法繁瑣，加工參數較多，一般很難得到齒面方程的解析表達式，多采用數值計算的方法，求得齒面上離散的點，即可用離散點構建有限元模型［4?5］。

將刀具切削面方程從刀具坐標系變換到輪坯坐標系下，聯立齒面的位矢方程和嚙合方程如下（文獻［6］）：

式中：M為刀具坐標系到輪坯坐標系的變換矩陣，r為刀盤方程，μ、θ分別為切削刃產生圓錐曲面坐標，φ為轉角，n、v分別為齒面點法矢和相對速度矢。由式（1）通過編程進行數值計算，可得出齒面上的離散點在輪坯坐標系下的坐標。

按照有限元軟件ABAQUS對單元網格以及節點的書寫語法要求，首先計算理論齒面上的點，將齒面點坐標繞軸線旋轉，對兩齒面之間采用點陣進行填充，然后對所有節點按照一定順序進行自然數編號；之后按順序編寫8節點單元，單元由8個節點自然序號構成，故點陣的密度決定了網格單元的密度。

裝配時，選取大輪接觸跡線與齒頂線的交點為嚙合參考點，即恰好進入嚙合的位置。利用輪齒接觸分析求得嚙合點的坐標以及小、大輪對應的轉角φ1、φ2（式（2）），采用坐標轉換矩陣，對大小輪分別旋轉φ1、φ2，即能通過理論計算完成一對弧齒錐齒輪的裝配。

式中：下標1、2分別表示大小齒輪，下標p、g分別表示大小輪的齒面參數，上標h表示在固定坐標系Sh下，分別表示齒面離散點位矢在旋轉投影面上的X、Y坐標，ha2表示大輪齒頂高，AFO表示外錐距，γ2為大輪根錐角。

通過編寫inp文件將裝配坐標系下大小輪的節點和單元數據導入ABAQUS軟件中，可直接得到有限元裝配模型（如圖3所示）。

圖3 有限元裝配模型Fig.3 Assembly model of finite element

2.2 有限元模型前處理

進行有限元齒輪接觸分析，單元一般采用六面體一次減縮積分單元（C3D8R），單元屬性在編寫inp文件生成有限元模型時定義（文獻［7］）；將已建好的弧齒錐齒輪的五齒有限元模型導入ABAQUS軟件中，需要進行以下設置（文獻［5，8］）：

1）定義材料：將大小輪賦予相同的材料屬性，其彈性模量為2.06×105MPa，泊松比為0.3，密度為7.8×10-9t／mm3。

2）定義接觸對：分別將大小輪的齒面、過渡曲面、齒頂等形成的連續面選中，并將此兩對連續面設置為接觸對。在接觸屬性中，采用無摩擦硬接觸，計算方法采用Kinematic接觸方法。

3）設置參考點：取軸線上任一點為參考點，并在參考點和小輪內圈和剖面間建立耦合剛體約束。

4）分析類型的定義：采用顯示動力學分析算法，設置分析時間并根據網格數量選擇適當的質量放大系數。顯式算法基于動力學方程，無需迭代，有較好的穩定性。

5）場變量和歷史變量設置：設置接觸壓強、接觸力、應力、位移等為場變量輸出量，動能、內能為歷史變量輸出量。

6）邊界條件和載荷：大小輪分別只釋放繞自身軸線旋轉的自由度，施加小輪的繞軸線的轉速為0.2 rad／s，大輪施加阻力矩為7.5×105N·mm，小輪轉速和大輪轉矩均采用平緩加載方式；從而模擬出以小輪驅動大輪，同時小輪受到大輪的阻力矩的作用，在齒面上產生壓力的工況。

同樣將以上設置步驟編寫入inp文件中，這樣inp文件就包括了所有模型信息和條件設置信息。直接將inp文件導入ABAQUS中，提交計算即可，無需再人工操作。

3 結果分析

3.1 有限元仿真結果

前處理結束后，用ABAQUS軟件分析可得到弧齒錐齒輪的有限元分析結果。圖4為某一時刻的小輪凹面和大輪凸面（工作面）的Mises應力云圖，圖中可以看到3對齒輪同時嚙合的狀態。

圖4 有限元計算Misse應力云圖Fig.4 Mises stress nephogram of finite element calculation

改變運動時間，當給小輪施加逆時針轉速時，小輪的凹面和大輪凸面發生接觸，小輪驅動大輪順時針旋轉，同時大輪的阻力矩模擬負載，使兩齒嚙合點處產生應力和變形。小輪凹面接觸橢圓從大端齒根逐漸向小端齒頂移動，大輪凸面接觸區域接觸橢圓從大端齒頂逐漸向小端齒根移動；最大主應力位于接觸橢圓區域的中心，接觸橢圓在齒面中間面積最大，表示在齒面中間輪齒承受載荷最大。

3.2 高重合度對齒面接觸力的影響

有限元計算結束生成結果文件，對運行編譯好的腳本文件，從計算結果模型中提取所需要的節點坐標，接觸壓強和Mises應力等，并通過數據處理軟件將應力曲線繪出，如圖5所示。

圖5 3種接觸路徑的接觸應力對比Fig.5 Contrast of contact stress with three different designs

通過對比容易發現，在計算載荷一定的情況下（大輪扭矩7.5×105N·mm），接觸路徑傾斜度越大，即傾斜角越小，則重合度越大，接觸應力越小。傾斜角80°設計中接觸應力最大值為1 111.84 MPa，而傾斜角20°設計中的最大值為1 164.98 MPa；接觸路徑傾斜前后相比，接觸應力下降4.56%，并且，隨著接觸路徑的傾斜，重合度的增大，齒面接觸應力變化趨于平緩，曲線上升下降都明顯變化緩慢，即加載明顯趨于平穩。

4 彎曲應力實驗對比

本實驗由調頻電機（30 kW）提供實驗驅動力，經減速器減速，利用轉矩轉速傳感器將動力輸給試驗齒輪小輪；試驗齒輪由磁粉制動器施加負荷，裝配在試驗箱內。通過測量齒根應變獲得齒根彎曲應力，通過轉矩轉速測試儀測得試驗齒輪的轉矩和轉速。電阻應變計貼在大輪凸面的齒根中點處，位于非嚙合區，在遠離嚙合齒面處的齒坯上貼2個溫度補償片，分別與凸面和凹面的應變計組成半橋連接。引線由集流環引出，經應變儀放大，通過多功能信號測試儀紀錄［9］。通過實驗測得齒根應力測試曲線如圖6所示。圖7為對應齒輪的有限元仿真計算所得到的大輪齒根中點彎曲應力過程，詳細數值比較如表2。通過圖6與圖7以及表2對比，大輪齒根中點應力變化過程中，實驗測得波形傾斜角20°設計最大電壓信號達1.053 V，而傾斜角80°設計最大電壓信號0.597 V，二者相比峰值下降43.3%；有限元分析結果中，傾斜角20°設計最大彎曲應力為126.08 MPa，而傾斜角80°設計最大彎曲應力值僅為83.28 MPa；接觸路徑傾斜前后相比，彎曲應力下降34%，且當接觸路徑越傾斜時，即重合度增大，彎曲應力數值明顯下降。隨著接觸路徑傾斜角加大，接觸路徑延伸增長，齒根彎曲應力逐漸降低，變化趨于平緩，仿真結果與實驗結果曲線變化規律一致，形狀近似。實驗數據與有限元分析結果有一定的誤差，誤差的主要來源有：實驗儀器條件，實驗系統誤差，應變計所貼位置，溫度補償片的選用；有限元所選用的單元類型、網格密度、計算方法、質量放大系數等。

圖6 三對齒輪齒根應力測定曲線Fig.6 Bending stress curves of experiment results

圖7 有限元大輪齒根中點彎曲應力過程Fig.7 Bending stress changes of midpoint on gear root

表2 大輪齒根應力實驗數據與有限元分析比較Table 2 Bending stress comparison between experi?ment results and finite element analysis

5 接觸印痕實驗對比

有限元分析結束后，編寫讀取abaqus結果文件（.odb文件）的腳本文件，求得某齒一個嚙合周期內每一瞬時的接觸區域內的點坐標，利用數據處理軟件編程，將所有點疊加到同一坐標系下，并將其投影到旋轉投影面，即得到有限元分析的接觸印痕區域。有限元計算和實驗得到接觸印痕區域分別如圖8和圖9所示。

圖8 有限元印痕提取Fig.8 Contact pattern of finite element analysis

圖9 實驗印痕Fig.9 Contact pattern of experiment

通過對比圖8和圖9，隨著傾斜度逐漸變大，接觸區域傾斜，接觸路徑增長，接觸橢圓長軸變短，有限元計算結果與實驗結果印痕區域一致，變化規律相同，由此證明此高精度建模方法及有限元分析方法的正確性。

6 結論

1）本文有限元建模采用數值計算的方法，得到節點坐標，并通過編號導入inp文件中直接建立有限元模型，省去了三維建模以及劃分網格的過程，通過與NURBS曲面重構結果對比，證明此方法有限元齒面精度及裝配精度極高，可重復性好，適用于高精度的齒輪接觸分析。

2）通過有限元計算與實驗結果對比，隨著接觸路徑傾斜度逐漸變大，有限元計算印痕接觸路徑逐漸增長，接觸橢圓長軸逐漸變短，印痕區域和變化規律與實驗一致；有限元計算齒根彎曲應力逐漸降低，變化趨于平緩，與實驗彎曲應力變化趨勢一致，形狀近似，驗證了本文方法的計算精度。

［1］鄧效忠，方宗德，張金良等.弧齒錐齒輪的高重合度設計［J］.中國機械工程，2002，14（19）：791?795.DENG Xiaozhong，FANG Zongde，ZHANG Jinliang，et al.Design of spiral bevel gears with high contact ratio［J］.Chi?na Mechanical Engineering，2002，14（19）：791?795.

［2］鄧效忠，方宗德，魏冰陽.高重合度弧齒錐齒輪的接觸應力與彎曲應力分析［C］／／中國工程機械學會2003年年會論文集，2003：107?111.DENG Xiaozhong，FANG Zongde，WEI Bingyang.Analysis of contact stress and bending stress on spiral bevel gears of high contact ratio［C］／／China Mechanical Engineering Soci?ety 2003 Annual Meeting Proceedings，2003：107?111.

［3］唐進元，蒲太平.弧齒錐齒輪動態嚙合有限元數值分析［J］.機械科學與技術，2009，8（8）：981?985.TANG Jinyuan，PU Taiping.FEM numerical analysis of the dynamic engagement of spiral bevel gears［J］.Mechanical Science and Technology for Aerospace Engineering，2009，8（8）：981?985.

［4］ASTOUL J，GENEIX J，MERMOZ E，et al.A simple and robust method for spiral bevel gear generation and tooth con?tact analysis［J］.International Journal on Interactive Design and Manufacturing，2013，7（1）：37?49.

［5］唐進元，蒲太平.基于有限元法的螺旋錐齒輪嚙合剛度計算［J］.機械工程學報，2011，47（11）：23?29.TANG Jinyuan，PU Taiping.Spiral bevel gear meshing stiff?ness calculations based on the finiteelement method［J］.Journal of Mechanical Engineering，2011，47（11）：23?29.

［6］LITVIN F L.Gear geometry and applied theory［M］.Cam?bridge University Press，2004：613?635.

［7］GONZALEZ?PEREZ I，ISERTE J L，FUENTES A.Imple?mentation of Hertz theory and validation of a finite element model for stress analysis of gear drives with localized bearing contact［J］.Mechanism and Machine Theory，2011，46（6）：765?783.

［8］SHEVELEVA G I，VOLKOV A E，MEDYEDEV V I.Algo?rithms for analysis of meshing and contact of spiral bevel gear［J］.Mechanism and Machine Theory，2007，42（2）：198?215.

［9］鄧效忠.高重合度弧齒錐齒輪的設計理論及實驗研究［D］.西安：西北工業大學，2002：49?77.DENG Xiaozhong.Research of design theory and experi?ments on spiral bevel gears with high contact ratio［D］.Xi'an：Northwestern Polytechnical University，2002：49?77.

［10］WANG P Y，FAN S C，HUANG Z G.Spiral bevel gear dy?namic contact and tooth impact analysis［J］.Journal of Me?chanical Design，2011，133：084501.

［11］GONZALEZ?PEREZ I，RODA?CASANOVA V，FUENTES A，et al.A finite element model for consideration of the tor?sional effect on the bearing contact of gear drives［J］.Jour?nal of Mechanical Design，2012，134：071007.

Contact analysis of spiral bevel gears based on finite element model

HOU Xiangying1，FANG Zongde1，DENG Xiaozhong2，CAI Xiangwei1，NING Chengfeng1

（1.Department of Mechanical Engineering，Northwestern Polytechnical University，Xi'an 710072，China；2.Department of Mechani?cal Engineering，Henan University of Science and Technology，Luoyang 471023，China）

To make analysis on gears whose tooth surfaces have micron level changes，the finite element method based on high?precision finite element model is proposed.The theoretical formulas of tooth surfaces are converted into numerical algorithm，the coordinates of nodes and finite mesh elements of the gears are obtained based on the finite element theory and program syntax of finite element software，and the high?precision finite element model of a pair of gears is built.Through analysis of spiral bevel gears with high contact ratio，the simulation results are in good agreement with the experimental ones.The bending stress curves show the same shape and accordant changing trend，but the contact patterns show the same location and changing trend.The finite element analysis with high?precision modeling could reflect the changes of stresses and contact pattern caused by modification to the surfaces of gears in micron level.

spiral bevel gear；finite element method；contact analysis；bending stress；contact pattern

10.3969／j.issn.1006?7043.201312055

TH132.41

：A

：1006?7043（2015）06?0826?06

http：／／www.cnki.net／kcms／detail／23.1390.u.20150428.0910.012.html

2013?12?27.網絡出版時間：2015?04?28.

國家自然科學基金資助項目（51175423，51375384，51205310）.

侯祥穎（1990?），男，博士生；方宗德（1948?），男，教授，博士生導師.

侯祥穎，E?mail：houxiangying＠126.com.