传动装置的高效、可靠、低振动、低噪声问题日益受到各行业的关注，齿轮传动装置作为主要的机械传动形式，其振动噪声的控制问题也不容忽视。在齿轮箱系统的动态特性传统计算中，常将齿轮传动系统与箱体分开考虑。国内外许多学者针对齿轮啮合副动力学模型和齿轮转子系统动力学模型做了大量研究^[1-4]，侧重分析齿轮副以及转子的响应；而箱体则单独计算，将轴承动载荷作为输入，利用有限元分析得到箱体的振动特性^[5]，多数情况下箱体的刚性较大，动态特性时不用考虑箱体，研究齿轮系统，忽略其柔性的影响。但是在现代船舶、航空等工程领域箱体的非刚性问题较为突出，大型风电齿轮箱、航空动力传动机匣及船用齿轮系统箱体结构较大，壁厚相对较薄，箱体薄壁结构带来的柔性支撑，会使齿轮系统与箱体之间存在动力耦合效应。因此，为获得更为准确的系统响应，更好地预测此类系统的振动噪声，有必要建立齿轮-轴-轴承-箱体全耦合模型来分析系统动力学特性。

目前对于齿轮系统耦合的研究主要集中在齿轮、传动轴、轴承耦合，而建立考虑箱体柔性的全耦合模型的较少。Choy^[6-7]采用模态综合法研究了齿轮箱振动对齿轮传动系统动态特性的影响，但是对箱体做了大量简化。梁明轩^[8]基于显式有限元法建立了箱体-轴承-齿轮副三维动态接触非线性模型，但是很难准确全面计入齿轮激励。林腾蛟^[9-10]在全耦合有限元模型中将齿轮简化为分度圆柱，在外部计算激励并施加到有限元模型中，对传动系统采用刚体模型而对箱体采用柔体模型的刚柔混合模型做简化计算。

本文对传动子系统建立广义有限元模型，分别用直接法和静态子结构法提取箱体子系统的动力学参数，通过轴承支承节点将二者相连。在考虑时变啮合刚度激励的情况下，对全耦合模型动态特性进行研究，为齿轮系统和箱体系统动态性能优化提供了理论依据。

1 齿轮箱系统动力学模型的构建

齿轮箱全耦合系统即齿轮-轴-轴承-箱体系统，分为齿轮传动系统和箱体系统两部分，如图 1所示，包括输入轴、输出轴、一对直齿轮副、4个滚动轴承以及与传动系统相连的箱体。全耦合系统的坐标系与齿轮传动系统的坐标系一致。

将系统各部分进行节点离散，形成轴节点l_ij(i=1，2；j=1, 2, …, n, n为第i根轴节点总数)、齿轮节点s₁和s₂、轴承节点p_i(i=1, 2, 3, 4)以及箱体节点b_i(i=1, 2, …, n, n为箱体节点总数)，并根据节点间的物理关系建立系统有限单元模型，如图 2所示。

1.1 轴段单元

采用Timoshenko梁单元来建立轴段单元方程，并得到其刚度矩阵和考虑到相邻节点间质量耦合作用的一致质量矩阵。如图 3所示为一段弹性轴段示意图，用q_l={x₁，y₁，z₁，θ_x1，θ_y1，θ_z1，x₂，y₂，z₂，θ_x2，θ_y2，θ_z2}来表示轴段节点的坐标，其中，x, y, z为节点沿坐标方向位移，θ_x、θ_y、θ_z为绕3个坐标的截面转角。

根据弹性力学理论，可得Timoshenko梁的单元刚度矩阵K_t和一致质量矩阵M_t为

$ {\mathit{\boldsymbol{K}}_t} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{K}}_{t11}}}&{{\mathit{\boldsymbol{K}}_{t12}}}\\ {{\mathit{\boldsymbol{K}}_{t21}}}&{{\mathit{\boldsymbol{K}}_{t22}}} \end{array}} \right] $

(1)

$ {\mathit{\boldsymbol{M}}_t} = \frac{{\rho Al}}{6}\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{M}}_{t11}}}&{{\mathit{\boldsymbol{M}}_{t12}}}\\ {{\mathit{\boldsymbol{M}}_{t21}}}&{{\mathit{\boldsymbol{M}}_{t22}}} \end{array}} \right] $

(2)

其中，

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{K}}_{t12}} = \mathit{\boldsymbol{K}}_{t21}^{\rm{T}}\\ {\mathit{\boldsymbol{M}}_{t21}} = {\mathit{\boldsymbol{M}}_{t12}}\\ {\mathit{\boldsymbol{M}}_{t22}} = {\mathit{\boldsymbol{M}}_{t11}} \end{array} \right. $

