CN105928707A - 一种滚动轴承-转子系统动态耦合建模方法 - Google Patents

Abstract

本发明涉及一种滚动轴承‑转子系统动态耦合建模方法，包括：1)收集系统的属性参数；2)计算每个部件的初始速度和加速度；3)建立转子的有限元模型；4)计算每个部件上的合力和合力矩；5)将合力代入到部件的动力学方程获得每个部件的刚体运动速度和加速度；6)求解转子的有限元模型，获得转子的弹性振动；7)将转子的弹性振动和转子的刚体运动进行叠加；8)计算整个轴承‑转子系统的时域响应。本发明考虑了每个部件的三维运动、滑动牵引力、轴承间隙、套圈局部故障和套圈波纹度，同时应用转子的有限元模型来求解转子的弹性变形，并与系统的动力学模型进行耦合，最终得到轴承‑转子系统的时域响应。

Description

一种滚动轴承-转子系统动态耦合建模方法

技术领域

本发明属于机械动力学领域，涉及一种滚动轴承-转子系统动态耦合建模方法。

背景技术

高速主轴、航空发动机、高速列车轮毂等滚动轴承支撑的旋转机械设备应用广泛，其关键部件包含一个由滚动轴承和转子组成的滚动轴承-转子系统。为了定量研究轴承损伤、间隙、转子碰磨等常见故障的振动机理，准确预测并分析转子的非线性振动响应、不平衡引起的同步回旋、高速状态下的失稳、固有频率随系统参数的变化等动力学行为，有必要建立整个滚动轴承-转子系统的动力学模型。通过模型求解和仿真研究，建立系统参数与动态响应特征的联系，从而为故障及振动机理研究提供有效的理论指导。总而言之，滚动轴承-转子系统动力学模型可以应用在故障诊断、机理研究、参数优化等诸多方面，是定量分析和研究系统动力学特性的关键步骤。

为了研究轴承-转子系统的动力学行为，须建立转子和滚动轴承的耦合模型。Zhang等(Zhang X,Han Q,Peng Z,et al.A new nonlinear dynamic model of the rotor-bearing systemconsidering preload and varying contact angle of the bearing[J].Commun.Nonlinear Sci.Numer.Simulat.2014,22(1-3):821-841.)基于Jones模型建立了一个滚动轴承-刚性转子动力学模型，每个安装在转子上的轴承都有5个自由度，施加在转子上的合力通过Jones模型来求解。吕运等(吕运，童大鹏，田野，赖亚辉.滚动轴承-转子系统动力学建模与仿真分析[J].机械强度，2005,37(6):1178-1185)着重考虑滚动轴承的接触非线性对滚动轴承-转子系统动力学特性的影响，将基座、轴承的振动纳入到研究范围内。张玉言等(张玉言,王晓力,闫晓亮.高速滚动轴承-转子系统的动力学特性研究[C].第十一届全国摩擦学大会论文集.2013.)考虑滚动体与滚道间的弹流润滑油膜，研究了高速滚动轴承平衡转子系统的动力学特性。袁茹等(袁茹,赵凌燕,王三民.滚动轴承转子系统的非线性动力学特性分析[J].机械科学与技术,2004,23(10):1175-1177.)在计及轴承接触非线性和径向间隙的条件下，建立了滚动轴承支承的水平刚性转子系统的非线性动力学模型。

从文献调研中可以发现，在众多轴承模型中，Gupta的动力学模型考虑了每个部件的三维运动、滑动牵引力和保持架效应，是目前考虑因素最全面、系统的动力学模型，然而Gupta模型到目前还没有应用到滚动轴承-转子系统动力学建模中。同时，转子的弹性变形通常被忽略，这在当转子受到重载或在高速状态下是不准确的。因此要想建立准确的滚动轴承-转子系统动力学模型，就必须考虑因素全面的轴承模型，且不能忽略转子的弹性变形。

发明内容

本发明的目的在于基于Gupta模型建立一个可以计算滚动轴承-转子系统刚体运动的动力学模型，并将其与转子有限元模型进行耦合来考虑转子的弹性振动，最终得到一种滚动轴承-转子系统动态耦合建模方法。

