CN103995940A - 一种考虑输入转矩变化的驱动桥动力学特性计算方法 - Google Patents

Abstract

本发明涉及一种考虑输入转矩变化的驱动桥动力学特性计算方法，包括以下步骤：包括以下步骤：1)定义驱动桥全局坐标系；2)建立轴系梁单元模型；3)建立非线性滚子轴承模型；4)建立齿轮模型；5)建立连接部件模型；6)建立壳体有限元模型及缩维模型；7)建立完整的驱动桥系统动力学模型；8)计算不同输入转矩工况下的轴承刚度；9)计算不同输入转矩工况下的驱动桥系统动力学特性。本发明以有限元方法和模态综合方法建立包含主减速器总成、差速器总成、轮毂总成和壳体等部件的完整驱动桥系统动力学模型，能够准确高效地计算驱动桥系统的动力学特性。

Description

一种考虑输入转矩变化的驱动桥动力学特性计算方法

技术领域

本发明涉及一种车辆传动领域中的零部件动力学特性计算方法，特别是关于一种考虑输入转矩变化的驱动桥动力学特性计算方法。

背景技术

驱动桥的振动噪声是汽车噪声的主要来源，其主减速器锥齿轮的传动误差是造成驱动桥振动噪声问题的根本原因。驱动桥是由轴、轴承、齿轮和壳体组成的动力学系统，齿轮的传动误差会产生动态激励，经过轴-轴承传递至壳体，形成整个系统的振动噪声问题。建立准确高效的驱动桥系统动力学模型是研究驱动桥系统动力学特性的有效方法，可以大幅缩短驱动桥产品的研发周期，节约成本。在驱动桥系统中，滚子轴承是连接轴系与壳体的关键部件，由于滚子和轴承内外圈通过接触作用传递载荷，轴承刚度具有各向耦合性和非线性特性，会随系统输入转矩的不同而变化，进而影响系统的动力学特性。

现有研究方法在进行驱动桥动力学计算时，大多采用简化的集中参数模型，将系统部件处理为集中质量，并将轴承处理为弹簧，一端与轴的集中质量点相连，一端接地，轴承刚度为常数，没有考虑轴承刚度受驱动桥输入转矩变化的影响，从而无法体现输入转矩变化对驱动桥系统动力学特性的影响。

发明内容

针对上述问题，本发明的目的是提供一种准确高效的考虑输入转矩变化的驱动桥动力学特性计算方法。

为实现上述目的，本发明采取以下技术方案：一种考虑输入转矩变化的驱动桥动力学特性计算方法，包括以下步骤：1)定义驱动桥全局坐标系：对驱动桥系统的全局坐标系进行定义，作为系统建模的基础；2)建立轴系梁单元模型：采用考虑剪切应变的欧拉伯努利空间梁单元对轴系部件进行模拟；3)建立非线性滚子轴承模型：采用具有耦合非线性刚度特性的轴承单元对滚子轴承进行模拟；4)建立齿轮模型：采用等效啮合模型对主减速器锥齿轮和轮间差速器齿轮进行建模，对齿轮的载荷传递进行模拟；5)建立连接部件模型：对驱动桥系统中的连接部件进行模拟；6)建立壳体有限元模型及缩维模型：建立驱动桥系统中所包含的壳体结构有限元模型，并对壳体的有限元模型进行缩维计算；7)建立完整的驱动桥系统动力学模型：将轴系模型与壳体缩维模型的相关刚度矩阵和质量矩阵按照节点自由度耦合关系组集，获得系统刚度矩阵和系统质量矩阵，建立包含主减速器总成、差速器总成、轮毂总成和桥壳在内的完整驱动桥系统动力学模型；8)计算不同输入转矩工况下的轴承刚度：根据驱动桥输入转矩大小的不同定义轻载、中载和重载工况，采用牛顿-拉普森方法迭代求解不同输入转矩下的系统静力学方程，得到静力平衡时的轴承刚度矩阵；9)计算不同输入转矩工况下的驱动桥系统动力学特性：在不同输入转矩工况下，采用模态叠加法计算单位谐波齿轮传动误差激励下系统的动力学特性。

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

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

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

上式中，取轴承的轴线方向为轴承局部坐标系z轴的方向；F_x和F_y分别为沿x方向和y方向的径向力；F_z为沿z方向的轴向力；M_x、M_y和M_z分别为绕x轴、y轴和z轴的力矩；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_{\mathrm{we}}^{8/9}}{{14.22}^{10/9}}$

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

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

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

上式中，δ_x和δ_y分别为轴承内外圈之间沿x轴和y轴的径向变形；δ_z为轴承内外圈之间沿z轴的轴向变形；θ_x和θ_y分别为轴承内外圈之间绕x和y轴的角变形；计算时，对轴承滚子单元的受力状态进行判断，若δ_j,k小于0，则取δ_j,k等于0，表示该轴承滚子单元未发生接触变形，不传递载荷。

在所述步骤5)中，所述连接部件在建模时均采用线性弹簧单元模拟，耦合相互作用的节点自由度，全局坐标系下线性弹簧单元的刚度矩阵K_c表示为：

K_c＝diag([k_x k_y k_z k_rx k_ry k_rz])

上式中，k_x为沿x轴方向的平动刚度；k_y为沿y轴方向的平动刚度；k_z为沿z轴方向的平动刚度；k_rx为绕x轴方向的转动刚度；k_ry为绕y轴方向的转动刚度；k_rz为绕z轴方向的转动刚度。