(3)

式中：K_t11、K_t12、K_t22、M_t11、M_t12对应的矩阵计算参考文献[11]。

阻尼矩阵的计算采用工程中常用的啮合单元阻尼矩阵计算公式为

$ \mathit{\boldsymbol{C}} = {\alpha _0}\mathit{\boldsymbol{M}} + {\alpha _1}\mathit{\boldsymbol{K}} $

(4)

式中：α₀为质量比例系数，α₁为刚度比例系数。

本齿轮传动系统中，将主动轴、从动轴分别划分11、12个轴段，共得到25个轴节点，形成维度为150×150的质量矩阵和刚度矩阵。

考虑弹性轴与齿轮、轴承之间的相互作用力，通过对Timoshenko梁单元进行有限元组装，可得如式(5)所示的整个轴段的运动微分方程。M_l、C_l、K_l为轴段单元的质量矩阵、阻尼矩阵、刚度矩阵，x_l、ẋ_l、ẍ_l为其位移、速度、加速度列向量。

$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}_l^s}&{{\mathit{\boldsymbol{M}}_{ls}}}&{}\\ {{\mathit{\boldsymbol{M}}_{sl}}}&{\mathit{\boldsymbol{M}}_l^l}&{{\mathit{\boldsymbol{M}}_{lp}}}\\ {}&{{\mathit{\boldsymbol{M}}_{pl}}}&{\mathit{\boldsymbol{M}}_l^p} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot x}}_l^s}\\ {\mathit{\boldsymbol{\ddot x}}_l^l}\\ {\mathit{\boldsymbol{\ddot x}}_l^p} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{C}}_l^s}&{{\mathit{\boldsymbol{C}}_{ls}}}&{}\\ {{\mathit{\boldsymbol{C}}_{sl}}}&{\mathit{\boldsymbol{C}}_l^l}&{{\mathit{\boldsymbol{C}}_{lp}}}\\ {}&{{\mathit{\boldsymbol{C}}_{pl}}}&{\mathit{\boldsymbol{C}}_l^p} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot x}}_l^s}\\ {\mathit{\boldsymbol{\dot x}}_l^l}\\ {\mathit{\boldsymbol{\dot x}}_l^p} \end{array}} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}}_l^s}&{{\mathit{\boldsymbol{K}}_{ls}}}&{}\\ {{\mathit{\boldsymbol{K}}_{sl}}}&{\mathit{\boldsymbol{K}}_l^l}&{{\mathit{\boldsymbol{K}}_{lp}}}\\ {}&{{\mathit{\boldsymbol{K}}_{pl}}}&{\mathit{\boldsymbol{K}}_l^p} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{x}}_l^s}\\ {\mathit{\boldsymbol{x}}_l^l}\\ {\mathit{\boldsymbol{x}}_l^p} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{F}}_l^s}\\ 0\\ {\mathit{\boldsymbol{F}}_l^p} \end{array}} \right]} \end{array} $

(5)

式中：M_l^s、C_l^s、K_l^s、x_l^s、ẋ_l^s、ẍ_l^s为与齿轮相连的轴节点矩阵参数，M_l^p、C_l^p、K_l^p、x_l^p、ẋ_l^p、ẍ_l^p表示与轴承相连的轴节点矩阵参数，M_l^l、C_l^l、K_l^l、x_l^l、ẋ_l^l、ẍ_l^l表示剩余轴节点矩阵参数，下角标为ls、sl、lb、bl的为3类节点之间的耦合参数矩阵，可通过Timoshenko梁单元之间的耦合项得到，F_l^s、F_l^b分别表示齿轮和轴承作用在弹性轴上的力。

1.2 啮合单元

图 4为直齿啮合单元动力学模型图，根据啮合齿轮位置关系，将齿轮振动位移向啮合线方向投影，齿轮副沿啮合线方向投影相对总变形为

$ \begin{array}{*{20}{c}} {{\delta _s} = {x_{s1}}\sin \phi + {y_{s1}}\cos \phi + {\theta _{zs1}}{r_{s1}} - }\\ {{x_{s2}}\sin \phi - {y_{s2}}\cos \phi + {\theta _{zs2}}{r_{s2}} - {e_m}} \end{array} $

(6)

啮合单元的运动微分方程为