为了实现上述目的，本发明采用如下技术方案：

一种滚动轴承-转子系统动态耦合建模方法，包括以下步骤：

1)收集滚动轴承-转子系统的属性参数；

2)根据Jones拟静力学模型，求解滚动轴承-转子系统中每个部件的速度和位移；

3)根据转子的物理属性参数来建立转子的有限元模型；

4)根据轴承各部件间的相互作用计算作用在转子上的合力和合力矩；

5)将计算得到的作用在转子上的合力代入到滚动轴承动力学模型获得滚动轴承-转子系统中每个部件的刚体运动速度和加速度；

6)使用Newmarkβ方法求解转子的有限元模型，获得转子的弹性振动；

7)将转子的弹性振动和刚体运动进行叠加，获得转子的耦合运动；

8)轴承-转子系统所有部件的速度和加速度求得后，使用变步长四阶龙格-库塔-费尔伯格方法得到下一时刻滚动轴承-转子系统中所有部件的速度和位移，重复进行直到达到预先设定的时刻，得到整个轴承-转子系统的时域响应。

进一步的，步骤4)中，滚球和套圈之间的合力为：

${F_k}^c={Q_k}^c+{f_k}^c-c_{br}{v}_{rbk}^c---\left(16\right)$

式中：c_br——由于接触区域的润滑引起的阻尼系数/N.s.m^-1；Q_k ^c——垂直于滚球和套圈接触面的接触力Q_k ^c，f_k ^c——平行于滚球和套圈接触面的牵引力；

作用在套圈上的力等于作用在滚球上的力，且两者方向相反；滚球质心的力矩和关于套圈质心的合力矩为

$\left\{\begin{array}{c}{M}_{bk}^c=r_{cpk}^c\×F_k^c\\ M_{rk}^{rk}=\underset{j=1}{\overset{z_k}{\Σ}}{r}_{prkj}^{rk}\×T_{c,rk}{F}_{kj}^c\end{array}\right.---\left(17\right)$

式中：r_cpk和r_prk——某接触点相对于滚球中心和套圈中心的位置向量；下标j——轴承的第j个滚球；上标rk——第k个套圈坐标系；上标c——接触坐标系；T_c,rk——从接触坐标系到套圈坐标系的变换矩阵；

作用在转子上的合力可以表示为

$F_r^i=\underset{k=1}{\overset{n}{\Σ}}{F}_{rk}^i+\underset{k=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{z_k}{\Σ}}{T}_{ci}{F}_{kj}^c+G_r^i---\left(18\right)$

式中：n——装配在转子上轴承的个数；T_ci——从接触坐标系到惯性坐标系的变换矩阵；——第k个轴承的滚球作用在内圈上的合力；——转子的重力；

作用在转子质心的合力矩为

$M_r^r=\underset{k=1}{\overset{n}{\Σ}}(M_{rk}^{rk}+M_{fk}^r)=\underset{k=1}{\overset{n}{\Σ}}(\underset{j=1}{\overset{z_k}{\Σ}}{r}_{prkj}^{rk}\×T_{cr}{F}_{kj}^c+\stackrel{\‾}{O_r{O}_{rk}}\×\underset{j=1}{\overset{z_k}{\Σ}}{T}_{cr}{F}_{kj}^c)---\left(19\right)$

式中：——由第k个轴承滚球和内圈之间的合力引起的力矩；T_cr——从接触坐标系到转子定体坐标系的转换矩阵。

进一步的，步骤5)中，在惯性坐标系中描述滚球的平动运动，在笛卡尔坐标系中描述套圈、保持架和转子的平动运动；将轴承外圈和轴承座孔的相互作用建模成N_p个均匀支撑在外圈上的弹簧和阻尼并列出外圈和轴承座的动力学方程，同时考虑外圈可活动时的动力学方程和旋转部件在各自定体坐标系中的欧拉方程，求出平动加速度和角加速度。

进一步的，步骤7)中，转子节点处的速度v_p和加速度a_p由下式获得