在所述步骤6)中，采用模态综合法对壳体的有限元模型进行缩维变换，根据壳体的连接特点，将有限元模型中的节点自由度分为边界自由度和内部自由度，将与轴系连接的节点自由度定义为边界自由度，其他节点自由度定义为内部自由度；将壳体有限元模型的刚度矩阵和质量矩阵通过模态综合法转换到边界自由度和模态自由度上，在保留原结构模态特性的同时实现对模型的缩维，采用固定截面模态综合方法求得壳体缩维刚度矩阵K_h和壳体缩维质量矩阵M_h。

在所述步骤7)中，所建立的驱动桥系统动力学模型为：

$M\stackrel{..}{\mathrm{\δ}}+C\stackrel{.}{\mathrm{\δ}}+\mathrm{K\δ}=f$

上式中，δ为节点自由度时域位移向量；f为动态载荷向量；M为系统质量矩阵，由梁单元质量矩阵M_s、壳体缩维质量矩阵M_h组集而成；C为系统阻尼矩阵；K为系统刚度矩阵，由梁单元刚度矩阵K_s、滚子轴承刚度矩阵K_b、准双曲面齿轮啮合刚度矩阵K_m、差速器齿轮刚度矩阵K_d和啮合刚度矩阵K_dm、壳体缩维刚度矩阵K_h和连接部件刚度矩阵K_c组集而成。

在所述步骤8)中，系统静力学非线性方程表示为：

K₀(δ₀,f₀(T))δ₀＝f₀(T)

上式中，K_b为求得的轴承刚度矩阵；T为输入转矩；f₀(T)为输入转矩工况对应的系统载荷向量，δ₀为系统位移向量。

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

F_mesh＝D_meshδ_mesh

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

D_mesh的计算公式为：

D_mesh＝[C_p+C_g]^-1

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

δ_mesh的求解方法如下：

采用振型叠加法计算系统动态响应：

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

上式中，δ为系统在单位谐波传动误差激励下的系统位移响应；f(t)为单位谐波传动误差引起的激振载荷；ω_i和φ_i为考虑准双曲面齿轮啮合刚度耦合时第i阶振动频率和正则振型；n为保留的模态阶数；λ_i＝ω/ω_i为第i阶频率比，ω为激振频率；ξ_i为第i阶模态阻尼比；

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

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

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

本发明由于采取以上技术方案，其具有以下优点：1、本发明采用了具有各向耦合性和非线性刚度特性的轴承单元，能够体现轴承刚度随驱动桥系统输入转矩变化的特性。2、本发明以有限元方法和模态综合方法建立包含主减速器总成、差速器总成、轮毂总成和壳体等部件的完整驱动桥系统动力学模型，能够准确高效地计算驱动桥系统的动力学特性。3、本发明所采用的方法基于经典的非线性轴承理论、有限元方法和模态综合方法，具有可靠的理论基础，且易于在各类常用的编程语言环境下编程实现，具有较高的计算效率。4、采用本发明中的方法可以求得驱动桥在对应工况下的动力学特性，能够为驱动桥的工况选取和减振降噪设计提供指导。

附图说明

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

图2是驱动桥系统传动系统示意图；

图3是小轮轴有限元模型示意图；

图4是半轴有限元模型示意图；

图5是十字轴有限元模型示意图；

图6是行星轮轴有限元模型示意图；

图7是太阳轮轴有限元模型示意图；

图8是圆锥滚子轴承示意图；

图9是不同输入转矩工况下的齿轮动态啮合力频响特性曲线图；

图2中各附图标记的含义如下：1、小轮轴；2、输入转矩；3、圆锥滚子轴承；4、圆锥滚子轴承；5、圆柱滚子轴承；6、半轴；7、轮毂；8、圆锥滚子轴承；9、圆锥滚子轴承；10、圆锥滚子轴承；11、太阳轮轴；12、行星轮轴；13、大轮；14、差速器壳；15、十字轴；16、圆锥滚子轴承。

具体实施方式

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

如图1所示，本发明提供的考虑输入转矩变化的驱动桥动力学特性计算方法包括以下步骤：

1)定义驱动桥全局坐标系：对驱动桥系统的全局坐标系进行定义，作为系统建模的基础。

在本实施例中，采用汽车的车辆坐标系作为系统的全局坐标系，即车辆前方为x轴正向，左侧为y轴正向，竖直向上为z轴正向，坐标原点为差速器十字轴中心位置。

2)建立轴系梁单元模型：采用考虑剪切应变的欧拉伯努利空间梁单元对轴系部件进行模拟。

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

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

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

在本实施例中，驱动桥的轴部件包含小轮轴、左右半轴、十字轴、四个行星轮轴和左右太阳轮轴，这些轴部件均采用考虑剪切变形的欧拉伯努利空间梁单元模拟。根据部件的尺寸特征和装配特点划分单元，参照部件的设计图纸定义梁单元几何参数，梁单元材料参数定义为部件的真实材料参数，确保轴模型准确。将小轮位置的梁单元几何参数定义为小轮节锥几何参数，以体现小轮的刚度和质量属性。将行星轮轴和太阳轮轴的梁单元几何参数定义为差速器齿轮的节锥参数，以体现差速器齿轮的刚度和质量属性。梁单元的质量矩阵采用考虑剪切变形的一致质量矩阵。

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

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

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

上式中，取轴承的轴线方向为轴承局部坐标系z轴的方向；F_x和F_y分别为沿x方向和y方向的径向力；F_z为沿z方向的轴向力；M_x、M_y和M_z分别为绕x轴、y轴和z轴的力矩；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_{\mathrm{we}}^{8/9}}{{14.22}^{10/9}}$

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

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

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