$ \left\{ \begin{array}{l} {m_{s1}}{{\ddot x}_{s1}} + \left( {{c_m}{{\dot \delta }_s} + {k_m}{\delta _s}} \right)\sin \phi = 0\\ {m_{s1}}{{\ddot y}_{s1}} + \left( {{c_m}{{\dot \delta }_s} + {k_m}{\delta _s}} \right)\cos \phi = 0\\ {I_{zs1}}{{\ddot \theta }_{zs1}} + \left( {{c_m}{{\dot \delta }_s} + {k_m}{\delta _s}} \right){r_{s1}} = 0\\ {m_{s2}}{{\ddot x}_{s2}} - \left( {{c_m}{{\dot \delta }_s} + {k_m}{\delta _s}} \right)\sin \phi = 0\\ {m_{s2}}{{\ddot y}_{s2}} - \left( {{c_m}{{\dot \delta }_s} + {k_m}{\delta _s}} \right)\cos \phi = 0\\ {I_{zs2}}{{\ddot \theta }_{zs2}} + \left( {{c_m}{{\dot \delta }_s} + {k_m}{\delta _s}} \right){r_{s2}} = {T_g} \end{array} \right. $

(7)

式中：x_s1、y_s1、θ_zs1、x_s2、y_s2、θ_zs2为啮合单元节点的位移自由度，r_s1, r_s2为主、从动轮的基圆半径，ϕ为安装相位角，m_i, I_zi(i=s₁, s₂)分别表示主、从动齿轮的质量和绕z轴的转动惯量，T_g表示齿轮副负载扭矩，e_m为齿轮副法线方向的综合啮合误差，k_m为齿轮副法向综合啮合刚度，c_m为齿轮副法向综合啮合阻尼，c_m的计算可参考文献[13]。

考虑齿轮与弹性轴段之间的相互作用力，整理后可以得到啮合单元的运动微分方程组为

$ \mathit{\boldsymbol{M}}_s^l\mathit{\boldsymbol{\ddot x}}_s^l + \mathit{\boldsymbol{C}}_s^l\left( {\mathit{\boldsymbol{\dot x}}_s^l - \mathit{\boldsymbol{\dot e}}} \right) + \mathit{\boldsymbol{K}}_s^l\left( {\mathit{\boldsymbol{x}}_s^l - \mathit{\boldsymbol{e}}} \right) = \mathit{\boldsymbol{F}}_s^l $

(8)

式中：M_s^l、K_s^l、C_s^l为啮合单元的质量矩阵、刚度矩阵、阻尼矩阵，x_s^l、ẋ_s^l、ẍ_s^l分别表示啮合副节点的位移向量、速度向量、加速度向量；e、ė为啮合单元综合啮合误差在各方向分解后的等效位移与等效速度的列向量，F_s^l表示弹性轴作用在齿轮上的力。

1.3 箱体单元

为建立箱体的动力学模型，需提取箱体动力学参数，包括质量矩阵、刚度矩阵和阻尼矩阵。直接法提取箱体参数即为提取箱体有限元模型上各个点的参数，而静态子结构法只提取箱体等效主节点的参数。本文分别采取直接法及静态子结构法提取箱体参数，箱体有限元模型如图 5所示。

采用静态子结构法提取箱体矩阵信息时，箱体节点分为内部节点和外部节点，将轴承孔面上的节点耦合到孔中心即得到外部节点，其余有限元节点为内部节点，箱体质量矩阵、刚度矩阵可通过子结构法将内部节点凝聚到外部节点来提取得到。

文献[12]对于无阻尼的子结构系统，可以写出其动力学方程为

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{M}}_{ww}}{{\mathit{\boldsymbol{\ddot X}}}_w} + {\mathit{\boldsymbol{M}}_{wn}}{{\mathit{\boldsymbol{\ddot X}}}_w} + {\mathit{\boldsymbol{K}}_{ww}}{\mathit{\boldsymbol{X}}_w} + {\mathit{\boldsymbol{K}}_{wn}}{\mathit{\boldsymbol{X}}_w} = {\mathit{\boldsymbol{F}}_w}\\ {\mathit{\boldsymbol{M}}_{nn}}{{\mathit{\boldsymbol{\ddot X}}}_n} + {\mathit{\boldsymbol{M}}_{nw}}{{\mathit{\boldsymbol{\ddot X}}}_w} + {\mathit{\boldsymbol{K}}_{nn}}{\mathit{\boldsymbol{X}}_n} + {\mathit{\boldsymbol{K}}_{nw}}{\mathit{\boldsymbol{X}}_w} = {\mathit{\boldsymbol{F}}_n} \end{array} \right. $

(9)

