CN103971006B - 一种考虑主减速器壳的驱动桥齿轮动力学特性确定方法 - Google Patents

Abstract

本发明涉及一种考虑主减速器壳的驱动桥齿轮动力学特性计算方法，包括以下步骤：1)定义全局坐标系；2)建立轴有限元模型；3)建立非线性轴承模型；4)建立齿轮力学模型；5)建立主减速器壳缩维有限元模型；6)建立驱动桥主减速器齿轮传动系统静力学模型；7)计算静力平衡时的轴承刚度；8)建立驱动桥主减速器齿轮传动系统动力学模型；9)驱动桥主减速器齿轮传动系统固有振动特性计算；10)驱动桥主减速器齿轮动力学特性计算。本发明以有限元方法和模态综合方法建立包含主减速器壳在内的齿轮传动系统动力学模型，考虑主减速器壳与传动系的相互影响，能够准确高效地对考虑主减速器壳影响的驱动桥齿轮动力学特性进行计算。本发明可广泛应用于各种包含壳体的齿轮传动结构的动力学特性计算分析。

Description

一种考虑主减速器壳的驱动桥齿轮动力学特性确定方法

技术领域

本发明涉及一种车辆传动领域中的零部件动力学特性计算方法，特别是关于一种考虑主减速器壳的驱动桥齿轮动力学特性计算方法。

背景技术

在驱动桥中，主减速器齿轮的传动误差是驱动桥系统振动噪声的主要激励源，一方面传动误差会引起齿轮自身的啸叫问题，另一方面齿轮传动误差产生的动态激励经过传动轴、轴承传递至壳体，传动系与壳体相互影响，形成整个系统的振动噪声问题。驱动桥主减速器齿轮传动系统动力学分析的难点在于，如何建立包含主减速器壳和传动系各个部件的准确完整的系统动力学分析模型，进行准确高效的计算分析。

现有研究方法在进行驱动桥主减速器齿轮传动系统动力学分析时，大多采用简化的集中参数模型，将系统部件处理为集中质量，并将轴承处理为弹簧，一端与轴的集中质量点相连，一端接地，无法考虑主减速器壳的影响。这种方法虽然计算效率高，但系统模型过于简化，无法体现传动系部件的尺寸特征，尤其没有考虑主减速器壳的影响，无法准确的体现系统的动力学特性。另有研究方法直接建立包含驱动桥所有零部件的实体单元有限元模型，采用有限元接触计算方法分析驱动桥系统的动力学特性，但是这种方法的建模过程十分复杂，系统模型的规模过大，接触分析需要消耗大量的计算资源，计算效率低。

发明内容

针对上述问题，本发明的目的是提供一种准确高效的考虑主减速器壳的驱动桥齿轮动力学特性计算方法。

为实现上述目的，本发明采取以下技术方案：一种考虑主减速器壳的驱动桥齿轮动力学特性计算方法，包括以下步骤：1)定义全局坐标系：对驱动桥主减速器齿轮传动系统的全局坐标系进行定义，作为系统建模的基础；2)建立轴有限元模型：采用考虑剪切应变的欧拉伯努利空间梁单元建立轴的有限元模型；3)建立非线性轴承模型：采用具有耦合非线性刚度特性的轴承单元对滚子轴承进行模拟；4)建立齿轮力学模型：在齿轮等效啮合节点之间建立沿等效啮合力作用线方向的空间弹簧单元，对准双曲面齿轮的啮合关系进行模拟；5)建立主减速器壳缩维有限元模型：采用四节点四面体单元对主减速器壳的几何模型进行网格划分，建立主减速器壳的有限元模型，采用模态综合法求得主减速器壳的缩维刚度矩阵和缩维质量矩阵；6)建立驱动桥主减速器齿轮传动系统静力学模型：用轴承非线性刚度矩阵组集传动系有限元模型刚度矩阵和主减 速器壳缩维刚度矩阵，获得完整传动系统的静力学有限元模型；7)计算静力平衡时的轴承刚度：采用牛顿-拉普森方法迭代求解对应输入转矩下的系统静力学方程，得到静力平衡时的轴承刚度矩阵；8)建立驱动桥主减速器齿轮传动系统动力学模型：以静力平衡时的轴承刚度矩阵作为轴承的线性刚度矩阵，组集传动系有限元模型刚度矩阵和主减速器壳缩维刚度矩阵，以得到整体系统的线性刚度矩阵，以轴的质量矩阵和主减速器壳的缩维质量矩阵，组集得到整体系统的质量矩阵，建立系统动力学模型；9)驱动桥主减速器齿轮传动系统固有振动特性计算：求解系统无阻尼自由振动方程的特征方程的特征根和特征向量，获得系统的振动频率和正则振型；10)驱动桥主减速器齿轮动力学特性计算：先计算齿轮的动态柔度和动态刚度，再计算单位谐波齿轮传动误差激励下系统的动态响应，获得齿轮啮合节点自由度的响应，由齿轮的动态刚度和齿轮啮合节点自由度的响应，求得齿轮的动态啮合力。

在所述步骤1)中，驱动桥主减速器齿轮传动系统采用标准汽车坐标系作为全局坐标系，即汽车前进方向为x轴正向，汽车左侧方向为y轴正向，竖直向上为z轴正向，坐标原点为差速器十字轴中心位置。

在所述步骤2)中，考虑剪切应变的欧拉伯努利空间梁单元是指在经典欧拉伯努利梁单元模型中，以剪切影响系数形式来引入剪切应变的影响，圆截面的剪切影响系数表示为：

上式中，E为轴材料的弹性模量；I为梁单元的截面惯性矩；G为轴材料的剪切模量；L为梁单元的长度；A为梁单元的截面面积。

在所述步骤3)中，滚子轴承的载荷计算公式表示为：

上式中，F_x和F_y分别为沿x方向和y方向的径向力；F_z为沿z方向的轴向力；M_x和M_y分别为绕x轴和y轴的力矩；K_n为滚子与内外圈的综合接触刚度；n_s为每个滚 子沿长度方向上划分的单元数；Z为滚子数；δ_j,k为第j个滚子的第k个单元的法向变形量；α为接触角，对于圆柱滚子轴承为零；ψ_j为第j个滚子的方位角；D_pw为滚子的节圆直径；x_k为每个滚子第k个单元中心与滚子中心的距离；

其中，K_n的计算公式表示为：

上式中，E为轴承材料的弹性模量；ν为轴承材料的泊松比；L_we为滚子的有效长度；

δ_j,k的计算公式表示为：