上式中，δ_x和δ_y分别为轴承内外圈之间沿x轴和y轴的径向变形；δ_z为轴承内外圈之间沿z轴的轴向变形；θ_x和θ_y分别为轴承内外圈之间绕x和y轴的角变形；计算时，对轴承滚子单元的受力状态进行判断，若δ_j,k小于0，则取δ_j,k等于0，表示该轴承滚子单元未发生接触变形，不传递载荷。

在本实施例中，计算得到的轴承刚度矩阵K_b表示为：

$K_b=\left[\begin{array}{cccccc}\frac{{\∂F}_x}{{\∂\mathrm{\δ}}_x}& \frac{{\∂F}_x}{{\∂\mathrm{\δ}}_y}& \frac{{\∂F}_x}{{\∂\mathrm{\δ}}_z}& \frac{{\∂F}_x}{{\∂\mathrm{\θ}}_x}& \frac{{\∂F}_x}{{\∂\mathrm{\θ}}_y}& 0\\ \frac{{\∂F}_y}{{\∂\mathrm{\δ}}_x}& \frac{{\∂F}_y}{{\∂\mathrm{\δ}}_y}& \frac{{\∂F}_y}{{\∂\mathrm{\δ}}_z}& \frac{{\∂F}_y}{{\∂\mathrm{\θ}}_x}& \frac{{\∂F}_y}{{\∂\mathrm{\θ}}_y}& 0\\ \frac{{\∂F}_z}{{\∂\mathrm{\δ}}_x}& \frac{{\∂F}_z}{{\∂\mathrm{\δ}}_y}& \frac{{\∂F}_z}{{\∂\mathrm{\δ}}_z}& \frac{{\∂F}_z}{{\∂\mathrm{\θ}}_x}& \frac{{\∂F}_z}{{\∂\mathrm{\θ}}_y}& 0\\ \frac{{\∂M}_x}{{\∂\mathrm{\δ}}_x}& \frac{{\∂M}_x}{{\∂\mathrm{\δ}}_y}& \frac{{\∂M}_x}{{\∂\mathrm{\δ}}_z}& \frac{{\∂M}_x}{{\∂\mathrm{\θ}}_x}& \frac{{\∂M}_x}{{\∂\mathrm{\θ}}_y}& 0\\ \frac{{\∂M}_y}{{\∂\mathrm{\δ}}_x}& \frac{{\∂M}_y}{{\∂\mathrm{\δ}}_y}& \frac{{\∂M}_y}{{\∂\mathrm{\δ}}_z}& \frac{{\∂M}_y}{{\∂\mathrm{\θ}}_x}& \frac{{\∂M}_y}{{\∂\mathrm{\θ}}_y}& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]$

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

4)建立齿轮模型：采用等效啮合模型对主减速器锥齿轮和轮间差速器齿轮进行建模，对齿轮的载荷传递进行模拟。

在本实施例中，小轮的尺寸特征通过上述小轮轴梁单元的几何参数体现，大轮采用后续的体单元缩维的有限元模型进行模拟。分别在小轮和大轮的理论啮合位置建立等效啮合节点，以模拟齿轮啮合的刚度耦合关系，等效啮合节点的坐标位置根据准双曲面齿轮参数求得。在小轮与等效啮合点之间建立刚性梁单元，用来实现载荷从小轮轴到啮合节点的传递。在小轮和大轮等效啮合点之间建立等效啮合力作用线方向的空间弹簧单元，等效啮合力作用线方向为n＝{n_x n_y n_z}，其中n_x、n_y和n_z分别为等效啮合力作用线方向向量n在全局坐标系下x方向、y方向和z方向的分量，同样根据准双曲面齿轮参数求得。准双曲面齿轮的等效啮合刚度系数k_m为常数，参考ISO标准求得。

在本实施例中，差速器直齿锥齿轮由4个行星轮和2个太阳轮组成，行星轮和太阳轮共有8个齿轮啮合对传递载荷，齿轮的尺寸特征通过上述齿轮轴梁单元的几何参数体现。同样采用等效啮合模型对差速器齿轮进行建模，分别在行星轮和太阳轮的理论啮合位置建立等效啮合节点，以模拟齿轮啮合的刚度耦合关系。在齿轮与等效啮合点之间建立刚性梁单元，用来实现载荷传递。在行星轮和太阳轮等效啮合点之间建立等效啮合力作用线方向的空间弹簧单元，计算得到差速器齿轮刚度矩阵K_d和啮合刚度矩阵K_dm。

5)建立连接部件模型：对驱动桥系统中的连接部件进行模拟。

在本实施例中，连接部件包括垫圈、垫片、花键、螺栓以及部件之间的接触配合关系，如十字轴与差速器壳十字轴孔的配合关系。驱动桥系统中的连接部件在建模时均采用线性弹簧单元模拟，耦合相互作用的节点自由度，全局坐标系下线性弹簧单元的刚度矩阵K_c表示为：

K_c＝diag([k_x k_y k_z k_rx k_ry k_rz])

上式中，k_x为沿x轴方向的平动刚度；k_y为沿y轴方向的平动刚度；k_z为沿z轴方向的平动刚度；k_rx为绕x轴方向的转动刚度；k_ry为绕y轴方向的转动刚度；k_rz为绕z轴方向的转动刚度。

6)建立壳体有限元模型及缩维模型：建立驱动桥系统中所包含的壳体结构有限元模型，并对壳体的有限元模型进行缩维计算。

