基于接触有限元的斜齿轮副啮合动态响应及频谱特性分析

陈彦齐，黄修长,2，华宏星,2

（1.上海交通大学 振动冲击噪声研究所，上海 200240；2.上海交通大学 机械系统与振动国家重点实验室，上海 200240）

齿轮传动装置是机械传动系统中的关键设备之一，被广泛运用于航空、船舶等领域。随着实际运用要求的提高，对齿轮的振动噪声水平以及传动系统的动态特性提出了更高的要求。齿轮啮合过程中由于啮合冲击、转速波动、时变啮合刚度、输入波动、基础振动等因素具有强烈的非线性特征。发展准确的齿轮系统的建模及计算方法，研究齿轮系统的动态特性及影响规律，对明确齿轮系统的振动机理及齿轮系统减振降噪具有重要意义。

针对齿轮系统的较为常用的建模方法是集中参数法。Tuplin[1]提出了采用质量-弹簧模型模拟齿轮模型，以啮合刚度等效齿轮间的啮合作用，成为集中参数法建模的理论基础。王峰[2]基于集中质量法建立了人字齿轮的12自由度模型，模型中考虑了时变刚度、啮合冲击以及齿面摩擦等因素，研究结果表明，啮合刚度激励占振动的主要成分。集中参数法将轴与轴承的刚度直接等效为弹簧，未能考虑包含轴、轴承在内的转子系统的动力学特性。常乐浩[3]采用子结构法推导了平行轴齿轮传动系统通用建模方法，较之于前人工作添加了箱体单元，采用Newmark积分法与傅里叶级数法分别求解，验证了该建模方法的有效性以及相较于集中质量法建模的优越性。Ouyang等[4]建立了考虑时变刚度与接触非线性特性的齿轮-滚子-轴承系统动力学模型，结果表明齿轮传动系统产生的动态激励是高速条件下的主要振动源。

随着有限元商业软件的发展，逐步采用有限元软件如ABAQUS、ANSYS 等建立了完整的齿轮传动系统有限元模型，具有更高的准确度，在建模过程中不需人为考虑非线性因素，求解更为简便，并且能得到齿轮副在啮合过程中的动态接触力等实验难以测得的量。Chen等[5]通过有限元方法研究了行星齿轮外圈厚度对啮合刚度的影响，并与集中参数法结合求解得到齿轮系统的动态响应，发现会出现调制现象，但未对调制现象做出详细解释。Stoyanov等[6]通过ABAQUS建立了行星齿轮的有限元模型，并计算得到各齿轮中心的响应，表明接触有限元方法对了解系统动态特性的有效性。吴勇军等[7－8]基于ANSYS软件建立了齿轮-轴-轴承耦合系统的动力学模型，求解得到了系统的动力学特性，结果表明转子系统响应中不仅包含啮合频率及其倍频，还包含齿轮系统平均等效频率、弯曲振动固有频率等成分。张涛等[9]在齿轮副的基础上，引入轴与轴承，采用ABAQUS软件对人字齿轮进行动力学仿真，研究了系统阻尼、输入转速等因素对齿轮冲击啮合特性的影响。Ericson等[10]通过接触有限元对行星齿轮弹性体进行仿真，仿真结果与实验结果较为吻合，结果表明采用有限元方法能够获得集中参数法没有的模态，且弹性体振动与啮合力激励高度耦合。曹茂鹏等[11]采用有限元方法对直齿面齿轮动态啮合力进行仿真，得到了动态啮合力的时域变化结果，研究了转速等因素对啮合力的影响。上述研究大都集中于验证有限元软件方法的计算正确性，在时域方面开展齿轮系统的动态响应规律分析，在频域方面未对齿轮系统的频域特征做出详细分析，通常频率分辨率较低、频谱特征不够丰富、频谱图中仅有啮合频率及其倍频的峰值、边频带等特征被忽略的研究，不能详细分析齿轮系统的振动机理。