上式中，δ_x和δ_y分别为轴承内外圈之间沿x轴和y轴的径向位移；δ_z为轴承内外圈之间沿z轴的轴向位移；θ_x和θ_y分别为轴承内外圈之间绕x和y轴的角位移；计算时，若δ_j,k小于0，则取δ_j,k等于0。

在所述步骤4)中，准双曲面齿轮中小轮和大轮的理论啮合位置在全局坐标系中的坐标表示为：

(x,y,z)＝(x₁,y₁,z₁)+(Δx₁,Δy₁,Δz₁)

上式中，(x₁,y₁,z₁)为小轮中心节点在全局坐标系中的坐标；Δx₁、Δy₁和Δz₁分别为理论啮合位置相对小轮中心节点坐标在全局坐标系x、y和z方向上的偏移量；

其中，(Δx₁,Δy₁,Δz₁)的计算公式表示为：

上式中，γ₁和γ₂分别为小轮和大轮的节锥角；E_g为准双曲面齿轮的偏置距；R_m2为大轮的平均节圆半径；k₁和k₂分别为小轮和大轮的朝向系数，当小轮朝前时，k₁取1，反之k₁取-1，当大轮朝右时，k₂取1，反之k₂取-1；k_h为小轮的旋向系数，当小轮右旋时，k_h取1，反之k_h取-1。

在所述步骤4)中，全局坐标系下小轮和大轮等效啮合节点之间的等效啮合力作用线方向向量(x_n,y_n,z_n)的计算公式表示为：

上式中，α₁为工作齿面的平均节圆压力角；β₁为小轮平均节圆螺旋角；k_l为齿轮工作状态系数，当小轮凹面和大轮凸面为工作齿面时，k_l取1，当小轮凸面和大轮凹面为工作齿面时，k_l取-1；θ₁为中间参数；

其中，θ₁的计算公式表示为：

上式中，γ₁和γ₂分别为小轮和大轮的节锥角。

在所述步骤6)中，主减速器齿轮传动系统的静力学方程表示为：

[K]{δ}＝{f}

上式中，[K]为系统刚度矩阵，由梁单元刚度矩阵[K_beam]、非线性轴承刚度矩阵[K_bearing]、齿轮刚度矩阵[K_gear]、齿轮啮合刚度矩阵[K_mesh]和主减速器壳缩维刚度矩阵[K_housing]组集而成；{δ}为节点自由度位移向量；{f}为外载荷向量。

在所述步骤8)中，驱动桥主减速器齿轮传动系统模型的动力学方程表示为：

上式中，{δ_d}为节点自由度时域位移向量；{f_d}为动态载荷向量；[K]为系统刚度矩阵，组集方式与静力学相同，其中轴承的刚度矩阵为对应工况静力平衡时的线性刚度；[M]为系统质量矩阵，由梁单元质量矩阵[M_beam]和主减速器壳缩维质量矩阵[M_housing]组集而成；[C]为系统阻尼矩阵。

在所述步骤10)中，主减速器齿轮的动力学特性通过齿轮的动态啮合力体现，齿轮的动态啮合力表示为：

{F_mesh}＝D_meshδ_mesh

上式中，δ_mesh为在激振力作用下小轮和大轮啮合节点沿齿轮啮合力作用线方向上的相对位移响应，D_mesh为齿轮沿啮合力作用线方向上的动态刚度，由小轮和大轮的动态柔度决定；

其中，D_mesh的计算公式表示为：

D_mesh＝[C_p+C_g]^-1

上式中，C_p和C_g分别为小轮和大轮的动态柔度，即单位谐波力激励下齿轮啮合节点在啮合力作用线方向上的位移响应幅值；

其中，小轮和大轮的动态柔度C_p和C_g分别表示为：

C_p＝{δ₀}_p{x_n,y_n,z_n}^T

C_g＝{δ₀}_g{x_n,y_n,z_n}^T

上式中，{δ₀}_p和{δ₀}_g分别为{δ₀}中对应的小轮和大轮啮合节点自由度的位移响应；{x_n,y_n,z_n}^T为步骤4)中小轮和大轮等效啮合节点之间的等效啮合力作用线方向向量；

上式中的{δ₀}可由下式求得：

上式中，ω_0i和{φ₀}_i为不考虑准双曲面齿轮啮合刚度耦合时第i阶振动频率和正则振型；n为保留的模态阶数；{F(t)}为小轮和大轮啮合节点上施加的沿啮合力作用线方向单位谐波力向量；λ_0i为第i阶频率比；ξ_0i为第i阶模态阻尼比；

其中，λ_0i的计算公式表示为：

上式中，ω为齿轮单位谐波传动误差的激振频率。

在所述步骤10)中，激振力作用下小轮和大轮啮合节点沿齿轮啮合力作用线方向上的相对位移响应δ_mesh的求解过程如下：

将齿轮单位谐波传动误差引起的激振力施加在齿轮耦合传动系统中，分别作用在小轮和大轮的啮合节点上，采用振型叠加法计算获得系统的位移响应表示为：

上式中，ω_i和{φ}_i为考虑准双曲面齿轮啮合刚度耦合时第i阶振动频率和正则振型；n为保留的模态阶数；{f(t)}为齿轮单位谐波传动误差引起的激振力向量；λ_i为第i阶频率比；ξ_i为第i阶模态阻尼比；

其中，λ_i的计算公式表示为：

由{δ_d}中对应的小轮和大轮啮合节点自由度的位移响应{δ_d}_p和{δ_d}_g，进一步计算出小轮和大轮啮合节点沿齿轮啮合力作用线方向上的相对位移响应δ_mesh：

{δ_mesh}＝|{δ_d}_p-{δ_d}_g|{x_n,y_n,z_n}^T

上式中，{δ_d}_p和{δ_d}_g分别为齿轮单位谐波传动误差激励下的小轮和大轮啮合节点自由度的位移响应；{x_n,y_n,z_n}^T为步骤4)中小轮和大轮等效啮合节点之间的等效啮合力作用线方向向量。

本发明由于采取以上技术方案，其具有以下优点：1、本发明采用了Harris(哈里斯)等人提出的考虑轴承刚度耦合性和非线性的轴承单元，保证了轴承计算模型的准确性。2、本发明以有限元方法和模态综合方法建立包含主减速器壳在内的齿轮传动系统动力学模型，考虑主减速器壳与传动系的相互影响，能够准确高效地对考虑主减速器壳影响的驱动桥齿轮动力学特性进行计算。3、本发明所采用的方法基于经典的非线性轴承理论、有限元方法和模态综合方法，具有可靠的理论基础，且易于在各类常用的编程语言环境下编程实现，具有较高的计算效率。本发明可广泛应用于各种包含壳体的齿轮传动结构的动力学特性计算分析。

