劉光磊，趙 寧，張瑞庭

(西北工業(yè)大學(xué)機電學(xué)院，陜西 西安 710072)

0 前言

在進行輪齒加載接觸分析時，如何獲得初始嚙合點處接觸橢圓長軸上離散點的柔度矩陣是關(guān)鍵技術(shù)之一。許多學(xué)者對此進行了研究［1-7］，所提出的方法大致可以分為兩種。一種是直接法［1-5］，即直接采用有限元法進行分析;另一種是間接法，即先采用有限元法形成齒面網(wǎng)格柔度矩陣，經(jīng)轉(zhuǎn)換，得到輪齒接觸橢圓長軸上離散點的柔度矩陣［6-7］。直接法計算結(jié)果精確，但計算過程復(fù)雜，花費時間長。間接法計算精度稍低，但計算速度快、適應(yīng)面廣。直接法已經(jīng)成熟，間接法尚有待改進。

間接法在操作中，文獻［6-7］提出了從齒面網(wǎng)格柔度矩陣向齒面接觸點柔度矩陣轉(zhuǎn)換的二元插值法，但并未給出具體插值的細節(jié)。為了實現(xiàn)從齒面網(wǎng)格柔度矩陣到齒面接觸點柔度矩陣的轉(zhuǎn)換，本文以有限元法中的數(shù)值計算原理為基礎(chǔ)，提出從齒面網(wǎng)格柔度矩陣到齒面接觸點柔度矩陣的混合插值法。在求解中發(fā)現(xiàn)，在齒面網(wǎng)格結(jié)點施加單位載荷時，相鄰齒面結(jié)點出現(xiàn)齒面畸變。這與彈性力學(xué)中的布辛涅斯克的計算結(jié)果相左。在網(wǎng)格單元包含兩個以上的接觸點時，從網(wǎng)格結(jié)點向著它所包圍的接觸點所進行的位移插值，有可能出現(xiàn)接觸點柔度矩陣主對角線不占優(yōu)的情況。為了解決這一問題，提出了齒面網(wǎng)格結(jié)點間距與接觸點間距必須協(xié)調(diào)，須保證一個網(wǎng)格結(jié)點單元只包含一個接觸點。所提出的方法，物理概念清楚，可操作性強。

1 混合插值方法

1.1 弧齒錐齒輪有限元分析模型

本文采用弧齒錐齒輪單齒模型，齒輪副的工作面為小輪凹面和大輪凸面。邊界條件為固定輪轂內(nèi)表面和側(cè)面，計算齒面網(wǎng)格結(jié)點的柔度矩陣時采用三維八節(jié)點非協(xié)調(diào)等參單元［8］。齒面加載接觸分析有限元網(wǎng)格模型［9-10］如圖1所示。

圖1 輪齒有限元模型Fig.1 Finite element models of gear teeth

1.2 齒面網(wǎng)格柔度矩陣和齒面接觸點柔度矩陣

齒面網(wǎng)格柔度矩陣用于反映齒面抵抗外加載荷變形的能力。依次對齒面上的每個網(wǎng)格結(jié)點施加單位法向載荷，求得每次加載時齒面所有結(jié)點沿該處法線方向的位移，也就是各結(jié)點的柔度系數(shù)，將它們按行排列，即可組成齒面網(wǎng)格柔度矩陣。圖1中的輪齒工作齒面上共有N個結(jié)點，則齒面網(wǎng)格柔度矩陣為

式中，腳標G為齒面網(wǎng)格;矩陣中的元素fij為在網(wǎng)格結(jié)點j作用單位載荷時網(wǎng)格結(jié)點i產(chǎn)生的位移。

齒面網(wǎng)格柔度矩陣反映了作用在輪齒齒面上各結(jié)點的載荷與變形關(guān)系，而在弧齒錐齒輪加載接觸分析過程中，需要用到反映接觸橢圓長軸上離散點的載荷與變形關(guān)系，即接觸點柔度矩陣。

接觸點柔度矩陣的含義與齒面網(wǎng)格柔度矩陣相同，即依次對各接觸點加載時，求得每次加載時所有接觸點沿法線方向的位移，也就是各接觸點的柔度系數(shù)，將它們按行排列，即可組成齒面接觸點柔度矩陣。

設(shè)齒面瞬時接觸橢圓長軸上的離散點有n個，按照接觸點柔度矩陣的定義，可得齒面瞬時接觸橢圓長軸上離散點的柔度矩陣為

式中，腳標C為齒面接觸點;矩陣中的元素fij為在接觸點j作用單位載荷時接觸點i產(chǎn)生的位移。

1.3 混合插值方法

從［F］G與［F］C的定義可知，盡管兩者物理含義相同，但不能通過簡單的網(wǎng)格點加密插值，實現(xiàn)從齒面網(wǎng)格柔度矩陣到接觸點柔度矩陣的轉(zhuǎn)換，而要依據(jù)柔度矩陣的力學(xué)含義，基于有限元數(shù)值計算原理，實現(xiàn)兩者的轉(zhuǎn)換。