基于ABAQUS商业有限元软件，采用三维接触有限元，建立了考虑轴及轴承刚度的一对斜齿轮副-轴-轴承系统动态啮合有限元模型，发展对齿轮系统的一般性有限元分析方法。从多种工况下分析齿轮系统的动态响应，验证该方法的有效性，并且相较于现有研究结果，得到较高频率分辨率、较为丰富的频域特征。

1 有限元模型建立

对斜齿轮副采用Solidworks软件进行三维模型建模，斜齿轮副参数如表1 所示。然后进行有限元网格划分，为得到高质量纯六面体网格，将斜齿轮三维模型导入Hypermesh 进行网格划分。

表1 齿轮副参数

将上述网格模型导入到ABAQUS，进行有限元仿真设置和计算。轴-轴承参数及运行仿真参数如表2所示。采用647 992个C3D8单元对齿轮进行建模、采28个B31梁单元对轴进行模拟、采用4个赋予刚度属性的Wire 单元对轴承进行模拟。采用Dynamic/Implicit求解器进行求解，设置分析时间总长为1 s，时间增量步长设为自动调整，最小增量步长设置为1×10-8。将主动轮与从动轮设为接触对，采用面-面接触类型，齿轮接触选用有限滑动接触，并设置其接触属性的法向行为为硬接触，在此不考虑齿面间的摩擦作用，因此切向设置为无摩擦。

表2 轴-轴承参数及运行仿真参数

在齿轮内孔面与轴的中点建立运动耦合约束，以对主动轮及从动轮的边界条件及载荷进行控制；在从动轮中心添加扭矩，在主动轮中心添加转速，以模拟齿轮系统运转工况。为保证计算模型分析收敛，共设置2 个分析步，第一个分析步固定从动轮，主动轮转动0.05 rad，保证两个齿轮充分接触；第二个分析步释放从动轮的自由度，如前所述添加主动轮角速度与从动轮扭矩。

根据释放从动轮自由度的个数，共划分3 种工况：

(1)从动轮轴与从动轮内孔面耦合的节点只释放绕轴旋转自由度；

(2)对从动轮轴与从动轮内孔面耦合的节点不作约束，此时靠轴端的弹簧单元约束系统的自由度；

(3)对从动轮轴与从动轮内孔面耦合的节点不作约束，并且在此节点上添加大小为1 000 N的轴向静推力。其中第一种工况齿轮旋转中心只有旋转自由度，是较为理想的工况，模拟轴、箱体、基座的刚度都较大的情况；第二种工况是一般性工况，考虑轴、轴承的刚度；第三种工况是研究轴向静推力负载下齿轮系统动力学特性。系统有限元模型如图1所示。

图1 斜齿轮副-轴-轴承系统有限元模型

2 动态响应结果及分析

计算完成后，导出大齿轮中心的角速度、角加速度和齿轮副的动态接触力。齿轮副在啮合点处的加速度可作为评估齿轮系统的动态性能指标之一[10]，其定义为：

式中：rp、分别为主动轮的基圆半径与角加速度；rg、分别为从动轮的基圆半径与角加速度。仿真中，主动轮的角加速度为0，因此根据从动轮加速度即可评估系统的动态性能。

主动轮转频为20.0 Hz，从动轮转频为7.7 Hz，系统的啮合频率为340 Hz，理论角速度为48.550 5 rad/s，理想平稳状态下齿轮的法向接触力可用式(2)进行计算：

式中：Tm为负载扭矩，mn为法向模数；zg为从动轮齿数；αt为端面压力角。将仿真模型参数代入公式进行计算，得到理论法向接触力为2 419 N。

2.1 从动轮中心仅释放旋转自由度工况

从动轮中心仅释放旋转自由度工况下从动轮的角速度、角加速度及齿轮副动态接触力时域结果如图2所示。由图2可知，由于时变啮合刚度等非线性因素存在，时域结果存在明显周期性波动与峰值，角速度的均值与均方根值都为48.550 8 rad/s。峰值的间隔周期约为0.05 s，与主动轮的转动周期较为接近，伴随明显的冲击现象。角加速度的均方根值为70.692 0 rad/s2，动态接触力的均方根值为2 453.6 N。相较平稳状态理论值增加1.4%，表明在从动轮中心仅释放旋转自由度工况下，齿轮系统的非线性因素对于齿轮动态接触力的影响不大。

