CN110162909B - 一种渐开线直齿轮传动系统动态特性求解方法 - Google Patents

Abstract

本发明提供一种渐开线直齿轮传动系统动态特性求解方法。本发明方法，包括：步骤S1：基于集中质量法建立渐开线直齿轮传动系统非线性动力学模型；步骤S2：考虑时变啮合刚度、动态传递误差、摩擦、偏心、修形、间隙、重力和非线性轴承力，基于拉格朗日方程推导渐开线直齿轮传动系统非线性动力学方程；步骤S3：基于Runger‑Kutta法求解渐开线直齿轮传动系统动态特性。本发明方法全面考虑时变啮合刚度和动态传递误差等非线性影响因素，将集中质量法、拉格朗日方程和Runger‑Kutta法相结合应用于渐开线直齿轮传动系统动态特性求解，提高了渐开线直齿轮传动系统动态特性的求解准确性和效率，对提高齿轮传动系统的啮合平稳性、承载能力和降低摩擦损耗等具有重要意义。

Description

一种渐开线直齿轮传动系统动态特性求解方法

技术领域

本发明涉及动力学技术领域，具体而言，尤其涉及一种渐开线直齿轮传动系统动态特性求解方法。

背景技术

作为主要传动系统，齿轮系统其传动性能的优劣直接影响着高端装备的性能和可靠性。齿轮系统传动过程中往往受多因素非线性耦合的影响，因而渐开线直齿轮传动系统非线性动力学的研究对改善高端装备机械性能和可靠性具有重要意义。

国内外许多学者对齿轮系统的动力学特性做了很多有意义的研究，并取得了一定的成果。Fernandez等建立了含有齿廓误差的非线性动力学模型，在此基础上分析了齿廓误差对传动性能的影响。张涛等分别以齿廓误差和齿距误差为对象，利用傅里叶变换量化分析了不同加工公差等级下的单项制造误差对齿轮副动态传递误差、角加速度特性的影响规律。Velex等研究了齿轮系统中传递误差对动态啮合力的影响。石照耀等建立了考虑时变啮合刚度和加工误差因素的单双齿交替的直齿轮动力学模型，研究了不同工况下啮合刚度和加工误差对振动响应的影响。Baguet和Jacquenot考虑时变啮合刚度和齿轮侧隙研究了齿轮传动系统的动态响应。李应刚等在同时考虑时变啮合刚度和齿侧间隙因素的影响下，建立了齿轮副系统扭转振动模型，并采用增量谐波平衡法研究了齿轮系统动态特性。胡鹏等建立了考虑静态传递误差、时变啮合刚度和齿侧间隙的惰轮系统非线性动力学模型，研究了时变刚度和载荷力矩对系统扭转振动响应的影响。Brethee等考虑摩擦因素研究了不同摩擦系数下的齿轮系统振动特性。郜浩冬等建立了考虑间隙和摩擦的多级齿轮扭转振动模型，研究了不同摩擦系数对系统动态响应的影响。Yu研究了齿轮偏心对动态传递误差的影响。Liu等考虑动态传递误差、齿侧间隙和偏心影响因素，推导了齿轮系统的运动微分方程，并通过Newmark方法对其进行了求解。张涛等基于ANSYS分析了齿廓误差和齿距误差对直齿圆柱齿轮动力学特性的影响。Wang等基于ANSYS/LS-DYNA研究了不同间隙下偏心对动态传递误差的影响。另外，一些学者通过有限元软件二次开发，分析了摩擦和偏心因素对齿轮系统振动响应的影响。

综上所述，虽然国内外学者在齿轮系统非线性动力学研究方面做了很多有意义的研究，但是大部分研究或者考虑影响因素不全面，或者计算效率不高，因而不能高效准确地反映齿轮系统的动态特性。

发明内容

根据上述技术问题，本发明提供了一种渐开线直齿轮传动系统动态特性求解方法，该方法能够高效准确地反映渐开线直齿轮传动系统的动态特性。

本发明采用的技术手段如下：

一种渐开线直齿轮传动系统动态特性求解方法，包括以下步骤：

步骤S1：基于集中质量法建立渐开线直齿轮传动系统非线性动力学模型；

步骤S2：考虑时变啮合刚度、动态传递误差、摩擦、偏心、修形、间隙、重力和非线性轴承力，基于拉格朗日方程推导渐开线直齿轮传动系统非线性动力学方程；