$\left\{\begin{array}{c}{v}_p=v_{ri}+v_{el}\\ a_p=a_{ri}+a_{el}\end{array}\right.---\left(27\right)$

式中：v_ri和a_ri——刚体运动的速度/m.s^-1和加速度/m.s^-2，通过动力学模型计算；v_el和a_el——弹性振动的速度/m.s^-1和加速度/m.s^-2，通过有限元模型求得；在每一个计算步长，轴承滚球和内圈之间的合力矩M_rk、以及惯性力F_el和惯性力矩M_el通过动力学模型求解，这些力施加在转子有限元模型上；获得转子的弹性振动后，轴承配合处节点的速度和加速度通过式(27)计算。

进一步的，步骤1)中所述属性参数包括：几何参数、材料属性、运行状态以及润滑模型。

进一步的，在步骤4)、5)中，基于Gupta模型建立了一个可以计算滚动轴承-转子系统刚体运动的动力学模型。Gupta轴承动力学模型考虑了轴承间隙、润滑效应、保持架效应等对轴承动力学行为有很大影响的因素，是目前考虑因素最为完备的轴承模型之一。因此，本发明考虑了每个部件(外圈、内圈、滚球和转子)的三维运动、滑动牵引力、轴承间隙、套圈局部故障和套圈波纹度，基于Gupta滚动轴承模型提出了一个滚动轴承-转子系统动力学模型。

相对于现有技术，本发明具有以下优点：Gupta轴承动力学模型考虑了轴承间隙、润滑效应、保持架效应等对轴承动力学行为有很大影响的因素，能够比较完全地描述轴承的振动特征，是目前最为完备的轴承模型之一。本发明基于Gupta模型建立了一个可以计算滚动轴承-转子系统刚体运动的动力学模型，并将其与转子有限元模型进行耦合来考虑转子的弹性振动，使得转子在受到重载或在高速状态下的时域响应更加准确。

附图说明

图1是本发明的转子有限元模型。

图2是本发明的轴承-转子系统几何相互作用图。

图3是本发明的轴承-转子试验台装配简图。

图4是本发明的转子有限元模型实例。

图5是本发明的节点8仿真和实验的振动响应比较。

图6是本发明的流程示意图。

具体实施方式

请参阅图1至图6所示，本发明一种滚动轴承-转子系统动态耦合建模方法，包括以下步骤：

1)收集滚动轴承-转子系统的属性参数，包括：几何参数、材料属性、运行状态以及润滑模型，为轴承-转子系统的建模提供数据支持；

2)根据Jones拟静力学模型，求解滚动轴承-转子系统中每个部件的速度和位移；

3)根据转子的物理属性参数来建立转子的有限元模型；

4)根据轴承各部件间的相互作用计算作用在转子上的合力和合力矩；

5)将计算得到的作用在转子上的合力代入到滚动轴承动力学模型获得滚动轴承-转子系统中每个部件的刚体运动速度和加速度；

6)使用Newmarkβ方法求解转子的有限元模型，获得转子的弹性振动；

7)将转子的弹性振动和刚体运动进行叠加，也就是将滚动轴承动力学模型与转子有限元模型进行耦合，获得转子的耦合运动，这是本发明最具创新性的一步；

8)轴承-转子系统所有部件的速度和加速度求得后，使用变步长四阶龙格-库塔-费尔伯格方法可以得到下一时刻滚动轴承-转子系统中所有部件的速度和位移，重复进行直到达到预先设定的时刻，便可以得到整个轴承-转子系统的时域响应。

步骤3)和步骤6)中，使用有限元方法求解转子的弹性振动。有限元模型只有当施加了边界条件之后才能进行求解。假设轴承仅仅限制了轴承配合处节点的平动自由度，而相应的转动自由度是自由的。如图1所示，第一、第二和第k个轴承分别装配在节点2、4和i，它们的平动自由度都被限制。由于轴承配合处节点的转动自由度没有被限制，滚球和内圈之间的合力矩M_r1、M_r2...M_rk将会对转子的弹性振动有贡献，如图1所示。

应用达朗贝尔原理可以获得由刚体运动产生的作用在转子上的惯性力。F_el和M_el分别是作用在转子第l个单元(图1只给出了作用在第一个单元上的惯性力)上的惯性力和惯性力矩，它们可以表示为

