一種齒面接觸點柔度矩陣的計算方法

劉光磊，趙 寧，張瑞庭

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

0 前言

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

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

1 混合插值方法

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

本文采用弧齒錐齒輪單齒模型，齒輪副的工作面為小輪凹面和大輪凸面。邊界條件為固定輪轂內(nèi)表面和側(cè)面，計算齒面網(wǎng)格結點的柔度矩陣時采用三維八節(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)格結點施加單位法向載荷，求得每次加載時齒面所有結點沿該處法線方向的位移，也就是各結點的柔度系數(shù)，將它們按行排列，即可組成齒面網(wǎng)格柔度矩陣。圖1中的輪齒工作齒面上共有N個結點，則齒面網(wǎng)格柔度矩陣為

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

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

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

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

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

1.3 混合插值方法

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

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

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

式中，F(xiàn)C為接觸點單位法向載荷;Ni(i=1，2，3，4)為四個網(wǎng)格結點的形函數(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)格結點相關的形函數(shù)定義于局部坐標系中，而實際齒面網(wǎng)格結點和接觸點都定義于總體坐標中，如圖2和圖3所示。為了計算形函數(shù)在接觸點處的值，需要建立如下坐標變換關系:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

FGi=［0.1406 0.2344 0.3906 0.2344］

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

由計算結果可以看出，加載接觸點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)格的尺寸應與接觸點距相協(xié)調(diào)。以一個網(wǎng)格單元正好包含一個接觸點，且相鄰接觸點落在相鄰網(wǎng)格單元為宜。

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

4 結論

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

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

(3)采用載荷－位移插值時，可能導致齒面接觸點柔度矩陣的主對角線不占優(yōu)。為此，應使齒面網(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］．機械科學與技術，2009，28(8):981－985．

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

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

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

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

［8］ 劉更．結構動力學有限元程序設計［M］．北京:國防科技大學出版社，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］．機械科學與技術，2010，29(12):1595－1601．