附图说明

图1是本发明方法的流程示意图；

图2是驱动桥主减速器齿轮传动系统示意图；

图3是齿轮柔度频响特性曲线图；

图4是齿轮动态啮合力频响特性曲线图。

图2中：1、大轮轴；2、差速器轴；3、主减速器壳轴承座；4、圆锥滚子轴承；5、主减速器壳轴承座；6、圆锥滚子轴承；7、主减速器壳轴承座；8、圆锥滚子轴承9、输入转矩10、小轮轴；11、小轮轮齿；12、主减速器壳轴承座；13、圆柱滚子轴承；14、圆锥滚子轴承；15、主减速器壳轴承座；16、输出转矩；17、大轮轮齿。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

如图1所示，本发明提供的考虑主减速器壳的驱动桥齿轮动力学特性计算方法包括以下步骤：

1)定义全局坐标系：对驱动桥主减速器齿轮传动系统的全局坐标系进行定义，作为系统建模的基础。

在本实施例中，驱动桥主减速器齿轮传动系统采用标准汽车坐标系作为全局坐标系，即汽车前进方向为x轴正向，汽车左侧方向为y轴正向，竖直向上为z轴正向，坐标原点为差速器十字轴中心位置。

2)建立轴有限元模型：采用考虑剪切应变的欧拉伯努利空间梁单元建立轴的有限元模型。

考虑剪切应变的欧拉伯努利空间梁单元是指在经典欧拉伯努利梁单元模型中，以剪切影响系数形式来引入剪切应变的影响，圆截面的剪切影响系数表示为：

上式中，E为轴材料的弹性模量；I为梁单元的截面惯性矩；G为轴材料的剪切模量；L为梁单元的长度；A为梁单元的截面面积。

3)建立非线性轴承模型：采用具有耦合非线性刚度特性的轴承单元对滚子轴承进行模拟。

在本实施例中，轴承受力与变形的计算关系参考了Harris(哈里斯)等人的研究(文献为罗继伟，马伟等译，滚动轴承分析，机械工业出版社，2010年)。滚子轴承的载荷计算公式表示为：

上式中，F_x和F_y分别为沿x方向和y方向的径向力；F_z为沿z方向的轴向力；M_x和M_y分别为绕x轴和y轴的力矩；K_n为滚子与内外圈的综合接触刚度；n_s为每个滚子沿长度方向上划分的单元数；Z为滚子数；δ_j,k为第j个滚子的第k个单元的法向变形量；α为接触角，对于圆柱滚子轴承为零；ψ_j为第j个滚子的方位角；D_pw为滚子的节圆直径；x_k为每个滚子第k个单元中心与滚子中心的距离；

其中，K_n的计算公式表示为：

上式中，E为轴承材料的弹性模量；ν为轴承材料的泊松比；L_we为滚子的有效长度；

δ_j,k的计算公式表示为：

上式中，δ_x和δ_y分别为轴承内外圈之间沿x轴和y轴的径向位移；δ_z为轴承内外圈之间沿z轴的轴向位移；θ_x和θ_y分别为轴承内外圈之间绕x和y轴的角位移；计算时，若δ_j,k小于0，则取δ_j,k等于0。

在本实施例中，计算得到的非线性轴承刚度矩阵[K_bearing]表示为：

上式中，轴承刚度矩阵对应的坐标系为轴承局部坐标系。

4)建立齿轮力学模型：在齿轮等效啮合节点之间建立沿等效啮合力作用线方向的空间弹簧单元，对准双曲面齿轮的啮合关系进行模拟。

在本实施例中，小轮和大轮的理论啮合位置在全局坐标系中的坐标表示为：

(x,y,z)＝(x₁,y₁,z₁)+(Δx₁,Δy₁,Δz₁)

上式中，(x₁,y₁,z₁)为小轮中心节点在全局坐标系中的坐标；Δx₁、Δy₁和Δz₁分别为理论啮合位置相对小轮中心节点坐标在全局坐标系x、y和z方向上的偏移量；

其中，(Δx₁,Δy₁,Δz₁)的计算公式表示为：

上式中，γ₁和γ₂分别为小轮和大轮的节锥角；E_g为准双曲面齿轮的偏置距；R_m2为大轮的平均节圆半径；k₁和k₂分别为小轮和大轮的朝向系数，当小轮朝前时，k₁取1，反之k₁取-1，当大轮朝右时，k₂取1，反之k₂取-1；k_h为小轮的旋向系数，当小轮右旋 时，k_h取1，反之k_h取-1。

在本实施例中，在小轮和大轮中心到各自的等效啮合节点之间分别建立的刚性梁单元同样为考虑剪切应变的欧拉伯努利梁单元，弹性模量取为一般钢材的10⁶倍，通过该齿轮刚性梁单元将齿轮啮合产生的载荷传递到轴上。

在本实施例中，全局坐标系下小轮和大轮等效啮合节点之间的等效啮合力作用线方向向量(x_n,y_n,z_n)的计算公式表示为：

上式中，α₁为工作齿面的平均节圆压力角；β₁为小轮平均节圆螺旋角；k_l为齿轮工作状态系数，当小轮凹面和大轮凸面为工作齿面时，k_l取1，当小轮凸面和大轮凹面为工作齿面时，k_l取-1；θ₁为中间参数；

其中，θ₁的计算公式表示为：

在本实施例中，全局坐标系下的小轮和大轮啮合节点之间空间弹簧单元的啮合刚度矩阵[K_mesh]表示为：

上式中，k_m为齿轮等效啮合刚度系数；{h}为方向向量；

其中，{h}的计算公式表示为：

{h}＝{x_n y_n z_n 0 0 0}

上式中，{x_n,y_n,z_n}为全局坐标系下小轮和大轮等效啮合节点之间的等效啮合力作用线方向向量。

5)建立主减速器壳缩维有限元模型：采用四节点四面体单元对主减速器壳的几何模型进行网格划分，建立主减速器壳的有限元模型，采用模态综合法求得主减速器壳的缩维刚度矩阵和缩维质量矩阵。

