斜齒輪動態嚙合及疲勞壽命分析

袁麗蕓, 楊曉濤, 黃院星, 付學中, 宋錦江

(廣西科技大學機械與汽車工程學院, 柳州 545006)

齒輪傳動廣泛應用于電動汽車領域,其動力學的研究在齒輪設計方面有著至關重要的作用,而齒輪由于具有傳動平穩、承載能力強、應用廣泛等特點,齒輪性能的改善有著深遠的意義與前景。但經過對齒輪的長期使用,開始出現一些問題,如齒輪傳動時會產生振動,這樣會對齒輪傳動效率產生很大影響,并且齒輪也會由于各種因素如重載、高速等產生疲勞失效,所以對齒輪還需更多的研究。

對于齒輪建模與仿真的研究有很多。杜進輔等[1]以高速斜齒輪為研究對象,建立六自由度的彎-扭-軸動力學模型,采用改進的承載接觸分析計算時變嚙合剛度、嚙合沖擊等,為高速斜齒輪傳動設計提供參考。周洋等[2]以圓柱齒輪傳動系統為研究對象,考慮齒面摩擦、齒側間隙和時變嚙合剛度等非線性因素,引入溫度變化的影響,建立六自由度的齒輪系統非線性動力學模型,并采用4～5階龍格-庫塔算法對模型進行求解,結合分岔圖、相圖和Poincare 映射圖,分析溫度變化和激勵頻率對齒輪系統動力學的影響。王旭偉[3]使用Rayleigh梁模型,考慮轉動慣量和陀螺效應,建立多自由度的轉子-軸承系統的有限元模型,使用Lagrange方程推導出系統的運動微分方程,得到系統的臨界轉速與不平衡響應,使用的梁單元沒有考慮剪切變形。劉浩[4]以高速齒輪箱中兩級斜齒輪為研究對象,采用集中質量法構建12自由度的系統模型,基于高速工況,采用4階龍格庫塔法求解振動響應。

劉明勇等[5]建立齒面接觸分析模型,研究了齒面摩擦等因素對嚙合性能的影響,驗證了模型的正確性。岳會軍等[6]以一對常嚙合斜齒輪為研究對象,采用有限元法分析斜齒輪靜態和動態接觸,得到齒面接觸應力的大小及分布,為斜齒輪傳動改善提供了基礎。閆博等[7]利用Palmgren-Miner線性損傷累計理論進行齒根疲勞損傷計算,確定齒根斷裂的主要原因,并通過增大齒輪螺旋角、齒根倒角以及齒面修形等方法對齒輪進行結構優化,優化后的斜齒輪未發生斷裂,滿足了疲勞壽命要求,為齒輪修形方式提供參考。陳興彬等[8]在不同載荷譜的情況下,基于材料的疲勞-壽命(S-N)曲線和Miner線性損傷累計理論,利用nCode軟件對直齒輪副進行疲勞可靠性分析,得到齒輪傳動在靜載和動載條件下的最大接觸應力和最小疲勞壽命的區域,為齒輪嚙合疲勞的分布提供參考。王玉等[9]以斜齒輪副為研究對象,利用ANSYS Workbench對不同工況下的齒輪進行應力分析,并基于載荷譜對齒輪副進行了強度校核與疲勞壽命計算,得到最小疲勞壽命出現在輪齒接觸面處,但沒有對修形齒輪進行進一步嚙合研究。韓嘯等[10]運用 ANSYS Workbench對單個輪齒進行靜力有限元分析,得到齒輪齒根應力云圖,對不同修形量之間的對比,得到隨著修形量的增加,齒根應力呈現減小的趨勢,但沒有進行瞬態分析研究。劉建剛等[11]利用KISSsoft對減速器齒輪進行計算,并將理論值與其進行比較分析,得到應力應變分布圖,為傳動設計提供參考。Li等[12]對傳統的疲勞壽命分析進行修正,通過考慮工況和材料的影響,得到反應面法對于電力動車組齒輪疲勞壽命預測比傳統更加準確和安全。Yang等[13]研究了齒輪的應力敏感性和疲勞壽命特性,利用ANSYS Workbench有限元軟件建立了有限元模型,并將有限元分析結果與理論計算結果進行了對比,證明了模型的正確性,但沒有針對齒輪修形進行對比分析。Feng等[14]利用ANSYS Workbench軟件對齒面進行靜態和動態接觸有限元模擬,得到失效斜齒輪的斷裂開始于齒根處,為齒輪動態分析提供參考。