步骤S3：基于Runger-Kutta法求解渐开线直齿轮传动系统动态特性。

进一步地，所述步骤S1建立的渐开线直齿轮传动系统非线性动力学模型中，在主、从动齿轮的理想中心A_i处建立固定坐标系A_i(x_i,y_i,z_i)(i＝1,2)，在轴承的理想中心B_i处建立固定坐标系B_i(x_bi,y_bi,z_bi)(i＝1～4)。坐标轴x_i，x_bi垂直于齿轮啮合线；坐标轴y_i，y_bi平行于齿轮啮合线；坐标轴z_i，z_bi穿过轴承理想中心。主、从动齿轮的旋转中心为O_i(x_i，y_i)(i＝1，2)，质心为G_i(x_gi,y_gi)(i＝1,2)。轴承的旋转中心为O_bi(x_bi,y_bi)(i＝1～4)。所述渐开线直齿轮传动系统非线性动力学模型的输入端转角为

输入端扭转振动角位移为θ_p；主动轮转角为

主动轮扭转振动角位移为θ₁；从动轮转角为

从动轮扭转振动角位移为θ₂；输出端转角为

输出端扭转振动角位移为θ_q，其关系可表示为：

式中，ω₁为主动轮角速度，ω₂为从动轮角速度，t为时间。

旋转中心O_i(x_i,y_i)与质心G_i(x_gi,y_gi)，(i＝1,2)之间的关系为：

式中，ρ₁为主动轮偏心量，ρ₂为从动轮偏心量。

主、从动轴的弹性变形δ_xi,δ_yi(i＝1,2)可表示为：

式中，齿轮位置系数η_i＝l_bi/l_j(i＝1,2,j＝1；i＝3,4,j＝2)；l₁，l₂为主、从动轴在两轴承间的长度；l_bi(i＝1～4)为主、从动轴上齿轮旋转中心到轴承质心的距离。

进一步地，所述步骤S2的具体步骤如下：

步骤S21：根据拉格朗日方程，分别建立非保守系统广义坐标X，系统的动能T，势能U，耗散函数R，系统除粘性耗散力以外的非保守广义力P和渐开线直齿轮传动系统非线性动力学方程，其表达式如下：

X＝[θ_p x₁ y₁ θ₁ x₂ y₂ θ₂ x_b1 y_b1 x_b2 y_b2 x_b3 y_b3 x_b4 y_b4 θ_q]^T (4)

上式中，J_p、J_q、J₁、J₂分别表示输入端转动惯量、输出端转动惯量、主动轮转动惯量、从动轮转动惯量，m₁、m₂、m_bi(i＝1～4)分别表示主动轮质量、从动轮质量和四个轴承的质量，r_b1、r_b2为主、从动齿轮基圆半径，L_p1、L_p2、L_g1、L_g2为齿对1、2摩擦力对主、从动齿轮力臂，ρ₁、ρ₂为主、从动齿轮偏心量，η_i(i＝1～4)为齿轮位置系数，δ_x1、δ_x2、δ_y1、δ_y2为主、从动轴分别沿x,y方向的弹性变形，ε为重合度，g为重力加速度。k_t1、k_t2分别表示中心轴1、2扭转刚度，c_t1、c_t2分别表示中心轴1、2扭转阻尼，k_s1、k_s2分别表示中心轴1、2弯曲刚度，c_s1、c_s2分别表示中心轴1、2弯曲阻尼，k_bxi,k_byi(i＝1～4)表示四个轴承处沿x,y方向的刚度，c_bxi,c_byi(i＝1～4)表示四个轴承处沿x,y方向的阻尼。M_p、M_q分别表示输入扭矩、输出扭矩，F_f1、F_f2为啮合齿对1、2的摩擦力，F_m1、F_m2为啮合齿对1、2的动态啮合力，F_bxi,F_byi(i＝1～4)表示四个轴承处沿x,y方向的非线性轴承力。

步骤S22：计算时变啮合刚度，其计算公式为：

式中，

分别为啮合齿对1、2的时变啮合刚度；平均啮合刚度k_d＝k₀+A₀/2，k₀＝k_max(ε-1)+k_min(2-ε)，A₀＝2Δ(2ε-3)；各谐波分量刚度幅值k_n＝4Δsin(nπ(ε-1))/nπ，Δ＝(k_max-k_min)/2，ε为重合度，一般取到9阶。其中，k_max、k_min分别为时变啮合刚度最大值、最小值；ω_e为啮合频率，ω_e＝2πn₁z₁/60，n₁为主动轮转速，z₁为主动轮齿数；啮合点经过单齿啮合区和双齿啮合区的时间