由于主减速器壳有限元模型包含大量节点自由度，无法直接与传动系模型连接，需要进行缩维处理。本实施例中采用模态综合法对壳体的有限元模型进行缩维，将主减速器壳有限元模型中的节点自由度分为边界自由度o和内部自由度i，主减速器壳有限元模型的边界自由度为轴承座与轴承外圈连接节点的自由度，分别在轴承中心位置建立节点，用刚性单元与轴承座的相关节点耦合，即将主减速器壳轴承座与轴承的连 接节点自由度体现在每个轴承中心节点上，定义每个轴承中心节点为边界节点，每个边界节点包含6个自由度，由于在实际驱动桥系统中，主减速器壳通过螺栓固定在桥壳上，因此在对主减速器壳有限元模型进行缩维计算时，约束主减速器壳与桥壳连接法兰面上的节点自由度，利用有限元软件采用模态综合法求得主减速器壳的缩维刚度矩阵[K_housing]和缩维质量矩阵[M_housing]。

6)建立驱动桥主减速器齿轮传动系统静力学模型：用轴承非线性刚度矩阵组集传动系有限元模型刚度矩阵和主减速器壳缩维刚度矩阵，获得完整传动系统的静力学有限元模型。

在本实施例中，将主减速器壳轴承中心连接位置的边界节点自由度与传动系梁单元对应的轴承内圈安装节点自由度之间用滚子轴承刚度矩阵[K_bearing]连接，即实现主减速器壳有限元模型与传动系有限元模型的组装。主减速器齿轮传动系统的静力学方程表示为：

[K]{δ}＝{f}

上式中，[K]为系统刚度矩阵，由梁单元刚度矩阵[K_beam]、非线性轴承刚度矩阵[K_bearing]、齿轮刚度矩阵[K_gear]、齿轮啮合刚度矩阵[K_mesh]和主减速器壳缩维刚度矩阵[K_housing]组集而成；{δ}为节点自由度位移向量；{f}为外载荷向量。

7)计算静力平衡时的轴承刚度：采用牛顿-拉普森方法迭代求解对应输入转矩下的系统静力学方程，得到静力平衡时的轴承刚度矩阵。

在本实施例中，采用牛顿-拉普森方法对系统静力学方程迭代求解时，采用相邻两次迭代所得的节点位移向量之差的模小于给定小量作为收敛准则。

8)建立驱动桥主减速器齿轮传动系统动力学模型：以静力平衡时的轴承刚度矩阵作为轴承的线性刚度矩阵，组集传动系有限元模型刚度矩阵和主减速器壳缩维刚度矩阵，以得到整体系统的线性刚度矩阵，以轴的质量矩阵和主减速器壳的缩维质量矩阵，组集得到整体系统的质量矩阵，建立系统动力学模型。

系统模型的动力学方程表示为：

上式中，{δ_d}为节点自由度时域位移向量；{f_d}为动态载荷向量；[K]为系统刚度矩阵，组集方式与静力学相同，其中轴承的刚度矩阵为对应工况静力平衡时的线性刚度；[M]为系统质量矩阵，由梁单元质量矩阵[M_beam]和主减速器壳缩维质量矩阵[M_housing]组集而成；[C]为系统阻尼矩阵。

9)驱动桥主减速器齿轮传动系统固有振动特性计算：求解系统无阻尼自由振动方程的特征方程的特征根和特征向量，获得系统的振动频率和正则振型。

其中，系统无阻尼自由振动方程表示为：

10)驱动桥主减速器齿轮动力学特性计算：先计算齿轮的动态柔度和动态刚度，再计算单位谐波齿轮传动误差激励下系统的动态响应，获得齿轮啮合节点自由度的响应，由齿轮的动态刚度和齿轮啮合节点自由度的响应，求得齿轮的动态啮合力。

在本实施例中，首先将齿轮啮合关系解耦，即齿轮的等效啮合刚度系数k_m取0，分别在小轮和大轮啮合节点上施加沿啮合力作用线方向的单位谐波力，由于系统动力学模型为线性模型，由模态叠加法可知，正则振型下的位移响应表示为：

上式中，ω_0i和{φ₀}_i为不考虑准双曲面齿轮啮合刚度耦合时第i阶振动频率和正则振型；n为保留的模态阶数；{F(t)}为小轮和大轮啮合节点上施加的沿啮合力作用线方向单位谐波力向量；λ_0i为第i阶频率比；ξ_0i为第i阶模态阻尼比；

其中，λ_0i的计算公式表示为：

上式中，ω为齿轮单位谐波传动误差的激振频率。

在本实施例中，由{δ₀}计算小轮和大轮的动态柔度C_p和C_g分别表示为：

C_p＝{δ₀}_p{x_n,y_n,z_n}^T

C_g＝{δ₀}_g{x_n,y_n,z_n}^T

上式中，{δ₀}_p和{δ₀}_g分别为{δ₀}中对应的小轮和大轮啮合节点自由度的位移响应；{x_n,y_n,z_n}^T为步骤4)中小轮和大轮等效啮合节点之间的等效啮合力作用线方向向量；

进一步求得齿轮副的动态刚度D_mesh表示为：

D_mesh＝[C_p+C_g]^-1

上式中，C_p和C_g分别为小轮和大轮的动态柔度；

在本实施例中，将齿轮单位谐波传动误差引起的激振力施加在齿轮耦合传动系统中，分别作用在小轮和大轮的啮合节点上，同样采用振型叠加法计算获得系统的位移响应表示为：

上式中，ω_i和{φ}_i为考虑准双曲面齿轮啮合刚度耦合时第i阶振动频率和正则振型；n为保留的模态阶数；{f(t)}为齿轮单位谐波传动误差引起的激振力向量；λ_i为第i阶频率比；ξ_i为第i阶模态阻尼比；

其中，λ_i的计算公式表示为：

在本实施例中，由{δ_d}中对应的小轮和大轮啮合节点自由度的位移响应{δ_d}_p和{δ_d}_g，进一步计算出小轮和大轮啮合节点沿齿轮啮合力作用线方向上的相对位移响应δ_mesh表示为：

{δ_mesh}＝|{δ_d}_p-{δ_d}_g|{x_n,y_n,z_n}^T

上式中，{δ_d}_p和{δ_d}_g分别为齿轮单位谐波传动误差激励下的小轮和大轮啮合节点自由度的位移响应；{x_n,y_n,z_n}^T为步骤4)中小轮和大轮等效啮合节点之间的等效啮合力作用线方向向量；

在本实施例中，主减速器齿轮的动力学特性通过齿轮的动态啮合力体现，齿轮的动态啮合力表示为：

{F_mesh}＝D_mesh{δ_mesh}

上式中，{δ_mesh}为在齿轮单位谐波传动误差激励下的小轮和大轮啮合节点沿齿轮啮合力作用线方向上的相对位移响应，D_mesh为齿轮副沿啮合力作用线方向上的动态刚度。