通過以上兩方面的研究對于模型建立都考慮多自由度的情況,但有許多構建的動力學方程還是不完善,而其中對于修形齒輪的動態嚙合研究較少,且許多研究對象是低速齒輪。所以針對這些問題,為了構建更完善的動力學方程,現對每個節點考慮6自由度來建立模型,進行瞬態分析且采用Newmarkβ法進行振動響應求解,之后使用KISSsoft 軟件進行一對斜齒輪副幾何模型的構建,并加以修形,在ANSYS Workbench中進行瞬態分析,且在高速工況下,對修形前、后齒輪的動態嚙合和疲勞壽命兩方面進行對比分析, 以證明該修形方式可行性,為高速斜齒輪研究提供參考。

1 轉子-齒輪-軸承系統振動研究

1.1 轉軸單元

轉軸單元采用Timoshenko梁單元進行轉軸建模,如圖1所示。

圖1 Timoshenko梁單元

包含彎-扭-軸的型函數,且每個節點考慮6個自由度(x,y,z,θx,θy,θz),表達式如式(1)所示。

(1)

(2)

式中:ξ=S/L,S為沿著單元的軸向位置,L為轉軸單元的長度;A、E、G、I、φ、μ、κ′分別為轉軸的橫截面積、彈性模量、剪切模量、截面慣性矩、剪切因子、泊松比、截面剪切因素;而考慮的6個自由度中,x、y、z分別為X、Y、Z方向的變形量,θx、θy、θz分別為繞X軸、Y軸、Z軸的扭轉角度。

根據積分公式求得梁單元的質量矩陣、陀螺矩陣及剛度矩陣[15]。

(3)

式(3)中:MT為彎曲質量矩陣;MR為旋轉慣量矩陣,Mθ為扭轉慣量矩陣;G為陀螺矩陣;K為剛度矩陣;KS考慮了彎曲和剪切變形;Kθ考慮了扭轉變形;Id為直徑截面慣性矩;JP、Jd為極轉動慣量、直徑轉動慣量;ρ為材料密度;上角標“T”表示該矩陣的轉置,“″”和“′”分別為矩陣的二階導數和一階導數。

矩陣表達式[15]及參數取值如式(4)～式(9)所示。

(4)

式(4)中:sym表示對稱矩陣。

(5)

(6)

式(6)中:

(7)

(8)

式(8)中:skew表示反對稱矩陣。

(9)

1.2 圓盤單元

圓盤單元視為剛性圓盤,在不考慮嚙合齒側間隙、齒面摩擦力,嚙合部分采用彈簧單元模擬,如圖2所示,其質量矩陣及陀螺矩陣如式(10)和式(11)所示。

ψ為兩齒輪嚙合面與豎直方向的夾角;α為齒輪的嚙合角;β為齒輪的螺旋角

式中:m1、m2分別為齒輪偏心質量;Id為直徑截面慣性矩;Ip為極截面慣性矩。

齒輪系統的動力學方程為

(10)

(11)

(12)

式(12)中:Ω為旋轉速度;F為外力向量,包括輸出轉矩T2、不平衡力(Fd1x、Fd2x、Fd1y、Fd1y);e為偏心距,兩齒輪不平衡力x方向相反,y方向相同,M為圓盤的質量矩陣;G為圓盤的陀螺矩陣;K為圓盤嚙合矩陣。

嚙合剛度計算是將單位接觸線長度的嚙合剛度視為常數,以嚙合線長度函數與前者的乘積作為時變嚙合剛度。

k(t)=k0L(τ)

(13)

式(13)中:k(t)為時變嚙合剛度;k0為單位接觸線長度的嚙合剛度;τ=t/T,t為時間,T為嚙合周期,τ為標準化時間;L(τ)為斜齒輪時變接 觸線長度函數,mm[6]。

(14)

式(14)中:εα、εβ分別為齒輪副端面重合度、軸向重合度;b為齒寬;βb為基圓螺旋角。

以輸入級轉速6 000 r/min為嚙合周期條件,結合所研究的齒輪傳動系統參數,以《漸開線圓柱齒輪承載能力計算方法》(GB/T 3480—1997)定義的嚙合剛度為單位接觸線長度的嚙合剛度,基于式(14)計算時變嚙合剛度。

齒輪嚙合剛度矩陣[16]如式(15)所示。

1.3 軸承單元

對于軸承只考慮了兩個自由度(x、y)軸承剛度矩陣及阻尼矩陣,即

(16)

根據等向同性原則,Kxx=Kyy=Kb,Kxy=-Kyv,Kxx=Kyy=Kb,Kxy=-Kyv。

1.4 轉子-軸承系統

通過有限元法建立轉子-軸承系統模型,包括轉軸單元、圓盤單元、偏心質量單元、軸承單元,整個系統總共12個軸單元,兩個圓盤單元,4個軸承單元,14個節點,軸長0.1 m,軸半徑為0.01 m,如圖3所示。