主动轮转速

mod()为取余函数，t为时间。

步骤S23：计算动态啮合力，根据黏弹性理论，其计算公式为：

式中，F_m1、F_m2为啮合齿对1、2的动态啮合力，c_m为啮合阻尼，

分别为啮合齿对1、2的时变啮合刚度，

mod()为取余函数，t为时间。则由几何关系可知，动态传递误差

其中，静态传递误差e(t)＝e₀+e_rsin(ω_et+φ_e)，e₀、e_r为静态传递误差的均值和波动幅值，φ_e为相位角；啮合点处修形量

分别为从动轮和主动轮啮合点处修形量，其中：

式中，Δ_max为最大修形量，L为修形长度，啮合周期t₁＝ε·t₀，t₀为啮合点经过单齿啮合区和双齿啮合区的时间，ω₁为主动轮角速度，r_b1为主动轮基圆半径，

mod()为取余函数，t为时间。

f(δ)为间隙函数：

式中，b为间隙，δ为动态传递误差。

步骤S24：计算摩擦力，其计算公式为：

式中，F_f1、F_f2为啮合齿对1、2的摩擦力，μ为摩擦系数，F_m1、F_m2为啮合齿对1、2的动态啮合力，方向系数

和

的计算公式为：

式中，齿轮从进入啮合到节圆处啮合时间

t₁为啮合周期，t₀为啮合点经过单齿啮合区和双齿啮合区的时间，ω₁为主动轮角速度，r_b1和r_b2为主、从动齿轮基圆半径，r₁为主动轮分度圆半径，r_a2为从动轮齿顶圆半径，α为节圆压力角，

mod()为取余函数，t为时间。

步骤S25：计算摩擦力臂，所述摩擦力臂包括齿对i啮合点距主动轮摩擦力臂L_pi和啮合点距从动轮摩擦力臂L_gi，其计算公式分别为：

式中，r_b1和r_b2为主、从动齿轮基圆半径，r_a2为从动轮齿顶圆半径，α为节圆压力角，ω₁为主动轮角速度，t₀为啮合点经过单齿啮合区和双齿啮合区的时间，

mod()为取余函数，t为时间。

步骤S26：计算非线性轴承力沿x方向分力F_bx和y方向分力F_by，其计算公式为：

式中，K_c为赫兹接触刚度系数；N_b为滚珠个数；

为δ_j的p次幂，当轴承为球轴承时，指数p取3/2，当轴承为滚子轴承时，指数p取10/9；H(δ_j)为Heaviside函数，δ_j≤0时，H(δ_j)＝0；δ_j＞0时，H(δ_j)＝1。

第j个滚动体与滚道的法向接触变形量δ_j为：

δ_j＝x_bjcosθ_j+y_bjsinθ_j-γ₀ (21)

式中，x_bj、y_bj为内圈中心振动位移，γ₀为轴承游隙。第j个滚动体在时间t内的转动角度为θ_j：

θ_j＝ω_b·t+2π(j-1)/N_b (22)

式中，保持架的角速度ω_b＝ω×r′/(R′+r′)，ω为轴承角速度，R′、r′为轴承内外圈半径。

与现有技术相比，本发明具有以下优点：

1、本发明方法综合考虑时变啮合刚度、动态传递误差、摩擦、偏心、修形、间隙、重力和非线性轴承力非线性影响因素，基于拉格朗日方程推导了渐开线直齿轮传动系统非线性动力学方程，为求解渐开线直齿轮传动系统动态特性提供了准确的数学模型。

2、本发明方法采用Runger-Kutta法在i5处理器和8GB内存计算机中求解渐开线直齿轮传动系统动态特性耗时为3min，效率高。

基于上述理由本发明可在动力学等领域广泛推广。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做以简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明方法的流程图。

图2为用本发明方法建立的渐开线直齿轮传动系统非线性动力学模型。

图3为用本发明方法求解的动态啮合力时域图。

图4为用本发明方法求解的摩擦力时域图。

图5为用本发明方法求解的主动轮x方向振动位移时域图。

图6为用本发明方法求解的主动轮y方向振动位移时域图。