$\left\{\begin{array}{c}{F}_{el}=-m_{el}{a}_{el}\\ M_{el}=-J_{el}{\mathrm{\β}}_{el}\end{array}\right.---\left(1\right)$

式中：m_el——第l个单元的质量/kg；J_el——第l个单元的转动惯量/kg.m²；a_el——第l个单元刚体运动的加速度/m.s^-2；β_el——第l个单元刚体运动的角加速度/rad.s^-2。

此外，也考虑了转子的重力

G_el＝m_elg (2)

式中：g——重力加速度/m.s^-2。

作用在转子上的合力F(t)可以表示为

$F\left(t\right)=F_o+\underset{l=1}{\overset{m}{\Σ}}(F_{el}+M_{el}+G_{el})+\underset{k=1}{\overset{n}{\Σ}}{M}_{rk}---\left(3\right)$

式中：F_o——外部力/N；m——单元数量；rk——套圈定体坐标系；k——第k个轴承的套圈。

所述步骤4)中，轴承滚球和转子几何相互作用如图2所示，位矢r_brk由滚球中心和套圈中心确定，使用上标i表示惯性坐标系，r_brk可表示为

$r_{brk}^i=r_{bk}^i-r_{rk}^i---\left(4\right)$

式中：表示惯性坐标系中滚球几何中心的位移矢量；表示惯性坐标系中套圈几何中心的位移矢量。

使用上标rk表示第k个轴承的套圈定体坐标系，则位矢r_brk从惯性坐标系和套圈定体坐标系的向量变换可表示为

$r_{brk}^{rk}=T_{ir}{r}_{brk}^i---\left(5\right)$

式中：T_ir——惯性坐标系和套圈定体坐标系之间的相互转换矩阵，可以根据套圈的姿态角确定。

然后在滚道定体坐标系中描述的滚球方位角可表示为

${\mathrm{\θ}}_{brk}^{rk}=arctan\left(\frac{-r_{brk2}^{rk}}{r_{brk3}^{rk}}\right)---\left(6\right)$

式中：下标2和3——向量的第二和第三分量。

则在套圈定体坐标系中向量可以表示为

$r_{crk}^{rk}=\left\{\begin{array}{c}0\\ -R_f{\mathrm{sin\θ}}_{brk}^{rk}\\ R_f{\mathrm{cos\θ}}_{brk}^{rk}\end{array}\right\}---\left(7\right)$

式中：R_f——套圈滚道中心轨迹半径/m，对于内圈

$R_f=\frac{d_m}{2}+(f-0.5)D{\mathrm{cos\α}}_0{\mathrm{\Δ}}_i-p_i---\left(8\right)$

对于外圈

$R_f=\frac{d_m}{2}-(f_0-0.5)D{\mathrm{cos\α}}_0+{\mathrm{\Δ}}_0+p_0---\left(9\right)$

式中：d_m——节圆半径/m；f_i和f_o——内圈和外圈的滚道曲率因子；D——滚球直径/m；α₀——初始接触角/rad；Δ_i和Δ_o——内外圈间隙/m；p_i和p_o——由于内圈滚道波纹度和外圈滚道波纹度引起的位移/m。

滚道曲率中心R_ck和滚球中心O_ak的相对位置为

$r_{bck}^i=r_{brk}^i-r_{crk}^i---\left(10\right)$

若滚球和套圈之间存在任何相互作用，则接触负荷将沿着r_bck作用，因此接触角α可以表示为

$\mathrm{\α}=arctan\left(\frac{r_{bck1}^{ak}}{r_{bck3}^{ak}}\right)---\left(11\right)$

式中：上标ak——滚球/套圈方位坐标系，可以通过旋转套圈定体坐标系的x轴角获得：下标1和3——向量的第一和第三分量。

为了方便，设定了一个称为接触坐标系c(x^c,y^c,z^c)的新坐标系，其原点和接触区域中心重合。一旦向量通过接触角α从滚球/套圈方位坐标系ak转换到接触坐标系c，滚球和套圈之间的弹性变形可表示为

$\mathrm{\δ}=r_{bck3}^c-(f-0.5)D-{\mathrm{\δ}}_d---\left(12\right)$

