基于弹性圆柱壳的柔性齿轮副建模及动力学分析

张 洋，关先磊，王青山

（中南大学 高性能复杂制造国家重点实验室，长沙 410083）

齿轮传动系统广泛应用于汽车、飞机等各种场合，被认为是最重要的机械部件之一。在工作过程中，传动齿轮存在横向与面内的弹性振动，由此会引起结构失效等动力学现象。这种结构失效现象在高速轻质重载齿轮中经常出现，例如高速动车牵引齿轮传动系统以及航空发动机齿轮系统。因此，建立考虑轮体柔性的齿轮传动动力学模型具有十分重要的意义。

为了获得高效齿轮，人们对齿轮系统的动态建模和振动特性分析进行了大量的研究。对于齿轮的轮体结构，通常有三种建模方法。其一是扭振模型即刚性模型，Han和Chu[1]研究了齿轮传动转子系统周期角角速度运动下的动力学特性。Chen等[2]研究了双螺旋齿轮的旋转动力学，考虑到了轴承、回转效果，其中采用铁木辛柯梁作为轴单元，齿轮则视为刚性体。He 等[3]提出了一种新颖的刚性齿轮啮合模型，并采用改进的势能法计算了齿轮偏心对时效啮合刚度的影响。Shi等[4]完成了考虑时变参数和侧隙的啮合圆柱齿轮对的建模和分析。其二是以弹性环作为齿轮的轮体模型，Christopher和Robert[5]研究以耦合弹性环为模型，分析研究高速柔性齿轮副的振动。Chen等[6]研究了带挠性正齿环行星齿轮组的动态仿真。Fan 等[7]根据壳理论和铁木辛柯梁理论建立了行星齿轮内环和转轴模型。其三是以柔性圆盘作为齿轮轮体模型。Vinayak 和Singh[8]建立了柔性齿轮的多体动力学模型，考虑了分布齿轮啮合刚度和齿胚模型，求解了一系列的模态问题和响应问题。

综合比较上述几种模型可知，刚性模型具有模型简单，计算量小的优点，但是刚性模型没有考虑到齿轮齿面变形，和实际齿轮贴合度存在差距，尤其是齿轮轮体厚度较薄的齿轮，刚性模型和柔性模型的差距较大。相比较于刚性模型，柔性轮体模型则更加符合工程实际，因此柔性模型的建模和分析是非常有必要的。基于弹性环的柔性齿轮副模型只适用于齿轮的齿宽较小的情况，而以柔性圆盘为齿轮轮体的齿轮模型虽然适用范围较广，但是存在模型收敛性不好的缺点。基于现有柔性齿轮副模型研究的不足，本文基于弹性圆柱壳单元，建立一种新颖的柔性齿轮轮体模型，相比于上述两种柔性模型，基于弹性圆柱壳的齿轮副模型，具有良好的收敛性，可适用于齿轮宽度较大的情况，还可以进行齿轮宽度方向的动力学特性研究。

综上所述，本文首先基于微分求积有限元法和旋转圆柱壳理论，推导旋转圆柱壳的动能和势能，完成柔性齿轮的轮体建模。两圆柱齿轮由空间固定的具有时变啮合刚度的弹簧耦合，表示齿轮之间的啮合效应。并且每个圆柱壳用弹性基础支承，来表示轴以及轴承的弹性。最后分别使用拉格朗日方程和Newmark 迭代来求解柔性齿轮副的模态和动态响应。本文以某航空航天齿轮副为例，通过收敛性分析，研究了本文模型的数值稳定性；然后与ANSYS的有限元模型进行模态特性的对比，验证本文模型的正确性。接下来，数值计算了大转速范围下柔性齿轮副系统固有频率随转速的变化规律。最后通过与传统的基于齿轮体刚性假设的扭振模型进行对比，阐述齿轮柔性对系统动力学特性带来的影响。

1 理论推导

1.1 柔性齿轮副的几何模型