圖3 轉子-軸承劃分圖

通過以上的矩陣采用有限元法進行組集,構建動力學方程為

(17)

式(17)中:C為總阻尼包括軸承阻尼,通過計算得到轉速為6 000 r/min轉子-軸承系統前六階的固有頻率及0～6 000 r/min的campbell圖,如表1和圖4所示。

(15)

表1 固有頻率

表1結果得到彎曲的固有頻率,而圖4為轉速與固有頻率的變化圖,且前四階變化幾乎無變化,五階時出現稍微下降,六階出現上升趨勢,但由于系統陀螺力矩較小導致總體變化不明顯,在6 000 r/min位置是該轉速的第五階固有頻率。

1.5 位移響應

考慮結構阻尼,使用瑞利阻尼來求系統的阻尼:C=αK+βM,其中α取0.02,β取0.03,K為系統的剛度矩陣,M為系統的質量矩陣,再通過Newmarkβ法計算位移響應,得到結果如圖5～圖8所示。

圖5 時域圖

圖5(a)為小齒輪的彎曲振動時域圖,而圖5(b)為0～0.4 s的放大圖,0.1 s之前從負方向開始偏移,逐漸向0變化,0.1 s之后圍繞0進行波動。圖6(a)為彎曲振動頻域圖,圖中出現的峰值對應的頻率為第一臨界轉速時的固有頻率,圖6(b)中的峰值對應的固有頻率為第二臨界轉速時的固有頻率。

圖6 頻域圖

圖7(a)為大齒輪的彎曲振動時域圖,而圖7(b)為其中0～0.4 s的放大圖,0.1 s之前也從負方向開始偏移,比小齒輪的偏移要大,逐漸向0變化,0.1 s之后圍繞0進行波動。圖8(a)為彎曲振動頻域圖,出現的峰值對應的頻率為第一臨界轉速時的固有頻率,圖8(b)中的峰值對應的固有頻率為第二臨界轉速時的固有頻率。

2 齒輪修型及嚙合仿真

2.1 齒輪修型

在所用的仿真模型為一對斜齒輪副。主要參數為Z1=21,Z2=65,a=65 mm,m=1.5 mm,b1=30 mm,b2=28 mm,Z1和Z2分別為小齒輪和大齒輪的齒數,a為兩齒輪間的中心矩,m為齒輪的模數,b1和b2分別為兩齒輪的齒寬。其中小齒輪為驅動輪,大齒輪為從動輪。使用KISSsoft 軟件進行以上參數的設置得到幾何模型。為了降低網格劃分的難度和對計算機硬件的要求,提高計算速度,在保證精度的前提下,省略了齒輪上的鍵槽和倒角,之后再進行齒輪的修形設置,基于湯兆平等[17]的修型方式比較單位長度載荷和傳遞誤差得到修形后的幾何模型。修形參數如表2所示,比較形式如圖9和圖10所示。

表2 修形參數

圖10 單位長度載荷

圖9中傳遞誤差從1.79 μm下降到1.51 μm,傳動平穩性得到提高,圖10中單位長度最大載荷從243 N/mm下降到224 N/mm,且修型后在中間位置載荷最小,使載荷分布更加均勻,有利于齒輪平穩嚙合,且齒輪的比滑率在(-1,1)表現出齒輪嚙合狀態良好。修形后斜齒輪的幾何模型如圖11所示。

圖11 齒輪幾何模型

2.2 齒輪嚙合仿真

將KISSsoft中的模型轉換成stp格式的文件,之后直接導入ANSYS Workbench中。將材料定義為20CrMnTiH,其中屈服強度為1 500 MPa,屈服極限為1 800 MPa,且彈性模2.07×1011MPa,泊松比0.29。考慮到齒輪的實際工作條件,接觸齒輪表面的接觸類型設置為摩擦,摩擦因數初始化為0.15,齒輪的動態接觸如圖12所示,齒面的接觸狀態如圖13所示。

圖12 動態接觸

圖13 接觸狀態

圖12中紅色部分為小齒輪所有接觸齒面,藍色部分為大齒輪所有接觸齒面,可以看到齒輪動態接觸的狀態,而設置的參數方法采用廣義的朗格朗日法(augmented Lagrange),關閉微小滑動(small sliding),方向為在高斯點上(on Gauss point),且滲透容差(penetration tolerance)和搜索半徑倍數(pinball factor)的分別為0.1 和0.2,齒輪接觸設為調整接觸(adjust to touch)。而圖13中為接觸部位的各種狀態分布,其中大部分都處于接近閉合狀態。