式中：下标w表示主节点的自由度，n表示从节点的自由度。wn与nw下标表示二者之间的耦合作用。X为节点位移矢量，F为作用力矢量。将箱体看作一个子结构，可以获得箱体转换到边界节点以后重新表示的质量矩阵M_b和刚度矩阵K_b,

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{K}}_b} = {\mathit{\boldsymbol{K}}_{ww}} - {\mathit{\boldsymbol{K}}_{wn}}\mathit{\boldsymbol{K}}_{nn}^{ - 1}{\mathit{\boldsymbol{K}}_{nw}}\\ {\mathit{\boldsymbol{M}}_b} = {\mathit{\boldsymbol{M}}_{ww}} - {\mathit{\boldsymbol{K}}_{wn}}\mathit{\boldsymbol{K}}_{nn}^{ - 1}{\mathit{\boldsymbol{M}}_{nw}} - {\mathit{\boldsymbol{M}}_{wn}}\mathit{\boldsymbol{K}}_{nn}^{ - 1}{\mathit{\boldsymbol{M}}_{nw}} + \\ \;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{M}}_{wn}}\mathit{\boldsymbol{K}}_{nn}^{ - 1}{\mathit{\boldsymbol{M}}_{nn}}\mathit{\boldsymbol{K}}_{nn}^{ - 1}{\mathit{\boldsymbol{M}}_{nw}} \end{array} \right. $

(10)

图 6所示为采用静态子结构法时箱体的有限元局部模型图，约束箱体底面固定，考虑4个轴承孔节点的x、y、z三个平移方向的自由度，则形成12×12的质量矩阵M_b和刚度矩阵K_b。

采用阻尼矩阵的计算公式C_b=α₀M_b+α₁K_b得到箱体阻尼矩阵。考虑轴承对箱体的作用力，将箱体的运动微分方程表示为

$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}_b^p}&{{\mathit{\boldsymbol{M}}_{bp}}}\\ {{\mathit{\boldsymbol{M}}_{pb}}}&{\mathit{\boldsymbol{M}}_b^b} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot x}}_b^p}\\ {\mathit{\boldsymbol{\ddot x}}_b^b} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{C}}_b^p}&{{\mathit{\boldsymbol{C}}_{bp}}}\\ {{\mathit{\boldsymbol{C}}_{pb}}}&{\mathit{\boldsymbol{C}}_b^b} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot x}}_b^p}\\ {\mathit{\boldsymbol{\dot x}}_b^b} \end{array}} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}}_b^p}&{{\mathit{\boldsymbol{K}}_{bp}}}\\ {{\mathit{\boldsymbol{K}}_{pb}}}&{\mathit{\boldsymbol{K}}_b^b} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{x}}_b^p}\\ {\mathit{\boldsymbol{x}}_b^b} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{F}}_b^p}\\ \mathit{\boldsymbol{0}} \end{array}} \right]} \end{array} $

(11)

式中：M_b、C_b、K_b为箱体单元的质量矩阵、阻尼矩阵、刚度矩阵，x_b、ẋ_b、ẍ_b为其位移、速度、加速度列向量。其中M_b^p、C_b^p、K_b^p、x_b^p、ẋ_b^p、ẍ_b^p表示与轴承相连的箱体节点矩阵参数，M_b^b、C_b^b、K_b^b、x_b^b、ẋ_b^b、ẍ_b^b表示箱体其余节点矩阵参数，下标为bp, pb的为两类节点之间的耦合参数矩阵，F_b^p表示轴承作用在箱体上的力。

1.4 轴承单元

将轴承等效为弹簧，则弹簧单元的两端节点分别为与轴承相连的轴节点和箱体节点。本文采用滚动轴承，保留K_xx、K_yy、K_z三项，K_xx、K_yy为径向刚度，K_z为轴向刚度。

采用直接法时，由于箱体轴承孔面上具有多个节点，则形成的每个轴承单元包含多个弹簧，其等效动力学模型如图 7所示。严格来讲，图 7中各弹簧的刚度大小不尽相同，应与轮轴系自重、啮合受力以及轴承的承载特性等相关，需要做进一步研究。本文主要考虑齿轮传动系统和箱体之间的振动耦合传递，作为等效的手段，将各个弹簧刚度的大小简化为一致，其中每个弹簧的刚度为轴承刚度的1/n，n为每个轴承孔面上节点个数，并约束这些节点的位移，形成每个弹簧单元的运动微分方程为

