唐進(jìn)元,彭方進(jìn)
(中南大學(xué) 現(xiàn)代復(fù)雜裝備設(shè)計(jì)與極端制造教育部重點(diǎn)實(shí)驗(yàn)室;中南大學(xué) 機(jī)電工程學(xué)院,長沙 410083)
螺旋錐齒輪動(dòng)態(tài)嚙合問題是一個(gè)邊界條件高度非線性的接觸動(dòng)力學(xué)問題,直接對(duì)其進(jìn)行解析分析難度很大,有限元方法成為了計(jì)算螺旋錐齒輪問題的最普遍最有效的方法。
結(jié)構(gòu)動(dòng)力學(xué)分析中,慣性載荷一般是指由于結(jié)構(gòu)的質(zhì)量引起的動(dòng)載荷。文獻(xiàn)表明,國內(nèi)外許多學(xué)者使用有限元方法對(duì)螺旋錐齒輪靜態(tài)或者準(zhǔn)靜態(tài)接觸特性進(jìn)行了大量的研究[1-5],但關(guān)于動(dòng)態(tài)嚙合特性的研究文獻(xiàn)很少。而進(jìn)行動(dòng)態(tài)嚙合有限元分析的必要前提之一就是建立準(zhǔn)確的有限元模型[6]。已有文獻(xiàn)中,螺旋錐齒輪的有限元模型可以概括為以下兩種:一種是以Litvin等[4-6]為代表的剛體參考點(diǎn)耦合約束模型,即在輪齒軸線上建立剛體參考點(diǎn),在參考點(diǎn)和小輪內(nèi)圈和剖面間建立耦合剛體約束。由于分析時(shí)采用的是部分齒輪模型,這種方法無法表達(dá)真實(shí)齒輪的質(zhì)量及其質(zhì)心位置,從而不能體現(xiàn)真實(shí)慣性載荷的影響,不適合進(jìn)行動(dòng)態(tài)嚙合分析;另一種則是李源等[7]采用的包含輪體的三齒嚙合模型,一方面該模型的網(wǎng)格數(shù)量多,耗費(fèi)了大量的機(jī)時(shí),另一方面,該模型沒有嚴(yán)格地反應(yīng)真實(shí)的齒輪質(zhì)量。
本文根據(jù)動(dòng)力學(xué)原理和Lagrange乘子約束,提出了考慮慣性載荷的螺旋錐齒輪動(dòng)態(tài)嚙合分析有限元模型,應(yīng)用到非線性接觸動(dòng)力學(xué)有限元分析中,并借助LS-DYNA軟件對(duì)準(zhǔn)雙曲面齒輪的動(dòng)態(tài)嚙合特性進(jìn)行了研究。
動(dòng)力學(xué)問題的基本運(yùn)動(dòng)方程為:
其中:M為質(zhì)量矩陣,P(t)為載荷矢量,K為剛度矩陣,C為阻尼矩陣,為節(jié)點(diǎn)加速度矢量,為節(jié)點(diǎn)速度矢量,{x}為節(jié)點(diǎn)位移矢量。M即為慣性載荷。
在LS-DYNA3D程序中,運(yùn)動(dòng)方程描述為:
慣性載荷是使結(jié)構(gòu)產(chǎn)生動(dòng)力響應(yīng)的本質(zhì)因素,而慣性載荷的產(chǎn)生又是由結(jié)構(gòu)的質(zhì)量引起的。動(dòng)態(tài)接觸問題屬于復(fù)雜的非線性動(dòng)力學(xué)問題,必須考慮慣性載荷的影響。因此,對(duì)其結(jié)構(gòu)的質(zhì)量位置及其運(yùn)動(dòng)的準(zhǔn)確描述成為了結(jié)構(gòu)動(dòng)力分析中的關(guān)鍵。
螺旋錐齒輪結(jié)構(gòu)復(fù)雜,大輪和小輪的齒數(shù)相差大,如果對(duì)其進(jìn)行完整一致的有限元網(wǎng)格劃分,將極大增加計(jì)算規(guī)模,尤其對(duì)于動(dòng)力學(xué)分析而言。
在有限元法中,應(yīng)用Lagrange乘子約束可以實(shí)現(xiàn)將剛性體直接粘附在柔性體上[9]。一個(gè)簡(jiǎn)單的二維剛-柔性體問題如圖1(a)所示,其中灰色為剛性體。圖1(b)表示在剛-柔界面上的一個(gè)單元之間的分解圖,這里需要強(qiáng)令單元的兩個(gè)界面節(jié)點(diǎn)的位置具有與相應(yīng)剛體上的點(diǎn)相同的變形位置。剛體的運(yùn)動(dòng)可以表示為:
其中,XI是位置,t是時(shí)間,RI是在未變形物體中的某個(gè)參考點(diǎn),ri是在變形物體中的同一個(gè)點(diǎn)的位置,ΛiI表示一個(gè)剛性轉(zhuǎn)動(dòng)。利用式(3)可以很容易寫出這種約束:
式中腳標(biāo)α代表節(jié)點(diǎn)數(shù)。應(yīng)用經(jīng)典Lagrange乘子法修改泛函使其包含該約束,經(jīng)過處理可以得到每一個(gè)剛體的方程組為:
式中:ψα和ψμ分別是在節(jié)點(diǎn)α和 μ處有限元計(jì)算的殘值,β是與節(jié)點(diǎn)α連接的任何其它剛體節(jié)點(diǎn),μ和ν是與節(jié)點(diǎn)α連接的柔性節(jié)點(diǎn),yα=xα-r是y的節(jié)點(diǎn)值。計(jì)算殘值的步驟可以在每個(gè)單元中進(jìn)行,于是,約束柔性體到剛性體的步驟仍是一個(gè)標(biāo)準(zhǔn)的有限元組裝過程,而且很容易得并入到一個(gè)求解系統(tǒng)中。
綜合考慮以上因素,針對(duì)螺旋錐齒輪動(dòng)態(tài)嚙合分析問題,本文提出一種新的有限元分析模型:建立大、小輪全齒模型,其中,對(duì)參與嚙合的部分輪齒網(wǎng)格細(xì)分,不參與嚙合的輪齒網(wǎng)格粗分,粗、細(xì)網(wǎng)格單元之間通過固連約束連接,這樣既體現(xiàn)了真實(shí)模型的質(zhì)量屬性,又在滿足分析精度的前提下盡量減小分析模型的單元規(guī)模。同時(shí),應(yīng)用Lagrange乘子約束,實(shí)現(xiàn)小輪與其內(nèi)圈和大輪與其底面之間的剛?cè)峒s束,建立考慮慣性載荷的螺旋錐齒輪動(dòng)態(tài)嚙合分析邊界加載模型,如圖2所示。其中,圖2(a)為完整的分析模型,圖2(b)為分解模型,大錐形圓環(huán)為大輪底面,小錐形圓環(huán)為小輪內(nèi)圈,大、小錐形圓環(huán)上的陰影曲面為固連約束的邊界。
圖1 剛性體和柔性體之間的Lagrange乘子約束Fig.1 Lagrange multiplier constrain between flexible and rigid bodies:(a)rigid-flexible body;(b)Lagrange multiplier
圖2 準(zhǔn)確的有限元模型Fig.2 Exact finite element model
文章通過數(shù)值計(jì)算得到齒面離散數(shù)據(jù)點(diǎn)(主要幾何參數(shù)見表1),由離散數(shù)據(jù)點(diǎn)擬合成精確的準(zhǔn)雙曲面齒輪展成曲面,并依此構(gòu)建出準(zhǔn)雙曲面齒輪三維實(shí)體模型,如圖3所示。詳細(xì)建模方法參照文獻(xiàn)[10],本文不再詳述。
表1 準(zhǔn)雙曲面齒輪幾何參數(shù)Tab.1 Geometric parameter of hypoid gear
圖3 準(zhǔn)雙曲面齒輪幾何模型Fig.3 Geometric model of hypoid gear
圖4 準(zhǔn)雙曲面齒輪有限元網(wǎng)格模型Fig.4 Finite element model of hypoid gear
有限元網(wǎng)格劃分是進(jìn)行有限元數(shù)值分析至關(guān)重要的一步,它直接影響著后續(xù)數(shù)值計(jì)算分析結(jié)果的精確性。六面體網(wǎng)格由于其求解精度高、抗變形能力強(qiáng)、單元數(shù)量少等優(yōu)點(diǎn)成為數(shù)值分析的首選網(wǎng)格。圖4為準(zhǔn)雙曲面齒輪的六面體有限元網(wǎng)格模型,其中圖4(b)為圖4(a)的局部放大。
大、小輪均采用彈性材料,彈性模量為2.1×1011Pa,泊松比為 0.3,密度為 7.8 × 103kg/m3,采用 solid164實(shí)體單元,單元積分形式為單點(diǎn)高斯積分,加上沙漏控制。大輪底面和小輪內(nèi)圈采用剛體材料,彈性模量為2.1 ×1011Pa,泊松比為0.3,采用 shell163 殼單元,殼單元屬性默認(rèn)。
定義大輪和小輪之間為自動(dòng)面面接觸,約束大、小輪除了沿其各自軸線的旋轉(zhuǎn)自由度以外的所有自由度,對(duì)小輪內(nèi)圈施加恒定轉(zhuǎn)速1200 r/min,同時(shí)對(duì)大輪底面施加恒定負(fù)載扭矩200 Nm。
基于上節(jié)所建立的有限元模型及確定的相關(guān)參數(shù),應(yīng)用LS-DYNA對(duì)螺旋錐齒輪的動(dòng)態(tài)嚙合特性進(jìn)行有限元分析。
由于小輪轉(zhuǎn)速及質(zhì)量屬性不變,可忽略其慣性載荷影響。因此,通過改變大輪的密度來改變大輪的質(zhì)量屬性,同時(shí)保持其他條件一致,以對(duì)比不同慣性載荷下準(zhǔn)雙曲面齒輪動(dòng)態(tài)嚙合的變化規(guī)律。分別建立三種對(duì)比模型,不同密度大輪與真實(shí)密度大輪的質(zhì)量比分別為:0.1∶1,1∶1,10∶1,分別用代號(hào)Ⅰ、Ⅱ和Ⅲ表示。其中,模型Ⅱ?yàn)檎鎸?shí)的齒輪模型。
為了避免強(qiáng)烈的初始嚙合沖擊,對(duì)小輪的轉(zhuǎn)速采用逐步施加,如圖5所示。
從大輪的轉(zhuǎn)速變化(圖6)可知,不同質(zhì)量比模型的大輪轉(zhuǎn)速均經(jīng)過波動(dòng)后近似于穩(wěn)定值22.56 rad/s,但仍然呈現(xiàn)小幅波動(dòng),這與螺旋錐齒輪的傳動(dòng)特性是吻合的。隨著慣性載荷的增大,大輪轉(zhuǎn)速波動(dòng)頻率減小。
圖5 模型Ⅱ的大小輪轉(zhuǎn)速對(duì)比圖Fig.5 Comparison diagram of the rotational speeds of gear and pinion in modelⅡ
圖6 不同質(zhì)量比模型的大輪轉(zhuǎn)速歷程Fig.6 The gear’s rotational speed time history of different model
由圖7、圖8和圖9可知,質(zhì)量比為0.1∶1時(shí),齒面最大接觸力為18.98 kN;質(zhì)量比為1∶1時(shí),最大接觸力為142.16 kN;質(zhì)量比為10∶1時(shí),輪齒最大沖擊接觸力為1194.20 kN。可見,隨著慣性載荷的增大,齒面接觸力波動(dòng)幅值顯著增大,波動(dòng)頻率明顯減小。
綜上,慣性載荷對(duì)螺旋錐齒輪運(yùn)動(dòng)特性、齒面接觸力波動(dòng)幅值和頻率均有顯著影響:隨著慣性載荷的增大,大輪轉(zhuǎn)速波動(dòng)頻率減小,齒面接觸力波動(dòng)的幅值增大,波動(dòng)頻率減小。
圖7 模型Ⅰ的齒面接觸力歷程Fig.7 The tooth contact force time history of modelⅠ
圖8 模型Ⅱ的齒面接觸力歷程Fig.8 The tooth contact force time history of modelⅡ
圖9 模型Ⅲ的齒面接觸力歷程Fig.8 The tooth contact force time history of modelⅢ
圖10、圖11 和圖12 分別表示質(zhì)量比為0.1∶1、1∶1和10∶1三種模型在相同時(shí)刻(0.0116 s)的小輪齒面應(yīng)力云圖。其中,圖中的H22431、H21123、H34854分別表示三種模型當(dāng)前時(shí)刻最大齒面應(yīng)力所在單元。
對(duì)比三種質(zhì)量比模型的齒面應(yīng)力云圖可知,齒面接觸區(qū)均為類橢圓形??傮w上看,隨著質(zhì)量比的增大,慣性載荷隨之增大,接觸橢圓的長軸半徑即接觸區(qū)寬度也逐漸增大。其中,質(zhì)量比為0.1∶1為單齒接觸,齒面最大接觸應(yīng)力位于齒面中間部位;1∶1模型為2齒接觸,最大接觸應(yīng)力位置偏向齒頂;而隨著慣性載荷的進(jìn)一步增大,質(zhì)量比為10∶1模型變?yōu)?齒接觸,接觸區(qū)更大,齒面最大接觸應(yīng)力出現(xiàn)在第三齒。
結(jié)合圖10、圖11和圖12,對(duì)比不同質(zhì)量比模型在相同位置(單元22311)的有效應(yīng)力歷程,分別如圖13、圖14及圖15所示。其中,質(zhì)量比為0.1∶1模型單元22311的最大應(yīng)力為204.84 MPa;質(zhì)量比為1∶1模型的最大應(yīng)力為512.18 MPa;質(zhì)量比為10∶1模型的最大應(yīng)力為2195.8 MPa??梢?,隨著質(zhì)量比即輪齒慣性載荷的增大,齒面接觸應(yīng)力增大。
圖10 模型Ⅰ的大輪齒面應(yīng)力云圖Fig.10 Pressure pattern of pinion in modelⅠ
圖11 模型Ⅱ的大輪齒面應(yīng)力云圖Fig.11 Pressure pattern of pinion in modelⅡ
圖12 模型Ⅲ的大輪齒面應(yīng)力云圖Fig.12 Pressure pattern of pinion in modelⅢ
圖13 模型Ⅰ齒面單元有效應(yīng)力歷程Fig.13 The effective stress time history of the surface element in modelⅠ
圖14 模型Ⅱ齒面單元有效應(yīng)力歷程Fig.14 The effective stress time history of the surface element in modelⅡ
圖15 模型Ⅲ齒面單元有效應(yīng)力歷程Fig.15 The effective stress time history of the surface element in modelⅢ
非線性動(dòng)力分析程序用于工程計(jì)算,最大的困難是耗費(fèi)機(jī)時(shí)過多,顯式積分的每一時(shí)步,單元計(jì)算的機(jī)時(shí)占總機(jī)時(shí)的主要部分。采用單點(diǎn)高斯積分的單元計(jì)算可以極大地節(jié)省數(shù)據(jù)存貯量和運(yùn)算次數(shù),但是單點(diǎn)積分可能引起零能模式(沙漏),需加以控制。根據(jù)LSTC公司的經(jīng)驗(yàn)推薦,沙漏能與內(nèi)能的比值小于10%是可以接受的。查看本算例的能量關(guān)系圖,如圖16所示??梢杂?jì)算得到沙漏能平均值是內(nèi)能平均值的2.91%,可見,該分析中單點(diǎn)積分引起的沙漏是可以忽略的。
圖16 沙漏能與內(nèi)能對(duì)比示意圖Fig.16 Comparison diagram of internal energy and hourglass energy
(1)根據(jù)動(dòng)力學(xué)原理和Lagrange乘子約束,提出了考慮慣性載荷的螺旋錐齒輪動(dòng)態(tài)嚙合分析有限元模型,為客觀、準(zhǔn)確地分析螺旋錐齒輪的動(dòng)態(tài)嚙合性能提供了必要的準(zhǔn)備。
(2)慣性載荷對(duì)齒輪運(yùn)動(dòng)特性和齒面接觸力的波動(dòng)幅值和頻率均有顯著影響:隨著慣性載荷的增大,齒面接觸力的波動(dòng)幅值增大,波動(dòng)頻率減小。
(3)慣性載荷對(duì)輪齒接觸的接觸區(qū)寬度、齒面接觸應(yīng)力以及最大應(yīng)力位置有明顯影響:隨著齒輪慣性載荷的增大,接觸區(qū)寬度增大,齒面接觸應(yīng)力增大。
[1]Wilcox L,Nowell G.Improved Finite Element Model for Calculating Stress in Bevel and Hypoid Gear Teeth[C].AGMA 97FTM5,Nov.1997.
[2]李盛鵬,方宗德.弧齒錐齒輪齒根彎曲應(yīng)力分析[J].航空動(dòng)力學(xué)報(bào),2007,22(5):843 -848.
[3]張金良,方宗德.弧齒錐齒輪齒面接觸應(yīng)力分析[J].機(jī)械科學(xué)與技術(shù),2007,26(10):1268 -1272.
[4]Argyris J,F(xiàn)uentes A,Litvin F L.Computerized integrated approach for design and stress analysis of spiral bevel gears[J].Computer methods applied in mechanics and engineering,2002,191:1057 -1095.
[5]Litvin F L,Vecchiato D.Computerized developments in design,generation,simulation of meshing,and stress analysis of gear drives[J].Meccanica,2005,40:291 -324.
[6]Robert F H,George D B.Comparison of Experimental and Analytical Tooth Bending Stress of Aerospace Spiral Bevel Gears[R].NASA/TM 208903,1999.
[7]李 源,袁杰紅.航空減速器準(zhǔn)雙曲面齒輪動(dòng)態(tài)嚙合仿真分析[J].機(jī)械傳動(dòng),2007,31(5):43 -44.
[8]北京理工大學(xué)ANSYS/LS-DYNA中國技術(shù)支持中心.ANSYS/LS-DYNA算法基礎(chǔ)與使用方法[M].北京;北京理工大學(xué),1999.
[9]Zienkiewicz O C,Taylor K.L.The Finite Element Method,F(xiàn)ifth Edition[M].Butterworth Heinemann,2000.
[10]唐進(jìn)元,曹 康.含過渡曲面的弧齒錐齒輪齒面精確建模[J].機(jī)械科學(xué)與技術(shù),2009,28(3):317 -321.