图7为用本发明方法求解的主动轮扭转振动位移时域图。

图8为用本发明方法求解的从动轮x方向振动位移时域图。

图9为用本发明方法求解的从动轮y方向振动位移时域图。

图10为用本发明方法求解的从动轮扭转振动位移时域图。

具体实施方式

为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。

实施例

本实施例以某渐开线直齿轮传动系统为研究对象，其具体参数如表1所示：

表1渐开线直齿轮传动系统基本参数

如图1所示，本发明提供了一种渐开线直齿轮传动系统动态特性求解方法，包括以下步骤：

步骤S1：基于集中质量法建立渐开线直齿轮传动系统非线性动力学模型；

图2所示的模型中，在主、从动齿轮的理想中心A_i处建立固定坐标系A_i(x_i,y_i,z_i)(i＝1,2)，在轴承的理想中心B_i处建立固定坐标系B_i(x_bi,y_bi,z_bi)(i＝1～4)。坐标轴x_i，x_bi垂直于齿轮啮合线；坐标轴y_i，y_bi平行于齿轮啮合线；坐标轴z_i，z_bi穿过轴承理想中心。主、从动齿轮的旋转中心为O_i(x_i,y_i)(i＝1,2)，质心为G_i(x_gi,y_gi)(i＝1,2)。轴承的旋转中心为O_bi(x_bi,y_bi)(i＝1～4)。所述渐开线直齿轮传动系统非线性动力学模型的输入端转角为

输入端扭转振动角位移为θ_p；主动轮转角为

主动轮扭转振动角位移为θ₁；从动轮转角为

从动轮扭转振动角位移为θ₂；输出端转角为

输出端扭转振动角位移为θ_q，其关系可表示为：

式中，ω₁为主动轮角速度，ω₂为从动轮角速度，t为时间。

旋转中心O_i(x_i,y_i)与质心G_i(x_gi,y_gi)，(i＝1,2)之间的关系为：

式中，ρ₁为主动轮偏心量，ρ₂为从动轮偏心量。

主、从动轴的弹性变形δ_xi,δ_yi(i＝1,2)可表示为：

式中，齿轮位置系数η_i＝l_bi/l_j(i＝1,2,j＝1；i＝3,4,j＝2)；l₁，l₂为主、从动轴在两轴承间的长度；l_bi(i＝1～4)为主、从动轴上齿轮旋转中心到轴承质心的距离。

步骤S2：考虑时变啮合刚度、动态传递误差、摩擦、偏心、修形、间隙、重力和非线性轴承力，基于拉格朗日方程推导渐开线直齿轮传动系统非线性动力学方程；

步骤S21：根据拉格朗日方程，分别建立非保守系统广义坐标X，系统的动能T，势能U，耗散函数R，系统除粘性耗散力以外的非保守广义力P和渐开线直齿轮传动系统非线性动力学方程，其表达式如下：

X＝[θ_p x₁ y₁ θ₁ x₂ y₂ θ₂ x_b1 y_b1 x_b2 y_b2 x_b3 y_b3 x_b4 y_b4 θ_q]^T (4)

上式中，J_p、J_q、J₁、J₂分别表示输入端转动惯量、输出端转动惯量、主动轮转动惯量、从动轮转动惯量，m₁、m₂、m_bi(i＝1～4)分别表示主动轮质量、从动轮质量和四个轴承的质量，r_b1、r_b2为主、从动齿轮基圆半径，L_p1、L_p2、L_g1、L_g2为齿对i摩擦力对主、从动齿轮力臂，ρ₁、ρ₂为主、从动齿轮偏心量，η_i(i＝1～4)为齿轮位置系数，δ_x1、δ_x2、δ_y1、δ_y2为主、从动轴分别沿x,y方向的弹性变形，ε为重合度，g为重力加速度。k_t1、k_t2分别表示中心轴1、2扭转刚度，c_t1、c_t2分别表示中心轴1、2扭转阻尼，k_s1、k_s2分别表示中心轴1、2弯曲刚度，c_s1、c_s2分别表示中心轴1、2弯曲阻尼，k_bxi,k_byi(i＝1～4)表示四个轴承处沿x,y方向的刚度，c_bxi,c_byi(i＝1～4)表示四个轴承处沿x,y方向的阻尼。M_p、M_q分别表示输入扭矩、输出扭矩，F_f1、F_f2为啮合齿对1、2的摩擦力，F_m1、F_m2为啮合齿对1、2的动态啮合力，F_bxi,F_byi(i＝1～4)表示四个轴承处沿x,y方向的非线性轴承力。