柔性齿轮副的模型如图1所示。齿轮简化为圆柱壳模型，齿轮啮合副简化为弹簧模型。分别在两个圆柱壳的中心层上定义正交圆柱坐标系来描述壳体的振动，坐标轴xi，θi，zi(i=1，2)分别沿着轴向，周向以及径向。以该坐标系为参考，壳体的位移用ui，vi，wi(i=1，2)来表示。齿轮的转速为Ωi，根据齿轮啮合原理可知Ω2=-R1Ω1/R2，齿轮的基圆半径为Ri，壳厚为hi，齿宽为li。假设齿轮的密度、泊松比、弹性模量以及剪切模量分别用ρi，μi，Ei，Gi表示。每个齿轮分别采用径向弹簧（kri）和切向弹簧（kθi）来代表轴的弹性，啮合弹簧的刚度为km。

图1 柔性齿轮副模型

1.2 齿轮副系统的能量分析

根据1 阶剪切变形理论，壳体的位移分量可由中间面的位移来表示[9]：

壳体上任意点的位置矢量和速度矢量分别为：

其中：i、j、k分别表示x、θ、z方向上的单位矢量。

旋转圆柱壳的动能为：

将式（2）代入式（3）可得：

圆柱壳的弹性势能公式为：

其中应力应变为：

其中：

E为弹性模量，μ为泊松比。

两齿轮对应的啮合点之间用弹簧连接，其弹簧的啮合势能为：

其中：v1和v2分别表示主、从动轮上啮合点位置的切向位移，km(t)为时变啮合弹簧的刚度。

在圆柱壳中间环内部，分别采用径向弹簧和切向弹簧模拟轴的弹性。支撑弹簧势能为：

kri和kθi分别表示径向和切向单位长度上的弹簧刚度。但是在有限元方法中，弹簧只能添加在节点上，因此，采用一下等效方式，得到节点上的等效弹簧刚度：

其中：C表示圆环周长，N表示有限元网格的节点数，keri表示施加在节点上的等效弹簧刚度。

因此系统的总应变能为：

1.3 系统的离散及求解

将系统的的总应变能和总动能代入拉格朗日方程可得齿轮啮合系统的动力学方程为：

即可得到圆柱壳齿轮啮合结构的整体质量矩阵M、刚度矩阵K、阻尼矩阵C分别为：

其中：

其中：Mk是对称矩阵。A和B为微分权系数矩阵；C为积分权系数矩阵[10]，E为单位矩阵。

最后使用Newmark 积分迭代方法进行动态响应求解[11]。Newmark积分迭代方法的算法流程图如图2所示。

2 算例分析

本文采用实际工程中的齿轮参数作为算例对象，主从动齿轮的尺寸参数如表1所示。材料参数为：弹性模量E=2.1×1011Pa；密度ρ=7 850 kg/m3；泊松比μ=0.3。需要特别说明支撑弹簧设置和载荷设置。在研究系统模态时，同时考虑径向和周向的支撑弹簧，模拟转轴和载荷作用下的工况。而在进行动力响应计算时，只在主动轮上施加周向载荷（即扭矩），同时去掉周向支撑弹簧，保留径向支撑弹簧；而从动轮上则不施加载荷，继续保留周向和径向支撑弹簧。

表1 齿轮副尺寸参数

根据参考文献得到的齿轮啮合刚度的解析模型[12]，可得到如图3所示的啮合刚度曲线，取算术平均值可得平均啮合刚度1.350 3×109N/m，这将作为后续进行模态特性分析的啮合刚度。齿轮的啮合刚度平均分配到对应的有限元啮合节点上。

图3 时变啮合刚度曲线

2.1 收敛性和模型验证

首先进行模态收敛性分析，通过计算主从动齿轮啮合模型的模态，来验证微分求积有限元法的收敛性。分别沿两圆柱壳的轴向和周向划分单元数为M和N，并且每个单元的微分求积节点数为n。由于计算量过大，先确定轴向单元数M=2；研究周向单元数和节点数对收敛性的影响。图中偏差计算公式为(fi-fend)/fend%，其中fi为某单元数N和单元节点数n时的频率，fend为N和n同时取最大值时得到的频率，由图4易得，模型的收敛性极好，根据收敛性结果，在后续的分析中选择单个圆柱壳轴向单元数量为2，周向单元数量为2×6，单元节点为6进行计算。