$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {\frac{{\mathit{\boldsymbol{M}}_p^l}}{n}}&{}\\ {}&{\frac{{\mathit{\boldsymbol{M}}_p^b}}{n}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot x}}_p^l}\\ {\mathit{\boldsymbol{\ddot x}}_p^b} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\frac{{\mathit{\boldsymbol{C}}_p^l}}{n}}&{\frac{{ - {\mathit{\boldsymbol{C}}_p}}}{n}}\\ {\frac{{ - {\mathit{\boldsymbol{C}}_p}}}{n}}&{\frac{{\mathit{\boldsymbol{C}}_p^b}}{n}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot x}}_p^l}\\ {\mathit{\boldsymbol{\dot x}}_p^b} \end{array}} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {\frac{{\mathit{\boldsymbol{K}}_p^l}}{n}}&{ - \frac{{{\mathit{\boldsymbol{K}}_p}}}{n}}\\ { - \frac{{{\mathit{\boldsymbol{K}}_p}}}{n}}&{\frac{{\mathit{\boldsymbol{K}}_p^b}}{n}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{x}}_p^l}\\ {\mathit{\boldsymbol{x}}_p^b} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\frac{{\mathit{\boldsymbol{F}}_p^l}}{n}}\\ {\frac{{\mathit{\boldsymbol{F}}_p^b}}{n}} \end{array}} \right]} \end{array} $

(12)

采用静态子结构法时，由于每个轴承孔面上的节点等效为一个外部节点，则形成的每个轴承单元包含一个弹簧，其单元运动微分方程为

$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}_p^l}&{}\\ {}&{\mathit{\boldsymbol{M}}_p^b} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot x'}}_p^l}\\ {\mathit{\boldsymbol{\ddot x'}}_p^b} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{C}}_p^l}&{ - {\mathit{\boldsymbol{C}}_p}}\\ { - {\mathit{\boldsymbol{C}}_p}}&{\mathit{\boldsymbol{C}}_p^b} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot x'}}_p^l}\\ {\mathit{\boldsymbol{\dot x'}}_p^b} \end{array}} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}}_p^l}&{ - {\mathit{\boldsymbol{K}}_p}}\\ { - {\mathit{\boldsymbol{K}}_p}}&{\mathit{\boldsymbol{K}}_p^b} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{x'}}_p^l}\\ {\mathit{\boldsymbol{x'}}_p^b} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{F}}_p^l}\\ {\mathit{\boldsymbol{F}}_p^b} \end{array}} \right]} \end{array} $

(13)

式中：M_p、C_p、K_p为轴承单元的质量矩阵、阻尼矩阵、刚度矩阵，x_p、ẋ_p、ẍ_p、x'_p、ẋ'_p、ẍ'_p分别直接法和子结构法时节点的位移、速度、加速度列向量。其中M_p^l、C_p^l、K_p^l、x_p^l、ẋ_p^l、ẍ_p^l、x'_p^l、ẋ'_p^l、ẍ'_p^l表示与弹性轴相连的轴承节点矩阵参数，M_p^b、C_p^b、K_p^b、x_p^b、ẋ_p^b、ẍ_p^b、x'_p^b、ẋ'_p^b、ẍ'_p^b表示与箱体相连的轴承节点矩阵参数，K_p^l=K_p^b= K_p, C_p^l= C_p^b=C_p, F_p^l, F_p^b分别表示弹性轴和箱体作用在轴承上的力。

1.5 系统总体动力学方程

根据齿轮传动系统各部分的位移连续关系和力平衡关系:

$ \begin{array}{*{20}{c}} {x_l^p = x_p^l\left( {x_l^p = x{'}_p^l} \right),x_l^s = x_s^l,x_p^b = x_b^p\left( {x{'}_p^b = x_b^p} \right),}\\ {F_l^p = - F_p^l,F_l^s = - F_s^l,F_p^b = - F_b^p} \end{array} $

(14)

将单元整合后得到齿轮传动系统整体动力学方程为

$ \mathit{\boldsymbol{M\ddot q}} + \mathit{\boldsymbol{C}}\left( {\mathit{\boldsymbol{\dot q}} - \mathit{\boldsymbol{\dot e}}} \right) + \mathit{\boldsymbol{K}}\left( {\mathit{\boldsymbol{q}} - \mathit{\boldsymbol{e}}} \right) = \mathit{\boldsymbol{F}} $

(15)

其中：

$ \mathit{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}_l^s + \mathit{\boldsymbol{M}}_s^l}&{{\mathit{\boldsymbol{M}}_{ls}}}&{}&{}&{}\\ {{\mathit{\boldsymbol{M}}_{sl}}}&{\mathit{\boldsymbol{M}}_l^l}&{{\mathit{\boldsymbol{M}}_{lp}}}&{}&{}\\ {}&{{\mathit{\boldsymbol{M}}_{pl}}}&{\mathit{\boldsymbol{M}}_l^p + \mathit{\boldsymbol{M}}_p^l}&{}&{}\\ {}&{}&{}&{\mathit{\boldsymbol{M}}_p^b + \mathit{\boldsymbol{M}}_b^p}&{{\mathit{\boldsymbol{M}}_{bp}}}\\ {}&{}&{}&{{\mathit{\boldsymbol{M}}_{pb}}}&{\mathit{\boldsymbol{M}}_b^b} \end{array}} \right] $