步骤S22：计算时变啮合刚度，其计算公式为：

式中，

分别为啮合齿对1、2的时变啮合刚度；平均啮合刚度k_d＝k₀+A₀/2，k₀＝k_max(ε-1)+k_min(2-ε)，A₀＝2Δ(2ε-3)；各谐波分量刚度幅值k_n＝4Δsin(nπ(ε-1))/nπ，Δ＝(k_max-k_min)/2，ε为重合度，一般取到9阶。其中，k_max、k_min分别为时变啮合刚度最大值、最小值；ω_e为啮合频率，ω_e＝2πn₁z₁/60，n₁为主动轮转速，z₁为主动轮齿数；啮合点经过单齿啮合区和双齿啮合区的时间

主动轮转速

mod()为取余函数，t为时间。

步骤S23：计算动态啮合力，根据黏弹性理论，其计算公式为：

式中，F_m1、F_m2为啮合齿对1、2的动态啮合力，c_m为啮合阻尼，

分别为啮合齿对1、2的时变啮合刚度，

mod()为取余函数，t为时间。则由几何关系可知，动态传递误差

其中，静态传递误差e(t)＝e₀+e_rsin(ω_et+φ_e)，e₀、e_r为静态传递误差的均值和波动幅值，φ_e为相位角；啮合点处修形量

分别为从动轮和主动轮啮合点处修形量，其中：

式中，Δ_max为最大修形量，L为修形长度，啮合周期t₁＝ε·t₀，t₀为啮合点经过单齿啮合区和双齿啮合区的时间，ω₁为主动轮角速度，r_b1为主动轮基圆半径，

mod()为取余函数，t为时间。

f(δ)为间隙函数：

式中，b为间隙，δ为动态传递误差。

步骤S24：计算摩擦力，其计算公式为：

式中，F_f1、F_f2为啮合齿对1、2的摩擦力，μ为摩擦系数，F_m1、F_m2为啮合齿对1、2的动态啮合力，方向系数

和

的计算公式为：

式中，齿轮从进入啮合到节圆处啮合时间

t₁为啮合周期，t₀为啮合点经过单齿啮合区和双齿啮合区的时间，ω₁为主动轮角速度，r_b1和r_b2为主、从动齿轮基圆半径，r₁为主动轮分度圆半径，r_a2为从动轮齿顶圆半径，α为节圆压力角，

mod()为取余函数，t为时间。

步骤S25：计算摩擦力臂，所述摩擦力臂包括齿对i啮合点距主动轮摩擦力臂L_pi和啮合点距从动轮摩擦力臂L_gi，其计算公式分别为：

式中，r_b1和r_b2为主、从动齿轮基圆半径，r_a2为从动轮齿顶圆半径，α为节圆压力角，ω₁为主动轮角速度，t₀为啮合点经过单齿啮合区和双齿啮合区的时间，

mod()为取余函数，t为时间。

步骤S26：计算非线性轴承力沿x方向分力F_bx和y方向分力F_by，其计算公式为：

式中，K_c为赫兹接触刚度系数；N_b为滚珠个数；

为δ_j的p次幂，当轴承为球轴承时，指数p取3/2，当轴承为滚子轴承时，指数p取10/9；H(δ_j)为Heaviside函数，δ_j≤1时，H(δ_j)＝0；δ_j＞时，H(δ_j)＝1。

第j个滚动体与滚道的法向接触变形量δ_j为：

δ_j＝x_bjcosθ_j+y_bjsinθ_j-γ₀ (21)

式中，x_bj、y_bj为内圈中心振动位移，γ₀为轴承游隙。第j个滚动体在时间t内的转动角度为θ_j：

θ_j＝ω_b·t+2π(j-1)/N_b (22)

式中，保持架的角速度ω_b＝ω×r′/(R′+r′)，ω为轴承角速度，R′、r′为轴承内外圈半径。