為了在保證精度的前提下簡化計算，將齒面網(wǎng)格和接觸橢圓長軸上的離散點，轉(zhuǎn)換到齒輪的旋轉(zhuǎn)軸截面中，如圖2所示。圖中○表示網(wǎng)格結(jié)點，□表示橢圓長軸上的接觸點。具體插值分兩步完成。

第一步，載荷等效即形成包圍接觸點的網(wǎng)格結(jié)點的等效結(jié)點載荷系數(shù)。以圖2中左邊第一個接觸點為例，設(shè)想在該點施加單位法向載荷，按照虛功原理，將該載荷等效轉(zhuǎn)換到包圍它的四個網(wǎng)格結(jié)點上，即

式中，F(xiàn)C為接觸點單位法向載荷;Ni(i=1，2，3，4)為四個網(wǎng)格結(jié)點的形函數(shù)在接觸點處的值;FGi為等效載荷。因為FC為單位載荷，所以FGi也就是等效載荷系數(shù)。

圖2 齒面網(wǎng)格旋轉(zhuǎn)投影示意圖Fig.2 Schematic of rotary prejection of tooth surface grid

這里的形函數(shù)采用雙線性等參數(shù)單元，即

式中，(ξ，η)為局部坐標，如圖3所示。

圖3 單元坐標系Fig.3 Element coordinate system

與網(wǎng)格結(jié)點相關(guān)的形函數(shù)定義于局部坐標系中，而實際齒面網(wǎng)格結(jié)點和接觸點都定義于總體坐標中，如圖2和圖3所示。為了計算形函數(shù)在接觸點處的值，需要建立如下坐標變換關(guān)系:

式中，(xi，yi)為網(wǎng)格結(jié)點坐標;(x，y)為接觸點坐標。將(xi，yi)和(x，y)代入上式，解之，得到與(x，y)對應(yīng)的(ξ，η)，再帶入式(4)，即可求得Ni(i=1，2，3，4)。

第二步，位移插值即從齒面網(wǎng)格柔度矩陣出發(fā)，計算接觸點的柔度系數(shù)，形成齒面柔度矩陣。在計算時假設(shè)輪齒為線彈性、小變形，因此可以采用疊加原理。

在接觸點施加單位載荷后，各接觸點均會產(chǎn)生變形。在得到包圍施加單位載荷的接觸點的網(wǎng)格結(jié)點的等效載荷系數(shù)之后，就可以應(yīng)用疊加原理，將網(wǎng)格結(jié)點的變形等效為接觸點處的變形，即

式中，對i的求和，表示在同一網(wǎng)格結(jié)點施加等效載荷時，對接觸點的插值;對j的求和，表示包圍同一接觸點的四個網(wǎng)格結(jié)點在接觸點的插值之和。

依次對各接觸點施加單位法向載荷，采用上述載荷－位移混合插值方法，即可得到單齒接觸點的柔度矩陣。將相互嚙合的輪齒接觸點柔度矩陣相加，最終得到輪齒加載接觸分析(LTCA)所需的接觸點柔度矩陣。

2 齒面網(wǎng)格柔度矩陣的畸變

得到齒面網(wǎng)格柔度矩陣以后，圖形顯示其計算結(jié)果。圖4和圖5為大輪凸面上靠近小端處的網(wǎng)格結(jié)點受到法向單位載荷作用時，整個齒面及其局部的變形情況。

從圖5可以看出，沿著齒長方向分布的5個結(jié)點1、2、3、4、5，加載點3的變形最大，結(jié)點2和4的變形比結(jié)點1和5的變形小，但它們卻比結(jié)點1和5更靠近加載點3，即齒面出現(xiàn)畸變。齒面畸變有可能導(dǎo)致接觸點柔度矩陣的非對角元出現(xiàn)負剛度，從而影響接觸點載荷的分布。這種現(xiàn)象與彈性力學(xué)中的布辛涅斯克的計算結(jié)果相左。為此，作者將有限元模型導(dǎo)入到Ansys中，計算齒面網(wǎng)格在某點施加法向單位載荷時的齒面變形，如圖6所示。

圖6 Ansys中分析大輪輪齒模型Fig.6 Ansys analysis model for bull wheel gear tooth

計算結(jié)果如下:加載點3與圖5中加載點對應(yīng)。結(jié)點1、2、3、4、5點的法向位移見表1。

表1 Ansys齒面變形計算結(jié)果Table 1 Results of tooth distortion calculated with Ansys ×10－6

由以上數(shù)據(jù)可以看出，Ansys有限元分析結(jié)果也存在齒面畸變的現(xiàn)象。