(16)

$ \mathit{\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}}_l^s + \mathit{\boldsymbol{K}}_s^l}&{{\mathit{\boldsymbol{K}}_{ls}}}&{}&{}&{}\\ {{\mathit{\boldsymbol{K}}_{sl}}}&{\mathit{\boldsymbol{K}}_l^l}&{{\mathit{\boldsymbol{K}}_{lp}}}&{}&{}\\ {}&{{\mathit{\boldsymbol{K}}_{pl}}}&{\mathit{\boldsymbol{K}}_l^p + \mathit{\boldsymbol{K}}_p^l}&{ - {\mathit{\boldsymbol{K}}_P}}&{}\\ {}&{}&{ - {\mathit{\boldsymbol{K}}_P}}&{\mathit{\boldsymbol{K}}_p^b + \mathit{\boldsymbol{K}}_b^p}&{{\mathit{\boldsymbol{K}}_{bp}}}\\ {}&{}&{}&{{\mathit{\boldsymbol{K}}_{pb}}}&{\mathit{\boldsymbol{K}}_b^b} \end{array}} \right] $

(17)

$ \mathit{\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{C}}_l^s + \mathit{\boldsymbol{C}}_s^l}&{{\mathit{\boldsymbol{C}}_{ls}}}&{}&{}&{}\\ {{\mathit{\boldsymbol{C}}_{sl}}}&{\mathit{\boldsymbol{C}}_l^l}&{{\mathit{\boldsymbol{C}}_{lp}}}&{}&{}\\ {}&{{\mathit{\boldsymbol{C}}_{pl}}}&{\mathit{\boldsymbol{C}}_l^p + \mathit{\boldsymbol{C}}_p^l}&{ - {\mathit{\boldsymbol{C}}_P}}&{}\\ {}&{}&{ - {\mathit{\boldsymbol{C}}_P}}&{\mathit{\boldsymbol{C}}_p^b + \mathit{\boldsymbol{C}}_b^p}&{{\mathit{\boldsymbol{C}}_{bp}}}\\ {}&{}&{}&{{\mathit{\boldsymbol{C}}_{pb}}}&{\mathit{\boldsymbol{C}}_b^b} \end{array}} \right] $

(18)

2 齿轮传动系统输入参数 2.1 齿轮传动系统输入参数

采用直齿轮进行研究，齿轮的齿数为24和79，模数为3 mm，压力角为20°，齿宽为60 mm。输入转速为2 000 r/min，输入功率为141.7 kW。不考虑轴承刚度。系统的内部激励仅考虑齿轮时变啮合刚度激励，时变啮合刚度计算参考文献[13]。即结合有限元法与解析接触力学理论，采用有限元法计算齿轮的弯曲-剪切变形，利用线接触理论计算各啮合点的接触变形，从而计算出齿轮时变静态传递误差x_s(t)和啮合刚度k_m，从而得到时变啮合刚度激励为F_k(t)=x_s(t)×k_m。图 8给出了该齿轮传动系统中齿轮的时变啮合刚度激励数值。

2.2 耦合箱体方法对比

由于该系统为多自由度系统，因此采用傅里叶级数法求解方程，得到系统各自由度的动态位移、速度、加速度。

求得系统各自由度动态位移后，沿啮合线方向的动态传递误差e_DTE为

$ \begin{array}{*{20}{c}} {{e_{{\rm{DTE}}}} = {x_{s1}}\sin \phi + {y_{s1}}\cos \phi + {\theta _{zs1}}{r_{s1}} - }\\ {{x_{s2}}\sin \phi + {y_{s2}}\cos \phi + {\theta _{zs2}}{r_{s2}}} \end{array} $

(19)

则齿轮副法向啮合动载荷可由下式得到：

$ {F_d} = {k_m}\left( {{e_{{\rm{DTE}}}} - {e_m}} \right) + {c_m}\left( {{{\dot e}_{{\rm{DTE}}}} - {{\dot e}_m}} \right) $

(20)

耦合箱体前轴承支反力为