在本实施例中，轴系均通过滚子轴承支撑在壳体结构上，在系统建模时必须准确建立壳体结构的模型。驱动桥系统的壳体结构主要包括整体桥壳、差速器壳和两侧轮毂。采用四节点四面体单元对壳体的几何模型进行网格划分，单元尺寸为4mm，各部件之间的焊接关系通过共用节点的方式模拟，螺栓连接通过梁单元和多点约束刚性单元模拟。

在本实施例中，整体桥壳包括桥壳、主减速器壳、加强板、轴承座、板簧座、推力杆座和半轴套管。在差速器壳模型中考虑了大轮体，大轮体的有限元模型中，采用大轮节锥尺寸对大轮进行模拟。

在本实施例中，采用模态综合法对壳体的有限元模型进行缩维变换，根据壳体的连接特点，将有限元模型中的节点自由度分为边界自由度和内部自由度，将与轴系连接的节点自由度定义为边界自由度，其他定义为内部自由度。将壳体有限元模型的刚度矩阵和质量矩阵通过模态综合法转换到边界自由度和模态自由度上，在保留原结构模态特性的同时实现对模型的缩维。采用固定截面模态综合方法求得壳体缩维刚度矩阵K_h和壳体缩维质量矩阵M_h。

7)建立完整的驱动桥系统动力学模型：将轴系模型与壳体缩维模型的相关刚度矩阵和质量矩阵按照节点自由度耦合关系组集，获得系统刚度矩阵和系统质量矩阵，建立包含主减速器总成、差速器总成、轮毂总成和桥壳在内的完整驱动桥系统动力学模型。

在本实施例中，将滚子轴承内圈节点自由度与外圈节点自由度用滚子轴承刚度矩阵K_b耦合，将小轮啮合点自由度与大轮啮合点自由度用齿轮啮合刚度K_m耦合，将差速器行星轮和太阳轮之间用差速器齿轮刚度矩阵K_d和啮合刚度矩阵K_dm耦合，将花键等连接关系对应的节点自由度用连接部件刚度矩阵K_c耦合，即获得完整的系统动力学模型，系统模型的动力学方程为：

$M\stackrel{..}{\mathrm{\δ}}+C\stackrel{.}{\mathrm{\δ}}+\mathrm{K\δ}=f$

上式中，δ为节点自由度时域位移向量；f为动态载荷向量；M为系统质量矩阵，由梁单元质量矩阵M_s、壳体缩维质量矩阵M_h组集而成(即将系统质量矩阵中对应节点自由度的元素用轴系梁单元质量矩阵和壳体缩维质量矩阵中对应节点自由度的元素表示)；C为系统阻尼矩阵；K为系统刚度矩阵，由梁单元刚度矩阵K_s、滚子轴承刚度矩阵K_b、准双曲面齿轮啮合刚度矩阵K_m、差速器齿轮刚度矩阵K_d和啮合刚度矩阵K_dm、壳体缩维刚度矩阵K_h和连接部件刚度矩阵K_c组集而成，即将系统刚度矩阵中对应节点自由度的元素用上述部件刚度矩阵中对应节点自由度的元素表示。

8)计算不同输入转矩工况下的轴承刚度：根据驱动桥输入转矩大小的不同定义轻载、中载和重载工况，采用牛顿-拉普森方法迭代求解不同输入转矩下的系统静力学方程，得到静力平衡时的轴承刚度矩阵。

载荷工况的定义方式如下：轻载工况对应的输入转矩T_min所在范围为0<T_min≤0.3T_M，中载工况对应的输入转矩T_mid所在范围为0.3T_M<T_mid<0.7T_M，重载工况对应的输入转矩T_max所在范围为0.7T_M≤T_max<T_M，其中T_M为发动机最大输入转矩。

在本实施例中，由于滚子轴承的刚度矩阵K_b具有各向耦合性和非线性特性，在动力学计算之前需要根据驱动桥系统的输入转矩求得静平衡状态下的轴承刚度，以该轴承刚度矩阵作为线性刚度，进行系统动力学特性计算。采用牛顿-拉夫逊方法对系统静力学方程进行非线性迭代计算，系统静力学非线性方程表示为：

K₀(δ₀,f₀(T))δ₀＝f₀(T)

上式中，K_b为求得的轴承刚度矩阵；T为输入转矩；f₀(T)为输入转矩工况对应的系统载荷向量，δ₀为系统位移向量。

在本实施例中，由于驱动桥系统施加不同大小的输入转矩T，静力平衡时的滚子轴承刚度矩阵不同，对应的系统动力学特性也不同。选取轻、中、重载三个典型工况进行计算，对应的输入转矩大小分别为T_min、T_mid和T_max。

9)计算不同输入转矩工况下的驱动桥系统动力学特性：在不同输入转矩工况下，采用模态叠加法计算单位谐波齿轮传动误差激励下系统的动力学特性。

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

F_mesh＝D_meshδ_mesh

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

其中，D_mesh的计算公式为：

D_mesh＝[C_p+C_g]^-1

上式中，C_p和C_g分别为小轮和大轮的动态柔度。求解小轮和大轮的动态柔度时，需要将齿轮啮合关系解耦，即齿轮的等效啮合刚度系数k_m＝0，分别在小轮和大轮啮合节点上施加沿啮合力作用线方向的单位谐波力e^-jωt，计算系统的动态响应表示为：

