斜齿轮动态啮合及疲劳寿命分析

2024-03-22袁丽芸杨晓涛黄院星付学中宋锦江

科学技术与工程 2024年5期

袁丽芸, 杨晓涛, 黄院星, 付学中, 宋锦江

(广西科技大学机械与汽车工程学院, 柳州 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,修形齿轮的安全性更加可靠。综上对于修形前后的等效应力、疲劳寿命及安全系数的分析,高速齿轮可以通过齿廓和齿向两方面的修形,降低了齿轮的应力集中,减少了齿轮之间的传递误差,增大了齿轮使用的疲劳寿命和安全系数。