步骤S3：基于Runger-Kutta法求解渐开线直齿轮传动系统动态特性。图3～图10分别给出了动态啮合力、摩擦力以及主、从动齿轮沿x，y轴横向振动位移和绕z轴扭转振动位移时域图。本发明提出的数值法(在i5处理器和8GB内存计算机中运行耗时为3min)比ANSYS和ADAMS联合仿真(在i5处理器和8GB内存计算机中运行共耗时为6.7h)效率高。文献《金属橡胶复合齿轮副振动特性分析及实验》的研究结果表明仿真解与实验值的偏差在5％，而采用该文献数据，应用本发明方法得出数值解与实验值的偏差在3％左右。因此，本发明提出的方法高效、更准确。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (2)

1.一种渐开线直齿轮传动系统动态特性求解方法，其特征在于，包括以下步骤：

步骤S1：基于集中质量法建立渐开线直齿轮传动系统非线性动力学模型；

步骤S2：考虑时变啮合刚度、动态传递误差、摩擦、偏心、修形、间隙、重力和非线性轴承力，基于拉格朗日方程推导渐开线直齿轮传动系统非线性动力学方程；

所述步骤S2的具体步骤如下：

步骤S21：根据拉格朗日方程，分别建立非保守系统广义坐标X，系统的动能T，势能U，耗散函数R，系统除粘性耗散力以外的非保守广义力P和渐开线直齿轮传动系统非线性动力学方程，其表达式如下：

X＝[θ_p x₁ y₁ θ₁ x₂ y₂ θ₂ x_b1 y_b1 x_b2 y_b2 x_b3 y_b3 x_b4 y_b4 θ_q]^T (4)

上式中，J_p、J_q、J₁、J₂分别表示输入端转动惯量、输出端转动惯量、主动轮转动惯量、从动轮转动惯量，m₁、m₂、m_bi(i＝1～4)分别表示主动轮质量、从动轮质量和四个轴承的质量，r_b1、r_b2为主、从动齿轮基圆半径，L_p1、L_p2、L_g1、L_g2为齿对1、2摩擦力对主、从动齿轮力臂，ρ₁、ρ₂为主、从动齿轮偏心量，η_i(i＝1～4)为齿轮位置系数，δ_x1、δ_x2、δ_y1、δ_y2为主、从动轴分别沿x，y方向的弹性变形，ε为重合度，g为重力加速度；k_t1、k_t2分别表示中心轴1、2扭转刚度，c_t1、c_t2分别表示中心轴1、2扭转阻尼，k_s1、k_s2分别表示中心轴1、2弯曲刚度，c_s1、c_s2分别表示中心轴1、2弯曲阻尼，k_bxi,k_byi(i＝1～4)表示四个轴承处沿x,y方向的刚度，c_bxi,c_byi(i＝1～4)表示四个轴承处沿x,y方向的阻尼；M_p、M_q分别表示输入扭矩、输出扭矩，F_f1、F_f2为啮合齿对1、2的摩擦力，F_m1、F_m2为啮合齿对1、2的动态啮合力，F_bxi,F_byi(i＝1～4)表示四个轴承处沿x,y方向的非线性轴承力；

步骤S22：计算时变啮合刚度，其计算公式为：

式中，

分别为啮合齿对1、2的时变啮合刚度；平均啮合刚度k_d＝k₀+A₀/2，k₀＝k_max(ε-1)+k_min(2-ε)，A₀＝2Δ(2ε-3)；各谐波分量刚度幅值k_n＝4Δsin(nπ(ε-1))/nπ，，Δ＝(k_max-k_min)/2，ε为重合度，取到9阶；其中，k_max、k_min分别为时变啮合刚度最大值、最小值；ω_e为啮合频率，ω_e＝2πn₁z₁/60，n₁为主动轮转速，z₁为主动轮齿数；啮合点经过单齿啮合区和双齿啮合区的时间

主动轮转速

mod()为取余函数，t为时间；

步骤S23：计算动态啮合力，根据黏弹性理论，其计算公式为：

式中，F_m1、F_m2为啮合齿对1、2的动态啮合力，c_m为啮合阻尼，

分别为啮合齿对1、2的时变啮合刚度，

mod()为取余函数，t为时间；则由几何关系可知，动态传递误差

其中，静态传递误差e(t)＝e₀+e_rsin(ω_et+φ_e)，e₀、e_r为静态传递误差的均值和波动幅值，φ_e为相位角；啮合点处修形量

分别为从动轮和主动轮啮合点处修形量，其中：

式中，Δ_max为最大修形量，L为修形长度，啮合周期t₁＝ε·t₀，t₀为啮合点经过单齿啮合区和双齿啮合区的时间，ω₁为主动轮角速度，r_b1为主动轮基圆半径，