下面通过一个具体的实施例，用以说明本发明的效果。如图2所示，某型驱动桥主减速器齿轮传动系统由一对准双曲面齿轮副、若干滚子轴承6、8、13、4、14、差速器轴2和主减速器壳组成。小轮轴10由一对圆锥滚子轴承6、8和一个圆柱滚子轴承13支撑，大轮轴1安装在差速器轴2上，差速器轴2由一对圆锥滚子轴承4、14支撑，所有滚子轴承的外圈均安装在主减速器壳轴承座3、5、7、12、15上。在小轮轴10端部施加输入转矩9，在差速器轴的差速器中心位置输出转矩16。

1)定义全局坐标系：主减速器齿轮传动系统采用标准汽车坐标系作为全局坐标系，汽车前进方向为x轴正向，汽车左侧方向为y轴正向，竖直向上为z轴正向，坐标原点为差速器十字轴中心位置。

2)建立轴有限元模型：如图2所示，传动系统包含3根轴：小轮轴10、大轮轴1和差速器轴2。小轮轴上共有23个节点和22个梁单元，小轮轴10材料的弹性模量为200GPa，泊松比为0.252，密度为7880kg/m³，小轮轴10的有限元模型节点编号自左 向右依次为1-23，将各梁单元的刚度矩阵组集得到小轮轴的整体刚度矩阵。大轮轴1上共有7个节点和6个梁单元，大轮轴1材料的弹性模量为207GPa，泊松比为0.29，密度为7800kg/m³，大轮轴1的有限元模型节点编号自左向右依次为24-30，将各梁单元的刚度矩阵组集得到大轮轴1的整体刚度矩阵。差速器轴2上共有21个节点和20个梁单元，差速器轴2材料的弹性模量为207GPa，泊松比为0.29，密度为7800kg/m³，差速器轴2的有限元模型节点编号自左向右依次为31-51，将各梁单元的刚度矩阵组集得到差速器轴2的整体刚度矩阵。实际中，大轮轴1是通过螺栓安装在差速器轴2上，在模型中将大轮轴1与差速器轴2螺栓连接位置对应的节点之间建立刚度较大的线性弹簧单元，用来模拟实际的螺栓连接关系。

3)建立非线性轴承模型：如图2所示，传动系统包含4个圆锥滚子轴承6、8、4、14和1个圆柱滚子轴承13。小轮轴10上的后圆柱滚子轴承13型号为FAG575867，轴承内径为40mm，外径为94mm，宽度为30mm，平均直径为67mm，滚子数为13，滚子直径为16mm，滚子有效长度为19mm，内圈对应连接在小轮轴10上的第2号节点。小轮轴10上的中圆锥滚子轴承6型号为FAG546439，轴承内径为70mm，外径为165mm，宽度为57mm，平均直径为117.5mm，滚子数为15，滚子直径为22.6mm，滚子有效长度为39.556mm，轴承接触角为25°，内圈对应连接在小轮轴上的第8号节点。小轮轴10上的前圆锥滚子轴承8型号为FAG31312，轴承内径为60mm，外径为130mm，宽度为33.5mm，平均直径为95mm，滚子数为16，滚子直径为17.18mm，滚子有效长度为19.8mm，轴承接触角为28.81°，内圈对应连接在小轮轴10上的第13号节点。差速器轴2上的左圆锥滚子轴承4型号为FAG32021，轴承内径为105mm，外径为160mm，宽度为35mm，平均直径为132.5mm，滚子数为28，滚子直径为13.4mm，滚子有效长度为23.48mm，轴承接触角为16.5°，内圈对应连接在差速器轴上的第33号节点。差速器轴2上的右圆锥滚子轴承14型号为SKF33021，轴承内径为105mm，外径为160mm，宽度为43mm，平均直径为132.5mm，滚子数为28，滚子直径为13.74mm，滚子有效长度为29.76mm，轴承接触角为10.67°，内圈对应连接在差速器轴2上的第50号节点。轴承材料的弹性模量均为210GPa，泊松比均为0.3。根据滚子轴承的非线性刚度计算公式，计算得到每个轴承各自的非线性刚度矩阵。

4)建立准双曲面齿轮力学模型：准双曲面齿轮参数如表1所示，小轮中心对应小轮轴10上的第5号节点，大轮中心对应大轮轴1上的第29号节点，小轮啮合点的节点编号为52，大轮啮合点的节点编号为53。计算得到的小轮和大轮理论啮合节点在全局坐标系中的坐标为(181.695mm，-5.054mm，-24.1mm)，计算得到的等效啮合力作用线方向向量为(0.7183，-0.1992，0.6666)。根据准双曲面齿轮的力学模型计算公式 计算得到齿轮单元刚度矩阵和等效啮合刚度矩阵。

表1准双曲面齿轮参数

5)建立主减速器壳缩维有限元模型：采用四节点四面体单元对主减速器壳的几何模型进行网格划分，单元尺寸为4mm，主减速器壳的材料为QT450-10，弹性模量为169GPa，泊松比为0.257，密度为7060kg/m³，整个主减速器壳有限元模型共包含641057个单元，154040个节点。主减速器壳有限元模型的边界自由度为轴承座与轴承外圈连接节点的自由度，分别在5个轴承中心位置建立主节点，用RBE2刚性单元与轴承座的相关节点耦合，即将主减速器壳轴承座3、5、7、12、15与轴承外圈的连接节点自由度体现在这5个轴承中心节点上，定义这5个轴承中心节点为边界节点，编号为54-58，每个边界节点包含6个自由度，因此整个主减速器壳有限元模型有30个边界自由度。在实际驱动桥系统中，主减速器壳通过螺栓固定在桥壳上，因此在主减速器壳有限元模型中，约束主减速器壳与桥壳连接法兰面上的节点自由度。利用有限元软件Nastran采用模态综合法进行模型缩维计算时保留60阶主模态，保证获得的主减速器壳缩维刚度矩阵和质量矩阵具有足够的精度，能够准确体现主减速器壳的模态特性。