為了保證結果的準確性,在網格劃分中,首先采用掃掠(sweep)的方法進行劃分,之后再對接觸的齒面再進行網格細化。得到非修形齒輪劃分的網格為146 979個節點和30 765個單元,修形齒輪劃分的網格為279 774個節點和61 275個單元, 且網格的質量均較高,且平均質量都達到0.7左右,網格劃分圖如圖14所示。

圖14 網格劃分

圖14中網格劃分采用六面體網格,它比四面體網格的計算精度更加準確,且根據Hanoca等[18]中得到網格劃分越小,結果越好,所以把網格尺寸設為1 mm,使得在齒厚方向劃分出兩個網格如圖14中右下角,這樣比Xing等[19]中齒厚方向上只有一列網格劃分進行瞬態分析的結果更加精確,可以看出齒厚方向的變形情況,使應力和疲勞計算更加精確,這樣更符合有限元劃分規范。

接下來進行分析設置為直接法(direct),打開弱彈簧(weak springs)和大變形(large deflection),非線性控制方法為完全法(full),并打開力收斂(force convergence),位移收斂(displacement convergence)及線搜索(line search),且穩定性(semi-implicit)設置為程序控制(program controlled)。速度及阻力矩設置如圖15所示。

圖15 速度與力矩的設置

圖中大齒輪上施加了2 000 N·mm的恒定阻力扭矩,小齒輪上施加了6 000 r/min的恒定轉速,使得齒輪處于高速的工作狀態,其中阻力矩值和轉速只是某一個試驗值,并不是根據準確的實際工作條件來設置的。

2.3 動力學仿真分析

根據有限元模型的計算,將計算結果導入后處理器,獲得了修形前、后的齒輪嚙合過程中的等效應力分布云圖,如圖16所示。

圖16 等效應力分布圖

圖16(a)為未修形齒輪的等效應力分布圖,圖16(b)為未修形齒輪等效應力分布圖,下左的區域圖為齒輪的最大等效應力區域放大圖,且處在兩齒輪接觸的大齒輪齒根處,最小在小齒輪齒根處,得到修形的斜齒輪最大等效應力從1 579 MPa減小到1 426 MPa,整體減小了6.6%,其中除了齒面接觸靠近齒根處有較大的等效應力產生,其他地方的等效應力相對于齒面接觸處都較小。

2.4 疲勞分析

將計算結果及材料S-N曲線圖導入疲勞工具中,S-N曲線圖如圖17所示。

圖17 20CrMnTiH材料S-N曲線圖

圖17中,橫坐標為旋轉圈數N且使用對數形式,縱坐標為應力S,在圈數1～6.3的應力隨圈數增大直線下降,之后應力趨于平穩。此外,疲勞載荷類型設置為(基于零)based-zero,疲勞應力系數也設置為0.8,經過后處理,得到修形前、后斜齒輪的疲勞壽命和安全系數分布圖,如圖18、圖19所示。

圖18 疲勞壽命分布圖

圖19 疲勞安全系數分布圖

圖18(a)為未修形齒輪的疲勞壽命分布圖,圖18(b)為修形齒輪疲勞壽命分布圖,而左下角為最小疲勞壽命區域放大圖,可知疲勞壽命最小區域為等效應力最大的區域,且修形齒輪的最小疲勞壽命增加。

圖19(a)為未修形齒輪安全系數分布圖,圖19(b)為修形齒輪安全系數分布圖,左下角為最小安全系數區域放大圖,最小安全系數從0.83增加到0.92,修形齒輪的安全性更加可靠。綜上對于修形前后的等效應力、疲勞壽命及安全系數的分析,高速齒輪可以通過齒廓和齒向兩方面的修形,降低了齒輪的應力集中,減少了齒輪之間的傳遞誤差,增大了齒輪使用的疲勞壽命和安全系數。

3 結論

對轉子-軸承模型進行了系統中齒輪的固有特性分析,在之后通過KISSsoft和ANSYS Workbench軟件以高速條件對修形前、后斜齒輪副進行瞬態仿真分析,并進行對比分析,為高速斜齒輪的研究提供了參考,最終得到以下結論。

(1)通過MATLAB對多自由度的轉子-軸承系統在考慮輸入輸出轉矩、不平衡力等外力,建立了更完善的動力學方程最終得到系統的固有頻率,齒輪彎曲的時域圖、頻域圖。

(2)通過應力分析可知斜齒輪在高速工作狀態下其隨時間變化最大應力隨周期性變化,最大應力處于兩齒輪接觸的齒根位置,且修形后齒輪等效應力減小。

(3)通過疲勞分析可知高速斜齒輪的最小疲勞壽命和最小安全系數處于最大應力處,且修形齒輪的疲勞壽命和安全系數相應增加,證明了該修形方式的正確性。