董昭榮,張 楚,董惠敏

(大連理工大學(xué)機(jī)械工程學(xué)院,大連 116024)

0 引言

齒輪傳動(dòng)在機(jī)械領(lǐng)域占有極其重要的地位,齒輪嚙合過(guò)程中的變形仿真是齒輪設(shè)計(jì)的關(guān)鍵步驟,然而由于齒輪三維動(dòng)態(tài)嚙合有限元接觸分析具有高度非線性以及大量的自由度特點(diǎn),其計(jì)算往往不能兼顧精度與效率。

為了提高計(jì)算效率,凝聚法是減少計(jì)算自由度的有效方法。GUYAN[1]提出靜凝聚方法,忽略慣性力的影響;CALLAHAN[2]改進(jìn)了Guyan縮聚法,計(jì)入了部分慣性力的影響;李偉明[3]利用逐級(jí)迭代縮減方法驗(yàn)證了迭代級(jí)數(shù)與計(jì)算精度的規(guī)律;FRISWELL等[4]迭代改進(jìn)方法將靜凝聚進(jìn)一步擴(kuò)展到動(dòng)力學(xué)領(lǐng)域;XIA等[5]重新推導(dǎo)計(jì)入質(zhì)量陣影響的迭代公式,使現(xiàn)有的動(dòng)凝聚方法收斂速度更快;WEI等[6]直接從簡(jiǎn)化的振動(dòng)方程推導(dǎo)出響應(yīng)靈敏度,提高基于響應(yīng)的大型結(jié)構(gòu)有限元模型修正的計(jì)算精度和效率。

齒輪的接觸分析對(duì)研究其傳動(dòng)性能有重要意義。VIDEKEY[7]提出應(yīng)用赫茲接觸理論得到齒輪應(yīng)力于嚙合線方向的變化規(guī)律;LITVIN等[8]利用微分幾何和解析法研究了齒輪嚙合的接觸條件;譚自然[9]基于解析法研究齒根裂紋對(duì)齒輪時(shí)變嚙合剛度及動(dòng)態(tài)接觸特性的的影響。ZHANG等[10]基于點(diǎn)面線性互補(bǔ)方法提出解決三維齒輪準(zhǔn)靜態(tài)接觸分析方法。

目前,現(xiàn)有的動(dòng)凝聚方法多數(shù)基于靜凝聚公式基礎(chǔ)上對(duì)其計(jì)入慣性力的影響,針對(duì)不同的計(jì)入方式出現(xiàn)許多動(dòng)凝聚的改進(jìn)方法,其中迭代改進(jìn)方法最為常見(jiàn)。物理型縮減技術(shù)方面,三級(jí)近似遞推縮減公式的精度最好[3],但是目前還未有人基于三級(jí)近似遞推縮減公式推導(dǎo)其迭代改進(jìn)的動(dòng)凝聚公式。據(jù)此針對(duì)大型結(jié)構(gòu)的有限元模型,本文推導(dǎo)了一種新的動(dòng)凝聚公式,并且應(yīng)用到齒輪的動(dòng)力學(xué)接觸分析研究中。有限元方法是研究齒輪的動(dòng)力學(xué)接觸分析的主要方法,但是其非線性運(yùn)算使其分析過(guò)程需要占據(jù)大量的資源,針對(duì)這一問(wèn)題,本文在有限元方法的基礎(chǔ)上,結(jié)合動(dòng)凝聚方法、Bozzak-Newmark方法以及文獻(xiàn)[11]中點(diǎn)面線性互補(bǔ)方法,推導(dǎo)出新的齒輪動(dòng)態(tài)接觸分析求解方程。并且利用本文研究方法計(jì)算了一對(duì)斜齒輪的接觸應(yīng)力,與ANSYS斜齒輪對(duì)的動(dòng)態(tài)接觸分析結(jié)果進(jìn)行對(duì)比,該方法計(jì)算效率與計(jì)算精度得到了驗(yàn)證。

1 新的動(dòng)凝聚方法分析

動(dòng)凝聚方法可以減少運(yùn)動(dòng)微分方程計(jì)算自由度,極大的提高結(jié)構(gòu)動(dòng)力分析的效率。對(duì)于具有N自由度模型的運(yùn)動(dòng)微分方程:

(1)

式中:[M]是質(zhì)量陣,[C]是阻尼陣,[K]是剛度陣,{Q(t)}為位移向量,{F}為系統(tǒng)所受外力向量。為減少計(jì)算量,構(gòu)造縮減矩陣[T],將整體自由度凝聚到主自由度上,即:

{Q(t)}=[T]{q(t)}

(2)

將式(2)代入式(1),并在式(1)兩邊同時(shí)乘以縮減矩陣的轉(zhuǎn)置可得:

(3)

式中:[MR]=[T]T[M][T]、[CR]=[T]T[C][T]、[KR]=[T]T[K][T]、{FR}=[T]T{F}分別為凝聚后的質(zhì)量陣、阻尼陣、剛度陣和外部力向量。凝聚方法的差別之處主要在于縮聚轉(zhuǎn)換矩陣[T]不同[1]。WEI等[6]針對(duì)大型結(jié)構(gòu)從簡(jiǎn)化的振動(dòng)方程推導(dǎo)一種動(dòng)凝聚方法提高了有限元模型的計(jì)算精度與計(jì)算效率。其最終縮減矩陣由以下公式計(jì)算:

(4)

(5)

為驗(yàn)證其與傳統(tǒng)動(dòng)凝聚關(guān)系,將上述公式[M]d反代入到迭代方程(4)得:

(6)

由此可見(jiàn),文獻(xiàn)[6]中利用二級(jí)近似遞推縮減公式迭代計(jì)算縮減矩陣,根據(jù)張德文推導(dǎo)縮減公式,三級(jí)近似遞推縮減公式計(jì)算精度最高,因此重新推導(dǎo)分塊后n自由度系統(tǒng)的運(yùn)動(dòng)微分方程:

(7)

根據(jù)上式的第二行可得:

(8)

(9)

(10)

(11)

(12)

(13)

(14)

三級(jí)近似模型縮聚公式精度高于二級(jí)近似模型縮聚[3],因此利用基于文獻(xiàn)[6]二級(jí)近似模型縮聚迭代改進(jìn)思想,利用迭代思想重新計(jì)入三級(jí)近似模型縮聚公式中的慣性力的影響,進(jìn)一步提高動(dòng)凝聚的計(jì)算精度。

動(dòng)態(tài)凝聚方法是處理大型有限元結(jié)構(gòu)問(wèn)題的有效方法。利用轉(zhuǎn)換矩陣[T],將整體計(jì)算自由度轉(zhuǎn)換為主自由度,從而使用變換矩陣形成更小計(jì)算自由度的簡(jiǎn)化系統(tǒng)。在凝聚之前的計(jì)算自由度達(dá)到108,凝聚之后的自由度僅僅為104。因此,對(duì)簡(jiǎn)化系統(tǒng)進(jìn)行接觸分析可節(jié)省大量計(jì)算時(shí)間。

2 斜齒輪系統(tǒng)的動(dòng)力學(xué)點(diǎn)面線性互補(bǔ)接觸分析

根據(jù)文獻(xiàn)[10]對(duì)于一對(duì)斜齒輪有限元模型,為建立主從動(dòng)齒輪之間的點(diǎn)面互補(bǔ)模型,假設(shè)主動(dòng)齒輪2為接觸體,從動(dòng)齒輪1為目標(biāo)體。接觸體上潛在接觸區(qū)域的點(diǎn)分別對(duì)應(yīng)于目標(biāo)體上潛在接觸區(qū)域面網(wǎng)格中的最小的面,如圖1所示。主從動(dòng)齒輪之間的法向接觸力{Fn}與接觸間隙{gn}形成互補(bǔ)關(guān)系,即為:

圖1 主從動(dòng)齒輪接觸單元(ok-pk為接觸面法線方向)

{Fn}·{gn}=0,{Fn}≥0,{gn}≥0

(15)

因此,本文第2節(jié)的內(nèi)容圍繞如何對(duì)上式建立懲罰方程進(jìn)行展開(kāi)。首先在目標(biāo)體建立全局坐標(biāo)系S{O:x,y,z},在接觸體上潛在接觸區(qū)域內(nèi)每一個(gè)節(jié)點(diǎn)為原點(diǎn)建立節(jié)點(diǎn)坐標(biāo)系Sk{ok;nk,αk,τk},k是接觸體上第k個(gè)節(jié)點(diǎn)ok。

本文假設(shè)邊界條件為初始速度、加速度和位移均為0。

斜齒輪的重合度一般為2～3,在嚙合過(guò)程中,最多有3個(gè)齒面同時(shí)接觸。因此分別選擇有限元模型上主動(dòng)輪3個(gè)接觸齒面節(jié)點(diǎn)和中心,以及從動(dòng)齒輪的3個(gè)接觸齒面的節(jié)點(diǎn)作為主自由度所在節(jié)點(diǎn),如圖2和圖3所示,并且在ANSYS驗(yàn)證時(shí)將3對(duì)接觸面設(shè)置接觸對(duì)。將凝聚后的坐標(biāo)代入兩個(gè)齒輪的運(yùn)動(dòng)微分方程,方程兩側(cè)同乘[T]T可得:

圖2 灰色圓點(diǎn)表示從動(dòng)齒輪1主自由度所在節(jié)點(diǎn) 圖3 灰色圓點(diǎn)表示主動(dòng)齒輪2主自由度所在節(jié)點(diǎn)

(16)

式中:[MR]=[T]T[M][T],[CR]=[T]T[C][T],[KR]=[T]T[K][T],{FR}=[T]T{F}。

由文獻(xiàn)[10]可知變形后間隙與初始間隙與目標(biāo)體位移與接觸體位移關(guān)系為:

(17)

式中:{gn}為變形后接觸體節(jié)點(diǎn)與目標(biāo)面間隙,{gn,0}為初始間隙,{q}t為目標(biāo)體節(jié)點(diǎn)位移向量,[S]t將目標(biāo)體上節(jié)點(diǎn)沿全局坐標(biāo)系下的位移變換為接觸體節(jié)點(diǎn)映射在目標(biāo)面上節(jié)點(diǎn)沿法向接觸方向的位移,{qn}c為接觸體節(jié)點(diǎn)沿法向接觸力方向位移。

根據(jù)換齒的過(guò)程,在整個(gè)齒輪嚙合的過(guò)程中,可視為換齒過(guò)程的循環(huán),因此可以只研究一個(gè)換齒時(shí)間段的應(yīng)力變化,將時(shí)間等距離散。

根據(jù)Bozzak-Newmark方法[11],令αB=0.2,βB=0.5,γB=0.25,接觸體凝聚后運(yùn)動(dòng)微分方程可以轉(zhuǎn)換為:

[K*]c{q}c,j+1={QR}c-{QRF}c,j+1+{F*}c,j+1

(18)

根據(jù)文獻(xiàn)[10]消去接觸體平衡方程中切向力影響,經(jīng)推導(dǎo)后接觸體沿法向接觸方向的平衡方程為:

(19)

同理,凝聚后的目標(biāo)體平衡方程經(jīng)Bozzak-Newmark方法轉(zhuǎn)化后,利用虛功原理可得到目標(biāo)體沿法向接觸方向的平衡方程[10]可以表示為:

(20)

將最終推導(dǎo)的接觸體和目標(biāo)體的運(yùn)動(dòng)微分方程代入到文獻(xiàn)[10]中推導(dǎo)所得幾何方程(17)消除節(jié)點(diǎn)位移,得到接觸間隙和力的關(guān)系式:

(21)

最后整理出接觸體和目標(biāo)體線性互補(bǔ)求解方程為:

(22)

式中:

根據(jù)接觸體平衡方程,目標(biāo)體平衡方程,間隙位移關(guān)系,消去接觸體和目標(biāo)體的節(jié)點(diǎn)位移得到最終的線性互補(bǔ)方程。此時(shí)結(jié)合lemke算法可以得到各個(gè)時(shí)刻的接觸壓力和接觸間隙,并反帶入到主從動(dòng)齒輪的運(yùn)動(dòng)微分方程,得到各時(shí)刻齒輪的變形位移。

3 驗(yàn)證示例

為驗(yàn)證推導(dǎo)得到動(dòng)凝聚方法的有效性,利用ANSYS中beam188單元建立軸-圓盤(pán)有限元模型,如圖4所示,內(nèi)孔直徑0.02 m,外徑從左到右依次為0.03 m、0.04 m、0.02 m、0.03 m,長(zhǎng)度從左到右依次為0.4 m、0.6 m、0.6 m、0.4 m,彈性模量2.01 GPa,泊松比0.3,密度7.8e-3 kg/m3,圓盤(pán)質(zhì)量0.35 kg,圓盤(pán)直徑0.04 m,設(shè)置在軸中心節(jié)點(diǎn)處。

圖4 軸-圓盤(pán)有限元模型

圖5 各方法計(jì)算中心節(jié)點(diǎn)位移結(jié)果對(duì)比圖

圖6 計(jì)算步驟示意圖

對(duì)軸中心節(jié)點(diǎn)處施加繞X軸轉(zhuǎn)矩M=10*sin(t),分別利用完全法、新凝聚方法計(jì)算得到中心節(jié)點(diǎn)位移變化。

