一種齒面接觸點(diǎn)柔度矩陣的計(jì)算方法

劉光磊，趙 寧，張瑞庭

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

0 前言

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

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

1 混合插值方法

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

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

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

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

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

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

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

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

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

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

1.3 混合插值方法

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

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

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

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

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

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

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

圖3 單元坐標(biāo)系Fig.3 Element coordinate system

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

把三個(gè)接觸點(diǎn)的局部坐標(biāo)帶入式(4)分別得到三個(gè)接觸點(diǎn)的形函數(shù)為

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

FGi=［0.1406 0.2344 0.3906 0.2344］

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

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

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

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

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

4 結(jié)論

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

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

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

［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］ 唐進(jìn)元，蒲太平．弧齒錐齒輪動(dòng)態(tài)嚙合有限元數(shù)值分析［J］．機(jī)械科學(xué)與技術(shù)，2009，28(8):981－985．

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

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

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

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

［8］ 劉更．結(jié)構(gòu)動(dòng)力學(xué)有限元程序設(shè)計(jì)［M］．北京:國(guó)防科技大學(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］．機(jī)械科學(xué)與技術(shù)，2010，29(12):1595－1601．