在齒面上其他網(wǎng)格結(jié)點加載時，也發(fā)現(xiàn)類似齒面畸變現(xiàn)象。

3 接觸點柔度矩陣主對角線不占優(yōu)及解決方法

當(dāng)網(wǎng)格結(jié)點間距過大時，一個單元內(nèi)可能包含多個接觸點。如圖2所示，接觸點1、2、3同時落在一個網(wǎng)格單元內(nèi)。按照前述載荷－位移插值法，計算接觸點柔度矩陣時，因載荷與位移插值分別進行，有可能導(dǎo)致接觸點柔度矩陣主對角線元素不占優(yōu)。這不符合柔度矩陣應(yīng)有的物理意義。

以圖2中的接觸點1、2、3為例。設(shè)它們的局部坐標分別為(－0.25，－0.25)、(0.25，025)和(0.75，0.75)。當(dāng)對第二個接觸點作用法向單位載荷時，從已知的齒面網(wǎng)格結(jié)點柔度矩陣中提取該網(wǎng)格單元結(jié)點的柔度系數(shù)為

把三個接觸點的局部坐標帶入式(4)分別得到三個接觸點的形函數(shù)為

把接觸點2的形函數(shù)帶入式(3)，得到等效載荷系數(shù)為

FGi=［0.1406 0.2344 0.3906 0.2344］

將上述數(shù)據(jù)分別代入式(6)得到三個接觸點的位移值分別為

由計算結(jié)果可以看出，加載接觸點2的變形已經(jīng)小于同一網(wǎng)格單元內(nèi)的接觸點1的變形。說明柔度矩陣的主對角線元素已經(jīng)不再占優(yōu)。

要解決接觸點柔度矩陣主對角線不占優(yōu)的問題，就要密化齒面網(wǎng)格，使一個單元內(nèi)最多包含一個接觸點，如圖7所示。從而避免出現(xiàn)柔度矩陣主對角線不占優(yōu)的現(xiàn)象。

為了保證計算符合物理意義，同時兼顧計算效率，密化網(wǎng)格的尺寸應(yīng)與接觸點距相協(xié)調(diào)。以一個網(wǎng)格單元正好包含一個接觸點，且相鄰接觸點落在相鄰網(wǎng)格單元為宜。

圖7 密化后的齒面網(wǎng)格Fig.7 Tooth surface grid after densification

4 結(jié)論

(1)本文提出了一種由齒面網(wǎng)格柔度矩陣計算齒面接觸點柔度矩陣的載荷－位移插值方法。

(2)在求解齒面網(wǎng)格柔度矩陣時，發(fā)現(xiàn)了齒面網(wǎng)格結(jié)點位移畸變現(xiàn)象。

(3)采用載荷－位移插值時，可能導(dǎo)致齒面接觸點柔度矩陣的主對角線不占優(yōu)。為此，應(yīng)使齒面網(wǎng)格單元的尺寸與接觸點距離相協(xié)調(diào)，保證一個網(wǎng)格單元僅包含一個接觸點。

［1］ BIBEL G D．Contact stress analysis of spiral bevel gears using finite element analysis［J］．ASME Journal of Mechanical Design，1995，117(2):235－240．

［2］ LITVIN F L，F(xiàn)UENTES A，HAYASAKA K．Design，manufacture，stress analysis，and experimental tests of low-noise high endurance spiral bevel gears［J］．Mechanism and Machine Theory，2006，41(1):83－118．

［3］ 唐進元，蒲太平．弧齒錐齒輪動態(tài)嚙合有限元數(shù)值分析［J］．機械科學(xué)與技術(shù)，2009，28(8):981－985．

［4］ 鄭昌啟．螺旋錐齒輪輪齒加載接觸分析計算原理［J］．機械工程學(xué)報，1993，29(4):50－54．

［5］ 李潤方．運轉(zhuǎn)中嚙合輪齒的三維應(yīng)力應(yīng)變數(shù)值分析及實驗研究［J］．機械工程學(xué)報，1994，30(2):38－44．

［6］ 方宗德．齒輪輪齒承載接觸分析(LTCA)的模型和方法［J］．機械傳動，1998，22(2):1－3．

［7］ 田行斌．弧齒錐齒輪的有摩擦承載接觸分析［J］．西北工業(yè)大學(xué)學(xué)報，2000，18(1):19－22．

［8］ 劉更．結(jié)構(gòu)動力學(xué)有限元程序設(shè)計［M］．北京:國防科技大學(xué)出版社，1993．

［9］ LITVIN F L．Computerized design，simulation of meshing，and contact and stress analysis of facemilled formate generated spiral bevel gears［J］．Mechanism and Machine Theory，2002，37:441－459．

［10］ 劉光磊，樊紅衛(wèi)．含齒根過渡曲線的弧齒錐齒輪有限元精確建模［J］．機械科學(xué)與技術(shù)，2010，29(12):1595－1601．