式中：f——滚道曲率因子；下标3——向量的第三个分量；δ_d——由于轴承套圈局部故障而引起的额外位移。

滚球和滚道之间的接触力可以通过赫兹点接触理论求解

$Q=\left\{\begin{array}{cc}{\mathrm{K\δ}}^{1.5}& \mathrm{\δ}>0\\ 0& \mathrm{\δ}\≤0\end{array}\right.---\left(13\right)$

式中：K——赫兹接触刚度系数/N.m^-1。

在图2中，当求得了转子在惯性坐标系中的速度和在转子定体坐标系的角速度在套圈中心O_rk处的速度和角速度可表示为

$\left\{\begin{array}{c}{v}_{rk}^i=v_r^i+T_{ri}{\mathrm{\ω}}_r^r\×\stackrel{\‾}{O_r{O}_{rk}}\\ {\mathrm{\ω}}_{rk}^{rk}={\mathrm{\ω}}_r^r\end{array}\right.---\left(14\right)$

式中：T_ri——从转子定体坐标系到惯性坐标系的转换矩阵，可根据转子的姿态角确定。

假设和分别是滚球和套圈在接触坐标系中的速度，则滚球和套圈在接触点处的相对速度可表示为

$v_{rbk}^c=v_{prk}^c-v_{pbk}^c---\left(15\right)$

滑动速度的牵引系数κ_x和κ_y可由轴承润滑模型确定，然后牵引力f^c等于径向力和牵引系数的乘积。当垂直于接触面的接触力和与接触面平行的牵引力求得后，滚球和套圈之间的合力便可以求得

$F_k^c=Q_k^c+f_k^c-c_{br}{v}_{rbk}^c---\left(16\right)$

式中：c_br——由于接触区域的润滑引起的阻尼系数/N.s.m^-1。

作用在套圈上的力等于作用在滚球上的力，且两者方向相反。此外，关于滚球质心的力矩和关于套圈质心的合力矩可表示为

$\left\{\begin{array}{c}{M}_{bk}^c=r_{cpk}^c\×F_k^c\\ M_{rk}^{rk}=\underset{j=1}{\overset{z_k}{\Σ}}{r}_{prkj}^{rk}\×T_{c,rk}{F}_{kj}^c\end{array}\right.---\left(17\right)$

式中：z_k——第k个轴承的滚球个数；r_cpk和r_prk——某接触点相对于滚球中心和套圈中心的位置向量；下标j——轴承的第j个滚球；上标rk——第k个套圈坐标系；上标c——接触坐标系；T_c,rk——从接触坐标系到套圈坐标系的变换矩阵。

作用在转子上的合力可以表示为

$F_r^i=\underset{k=1}{\overset{n}{\Σ}}{F}_{\mathrm{rk}}^i+\underset{k=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{z_k}{\Σ}}{T}_{ci}{F}_{kj}^c+G_r^i---\left(18\right)$

式中：n——装配在转子上轴承的个数；T_ci——从接触坐标系到惯性坐标系的变换矩阵；——第k个轴承的滚球作用在内圈上的合力；——转子的重力。

作用在转子质心的合力矩可表示为

$M_r^r=\underset{k=1}{\overset{n}{\Σ}}(M_{rk}^{rk}+M_{fk}^r)=\underset{k=1}{\overset{n}{\Σ}}(\underset{j=1}{\overset{z_k}{\Σ}}{r}_{prkj}^{rk}\×T_{cr}{F}_{kj}^c+\stackrel{\‾}{O_r{O}_{rk}}\×\underset{j=1}{\overset{z_k}{\Σ}}{T}_{cr}{F}_{kj}^c)---\left(19\right)$

式中：——由第k个轴承滚球和内圈之间的合力引起的力矩；T_cr——从接触坐标系到转子定体坐标系的转换矩阵。

所述步骤5)中，轴承-转子系统部件的运动可以分解成质心的平动运动和绕质心的转动。不同部件的动力学方程可以根据其运动特点在不同的坐标系中描述。滚球的平动运动在惯性柱坐标系中描述比较方便

$\left\{\begin{array}{c}{m}_b\stackrel{\·\·}{x}=F_x\\ m_b\stackrel{\·\·}{r}-m_{b}r{\stackrel{\·}{\mathrm{\θ}}}^2=F_r\\ m_{b}r\stackrel{\·\·}{\mathrm{\θ}}+2m_b\stackrel{\·}{r}\stackrel{\·}{\mathrm{\θ}}=F_{\mathrm{\θ}}\end{array}\right.---\left(20\right)$

式中：m_b——滚球的质量/kg；(F_x,F_r,F_θ)——在柱坐标系三个方向上的合力分量/N。

套圈、保持架和转子的平动运动在笛卡尔坐标系中描述比较方便

$\left\{\begin{array}{c}{m}_r\stackrel{\·\·}{x}=F_x\\ m_r\stackrel{\·\·}{y}=F_y\\ m_r\stackrel{\·\·}{z}=F_z\end{array}\right.---\left(21\right)$

式中：m_r——套圈、保持架或者转子的质量/kg；(F_x,F_y,F_z)——笛卡尔坐标系三个方向上的合力分量/N。

轴承外圈和轴承座孔的相互作用被建模成N_p个均匀支承在外圈上的弹簧和阻尼，外圈y和z方向的动力学方程为

$\{\begin{array}{c}{m}_{or}\stackrel{\·\·}{y}+\underset{i=1}{\overset{N_p}{\Σ}}\[k_{pi}{(y{\mathrm{cos\θ}}_{pi}+z{\mathrm{sin\θ}}_{pi}-r_{rp})}_{+}+c_{pi}(\stackrel{\·}{y}{\mathrm{cos\θ}}_{pi}+\stackrel{\·}{z}{\mathrm{sin\θ}}_{pi})\]{\mathrm{cos\θ}}_{pi}=F_y\\ m_{or}\stackrel{\·\·}{z}+\underset{i=1}{\overset{N_p}{\Σ}}{\[k_{pi}(y{\mathrm{cos\θ}}_{pi}+z{\mathrm{sin\θ}}_{pi}-r_{rp})\]}_{+}+c_{pi}(\stackrel{\·}{y}{\mathrm{cos\θ}}_{pi}+\stackrel{\·}{z}{\mathrm{sin\θ}}_{pi})\]{\mathrm{sin\θ}}_{pi}=F_z\end{array}---\left(22\right)$

式中：m_or——轴承外圈质量/kg；k_pi、c_pi和θ_pi——弹簧刚度/N.m^-1、阻尼/N.s.m^-1和方位角/rad；r_rp——轴承座孔和外圈之间的间隙/m；右端括号的下标“+”表示当括号中的值大于0时，弹簧是压缩的，有恢复力存在；而当括号中的值小于0时，括号中的值取0，此时弹簧未压缩，不存在恢复力。而轴承座的动力学方程可表示为

$\left\{\begin{array}{c}{m}_p\stackrel{\·\·}{y}+c_{py}\stackrel{\·}{y}+k_{py}y=\\ \underset{i=1}{\overset{N_p}{\Σ}}\[k_{pi}{(y{\mathrm{cos\θ}}_{pi}+z{\mathrm{sin\θ}}_{pi}-r_{rp})}_{+}+c_{pi}(\stackrel{\·}{y}{\mathrm{cos\θ}}_{pi}+\stackrel{\·}{z}{\mathrm{sin\θ}}_{pi})\]{\mathrm{cos\θ}}_{pi}\\ m_p\stackrel{\·\·}{z}+c_{pz}\stackrel{\·}{z}+k_{pz}z=\\ \underset{i=1}{\overset{N_p}{\Σ}}\[k_{pi}{(y{\mathrm{cos\θ}}_{pi}+z{\mathrm{sin\θ}}_{pi}-r_{rp})}_{+}+c_{pi}(\stackrel{\·}{y}{\mathrm{cos\θ}}_{pi}+\stackrel{\·}{z}{\mathrm{sin\θ}}_{pi})\]{\mathrm{sin\θ}}_{pi}\end{array}\right.---\left(23\right)$

式中：m_p——轴承座质量/kg；c_py和c_pz——y和z方向的阻尼系数/N.s.m^-1；k_py和k_pz——y和z方向的刚度/N.m^-1。

若外圈x方向是活动的，则相应的动力学方程为

$m_{or}\stackrel{\·\·}{x}+c_{px}\stackrel{\·}{x}=F_x---\left(24\right)$

式中：c_px——x方向阻尼/N.s.m^-1。

对于任何部件的旋转运动可在各自的定体坐标系中用欧拉方程来描述

$\left\{\begin{array}{c}{I}_1{\stackrel{\·}{\mathrm{\ω}}}_1-(I_2-I_3){\mathrm{\ω}}_2{\mathrm{\ω}}_3=M_1\\ I_2{\stackrel{\·}{\mathrm{\ω}}}_2-(I_3-I_1){\mathrm{\ω}}_3{\mathrm{\ω}}_1=M_2\\ I_3{\stackrel{\·}{\mathrm{\ω}}}_3-(I_1-I_2){\mathrm{\ω}}_1{\mathrm{\ω}}_2=M_3\end{array}\right.---\left(25\right)$

式中：(I₁,I₂,I₃)——主惯性质量/kg.m²；(ω₁,ω₂,ω₃)——角速度分量/rad.s^-1；(M₁,M₂,M₃)——合力矩分量/N.m。

当速度、角速度、合力和合力矩求得后，便可以通过(20)～(25)来求得平动加速度和角加速度。

由于轴承内圈固定在转子上，套圈中心的加速度由转子质心的加速度确定。假设O_rk处的平动加速度和角加速度分别是α_rk和β_rk，转子质心O_r处的平动加速度和角加速度分别是α_r和β_r，则O_rk处的平动加速度和角加速度为

$\{\begin{array}{c}{\mathrm{\α}}_{rk}^i={\mathrm{\α}}_r^i+T_{ri}{\mathrm{\β}}_r^r\×\stackrel{\‾}{O_r{O}_{rk}}+T_{ri}{\mathrm{\ω}}_r^r\×({\mathrm{\ω}}_r^r\×\stackrel{\‾}{O_r{O}_{rk}})\\ {\mathrm{\β}}_{rk}^{rk}={\mathrm{\β}}_r^r\end{array}--\left(26\right)$

所述步骤7)中，应用转子的有限元模型与系统的动力学模型进行了耦合。转子节点处的速度v_p和加速度a_p可以由下式获得

$\left\{\begin{array}{c}{v}_p=v_{ri}+v_{el}\\ a_p=a_{ri}+a_{el}\end{array}\right.---\left(27\right)$

式中：v_ri和a_ri——刚体运动的速度/m.s^-1和加速度/m.s^-2，可以通过动力学模型计算；而v_el和a_el——弹性振动的速度/m.s^-1和加速度/m.s^-2，可以通过有限元模型求得。这两个模型可以同时进行求解。在每一个计算步长，轴承滚球和内圈之间的合力矩M_rk、以及惯性力F_el和惯性力矩M_el可以通过动力学模型求解，这些力是施加在转子有限元模型上的。另一方面，一旦获得转子的弹性振动，轴承配合处节点的速度和加速度可以通过式(27)计算，这些值被代入到动力学模型来计算下一时刻的刚体运动。

下面结合一个实例对本发明的一种基于滚动轴承动力学模型与转子有限元模型耦合建模的方法作进一步详细说明，同时验证本发明在工程应用中的有效性，但本实例并不用于限制本发明。

轴承-转子试验台的装配简图如图3所示。转子由两个角接触轴承(Timken LM11749)支承。装有不平衡质量的转盘固定在转子上。当转子旋转时，偏心质量将会使转子产生振动。转子在两个位置上的振动响应分别由两个Lion位移传感器(灵敏度：80mV/μm)来测量。

Timken LM11749轴承和转子的参数如表1所示。图4给出了转子的有限元模型，由15个单元组成。两个轴承分别安装在节点3和节点14上。因此节点3和节点14的平动自由度是固定的。由偏心质量产生的离心力施加在节点9上。传感器测量节点8和节点11的振动信号。

在仿真中，外圈间隙和内圈间隙均被设为7.5μm。图5给出了当转速为1800r/min时节点8仿真的振动响应和实验所测值的比较。在图5中，F-M表示考虑了转子弹性变形的耦合模型，而R-M表示刚体运动模型。从图5可以看出F-M的仿真结果和测量到振动响应匹配得比较好。由于R-M只可以计算转子的刚体运动，其仿真的幅值要比测量的振动响应小。

表1Timken LM11749轴承和转子的参数

Claims (5)

1.一种滚动轴承-转子系统动态耦合建模方法，其特征在于，包括以下步骤：

1)收集滚动轴承-转子系统的属性参数；

2)根据Jones拟静力学模型，求解滚动轴承-转子系统中每个部件的速度和位移；

3)根据转子的物理属性参数来建立转子的有限元模型；

4)根据轴承各部件间的相互作用计算作用在转子上的合力和合力矩；

5)将计算得到的作用在转子上的合力代入到滚动轴承动力学模型获得滚动轴承-转子系统中每个部件的刚体运动速度和加速度；

6)使用Newmarkβ方法求解转子的有限元模型，获得转子的弹性振动；

7)将转子的弹性振动和刚体运动进行叠加，获得转子的耦合运动；