图4 柔性啮合副模型的前4阶频率收敛性

为进一步验证模型结果的正确性，接下来进行模态对比验证。在ANSYS中建立有限元模型，圆柱壳采用Shell单元，每个圆柱壳40×80个单元；啮合弹簧和支撑弹簧使用COMBIN14 单元有限元模型如图5所示。

图5 ANSYS模型

表2展示了本模型和ANSYS 模型的前12 阶自然频率对比，可以发现，两者结果吻合很好，图6展示了齿轮啮合模型的前4 阶模态振型，其结果也与ANSYS 结果吻合。上述结论充分验证了模型的正确性。

图6 柔性齿轮副前4阶振型结果对比

表2 两种模型自然频率对比

2.2 坎贝尔图分析

在进行动态响应分析之前，首先计算系统的自然频率并得坎贝尔图。在啮合刚度为1.350 3×107N/m的情况下，计算不同转速时系统的自然频率，结果如图7所示。在零速时，由于啮合刚度的存在破坏了圆柱壳的轴对称性，导致固有频率不同。

图7 柔性齿轮副坎贝尔图

当速度从零开始增加时，一些固有频率增加而另一些固有频率减少。当自然频率减小到0 时，出现临界转速。由于此模型的尺寸较小，导致在转速增加到3 000 rad/s时，还未出现第一临界转速。

2.3 时间步长收敛性研究

从1.3的模型求解部分，可以知道进行动态响应计算时，需要选择合适的时间步长，时间步长Δt的选择会影响到动态响应结果的收敛性。时间步长的选取依据为：所关注频段内的结构主要贡献的若干振型的最小周期（最大频率）的1/10～1/20[11]。设取第20 阶频率为最大频率，取其倒数的1/c2的值作为一种时间步。对应的时间步长区间为：Δt=4.518×10-6～9.036×10-6s；

同时，还需考虑齿轮啮合的频率，不同转速对应的时间步长为：

需要注意的是，取两者时间步长的较小值。现在转速为100 rad/s 的工况下，选取主动轮上的中间啮合点作为研究对象，c1取200，研究c2取值对收敛性的影响。表3为c2的取值对应的时间步长。

表3 时间步长取值

图8是时间步长收敛性结果示意图，其中局部放大图的横纵坐标单位与整体图的单位一致。从图中可以发现，在合理区间内，不同时间步长对模型的动态响应曲线影响很小，所得响应曲线几乎完全重合。说明在时间步长足够小的情况下，时间步长的选择已经不影响动态响应分析结果，从而验证了本模型在计算动态响应时时间步的收敛性很好。

图8 动态响应时间步长收敛性

2.4 刚柔模型动态特性对比

接下来，通过和传统的刚性扭振模型[13]进行动态特性结果对比，通过两种模型的差异性，阐述轮体柔性对齿轮副系统的影响。两种模型处于相同的工况之下，即载荷以及支承弹簧，时变啮合刚度，模型尺寸都相同。图9为不同转速下刚柔模型的稳定之后的动态响应曲线对比，其中柔性模型的响应曲线为主动轮上中间啮合点处的响应曲线，刚性模型的响应曲线为主动轮基圆半径上的位移响应曲线。

由图可知，两种模型的动态响应周期完全一致，这是由于两者的时变啮合刚度一致。刚性模型的动态响应曲线波动范围明显大于柔性模型，振动平衡位置也大于柔性模型。这是由于柔性模型中圆柱壳模型的柔性引起。在刚性模型中，仅考虑了啮合弹簧和支撑弹簧的柔性，因此在相同载荷下，外力做功全部转化成弹簧的势能；而在柔性模型中，外力做功一部分转化成为弹簧的势能，还有一部分转化为柔性圆柱壳的应变能。换而言之，圆柱壳的柔性对于齿轮副系统有减振吸能效果，因此柔性模型中啮合弹簧的势能小于刚性模型中啮合弹簧的势能，从而表现为刚性结果的平衡位置和振动波动范围都大于柔性模型。