${\mathrm{\&delta;}}_0=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&phi;}}_{0i}^{T}F\left(t\right){\mathrm{\&phi;}}_{0i}}{{\mathrm{\&omega;}}_{0i}^2[1-{\mathrm{\&lambda;}}_{0i}^2+2j{\mathrm{\&xi;}}_{0i}{\mathrm{\&lambda;}}_{0i}]}$

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

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

C_p＝δ_0p{x_n,y_n,z_n}^T

C_g＝δ_0g{x_n,y_n,z_n}^T

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

在本实施例中，采用振型叠加法计算系统动态响应：

$\mathrm{\&delta;}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&phi;}}_i^{T}f\left(t\right){\mathrm{\&phi;}}_i}{{\mathrm{\&omega;}}_i^2[1-{\mathrm{\&lambda;}}_i^2+2j{\mathrm{\&xi;}}_i{\mathrm{\&lambda;}}_i]}$

上式中，δ为系统在单位谐波传动误差激励下的系统位移响应；f(t)为单位谐波传动误差引起的激振载荷；ω_i和φ_i为考虑准双曲面齿轮啮合刚度耦合时第i阶振动频率和正则振型；n为保留的模态阶数；λ_i＝ω/ω_i为第i阶频率比，ω为激振频率；ξ_i为第i阶模态阻尼比。

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

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

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

下面通过一个具体的实施例，用以说明本发明的效果。

如图2所示，以某后驱动桥系统为例，该系统包括主减速器总成、差速器总成、轮毂总成和桥壳。主减速器齿轮为准双曲面齿轮，小轮轴1由一个圆柱滚子轴承5和一对圆锥滚子轴承3、4支撑，大轮13与差速器壳14通过螺栓连接，差速器壳14由一对圆锥滚子轴承10、16支撑。输入转矩2作用在小轮轴1上，经由大轮13、差速器壳14、十字轴15、行星轮轴12、太阳轮轴11、半轴6，最终传递至轮毂7。主减速器滚子轴承3、4、5、10、16的内圈与轴系连接，外圈安装在主减速器壳轴承座上，轮毂圆锥滚子轴承8、9的内圈安装在半轴套管上，外圈与轮毂7连接。

1)定义驱动桥全局坐标系：采用汽车的车辆坐标系作为系统的全局坐标系，即车辆前方为x轴正向，左侧为y轴正向，竖直向上为z轴正向，坐标原点为差速器十字轴15的中心位置。

2)建立轴系梁单元模型：如图2所示，驱动桥的轴部件包含小轮轴1、半轴6、十字轴15、行星轮轴12和太阳轮轴11。

小轮轴1有限元模型的示意图如图3所示，小轮轴上共有23个节点和22个梁单元，小轮轴材料的弹性模量为200GPa，泊松比为0.252，密度为7880kg/m³，小轮轴的有限元模型节点编号自左向右依次为1-23，将各梁单元的刚度矩阵组集得到小轮轴的整体刚度矩阵。

半轴6有限元模型的示意图如图4所示，左右两侧半轴的模型相同，左右半轴上共有20个节点和18个梁单元，半轴材料的弹性模量为207GPa，泊松比为0.29，密度为7800kg/m³，半轴有限元模型节点编号依次为24-43，将各梁单元的刚度矩阵组集得到半轴的整体刚度矩阵。

十字轴15共有4个沿全局坐标系Y轴对称分布的分支轴，其中一个分支轴的有限元模型示意图如图5所示，整个十字轴上共有28个节点和24个梁单元，4个分支轴在十字轴中心位置的节点自由度耦合。十字轴材料的弹性模量为213GPa，泊松比为0.286，密度为7870kg/m³，十字轴的有限元模型节点编号依次为44-71，将各分支梁单元的刚度矩阵组集得到十字轴的整体刚度矩阵。

行星轮轴12有限元模型示意图如图6所示，用来模拟行星轮的尺寸特征，一共有四个相同的行星轮轴，沿全局坐标系Y轴对称分布，行星轮轴共有20个节点和16个梁单元。行星轮轴材料的弹性模量为207GPa，泊松比为0.29，密度为7800kg/m³，行星轮轴的有限元模型节点编号依次为72-91，将各梁单元的刚度矩阵组集得到行星轮轴的整体刚度矩阵。

太阳轮轴11有限元模型的示意图如图7所示，用来模拟太阳轮的尺寸特征，左右两侧太阳轮轴的模型相同，太阳轮轴共有14个节点和12个梁单元，太阳轮轴材料的弹性模量为207GPa，泊松比为0.29，密度为7800kg/m³，太阳轮轴有限元模型节点编号依次为92-105，将各梁单元的刚度矩阵组集得到太阳轮轴的整体刚度矩阵。