8)轴承-转子系统所有部件的速度和加速度求得后，使用变步长四阶龙格-库塔-费尔伯格方法得到下一时刻滚动轴承-转子系统中所有部件的速度和位移，重复进行直到达到预先设定的时刻，得到整个轴承-转子系统的时域响应。

2.根据权利要求1所述的一种滚动轴承-转子系统动态耦合建模方法，其特征在于，步骤4)中，滚球和套圈之间的合力为：

式中：c_br——由于接触区域的润滑引起的阻尼系数/N.s.m^-1；Q_k ^c——垂直于滚球和套圈接触面的接触力Q_k ^c，f_k ^c——平行于滚球和套圈接触面的牵引力；

作用在套圈上的力等于作用在滚球上的力，且两者方向相反；滚球质心的力矩和关于套圈质心的合力矩为

式中：r_cpk和r_prk——某接触点相对于滚球中心和套圈中心的位置向量；下标j——轴承的 第j个滚球；上标rk——第k个套圈坐标系；上标c——接触坐标系；T_c,rk——从接触坐标系到套圈坐标系的变换矩阵；

作用在转子上的合力可以表示为

式中：n——装配在转子上轴承的个数；T_ci——从接触坐标系到惯性坐标系的变换矩阵； ——第k个轴承的滚球作用在内圈上的合力；——转子的重力；