图10展示了不同转速下刚柔模型的频谱成分，图10(a)为刚性模型，可以发现图中频率成分以一倍啮合频率和多倍啮合频率成分为主，当一倍频率和多倍频率的幅值线与一条平行于转速轴的直线相交时，频谱图的幅值将会达到峰值，该条直线对应的频率为该模型下的第1阶固有频率f1（扭转模型的第1、2 阶固有频率分别为1 156.5 Hz 和15 368 Hz），这种现象叫做参数共振，这是由于啮合频率是该模型的激振频率，当激励频率接近固有频率或接近固有频率的两倍时，就会发生这种共振。由于频谱图的关注范围为0～10 000 Hz，所以只出现了一个共振区域。除去参数共振区域附近，一倍频和高倍频对应的峰值随转速增加，总体上也呈现线性增加的趋势，且一倍频率的峰值明显大于多倍频率的峰值。

图10(b)为柔性模型，可以发现图中频率成分仍然以一倍啮合频率和高倍啮合频率成分为主，但是峰值随转速的变化波动很大。和刚性模型类似，柔性模型的频谱图的幅值也会出现参数共振现象，但是由于柔性模型使用的圆柱壳单元，该模型的第一阶固有频率为1 230.5 Hz。除去峰值附近区域，一倍频和高倍频对应的峰值随转速增加呈现波动上升的趋势，但是在高转速下，一倍频率的峰值并不占有明显优势。

图10 刚柔模型的频谱图

从两种模型频谱图的整体幅值量级上看，刚性模型的频谱幅值大于柔性模型的频谱幅值。图9的动态响应曲线也可以很好说明这个现象，在转速为10 rad/s和500 rad/s的转速下，刚性的振动范围都远大于柔性的波动范围，因此将时域曲线转换为频谱图时，就会出现刚性频谱幅值大于柔性频谱幅值。

3 结语

(1)本文基于圆柱壳理论，考虑时变啮合刚度以及弹性支撑的影响，建立了柔性齿轮副的等效简化模型。该模型具有适用范围广以及收敛性好的优点，丰富了现有的理论知识。

(2)基于该模型，开展了柔性齿轮副动力学特性验证和研究。推导了基于微分求积有限单元法的圆柱壳齿轮啮合模型的线性微分方程组。在确定轴向单元数的前提下，分别研究不同周向单元数以及单元节点数下模型的自然频率，讨论了模型的数值收敛性。接下来，与ANSYS的有限元模型对比自然频率，来验证本文模型的正确性。

(3)建立了和圆柱壳齿轮啮合副模型等效的集中质量刚性模型。通过两种模型的动态响应曲线、位移响应频谱图的对比，来研究齿轮柔性的影响。

(4)在相同的工况下，柔性模型和刚性模型的动态响应曲线周期相同，并且两种模型在同一转速下的响应曲线形状相似，但是柔性模型结果的平衡位置和振动波动范围都小于刚性模型。这是由于柔性模型中圆柱壳模型的柔性引起。在刚性模型中，仅考虑了啮合弹簧和支撑弹簧的柔性，因此在相同载荷下，外力做功全部转化成弹簧的势能；而在柔性模型中，外力做功一部分转化成为弹簧的势能，还有一部分转化为柔性圆柱壳的应变能。

(5)从位移响应频谱图可以看出，两种模型的频率成分都以一倍啮合频率和高倍啮合频率成分为主。但是由于圆柱壳模型的柔性，导致了其固有频率与刚性模型不一致，从而导致了参数共振区域的不一致，进而影响了一倍啮合频率和高倍啮合频率的幅值。并且由于齿轮柔性的影响，柔性模型频谱图的幅值整体上小于刚性模型的幅值。