利用該凝聚方法計(jì)算中心節(jié)點(diǎn)Y方向平穩(wěn)后一個(gè)變化周期內(nèi)最大位移為3.583×10-5,Wei Tian所提出凝聚法的計(jì)算誤差0.885%,完全法計(jì)算得到3.615×10-5,計(jì)算誤差為0.58%,計(jì)算時(shí)間由完全法的9.23′′縮短為1.201′′,計(jì)算精度約為Wei Tian法的1.53倍,由此驗(yàn)證了所推導(dǎo)得到的新凝聚公式的有效性。

為驗(yàn)證動(dòng)態(tài)線性互補(bǔ)接觸分析方法的有效性,以本文研究方法和ANSYS瞬態(tài)動(dòng)力學(xué)完全法對(duì)斜齒輪對(duì)在固定轉(zhuǎn)矩和固定轉(zhuǎn)速條件下進(jìn)行接觸分析。

兩個(gè)齒輪的法向模數(shù)為m=3.0,齒數(shù)為36,壓力角20°,齒寬0.03 m,螺旋角15°/-15°,彈性模量2.06 GPa,泊松比0.3,密度7.85e3 kg/m3,材料阻尼0.004 N·s/m。對(duì)主動(dòng)齒輪中心施加100 rad/s的轉(zhuǎn)速,對(duì)從動(dòng)輪施加100 N·m的邊界條件,為防止產(chǎn)生較大沖擊,采用斜坡函數(shù)形式加載轉(zhuǎn)速和轉(zhuǎn)矩。設(shè)置收斂條件為:

|ωk-ωk-1|/ωk≤ε

(23)

式中:k為迭代次數(shù),ω為固有頻率,設(shè)誤差容差ε為10-4。以下是利用ANSYS和本研究方法計(jì)算的在各時(shí)刻主動(dòng)齒輪最大接觸壓力和傳動(dòng)誤差隨時(shí)間變化圖形。

根據(jù)圖7和圖8可知,ANSYSY瞬態(tài)動(dòng)力學(xué)與所提出的線性互補(bǔ)接觸分析方法結(jié)算結(jié)果中,最大接觸壓力相對(duì)誤差為4.01%,最大傳動(dòng)誤差差值相差0.94%,將圖2和圖3中的主自由度對(duì)應(yīng)的接觸對(duì)分別設(shè)置為1、2、3。兩種齒輪動(dòng)態(tài)接觸分析方法在0.006 s和0.007 2 s齒面接觸壓力分布圖對(duì)比情況為如圖9和圖10所示。

圖7 最大接觸壓力時(shí)域圖 圖8 傳動(dòng)誤差時(shí)域圖

圖9 0.006 s各齒面接觸壓力分布情況對(duì)比

圖10 0.007 2 s各齒面接觸壓力分布情況對(duì)比

0.006 s時(shí)刻ANSYS計(jì)算最大壓力為473.1 MPa,線性互補(bǔ)方法計(jì)算值為455.5 MPa,相對(duì)誤差為3.72%;0.007 2 s時(shí)刻ANSYS計(jì)算最大壓力為536.3 MPa,線性互補(bǔ)方法計(jì)算值為550.7 MPa,相對(duì)誤差為2.69%。

利用縮減轉(zhuǎn)換矩陣將運(yùn)動(dòng)微分方程中的全階自由度凝聚到更少的主自由度上,在凝聚前自由度大約為109,凝聚后的自由度約為104,極大的減少計(jì)算量,利用ANSYS瞬態(tài)動(dòng)力學(xué)分析計(jì)算需要總時(shí)間為3060′58′′,而利用線性互補(bǔ)方法求解需要總時(shí)間為395′59′′。因此該方法具有較好的精度與較高的計(jì)算效率。

4 結(jié)論

針對(duì)大型結(jié)構(gòu)有限元模型,本文利用迭代方法,對(duì)三級(jí)近似迭代動(dòng)凝聚公式重新計(jì)入慣性力的影響,提高動(dòng)凝聚公式計(jì)算精度。利用動(dòng)凝聚方法可將計(jì)算自由度縮減到所選主自由度上,減少分析過(guò)程計(jì)算量,然后結(jié)合點(diǎn)-面線性互補(bǔ)方法完成了齒輪的動(dòng)態(tài)接觸分析。由算例可知該方法在一定的轉(zhuǎn)速和扭矩下具有較高的計(jì)算精度和計(jì)算效率。本文為動(dòng)凝聚方法以及齒輪動(dòng)態(tài)有限元接觸分析提供了新的有效的思路。