3)建立非线性滚子轴承模型：如图8所示，传动系统包含8个圆锥滚子轴承和1个圆柱滚子轴承，轴承的局部坐标示意图如图4所示。小轮轴上的后圆柱滚子轴承型号为FAG575867，轴承内径为40mm，外径为94mm，宽度为30mm，滚子数为13，滚子直径为16mm，滚子有效长度为19mm。小轮轴上的中圆锥滚子轴承型号为FAG546439，轴承内径为70mm，外径为165mm，宽度为57mm，滚子数为15，滚子直径为22.6mm，滚子有效长度为39.556mm，接触角为25°。小轮轴上的前圆锥滚子轴承型号为FAG31312，轴承内径为60mm，外径为130mm，宽度为33.5mm，滚子数为16，滚子直径为17.18mm，滚子有效长度为19.8mm，接触角为28.81°。差速器壳上的左圆锥滚子轴承型号为FAG32021，轴承内径为105mm，外径为160mm，宽度为35mm，滚子数为28，滚子直径为13.4mm，滚子有效长度为23.48mm，接触角为16.5°。差速器壳上的右圆锥滚子轴承型号为SKF33021，轴承内径为105mm，外径为160mm，宽度为43mm，滚子数为28，滚子直径为13.74mm，滚子有效长度为29.76mm，接触角为10.67°。两侧轮毂各有一对圆锥滚子轴承，对应的型号相同，外侧轮毂轴承的型号为FAG33019，轴承内径为95mm，外径为145mm，宽度为39mm，滚子数为28，滚子直径为12.617mm，滚子有效长度为29.3mm，接触角为10.5°。内侧轮毂轴承的型号为FAG30222，轴承内径为110mm，外径为200mm，宽度为41mm，滚子数为20，滚子直径为22.944，滚子有效长度为26.6mm，接触角为16°。轴承材料的弹性模量均为210GPa，泊松比均为0.3。根据滚子轴承的非线性刚度计算公式，求得每个轴承各自的非线性刚度矩阵。

4)建立齿轮模型：主减速器准双曲面齿轮参数如表1所示，求得的小轮和大轮理论啮合节点在全局坐标系中的坐标为(181.695mm，-5.054mm，-24.1mm)，求得的等效啮合力作用线方向向量为(0.7183，-0.1992，0.6666)。根据准双曲面齿轮的力学模型计算公式求得齿轮单元刚度矩阵和等效啮合刚度矩阵。差速器齿轮参数如表2所示，由4个行星轮和2个太阳轮组成，行星轮和太阳轮共有8个齿轮啮合对，同样采用等效啮合模型对差速器齿轮进行建模，建模方法与准双曲面齿轮相同。

表1 主减速器准双曲面齿轮参数

表2 差速器齿轮参数

5)建立连接部件模型：驱动桥系统中的连接部件包括：行星轮与差速器壳之间的垫片、太阳轮与差速器壳之间的垫片、十字轴与差速器壳十字轴孔的配合关系、行星轮与十字轴之间的垫圈、半轴与太阳轮之间的花键及半轴与轮毂之间的螺栓连接，将上述连接部件和连接关系处理为空间弹簧，在实际作用力方向上施加线性刚度，以模拟部件之间的连接关系。

6)建立壳体有限元模型及缩维模型：采用四节点四面体单元对相关壳体的几何模型进行网格划分，单元尺寸均为4mm。整体桥壳的有限元模型包括桥壳、主减速器壳、半轴套管、上推力杆支座、板簧支座、下推力杆支座、等与内部轴系相关的完整外部壳体模型，建立的整体桥壳有限元模型共包含306万个单元，71万个节点。差速器壳有限元模型包括差速器壳和大轮体，采用共用节点的方式模拟差速器壳和大轮体之间的螺栓连接关系，建立的差速器壳有限元模型共包含112万个单元，24万个节点。两侧的轮毂模型相同，单侧轮毂有限元模型包含117万个单元，26万个节点。由于壳体有限元模型包含大量的节点自由度，无法直接与轴系模型连接，采用模态综合法对壳体的有限元模型进行缩维变换，根据壳体的连接特点，将有限元模型中的节点自由度分为边界自由度和内部自由度，将与轴系连接的节点自由度定义为边界自由度，其他定义为内部自由度。将原结构的刚度矩阵和质量矩阵通过模态综合法转换到边界自由度和模态自由度上，在保留原结构模态特性的同时实现对模型的缩维。在整体桥壳与轴系连接位置建立边界节点。整体桥壳模型通过9个滚子轴承与轴系连接，在每个滚子轴承中心位置建立边界节点，每个边界节点包含6个边界自由度，用多点约束刚性单元将轴承中心节点自由度与轴承座表面节点自由度耦合。由于整体桥壳受到来自推力杆和板簧的约束作用，需要约束推力杆座和板簧座的相应节点自由度，以消除系统的刚体位移。整体桥壳缩维模型共有9个边界节点，计算时保留100阶主模态。在差速器壳与其他部件连接位置建立边界节点。差速器壳模型与十字轴、差速器齿轮有接触关系，在对应的连接位置建立边界节点，边界节点与轴系部件的接触关系用刚度较大的线性弹簧模拟，差速器壳还与差速器滚子轴承的内圈连接，在滚子轴承中心位置建立边界节点，用多点约束刚性单元将轴承中心节点自由度与差速器壳表面轴承内圈节点耦合，该内圈边界节点与整体桥壳中的外圈边界节点通过轴承刚度连接。由于考虑了大轮的轮坯有限元模型，在大轮节锥的理论啮合位置建立边界节点，作为大轮的啮合节点。差速器壳缩维模型共有9个边界节点，计算时保留40阶主模态。在轮毂与其他部件连接位置建立边界节点，分别为轮毂模型与半轴的螺栓连接中心、轮毂轴承中心和轮距加载中心，因此轮毂缩维模型共有4个边界节点，计算时保留40阶主模态。选取壳体的模态综合截断频率足够大于系统关心频率，保证获得的缩维刚度矩阵和质量矩阵具有足够精度，以准确体现壳体的动态特性。

7)建立完整的驱动桥系统动力学模型：将滚子轴承内圈节点自由度与外圈节点自由度用滚子轴承刚度矩阵耦合，将小轮啮合点自由度与大轮啮合点自由度用齿轮啮合刚度耦合，将差速器行星轮和太阳轮之间用差速器齿轮刚度矩阵和啮合刚度矩阵耦合，将花键等连接关系对应的节点自由度用连接部件刚度矩阵耦合，即获得完整的系统动力学模型。最终建立的由轴系和壳体组成的驱动桥系统动力学模型共有228个节点和220阶壳体主模态自由度，共1588个自由度。