图2 从动轮中心仅释放旋转自由度工况下从动轮时域结果

对信号进行去均值处理之后采用傅里叶变换将时域结果变换至频域，频率分辨率为1 Hz，得到角速度、角加速度、动态接触力的脉动量频谱结果如图3所示。可见，主要成分为主动轮转频20 Hz 及其倍频；频谱的突出峰值为340 Hz、680 Hz、640 Hz、720 Hz 等频率，其中340 Hz 是该齿轮系统的啮合频率，680 Hz 为啮合频率二倍频，其余频率皆为主动轮转动频率的多倍频。在从动轮中心仅释放旋转自由度工况下，响应频率主要由主动轮转频及其多倍频与啮合频率及其倍频组成。

2.2 从动轮中心不作约束工况

从动轮中心不作约束工况下的角速度、角加速度、齿轮副法向接触力时域结果如图4所示。

将图4 与图2 比较，当齿轮中心浮动，引入轴的刚度与轴承的刚度，角速度、角加速度及动态接触力没有呈现明显的冲击现象，依旧保持明显的周期特性伴随着波动幅值增加。在此工况下，角速度的均方根值为48.556 3 rad/s，与理论角速度相差0.005 8 rad/s，但中心固定时的角速度均方根值仅为0.000 3 rad/s。角加速度的均方根值增至143.695 4 rad/s2，反映系统的振动加剧。动态接触力均方根值增至2 677.5 N，动态接触力相较于静态接触力或中心固定工况下的接触力增加了约258 N，增加了10%。3个物理量均能反映在引入轴刚度与轴承刚度之后，齿轮系统的振动加剧。

图3 从动轮中心只释放旋转自由度工况下从动轮频域结果

图4 从动轮中心不作约束工况下从动轮响应时域结果

将时域结果进行傅里叶变换，频域结果如图5所示。为与中心只有旋转自由度的工况进行对比，将幅值范围与图3 统一，其中啮合频率附近的幅值较大，此部分的频谱图如框线中的频谱所示。由图5(a)可知，300 Hz 以内存在主动轮转频20 Hz 及其多倍频。与图3(a)相比，400 Hz 及以上的频率成分减少，啮合频率及其附近的幅值显著增大。图3(a)中啮合频率处的幅值大小为0.005，图5(a)中啮合频率处的幅值增至0.03。此现象也同样存在于角加速度与动态接触力的频谱中，较高频段处的幅值减小，啮合频率及其附近频段的幅值显著增大。

除此之外，对比图3 与图5 可知，考虑系统的轴刚度及轴承刚度时，在啮合频率处的幅值大幅增加同时，会出现以从动轮转频为调制频率的边频带。如图5(a)中的332 Hz是由啮合频率340 Hz减去从动轮转频所得，其幅值为0.01；353 Hz是由啮合频率加上从动轮二倍转频，其幅值为0.008；此现象同样也存在于角加速度与法向接触力的频谱中。

图5 从动轮中心不作约束工况下从动轮响应频域结果

丁康等[12]发现齿轮啮合过程中，会出现以啮合频率及其倍频为载波频率、齿轮的转动频率为调制频率的调制现象并对其进行定性分析。可从幅值调制与频率调制两方面进行分析。从幅值调制的角度认为齿轮系统的动态响应如式(3)所示：

其中：Am是第m阶幅值，fz是啮合频率，[1+am(t)]是与齿轮转动频率有关的调制信号，可以被表示为am(t)=βmcos(2πmfnt)，其中fn为从动轮的转频。利用积化和差公式将式(3)的第m阶展开可得：

式(4)表明，除了啮合频率成分外，还会出现频率为fz+fn与fz-fn的成分。