作用在转子质心的合力矩为

式中：——由第k个轴承滚球和内圈之间的合力引起的力矩；T_cr——从接触坐标系到转子定体坐标系的转换矩阵。

3.根据权利要求1所述的一种滚动轴承-转子系统动态耦合建模方法，其特征在于，步骤5)中，在惯性坐标系中描述滚球的平动运动，在笛卡尔坐标系中描述套圈、保持架和转子的平动运动；将轴承外圈和轴承座孔的相互作用建模成N_p个均匀支撑在外圈上的弹簧和阻尼并列出外圈和轴承座的动力学方程，同时考虑外圈可活动时的动力学方程和旋转部件在各自定体坐标系中的欧拉方程，求出平动加速度和角加速度。

4.根据权利要求1所述的一种滚动轴承-转子系统动态耦合建模方法，其特征在于，步骤7)中，转子节点处的速度v_p和加速度a_p由下式获得

式中：v_ri和a_ri——刚体运动的速度/m.s^-1和加速度/m.s^-2，通过动力学模型计算；v_el和a_el——弹性振动的速度/m.s^-1和加速度/m.s^-2，通过有限元模型求得；在每一个计算步长，轴承滚球和内圈之间的合力矩M_rk、以及惯性力F_el和惯性力矩M_el通过动力学模型求解，这些力 施加在转子有限元模型上；获得转子的弹性振动后，轴承配合处节点的速度和加速度通过式(27)计算。

5.根据权利要求1所述的一种滚动轴承-转子系统动态耦合建模方法，其特征在于，步骤1)中所述属性参数包括：几何参数、材料属性、运行状态以及润滑模型。