6)建立驱动桥主减速器齿轮传动系统静力学模型：基于轴的刚度矩阵和主减速器壳的缩维刚度矩阵，用小轮刚性梁单元刚度矩阵连接5号小轮中心节点和52号小轮啮合节点，用大轮刚性梁单元刚度矩阵连接29号大轮中心节点和53号大轮啮合节点，用齿轮啮合刚度矩阵连接52号小轮啮合节点和53号大轮啮合节点，用后圆柱滚子轴承刚度矩阵连接小轮轴2号节点和主减速器壳54号边界节点，用中圆锥滚子轴承刚度矩阵连接小轮轴8号节点和主减速器壳55号边界节点，用前圆锥滚子轴承刚度矩阵连 接小轮轴13号节点和主减速器壳56号边界节点，用左圆锥滚子轴承刚度矩阵连接差速器轴33号节点和主减速器壳57号边界节点，用右圆锥滚子轴承刚度矩阵连接差速器轴50号节点和主减速器壳58号边界节点，最终组集得到完整的系统刚度矩阵。在小轮轴上的第15号节点位置施加输入转矩，大小为7297.3Nm，约束差速器轴上的第43号节点的轴向转动自由度。最终建立得到整个系统的静力学有限元模型，共包含58个节点，每个节点有6个自由度，同时还包含60阶主减速器壳模态自由度，因此系统总自由度数为408。

7)计算静力平衡时的轴承刚度：采用牛顿-拉普森方法迭代求解系统静力学方程，以相邻两次迭代所得的节点位移向量之差的模小于给定小量作为收敛准则，收敛容差取为10^-6，不满足则继续进行迭代，满足则判断迭代过程收敛。实际计算时迭代7次计算收敛，获得该静力平衡时的轴承刚度矩阵。

8)建立驱动桥主减速器齿轮传动系统动力学模型：系统动力学有限元模型的刚度矩阵与静力学模型中的刚度矩阵相同，其中轴承的刚度矩阵为对应载荷静力平衡时的线性刚度矩阵，质量矩阵为轴的质量矩阵与壳体的缩维质量矩阵的组集，组集方法与刚度矩阵的组集方法相同，阻尼以模态阻尼系数的形式体现，在计算时模态阻尼系数取为0.01，动态载荷为齿轮啮合传动误差产生的激励载荷，最终建立得到系统动力学方程。

9)驱动桥主减速器齿轮传动系统的振动模态分析：采用Cholesky(乔里斯基)分解法求解系统无阻尼自由振动方程的特征方程的特征根和特征向量，获得系统的振动频率和正则振型。通常关心2000Hz以内的齿轮激振，当考虑主减速器壳时，2000Hz以内的系统耦合振动模态频率如表2所示。当不考虑主减速器壳时，2000Hz以内的系统耦合振动模态频率如表3所示。

表2考虑主减速器壳时2000Hz以内的系统耦合振动模态频率

模态阶数	频率/Hz
		1	320.3
2	798.48
		3	800.71
4	889.93
		5	945.06
6	1029.2
		7	1102
8	1470.7
		9	1578.5
10	1967.4
		11	1993.5

表3不考虑主减速器壳时2000Hz以内的系统耦合振动模态频率

模态阶数	频率/Hz
		1	696.56
2	1395.9
		3	1406.7
4	1681.9
		5	1684.5
6	1929.6

10)驱动桥主减速器齿轮动力学特性计算：首先计算齿轮在非啮合状态下在单位谐波力激励下的系统动态响应，获得齿轮的动态柔度，频响特性曲线如图3所示。再计算齿轮啮合状态下在单位谐波传动误差激励下的系统动态响应，最终计算得到齿轮动态啮合力的频响特性曲线如图4所示。整个计算过程可通过Matlab程序实现，通用性强，采用双核CPU主频2.6GHz，4G内存的计算机进行求解的时间不超过1分钟，计算效率高。

通过上述方法，最终得到了驱动桥主减速器齿轮动力学特性，图4给出了采用本发明方法获得的准双曲面齿轮动态啮合力频响特性曲线(实线)与不考虑主减速器壳时的准双曲面齿轮动态啮合力频响特性曲线(虚线)的对比情况，通过对比可知如果不考虑主减速器壳的影响，系统固有振动特性计算结果误差很大，会丢失大量的动态 特性信息，分析得到的动力学特性与考虑主减速器壳时的偏差很大，峰值频率和峰值力均不准确。基于齿轮动态啮合力频响特性曲线，可以在驱动桥主减速器齿轮传动系统的设计和使用过程中避免产生啮合力峰值，以减小系统的振动，同时可以对系统结构的改进提供指导依据。

综上所述，本发明适合于驱动桥主减速器齿轮动力学特性计算，由于考虑了主减速器壳的影响，避免了未考虑壳体的集中参数模型无法真实反映系统的动力学特性的缺点和建立整个系统的实体单元有限元模型进行接触分析时计算量过大，效率低的问题，本发明能够准确高效地实现对驱动桥主减速器齿轮传动系统的动力学建模和主减速器齿轮动力学特性的计算分析。

上述各实施例仅用于对本发明的目的、技术方案和有益效果进行了进一步详细说明，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种考虑主减速器壳的驱动桥齿轮动力学特性确定方法，包括以下步骤：

1)定义全局坐标系：对驱动桥主减速器齿轮传动系统的全局坐标系进行定义，作为系统建模的基础；

2)建立轴有限元模型：采用考虑剪切应变的欧拉伯努利空间梁单元建立轴的有限元模型；

3)建立非线性轴承模型：采用具有耦合非线性刚度特性的轴承单元对滚子轴承进行模拟；

4)建立齿轮力学模型：在齿轮等效啮合节点之间建立沿等效啮合力作用线方向的空间弹簧单元，对准双曲面齿轮的啮合关系进行模拟；

5)建立主减速器壳缩维有限元模型：采用四节点四面体单元对主减速器壳的几何模型进行网格划分，建立主减速器壳的有限元模型，采用模态综合法获得主减速器壳的缩维刚度矩阵和缩维质量矩阵；

6)建立驱动桥主减速器齿轮传动系统静力学模型：用轴承非线性刚度矩阵组集传动系有限元模型刚度矩阵和主减速器壳缩维刚度矩阵，获得完整传动系统的静力学有限元模型；

7)计算静力平衡时的轴承刚度：采用牛顿-拉普森方法迭代求解对应输入转矩下的系统静力学方程，得到静力平衡时的轴承刚度矩阵；

8)建立驱动桥主减速器齿轮传动系统动力学模型：以静力平衡时的轴承刚度矩阵作为轴承的线性刚度矩阵，组集传动系有限元模型刚度矩阵和主减速器壳缩维刚度矩阵，以得到整体系统的线性刚度矩阵，以轴的质量矩阵和主减速器壳的缩维质量矩阵，组集得到整体系统的质量矩阵，建立系统动力学模型；

9)驱动桥主减速器齿轮传动系统固有振动特性计算：求解系统无阻尼自由振动方程的特征方程的特征根和特征向量，获得系统的振动频率和正则振型；