杨小青[13]提出齿轮动态响应的调频现象可从齿轮的转速波动对啮合频率的影响进行分析，将齿轮副的啮合刚度以啮合频率为基频进行傅里叶展开：

其中：K0为平均啮合刚度，fz为啮合频率，ki和αi为第i阶分量的幅值与相位。较为普遍的集中参数模型建模方法中假设啮合频率是一个在稳定转速下计算出来的恒定值，但是由于齿轮实际上为弹性体，齿轮在转动过程中包含摩擦、冲击等非线性因素，如图2、图4 所示，从动轮转速会存在波动。将转速的波动成分以傅里叶级数形式表示，波动转速表示如下：

其中：n0为理论转速，fn为从动轮转频，考虑转速波动时啮合频率计算公式为：

将式(6)与式(7)代入式(5)，可得到考虑波动转速的啮合频率表达式：

将式(8)的第i阶分量的相位记为Φi(t)，则Φi(t)表达如式(9)所示，对其相位进行求导，如式(10)所示，即可得到频率表达式：

由式(10)发现考虑转速波动时，除理论的啮合频率以外，啮合刚度会出现与转频相关的频率成分。通过对调制现象进行定性分析，说明了仿真结果中出现调制现象的合理性，同时也证明采用有限元方法能够得到更详细的频率特征。

2.3 从动轮中心不作约束且添加轴向静推力工况

从动轮中心浮动且在从动轮中心添加轴向1 000 N 静推力工况下的角速度、角加速度、齿轮副法向接触力时域结果如图6所示。添加轴向静推力之后的时域响应相较于没有添加轴向静推力的结果没有明显差异，呈现周期性波动特征。其角速度均方根值为48.557 4 rad/s，与理论角速度相差0.006 9 rad/s；其角加速度均方根值为153.968 5 rad/s2，相较于不添加轴向静推力的工况增加10.273 1 rad/s2；动态接触力均方根值增至2 680.8 N，相较于不添加轴向静推力工况增加约3 N，可以忽略不计。

图6 从动轮中心不作约束且添加轴向静推力工况下从动轮响应时域结果

3种工况下各指标的时域结果比较如表3所示，从表中可以看出，引入轴及轴承的刚度之后，系统的振动加剧，各指标的均方根值都有所增加，其中角加速度与动态接触力的增加较为明显；添加1 000 N轴向静推力对系统的动态响应特性影响较小。

表3 3种工况下各指标均值、均方值比较

将时域结果通过傅里叶变换至频域，如图7 所示。从图7 可知，添加轴向力的频谱图与未添加轴向力的频谱图较为接近。与未添加轴向力的频谱相比较，3个指标的频率特征没有明显变化，保持以啮合频率为主，啮合频率附近出现以从动轮转频及其倍频为调制频率的调制现象。

图7 从动轮中心不作约束且添加轴向静推力工况下从动轮响应频域结果

时域结果以及频域结果都能表明，在从动轮中心添加大小为1 000 N 的静推力对于齿轮系统的动态响应影响很小。

3 结语

基于接触有限元方法，建立了齿轮副-轴-轴承耦合系统，发展了更精确地分析齿轮系统动态响应及特性的方法，研究了从动轮中心仅释放旋转自由度、从动轮中心不作约束、从动轮中心不作约束且添加轴向静推力3 种工况下的动态响应，得出了以下结论：

(1)从动轮中心仅释放旋转自由度工况下的系统角速度、角加速度、动态接触力的波动要小于其余两种工况，反映了在考虑轴、轴承的刚度之后系统的振动加剧，且存在较为明显的周期性冲击现象；添加1 000 N 轴向静推力与未添加轴向静推力相比对系统的动态响应影响不大；

(2)从动轮中心仅释放旋转自由度工况下的频谱特征以主动轮转频及其倍频为主；对从动轮中心不作约束时，频谱特征呈现以啮合频率载波频率，以从动轮转频及其倍频为调制频率的调制现象，并从幅值调制和频率调制两方面定性分析了频谱出现边频带的原因。