mod()为取余函数，t为时间；

f(δ)为间隙函数：

式中，b为间隙，δ为动态传递误差；

步骤S24：计算摩擦力，其计算公式为：

式中，F_f1、F_f2为啮合齿对1、2的摩擦力，μ为摩擦系数，F_m1、F_m2为啮合齿对1、2的动态啮合力，方向系数

和

的计算公式为：

式中，齿轮从进入啮合到节圆处啮合时间

t₁为啮合周期，t₀为啮合点经过单齿啮合区和双齿啮合区的时间，ω₁为主动轮角速度，r_b1和r_b2为主、从动齿轮基圆半径，r₁为主动轮分度圆半径，r_a2为从动轮齿顶圆半径，α为节圆压力角，

mod()为取余函数，t为时间；

步骤S25：计算摩擦力臂，所述摩擦力臂包括齿对i的啮合点距主动轮摩擦力臂L_pi和啮合点距从动轮摩擦力臂L_gi，其计算公式分别为：

式中，r_b1和r_b2为主、从动齿轮基圆半径，r_a2为从动轮齿顶圆半径，α为节圆压力角，ω₁为主动轮角速度，t₀为啮合点经过单齿啮合区和双齿啮合区的时间，

mod()为取余函数，t为时间；

步骤S26：计算非线性轴承力沿x方向分力F_bx和y方向分力F_by，其计算公式为：

式中，K_c为赫兹接触刚度系数；N_b为滚珠个数；

为δ_j的p次幂，当轴承为球轴承时，指数p取3/2，当轴承为滚子轴承时，指数p取10/9；H(δ_j)为Heaviside函数，δ_j≤0时，H(δ_j)＝0；δ_j＞0时，H(δ_j)＝1；

第j个滚动体与滚道的法向接触变形量δ_j为：

δ_j＝x_bjcosθ_j+y_bjsinθ_j-γ₀ (21)

式中，x_bj、y_bj为内圈中心振动位移，γ₀为轴承游隙；第j个滚动体在时间t内的转动角度为θ_j：

θ_j＝ω_b·t+2π(j-1)/N_b (22)

式中，保持架的角速度ω_b＝ω×r′/(R′+r′)，ω轴承角速度，R′、r′为轴承内外圈半径；

步骤S3：基于Runger-Kutta法求解渐开线直齿轮传动系统动态特性。

2.根据权利要求1所述的渐开线直齿轮传动系统动态特性求解方法，其特征在于，所述步骤S1建立的渐开线直齿轮传动系统非线性动力学模型中，在主、从动齿轮的理想中心A_i处建立固定坐标系A_i(x_i,y_i,z_i)(i＝1,2)，在轴承的理想中心B_i处建立固定坐标系B_i(x_bi,y_bi,z_bi)(i＝1～4)；坐标轴x_i，x_bi垂直于齿轮啮合线；坐标轴y_i，y_bi平行于齿轮啮合线；坐标轴z_i，z_bi穿过轴承理想中心；主、从动齿轮的旋转中心为O_i(x_i,y_i)(i＝1,2)，质心为G_i(x_gi,y_gi)(i＝1,2)；轴承的旋转中心为O_bi(x_bi,y_bi)(i＝1～4)；所述渐开线直齿轮传动系统非线性动力学模型的输入端转角为

输入端扭转振动角位移为θ_p；主动轮转角为

主动轮扭转振动角位移为θ₁；从动轮转角为

从动轮扭转振动角位移为θ₂；输出端转角为

输出端扭转振动角位移为θ_q，其关系可表示为：

式中，ω₁为主动轮角速度，ω₂为从动轮角速度，t为时间；

旋转中心O_i(x_i,y_i)(i＝1,2)与质心G_i(x_gi,y_gi)(i＝1,2)之间的关系为：

式中，ρ₁为主动轮偏心量，ρ₂为从动轮偏心量；

主、从动轴沿x,y方向的弹性变形δ_xi,δ_yi(i＝1,2)可表示为：

式中，齿轮位置系数η_i＝l_bi/l_j(i＝1,2,j＝1；i＝3,4,j＝2)；l₁，l₂为主、从动轴在两轴承间的长度；l_bi(i＝1～4)为主、从动轴上齿轮旋转中心到轴承质心的距离。

Images