10)驱动桥主减速器齿轮动力学特性计算：先计算齿轮的动态柔度和动态刚度，再计算单位谐波齿轮传动误差激励下系统的动态响应，获得齿轮啮合节点自由度的响应，由齿轮的动态刚度和齿轮啮合节点自由度的响应，求得齿轮的动态啮合力。

2.如权利要求1所述的一种考虑主减速器壳的驱动桥齿轮动力学特性确定方法，其特征在于：在所述步骤1)中，驱动桥主减速器齿轮传动系统采用标准汽车坐标系作为全局坐标系，即汽车前进方向为x轴正向，汽车左侧方向为y轴正向，竖直向上为z轴正向，坐标原点为差速器十字轴中心位置。

3.如权利要求1所述的一种考虑主减速器壳的驱动桥齿轮动力学特性确定方法，其特征在于：在所述步骤2)中，考虑剪切应变的欧拉伯努利空间梁单元是指在经典欧拉伯努利梁单元模型中，以剪切影响系数形式来引入剪切应变的影响，圆截面的剪切影响系数表示为：

$\mathrm{\φ}=\frac{40EI}{3{\mathrm{GL}}^{2}A}$

上式中，E为轴材料的弹性模量；I为梁单元的截面惯性矩；G为轴材料的剪切模量；L为梁单元的长度；A为梁单元的截面面积。

4.如权利要求1所述的一种考虑主减速器壳的驱动桥齿轮动力学特性确定方法，其特征在于：在所述步骤3)中，滚子轴承的载荷计算公式表示为：

$\left\{\begin{array}{c}{F}_x=-\frac{K_n}{n_s}\underset{j=1}{\overset{Z}{\Σ}}\[\underset{k=1}{\overset{n_s}{\Σ}}\begin{array}{c}\left({\mathrm{\δ}}_{j,k}^{10/9}cos\mathrm{\α}sin{\mathrm{\ψ}}_j\right)\end{array}\]\\ F_y=\frac{K_n}{n_s}\underset{j=1}{\overset{Z}{\Σ}}\[\underset{k=1}{\overset{n_s}{\Σ}}\begin{array}{c}\left({\mathrm{\δ}}_{j,k}^{10/9}cos\mathrm{\α}\cos {\mathrm{\ψ}}_j\right)\end{array}\]\\ F_z=\frac{K_n}{n_s}\underset{j=1}{\overset{Z}{\Σ}}\[\underset{k=1}{\overset{n_s}{\Σ}}\begin{array}{c}\left({\mathrm{\δ}}_{j,k}^{10/9}\sin \mathrm{\α}\right)\end{array}\]\\ M_x=\frac{K_n}{n_s}\underset{j=1}{\overset{Z}{\Σ}}\[\underset{k=1}{\overset{n_s}{\Σ}}\left(\right(\frac{D_{pw}}{2}\sin \mathrm{\α}-x_k\left){\mathrm{\δ}}_{j,k}^{10/9}{\mathrm{cos\ψ}}_j\right)\]\\ M_y=\frac{K_n}{n_s}\underset{j=1}{\overset{Z}{\Σ}}\[\underset{k=1}{\overset{n_s}{\Σ}}\left(\right(\frac{D_{pw}}{2}\sin \mathrm{\α}-x_k\left){\mathrm{\δ}}_{j,k}^{10/9}{\mathrm{sin\ψ}}_j\right)\]\end{array}\right.$

上式中，F_x和F_y分别为沿x方向和y方向的径向力；F_z为沿z方向的轴向力；M_x和M_y分别为绕x轴和y轴的力矩；K_n为滚子与内外圈的综合接触刚度；n_s为每个滚子沿长度方向上划分的单元数；Z为滚子数；δ_j,k为第j个滚子的第k个单元的法向变形量；α为接触角，对于圆柱滚子轴承为零；ψ_j为第j个滚子的方位角；D_pw为滚子的节圆直径；x_k为每个滚子第k个单元中心与滚子中心的距离；

其中，K_n的计算公式表示为：

$K_n=\mathrm{\π}\frac{E}{1-v^2}\frac{L_{we}^{8/9}}{{14.22}^{10/9}}$

上式中，E为轴承材料的弹性模量；ν为轴承材料的泊松比；L_we为滚子的有效长度；

δ_j,k的计算公式表示为：

$\begin{array}{c}{\mathrm{\δ}}_{j,k}=\[{\mathrm{\δ}}_z+\frac{D_{pw}}{2}({\mathrm{\θ}}_x{\mathrm{cos\ψ}}_j+{\mathrm{\θ}}_y{\mathrm{sin\ψ}}_j)\]\sin \mathrm{\α}+(-{\mathrm{\δ}}_x{\mathrm{sin\ψ}}_j+{\mathrm{\δ}}_y{\mathrm{cos\ψ}}_j)\cos \mathrm{\α}\\ +x_k(-{\mathrm{\θ}}_x{\mathrm{cos\ψ}}_j-{\mathrm{\θ}}_y{\mathrm{sin\ψ}}_j)\end{array}$

上式中，δ_x和δ_y分别为轴承内外圈之间沿x轴和y轴的径向位移；δ_z为轴承内外圈之间沿z轴的轴向位移；θ_x和θ_y分别为轴承内外圈之间绕x和y轴的角位移；计算时，若δ_j,k小于0，则取δ_j,k等于0。

5.如权利要求1所述的一种考虑主减速器壳的驱动桥齿轮动力学特性确定方法，其特征在于：在所述步骤4)中，准双曲面齿轮中小轮和大轮的理论啮合位置在全局坐标系中的坐标表示为：

(x,y,z)＝(x₁,y₁,z₁)+(Δx₁,Δy₁,Δz₁)

上式中，(x₁,y₁,z₁)为小轮中心节点在全局坐标系中的坐标；Δx₁、Δy₁和Δz₁分别为理论啮合位置相对小轮中心节点坐标在全局坐标系x、y和z方向上的偏移量；

其中，(Δx₁,Δy₁,Δz₁)的计算公式表示为：