$ \left\{ \begin{array}{l} {F_{px}} = {K_{xx}}x_{px}^l + {C_{xx}}\dot x_{px}^l\\ {F_{py}} = {K_{yy}}x_{py}^l + {C_{yy}}\dot x_{py}^l\\ {F_{pz}} = {K_z}x_{pz}^l + {C_z}\dot x_{pz}^l \end{array} \right. $

(21)

静态子结构法耦合箱体的轴承支反力为

$ \left\{ \begin{array}{l} {{F'}_{px}} = {K_{xx}}x_{px}^l + {C_{xx}}\dot x_{px}^l - {K_{xx}}x_{px}^b - {C_{xx}}\dot x_{px}^b\\ {{F'}_{py}} = {K_{yy}}x_{py}^l + {C_{yy}}\dot x_{py}^l - {K_{yy}}x_{py}^b - {C_{yy}}\dot x_{py}^b\\ {{F'}_{pz}} = {K_z}x_{pz}^l + {C_z}\dot x_{pz}^l - {K_y}x_{py}^b - {C_y}\dot x_{py}^b \end{array} \right. $

(22)

直接法耦合箱体的轴承支反力为

$ \left\{ \begin{array}{l} {{F''}_{px}} = {K_{xx}}x_{px}^l + {C_{xx}}\dot x_{px}^l - \sum\limits_{i = 1}^n {\left( {\frac{{{K_{xx}}}}{n}x_{pxi}^b + \frac{{{C_{xx}}}}{n}\dot x_{pxi}^b} \right)} \\ {{F''}_{py}} = {K_{yy}}x_{py}^l + {C_{yy}}\dot x_{py}^l - \sum\limits_{i = 1}^n {\left( {\frac{{{K_{yy}}}}{n}x_{pyi}^b + \frac{{{C_{yy}}}}{n}\dot x_{pyi}^b} \right)} \\ {{F''}_{pz}} = {K_z}x_{pz}^l + {C_z}\dot x_{pz}^l - \sum\limits_{i = 1}^n {\left( {\frac{{{K_z}}}{n}x_{pzi}^b + \frac{{{C_z}}}{n}\dot x_{pzi}^b} \right)} \end{array} \right. $

(23)

式中：F_px、F_py、F_pz、F'_px、F'_py、F'_pz、F″_px、F″_py、F″_pz为轴承支反力在坐标轴方向的分量，x_px^l、x_py^l、x_pz^l、ẋ_px^l、ẋ_py^l、ẋ_pz^l分别为轴承单元中与弹性轴段相连节点的动态位移、动态速度在坐标轴方向的分量，x_px^b、x_py^b、x_pz^b、ẋ_px^b、ẋ_py^b、ẋ_pz^b分别为轴承单元中与箱体相连节点的动态位移、动态速度在坐标轴方向的分量，n为每个轴承孔面上节点个数。

为了比较本文采用的两种耦合箱体的方法，且由于直接法计算耗时较长，很难计算全转速动态特性，在3.4、3.5节全转速动态响应分析中，得到在转速2 880 r/min时耦合箱体前后齿轮动态啮合力均方差差距最大，在转速3 720 r/min时耦合箱体前后轴承动载荷均方差差距最大，因此针对图 1的模型，分别采用直接法和子结构法计算得到耦合前后2 880 r/min时齿轮动态啮合力与3 720 r/min时3号轴承p₃的轴承支反力，来比较两种耦合方法的差异性。

图 9所示为2 880 r/min时耦合前后齿轮动态啮合力对比图。两种方法耦合箱体后，动态啮合力两条曲线基本重合。子结构法耦合箱体后，相对于未耦合箱体，动态啮合力变化范围为-5.035%~5.038%，峰值减小了1.634%；直接法耦合箱体后，动态啮合力变化范围为-5.425%~4.689%，峰值减小了1.703%。直接法相对于子结构法，动态啮合力的变化范围为-0.381%~0.421%，峰值减小了0.407%。无论是相对于未耦合箱体时的动态数据，还是两种方法的直接对比，都反映出两种方法耦合箱体后动态啮合力差异不大。

图 10所示为3 720 r/min时轴承支反力对比图。子结构法耦合箱体后，相对于未耦合箱体，支反力变化率范围为-4.129%~4.353%，峰值增大了4.327%；直接法耦合箱体后，支反力变化率范围为-3.726%~3.835%，峰值增大了3.508%。直接法相对于子结构法，动态啮合力的变化范围为-0.785%~0.639%，峰值减小了0.785%。对比数据可以得到，相对于耦合前后支反力的变化，两种耦合方法的差异性对轴承支反力的影响较小。

在该计算工况下，两种方法耦合箱体计算结果有一定的差别，但差别不大，都远小于箱体柔性对动态响应的影响。理论上讲直接法耦合箱体得到的动态响应更加真实，且采用直接法可得到箱体上各个点的动态响应。但本文模型中直接法耦合箱体共4 023个节点，计算用时6 290 s，静态子结构法耦合箱体共162个节点，用时99 s，仅为直接法的1.57%，大大缩短了计算时间。因此，静态子结构法更适合计算大型多自由度耦合模型，本文以下计算均采用静态子结构法。

2.3 耦合箱体对系统固有频率的影响

系统固有特性问题可转化为求解如下特征值问题：

$ \omega _i^2\mathit{\boldsymbol{M}}{\varphi _i} = \mathit{\boldsymbol{K}}{\varphi _i} $

(24)

式中：M为系统质量矩阵，K为系统刚度矩阵，ω_i为第i阶固有频率，φ_i为第i阶固有频率对应的振型。

根据式(24)计算得到耦合前后箱体的固有频率。耦合箱体增加了传动系统的柔性，也增加了系统自由度，引入更多阶数的固有频率，对比系统前150阶固有频率，由图 11可知，耦合箱体后系统的固有频率明显降低了，降低幅度为1.5%~67.8%。

2.4 耦合箱体对齿轮啮合力的影响

在系统基本参数不变的情况下，改变输入转速，分析箱体柔性对齿轮动态啮合力的影响。

考察全转速下箱体柔性对齿轮动态啮合力的影响, 图 12为0~16 000 r/min时齿轮动态啮合力随转速变化的均方差。耦合箱体后两个主共振峰未发生偏移，对应转速下的啮合频率分别为2 328.4 Hz和4 657.2 Hz，对应耦合箱体前的第16阶固有频率及1/2倍频，耦合后的第25阶固有频率及1/2倍频。耦合箱体后，主共振峰的峰值分别由9 043.7和17 005.0降低到9 007.2和16 997.4，分别降低了0.4%和0.047%。在全转速下，耦合前后动态啮合力均方差的变化范围为-11.7%~8.6%，其中，在2 880 r/min时耦合前后啮合力均方差变化最大，降低了11.7%。

选择输入转速2 000、2 880 r/min及低阶共振转速5 821 r/min为代表转速，计算该3种转速时的齿轮动态啮合力，得到耦合箱体对动态啮合力的影响如图 13所示。

当转速为2 000 r/min及5 821 r/min时，耦合箱体前后动态啮合力的曲线基本重合，动态啮合力的变化范围分别为-0.886%~1.00%、-1.01%~1.15%，峰值和波动量的变化率都小于0.7%。这两种转速下的动态啮合力变化范围以及峰值、波动量的变化都小于2 880 r/min时，这与2 880 r/min时啮合力均方差变化大也是对应的。

2.5 耦合箱体对轴承支反力的影响

考察全转速下箱体柔性对轴承支反力的影响, 图 14为0~16 000 r/min时3号轴承的轴承支反力均方差随转速变化图。可以看出耦合箱体前主共振峰分别对应系统第5阶和第9阶固有频率, 耦合箱体后对应系统第6阶和第10阶固有频率, 主共振峰的峰值分别降低了6.4%和39.6%。在全转速下，耦合前后动态啮合力均方差的变化范围为-64.5%~97.2%，其中，在3 720 r/min时耦合前后啮合力均方差变化最大，增大了97.2%。可以看出，在全转速范围内，相对于动态啮合力，耦合箱体对轴承支反力的影响更大。

图 15所示为2 000 r/min、3 720 r/min及5 821 r/min时耦合箱体前后两种模型中耦合箱体对轴承支反力的影响。2 000 r/min、5 821 r/min转速下，耦合箱体前后轴承支反力的变化范围分别为-0.750%~0.941%、-3.78%~3.77%，峰值分别增大了0.239%、1.849%，波动量分别增大了17.66%、和121.1%。这两种转速下的轴承支反力和变化范围以及峰值、波动量的变化都小于3 720 r/min时，这也与3 720 r/min时支反力均方差变化大是对应的。

3 结论

1) 分别采用直接法和静态子结构法耦合箱体，两种方法耦合后齿轮动态啮合力和轴承支反力差别不大，且采用静态子结构法计算时间显著减少。因此，在大型多自由度计算模型中，采用静态子结构法更加合适。

2) 耦合箱体增加了传动系统的柔性，使得系统的固有频率明显降低。

3) 箱体柔性对齿轮动态啮合力和轴承支反力的影响在不同转速下表现不同，分析箱体柔性对齿轮传动系统动态特性的影响需在给定工况下具体分析。考察全转速下啮合力均方差和轴承支反力均方差，以耦合前后共振峰的变化率和均方差变化范围为衡量指标，可以得出，箱体柔性对轴承支反力的影响更大，该结论与文献[14]中采用阻抗综合法建立全耦合模型计算出的结果一致。