8)计算不同输入转矩工况下的轴承刚度：选取轻、中、重载三个典型工况进行计算，对应的输入转矩大小分别为T_min＝200Nm、T_mid＝3000Nm和T_max＝7000Nm。采用牛顿-拉普森方法迭代求解系统静力学方程，以相邻两次迭代所得的节点位移向量之差的模小于给定小量作为收敛准则，收敛容差取为10^-6，不满足则继续进行迭代，满足则判断迭代过程收敛，求得不同输入转矩工况静力平衡时的轴承刚度矩阵。

9)不同输入转矩工况下的驱动桥动力学特性计算：首先计算齿轮在非啮合状态下在单位谐波力激励下的系统动力学响应，获得齿轮的动态柔度，再计算齿轮啮合状态下在单位谐波传动误差激励下的系统动力学响应，最终求得不同输入转矩工况下齿轮动态啮合力的频响特性曲线如图9所示。由系统节点的位移响应可以进一步计算得到各个节点的速度和加速度响应，以及轴承的载荷响应，从而求得不同转矩工况下的驱动桥系统的动力学特性。

综上所述，本发明适合于不同输入转矩影响下的驱动桥动力学特性计算和分析，由图9不同输入转矩工况下的齿轮动态啮合力频响特性曲线可知，不同输入转矩下驱动桥具有不同的动力学特性，本方法由于考虑了不同输入转矩的影响，能够为驱动桥使用工况的选取提供指导，同时动力学计算和分析结果可以为驱动桥的减振降噪设计提供依据，本发明能够准确高效地实现对完整驱动桥系统的动力学建模和在任意输入转矩工况下的动力学特性计算。

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

Claims (10)

1.一种考虑输入转矩变化的驱动桥动力学特性计算方法，包括以下步骤：

1)定义驱动桥全局坐标系：对驱动桥系统的全局坐标系进行定义，作为系统建模的基础；

2)建立轴系梁单元模型：采用考虑剪切应变的欧拉伯努利空间梁单元对轴系部件进行模拟；

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

4)建立齿轮模型：采用等效啮合模型对主减速器锥齿轮和轮间差速器齿轮进行建模，对齿轮的载荷传递进行模拟；

5)建立连接部件模型：对驱动桥系统中的连接部件进行模拟；

6)建立壳体有限元模型及缩维模型：建立驱动桥系统中所包含的壳体结构有限元模型，并对壳体的有限元模型进行缩维计算；

7)建立完整的驱动桥系统动力学模型：将轴系模型与壳体缩维模型的相关刚度矩阵和质量矩阵按照节点自由度耦合关系组集，获得系统刚度矩阵和系统质量矩阵，建立包含主减速器总成、差速器总成、轮毂总成和桥壳在内的完整驱动桥系统动力学模型；

8)计算不同输入转矩工况下的轴承刚度：根据驱动桥输入转矩大小的不同定义轻载、中载和重载工况，采用牛顿-拉普森方法迭代求解不同输入转矩下的系统静力学方程，得到静力平衡时的轴承刚度矩阵；

9)计算不同输入转矩工况下的驱动桥系统动力学特性：在不同输入转矩工况下，采用模态叠加法计算单位谐波齿轮传动误差激励下系统的动力学特性。

2.如权利要求1所述的一种考虑输入转矩变化的驱动桥动力学特性计算方法，其特征在于：在所述步骤1)中，采用汽车的车辆坐标系作为驱动桥系统的全局坐标系，车辆前方为x轴正向，左侧为y轴正向，竖直向上为z轴正向，坐标原点为差速器十字轴中心位置。

3.如权利要求1所述的一种考虑输入转矩变化的驱动桥动力学特性计算方法，其特征在于：在所述步骤3)中，滚子轴承的载荷计算公式表示为：

$\left\{\begin{array}{c}{F}_x=-\frac{K_n}{n_s}\underset{j=1}{\overset{Z}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n_s}{\mathrm{\&Sigma;}}}{\mathrm{\&delta;}}_{j,k}^{10/9}\cos \mathrm{\&alpha;}\sin {\mathrm{\&psi;}}_j\\ F_y=\frac{K_n}{n_s}\underset{j=1}{\overset{Z}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n_s}{\mathrm{\&Sigma;}}}{\mathrm{\&delta;}}_{j,k}^{10/9}\cos \mathrm{\&alpha;}\cos {\mathrm{\&psi;}}_j\\ F_z=\frac{K_n}{n_s}\underset{j=1}{\overset{Z}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n_s}{\mathrm{\&Sigma;}}}{\mathrm{\&delta;}}_{j,k}^{10/9}\sin \mathrm{\&alpha;}\\ M_x=\frac{K_n}{n_s}\underset{j=1}{\overset{Z}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n_s}{\mathrm{\&Sigma;}}}(\frac{D_{\mathrm{pw}}}{2}\sin \mathrm{\&alpha;}-x_k){\mathrm{\&delta;}}_{j,k}^{10/9}\cos {\mathrm{\&psi;}}_j\\ M_y=\frac{K_n}{n_s}\underset{j=1}{\overset{Z}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n_s}{\mathrm{\&Sigma;}}}(\frac{D_{\mathrm{pw}}}{2}\sin \mathrm{\&alpha;}-x_k){\mathrm{\&delta;}}_{j,k}^{10/9}\sin {\mathrm{\&psi;}}_j\\ M_z=0\end{array}\right.$

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

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

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

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

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