$\left\{\begin{array}{c}\mathrm{\Δ}{x}_1=0\\ \mathrm{\Δ}{y}_1=(\frac{{\mathrm{sin\γ}}_2}{\sqrt{{\cos }^2{\mathrm{\γ}}_1-{\sin }^2{\mathrm{\γ}}_2}}{E}_g-R_{m2}tan{\mathrm{\γ}}_2)k_2\\ {\mathrm{\Δz}}_1=(E_g-\frac{\sqrt{{\cos }^2{\mathrm{\γ}}_1-{\sin }^2{\mathrm{\γ}}_2}}{{\mathrm{cos\γ}}_2}{R}_{m2})k_h{k}_1{k}_2\end{array}\right.$

上式中，γ₁和γ₂分别为小轮和大轮的节锥角；E_g为准双曲面齿轮的偏置距；R_m2为大轮的平均节圆半径；k₁和k₂分别为小轮和大轮的朝向系数，当小轮朝前时，k₁取1，反之k₁取-1，当大轮朝右时，k₂取1，反之k₂取-1；k_h为小轮的旋向系数，当小轮右旋时，k_h取1，反之k_h取-1。

6.如权利要求1所述的一种考虑主减速器壳的驱动桥齿轮动力学特性确定方法，其特征在于：在所述步骤6)中，主减速器齿轮传动系统的静力学方程表示为：

[K]{δ}＝{f}

上式中，[K]为系统刚度矩阵，由梁单元刚度矩阵[K_beam]、非线性轴承刚度矩阵[K_bearing]、齿轮刚度矩阵[K_gear]、齿轮啮合刚度矩阵[K_mesh]和主减速器壳缩维刚度矩阵[K_housing]组集而成；{δ}为节点自由度位移向量；{f}为外载荷向量。

7.如权利要求1所述的一种考虑主减速器壳的驱动桥齿轮动力学特性确定方法，其特征在于：在所述步骤8)中，驱动桥主减速器齿轮传动系统模型的动力学方程表示为：

$\[M\]\left\{{\stackrel{\·\·}{\mathrm{\δ}}}_d\right\}+\[C\]\left\{{\stackrel{\·}{\mathrm{\δ}}}_d\right\}+\[K\]\left\{{\mathrm{\δ}}_d\right\}=\left\{f_d\right\}$

上式中，{δ_d}为节点自由度时域位移向量；{f_d}为动态载荷向量；[K]为系统刚度矩阵，组集方式与静力学相同，其中轴承的刚度矩阵为对应工况静力平衡时的线性刚度；[M]为系统质量矩阵，由梁单元质量矩阵[M_beam]和主减速器壳缩维质量矩阵[M_housing]组集而成；[C]为系统阻尼矩阵。

8.如权利要求1所述的一种考虑主减速器壳的驱动桥齿轮动力学特性确定方法，其特征在于：在所述步骤10)中，主减速器齿轮的动力学特性通过齿轮的动态啮合力体现，齿轮的动态啮合力表示为：

{F_mesh}＝D_meshδ_mesh

上式中，δ_mesh为在激振力作用下小轮和大轮啮合节点沿齿轮啮合力作用线方向上的相对位移响应，D_mesh为齿轮沿啮合力作用线方向上的动态刚度，由小轮和大轮的动态柔度决定；

其中，D_mesh的计算公式表示为：

D_mesh＝[C_p+C_g]^-1

上式中，C_p和C_g分别为小轮和大轮的动态柔度，即单位谐波力激励下齿轮啮合节点在啮合力作用线方向上的位移响应幅值；

其中，小轮和大轮的动态柔度C_p和C_g分别表示为：

C_p＝{δ₀}_p{x_n,y_n,z_n}^T

C_g＝{δ₀}_g{x_n,y_n,z_n}^T

上式中，{δ₀}_p和{δ₀}_g分别为{δ₀}中对应的小轮和大轮啮合节点自由度的位移响应；{x_n,y_n,z_n}^T为步骤4)中小轮和大轮等效啮合节点之间的等效啮合力作用线方向向量；

上式中的{δ₀}可由下式求得：

$\left\{{\mathrm{\δ}}_0\right\}=\underset{i=1}{\overset{n}{\Σ}}\frac{{\left\{{\mathrm{\φ}}_0\right\}}_i^T\left\{F\left(t\right)\right\}{\left\{{\mathrm{\φ}}_0\right\}}_i}{{\mathrm{\ω}}_{0i}^2\[1-{\mathrm{\λ}}_{0i}^2+2{\mathrm{j\ξ}}_{0i}{\mathrm{\λ}}_{0i}\]}$

上式中，ω_0i和{φ₀}_i为不考虑准双曲面齿轮啮合刚度耦合时第i阶振动频率和正则振型；n为保留的模态阶数；{F(t)}为小轮和大轮啮合节点上施加的沿啮合力作用线方向单位谐波力向量；λ_0i为第i阶频率比；ξ_0i为第i阶模态阻尼比；

其中，λ_0i的计算公式表示为：

${\mathrm{\λ}}_{0i}=\frac{\mathrm{\ω}}{{\mathrm{\ω}}_{0i}}$

上式中，ω为齿轮单位谐波传动误差的激振频率。

9.如权利要求8所述的一种考虑主减速器壳的驱动桥齿轮动力学特性确定方法，其特征在于：在所述步骤10)中，激振力作用下小轮和大轮啮合节点沿齿轮啮合力作用线方向上的相对位移响应δ_mesh的求解过程如下：

将齿轮单位谐波传动误差引起的激振力施加在齿轮耦合传动系统中，分别作用在小轮和大轮的啮合节点上，采用振型叠加法计算获得系统的位移响应表示为：

$\left\{{\mathrm{\δ}}_d\right\}=\underset{i=1}{\overset{n}{\Σ}}\frac{{\left\{\mathrm{\φ}\right\}}_i^T\left\{f\left(t\right)\right\}{\left\{\mathrm{\φ}\right\}}_i}{{\mathrm{\ω}}_i^2\[1-{\mathrm{\λ}}_i^2+2{\mathrm{j\ξ}}_i{\mathrm{\λ}}_i\]}$

上式中，ω_i和{φ}_i为考虑准双曲面齿轮啮合刚度耦合时第i阶振动频率和正则振型；n为保留的模态阶数；{f(t)}为齿轮单位谐波传动误差引起的激振力向量；λ_i为第i阶频率比；ξ_i为第i阶模态阻尼比；

其中，λ_i的计算公式表示为：

${\mathrm{\λ}}_i=\frac{\mathrm{\ω}}{{\mathrm{\ω}}_i};$

由{δ_d}中对应的小轮和大轮啮合节点自由度的位移响应{δ_d}_p和{δ_d}_g，进一步计算出小轮和大轮啮合节点沿齿轮啮合力作用线方向上的相对位移响应δ_mesh：

{δ_mesh}＝|{δ_d}_p-{δ_d}_g|{x_n,y_n,z_n}^T

上式中，{δ_d}_p和{δ_d}_g分别为齿轮单位谐波传动误差激励下的小轮和大轮啮合节点自由度的位移响应；{x_n,y_n,z_n}^T为步骤4)中小轮和大轮等效啮合节点之间的等效啮合力作用线方向向量。