$\begin{array}{c}{\mathrm{\&delta;}}_{j,k}=[{\mathrm{\&delta;}}_z+\frac{D_{\mathrm{pw}}}{2}({\mathrm{\&theta;}}_x\cos {\mathrm{\&psi;}}_j+{\mathrm{\&theta;}}_y\sin {\mathrm{\&psi;}}_j)]\sin \mathrm{\&alpha;}+(-{\mathrm{\&delta;}}_x\sin {\mathrm{\&psi;}}_j+{\mathrm{\&delta;}}_y\cos {\mathrm{\&psi;}}_j)\cos \mathrm{\&alpha;}\\ +x_k\left({-\mathrm{\&theta;}}_x\cos {\mathrm{\&psi;}}_j{-\mathrm{\&theta;}}_y\sin {\mathrm{\&psi;}}_j\right)\end{array}$

上式中，δ_x和δ_y分别为轴承内外圈之间沿x轴和y轴的径向变形；δ_z为轴承内外圈之间沿z轴的轴向变形；θ_x和θ_y分别为轴承内外圈之间绕x和y轴的角变形；计算时，对轴承滚子单元的受力状态进行判断，若δ_j,k小于0，则取δ_j,k等于0，表示该轴承滚子单元未发生接触变形，不传递载荷。

4.如权利要求1所述的一种考虑输入转矩变化的驱动桥动力学特性计算方法，其特征在于：在所述步骤5)中，所述连接部件在建模时均采用线性弹簧单元模拟，耦合相互作用的节点自由度，全局坐标系下线性弹簧单元的刚度矩阵K_c表示为：

K_c＝diag([k_x k_y k_z k_rx k_ry k_rz])

上式中，k_x为沿x轴方向的平动刚度；k_y为沿y轴方向的平动刚度；k_z为沿z轴方向的平动刚度；k_rx为绕x轴方向的转动刚度；k_ry为绕y轴方向的转动刚度；k_rz为绕z轴方向的转动刚度。

5.如权利要求1所述的一种考虑输入转矩变化的驱动桥动力学特性计算方法，其特征在于：在所述步骤6)中，采用模态综合法对壳体的有限元模型进行缩维变换，根据壳体的连接特点，将有限元模型中的节点自由度分为边界自由度和内部自由度，将与轴系连接的节点自由度定义为边界自由度，其他节点自由度定义为内部自由度；将壳体有限元模型的刚度矩阵和质量矩阵通过模态综合法转换到边界自由度和模态自由度上，在保留原结构模态特性的同时实现对模型的缩维，采用固定截面模态综合方法求得壳体缩维刚度矩阵K_h和壳体缩维质量矩阵M_h。

6.如权利要求1所述的一种考虑输入转矩变化的驱动桥动力学特性计算方法，其特征在于：在所述步骤7)中，所建立的驱动桥系统动力学模型为：

$M\stackrel{..}{\mathrm{\&delta;}}+C\stackrel{.}{\mathrm{\&delta;}}+\mathrm{K\&delta;}=f$

上式中，δ为节点自由度时域位移向量；f为动态载荷向量；M为系统质量矩阵，由梁单元质量矩阵M_s、壳体缩维质量矩阵M_h组集而成；C为系统阻尼矩阵；K为系统刚度矩阵，由梁单元刚度矩阵K_s、滚子轴承刚度矩阵K_b、准双曲面齿轮啮合刚度矩阵K_m、差速器齿轮刚度矩阵K_d和啮合刚度矩阵K_dm、壳体缩维刚度矩阵K_h和连接部件刚度矩阵K_c组集而成。

7.如权利要求1所述的一种考虑输入转矩变化的驱动桥动力学特性计算方法，其特征在于：在所述步骤8)中，系统静力学非线性方程表示为：

K₀(δ₀,f₀(T))δ₀＝f₀(T)

上式中，K_b为求得的轴承刚度矩阵；T为输入转矩；f₀(T)为输入转矩工况对应的系统载荷向量，δ₀为系统位移向量。

8.如权利要求1所述的一种考虑输入转矩变化的驱动桥动力学特性计算方法，其特征在于：在所述步骤9)中，驱动桥的动力学特性通过主减速器齿轮的动态啮合力体现，主减速器齿轮的动态啮合力为F_mesh：

F_mesh＝D_meshδ_mesh

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

9.如权利要求8所述的一种考虑输入转矩变化的驱动桥动力学特性计算方法，其特征在于：D_mesh的计算公式为：

D_mesh＝[C_p+C_g]^-1

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

10.如权利要求8所述的一种考虑输入转矩变化的驱动桥动力学特性计算方法，其特征在于：δ_mesh的求解方法如下：

采用振型叠加法计算系统动态响应：

$\mathrm{\&delta;}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&phi;}}_i^{T}f\left(t\right){\mathrm{\&phi;}}_i}{{\mathrm{\&omega;}}_i^2[1-{\mathrm{\&lambda;}}_i^2+2j{\mathrm{\&xi;}}_i{\mathrm{\&lambda;}}_i]}$

上式中，δ为系统在单位谐波传动误差激励下的系统位移响应；f(t)为单位谐波传动误差引起的激振载荷；ω_i和φ_i为考虑准双曲面齿轮啮合刚度耦合时第i阶振动频率和正则振型；n为保留的模态阶数；λ_i＝ω/ω_i为第i阶频率比，ω为激振频率；ξ_i为第i阶模态阻尼比；

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

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

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