一种用于转子系统动力学建模的刚体单元法.pdf

本发明公开一种用于转子系统动力学建模的刚体单元法，包括：1)收集转子的几何参数、材料属性和运行状态；2)将转子划分为若干个离散的刚体单元，任意相邻的两个刚体单元之间通过四个假想的弹簧相连接；3)根据转子相邻刚体单元之间相互作用模型计算作用在每个刚体单元上的合力和合力矩，并建立每个刚体单元的动力学微分方程；4)使用变步长的四阶龙格‑库塔‑费尔伯格方法求解各个刚体单元的动力学方程来获得刚体单元的振动响应。利用本发明建立的转子系统动力学模型，不需要任何假设，就能够与Gupta轴承动力学模型实现耦合，建立整个滚动轴承‑转子系统的完全动力学模型。

1.一种用于转子系统动力学建模的刚体单元法，其特征在于，包括以下步骤：
1)收集转子的几何参数、材料属性和运行状态；
2)将转子划分为若干个离散的刚体单元，任意相邻的两个刚体单元之间通过四个假想
的弹簧相连接，所述四个假想的弹簧包括一个约束平移运动的拉伸弹簧和三个约束旋转运
动的扭转弹簧；
3)根据转子各个相邻刚体单元之间相互作用模型，计算作用在每个转子刚体单元上的
合力和合力矩，并建立每个刚体单元的动力学微分方程；
4)使用变步长的四阶龙格-库塔-费尔伯格方法求解各个刚体单元的动力学微分方程
来获得每个刚体单元的振动响应。
2.根据权利要求1所述的一种用于转子系统动力学建模的刚体单元法，其特征在于，所
述步骤3)中，定义惯性坐标系i(xi,yi,zi)和第j个刚体单元的定体坐标系rj(xrj,yrj,zrj)，
每个刚体单元有两个端面，左端被命名为α平面，右端被命名为β平面，Mj和Nj分别是第j个刚
体单元α平面和β平面的中心，同样，Mj+1和Nj+1分别是第j+1个刚体单元α平面和β平面的中
心，第j个刚体单元和第j+1个刚体单元之间的几何作用关系为：当转子在自由状态时，第j
个刚体单元的β平面与第j+1个刚体单元的α平面重合，此时它们之间的弹簧处于松弛状态；
当转子受到外力时，由于相邻刚体单元运动的不一致，连接相邻刚体单元平面之间的弹簧
将发生压缩或者拉伸，从而产生了相互作用力；使用上标i表示惯性坐标系，在惯性坐标系
中，Nj和Mj+1分别由向量和确定，由拉伸弹簧产生的作用在第j个刚体单元上的相互作
用力和作用在第j+1个刚体单元上的相互作用力表示
为：
$\left\{\begin{array}{c}{F}_{{\mathrm{\β}}_j}^i=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^s\left(r_{M_{j+1}}^i-r_{N_j}^i\right)=\frac{E\·\min \left(A_j,A_{j+1}\right)}{\left(\frac{l_j+l_{j+1}}{2}\right)}\left(r_{M_{j+1}}^i-r_{N_j}^i\right)\\ F_{{\mathrm{\α}}_{j+1}}^i=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^s\left(r_{N_j}^i-r_{M_{j+1}}^i\right)=\frac{E\·\min \left(A_j,A_{j+1}\right)}{\left(\frac{l_j+l_{j+1}}{2}\right)}\left(r_{N_j}^i-r_{M_{j+1}}^i\right)\end{array}\right.---\left(1\right)$
式中：——第j个刚体单元β面与第j+1个刚体单元α面之间的拉伸弹簧的刚度/
N.m-1；E——转子的弹性模量/Pa；lj和lj+1——第j个和第j+1个刚体单元的长度/m；Aj和
Aj+1——第j个和第j+1个刚体单元端面的面积/m2；函数min(a,b)计算a和b的最小值；
除了拉伸弹簧的相互作用力，同时还有由三个方向上的扭转弹簧产生的作用在第j个
刚体单元的作用力矩和作用在第j+1个刚体单元的作用力矩
相邻刚体单元的姿态角为(ηj,ξj,λj)和(ηj+1,ξj+1,λj+1)，相互作用
力矩和由下式确定：
$\left\{\begin{array}{c}{M}_j^{\mathrm{\β}x}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rx}\frac{{\mathrm{\η}}_{j+1}-{\mathrm{\η}}_j}{l_j+l_{j+1}}=2G\·min(I_{pj},I_{pj+1})\frac{{\mathrm{\η}}_{j+1}-{\mathrm{\η}}_j}{l_j+l_{j+1}}\\ M_{j+1}^{\mathrm{\α}x}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rx}\frac{{\mathrm{\η}}_j-{\mathrm{\η}}_{j+1}}{l_{j+1}+l_j}=2G\·min(I_{pj},I_{pj+1})\frac{{\mathrm{\η}}_j-{\mathrm{\η}}_{j+1}}{l_{j+1}+l_j}\end{array}\right.---\left(2\right)$
式中：——第j个刚体单元β面与第j+1个刚体单元α面之间x方向扭转弹簧的刚度/
N.m.rad-1，G——转子的剪切模量/Pa；Ipj和Ipj+1——第j个和第j+1个刚体单元的极惯性矩/
m4；
力矩分量和通过下式确定：
$\left\{\begin{array}{c}{M}_j^{\mathrm{\β}y}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{ry}\frac{{\mathrm{\ξ}}_j-{\mathrm{\ξ}}_{j+1}}{l_j+l_{j+1}}=2E\·\min \left(I_{yj},I_{yj+1}\right)\frac{{\mathrm{\ξ}}_j-{\mathrm{\ξ}}_{j+1}}{l_j+l_{j+1}}\\ M_{j+1}^{\mathrm{\α}y}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{ry}\frac{{\mathrm{\ξ}}_{j+1}-{\mathrm{\ξ}}_j}{l_j+l_{j+1}}=2E\·\min \left(I_{yj},I_{yj+1}\right)\frac{{\mathrm{\ξ}}_{j+1}-{\mathrm{\ξ}}_j}{l_j+l_{j+1}}\end{array}\right.---\left(3\right)$
式中：——第j个刚体单元β面与第j+1个刚体单元α面之间y方向扭转弹簧的刚度/
N.m.rad-1；Iyj和Iyj+1——第j个刚体单元和第j+1个刚体单元相对于y轴的惯性矩/m4；
力矩分量和表示为：
$\left\{\begin{array}{c}{M}_j^{\mathrm{\β}z}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rz}\frac{{\mathrm{\λ}}_j-{\mathrm{\λ}}_{j+1}}{l_j+l_{j+1}}=2E\·\min \left(I_{zj},I_{zj+1}\right)\frac{{\mathrm{\λ}}_j-{\mathrm{\λ}}_{j+1}}{l_j+l_{j+1}}\\ M_{j+1}^{\mathrm{\α}z}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rz}\frac{{\mathrm{\λ}}_{j+1}-{\mathrm{\λ}}_j}{l_j+l_{j+1}}=2E\·\min \left(I_{zj},I_{zj+1}\right)\frac{{\mathrm{\λ}}_{j+1}-{\mathrm{\λ}}_j}{l_j+l_{j+1}}\end{array}\right.---\left(4\right)$
式中：——第j个刚体单元β面与第j+1个刚体单元α面之间z方向扭转弹簧的刚度/
N.m.rad-1；Izj和Izj+1——第j个刚体单元和第j+1个刚体单元相对于z轴的惯性矩/m4。
3.根据权利要求1所述的一种用于转子系统动力学建模的刚体单元法，其特征在于，步
骤3)中，将每个刚体单元的运动分解成为质心的平动运动和绕其质心的转动；
和分别是由不平衡质量产生的不平衡力和施加在第j个刚
体单元上的外力，则第j个刚体单元的平动运动方程为：
$\left\{\begin{array}{c}{m}_j{\stackrel{\·\·}{x}}_j=F_{{\mathrm{\α}}_j}^x+F_{{\mathrm{\β}}_j}^x+F_{cj}^x+F_{ej}^x\\ m_j{\stackrel{\·\·}{y}}_j=F_{{\mathrm{\α}}_j}^y+F_{{\mathrm{\β}}_j}^y+F_{cj}^y+F_{ej}^y\\ m_j{\stackrel{\·\·}{z}}_j=F_{{\mathrm{\α}}_j}^z+F_{{\mathrm{\β}}_j}^z+F_{cj}^z+F_{ej}^z+G_j\end{array}\right.---\left(5\right)$
式中：mj——第j个刚体单元的质量/kg；(xj,yj,zj)——质心Orj在惯性坐标系中的坐标
分量；Gj——刚体单元的重力/N；不平衡力表示为：
$F_{cj}^i=T_{rj,i}\left\{\begin{array}{c}0\\ m_{uj}{r}_{uj}{\stackrel{\·}{\mathrm{\η}}}_j^2{\mathrm{cos\α}}_j\\ m_{uj}{r}_{uj}{\stackrel{\·}{\mathrm{\η}}}_j^2{\mathrm{sin\α}}_j\end{array}\right\}---\left(6\right)$
式中：mujruj——不平衡质量和其半径的乘积/kg.m；αj——不平衡质量的初始方位角/
rad；Trj,i——从刚体单元定体坐标系到惯性坐标系的变换矩阵；
第j个刚体单元的角速度为ωj(ωjx,ωjy,ωjz)，则刚体单元的旋转运动方程表示为：
$\left\{\begin{array}{c}{I}_{jx}{\stackrel{\·}{\mathrm{\ω}}}_{jx}-\left(I_{jy}-I_{jz}\right){\mathrm{\ω}}_{jy}{\mathrm{\ω}}_{jz}=M_j^{ax}+M_j^{\mathrm{\β}x}+M_{Fj}^{\mathrm{\α}x}+M_{Fj}^{\mathrm{\β}x}+M_{ej}^x\\ I_{jy}{\stackrel{\·}{\mathrm{\ω}}}_{jy}-\left(I_{jz}-I_{jx}\right){\mathrm{\ω}}_{jz}{\mathrm{\ω}}_{jx}=M_j^{ay}+M_j^{\mathrm{\β}y}+M_{Fj}^{\mathrm{\α}y}+M_{Fj}^{\mathrm{\β}y}+M_{ej}^y\\ I_{jz}{\stackrel{\·}{\mathrm{\ω}}}_{jz}-\left(I_{jx}-I_{jy}\right){\mathrm{\ω}}_{jx}{\mathrm{\ω}}_{jy}=M_j^{az}+M_j^{\mathrm{\β}z}+M_{Fj}^{\mathrm{\α}z}+M_{Fj}^{\mathrm{\β}z}+M_{ej}^z\end{array}\right.---\left(7\right)$
式中施加在第j个刚体单元上的外部力矩和
是由和产生的力矩，表示为：

式中：Ti,rj——从惯性坐标系到刚体单元定体坐标系的变换矩阵。
4.根据权利要求1所述的一种用于转子系统动力学建模的刚体单元法，其特征在于，步
骤4)中，使用变步长四阶龙格-库塔-费尔伯格方法进行积分，得到下一时刻各刚体单元的
速度和位移；重复此步骤求得每个刚体单元的振动响应，对应于转子不同位置的振动响应。
5.根据权利要求1所述的一种用于转子系统动力学建模的刚体单元法，其特征在于，步
骤2)中，采用如下标准将转子划分为若干个离散的刚体单元：
(1)同一个刚体单元不能具有不同的内径，同时也不能具有不同的外径；
(2)划分后的每一个刚体单元的长径比不超过5；
(3)每个轴承在转子上的支承位置的中心必须与它所对应的刚体单元的节点对齐，而
不能与两个相邻刚体单元之间的面对齐。

一种用于转子系统动力学建模的刚体单元法

技术领域

本发明涉及一种用于转子系统动力学建模的方法。

背景技术

在转子动力学研究中，合理地对转子系统进行动力学建模是进行后续分析和计算
的基础。从现有文献检索发现，转子动力学建模方法主要有理论法、传递矩阵法和有限元
法，其中基于理论力学的建模方法一般只能用于简单的转子系统建模，对于复杂的转子系
统，较常用的是传递矩阵法和有限元法。传递矩阵法的主要特点是矩阵的阶数不随系统的
自由度数增大而增加，因此编程简单，所需内存小，运算速度快，但是该方法在考虑转子周
围结构时分析计算比较困难，而且传递矩阵法对转子系统进行了简化和等效，无法保证模
型的完整性和分析的准确度。

有限元法基于转子的几何结构，将模型离散为一系列单元进行分析，其表达式简
单规范，能够同时保证模型的完整性和计算准确度，也很容易借助有限元软件实现，很好地
兼顾到模型的完整性和计算效率。另外，利用有限元分析软件可以很好地解决转子动力特
性分析中的“陀螺效应”问题，可以提高准确度和分析效率。申请号为201510816868.3的中
国发明专利申请公开了一种基于双转子简化动力学模型设计的航空发动机仿真试验台，其
特点在于首先采用有限元法建立多轮盘的复杂离散动力学模型，然后基于质心集中方法对
该模型进行了简化。可以看出，利用有限元方法建立的模型结构复杂，适合用于线性振动特
性的定量计算，但对于非线性问题的定量分析计算效率较低，定性分析难以实现；且存在自
由度数特别大，耗费计算机时等问题。

此外，在研究轴承-转子系统的动力学行为时，须建立转子和滚动轴承的耦合模
型。为了考虑转子的弹性振动和柔性，目前较为理想的方法是采用转子有限元模型与Gupta
滚动轴承动力学模型实现耦合，但两者进行耦合时，为了求解有限元模型，需要采用一定的
假设，而这些假设的合理性还有待于进一步验证，因此采用转子有限元模型与Gupta轴承模
型耦合的方法具有一定的局限性。

发明内容

本发明的目的在于提供一种用于转子系统动力学建模的刚体单元法：通过将转子
分隔成由弹簧相连接的离散的刚体单元，根据转子相邻刚体单元之间相互作用模型计算作
用在每个刚体单元上的合力和合力矩，并建立每个刚体单元的动力学微分方程，然后使用
变步长的四阶龙格-库塔-费尔伯格方法求解各个刚体单元的动力学方程来获得刚体单元
的振动响应。本发明提出了一种全新的转子系统动力学建模方法，与传统方法相比较，其显
著优点是：考虑了转子的弹性振动和柔性，操作方法稳定可靠、效率较高；并且利用该方法
建立的转子系统动力学模型，不需要任何假设，就能够与Gupta轴承动力学模型实现耦合，
建立整个滚动轴承-转子系统的完全动力学模型。

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

一种用于转子系统动力学建模的刚体单元法，包括以下步骤：

1)收集转子的几何参数、材料属性和运行状态；

2)将转子划分为若干个离散的刚体单元，任意相邻的两个刚体单元之间通过四个
假想的弹簧相连接，所述四个假想的弹簧包括一个约束平移运动的拉伸弹簧和三个约束旋
转运动的扭转弹簧；

3)根据转子各个相邻刚体单元之间相互作用模型，计算作用在每个转子刚体单元
上的合力和合力矩，并建立每个刚体单元的动力学微分方程；

4)使用变步长的四阶龙格-库塔-费尔伯格方法求解各个刚体单元的动力学微分
方程来获得每个刚体单元的振动响应。

所述步骤2)中，在将转子划分为若干个离散的刚体单元时，为了使建立的转子动
力学模型足够合理，并方便计算作用在刚体单元上的合力和合力矩，需要遵循以下标准：

(1)划分的每一个刚体单元的内径或外径必须一致，同一个刚体单元不能具有不
同的内径或外径；

(2)为了保证各个节点振幅相对误差在合理范围内，建议划分后的每一个刚体单
元的长径比不超过5；

(3)每个轴承在转子上的支承位置的中心必须与它所对应的刚体单元的节点对
齐，而不能与两个相邻刚体单元之间的面对齐。

进一步的，所述步骤3)中，定义惯性坐标系i(xi,yi,zi)和第j个刚体单元的定体坐
标系rj(xrj,yrj,zrj)，可通过刚体单元的姿态角确定。每个刚体单元有两个端面，左端被命
名为α平面，右端被命名为β平面。Mj和Nj分别是第j个刚体单元α平面和β平面的中心。同样，
Mj+1和Nj+1分别是第j+1个刚体单元α平面和β平面的中心。第j个刚体单元和第j+1个刚体单
元之间的几何作用关系为：当转子在自由状态时，第j个刚体单元的β平面与第j+1个刚体单
元的α平面重合，此时它们之间的弹簧处于松弛状态；当转子受到外力时，由于相邻刚体单
元运动的不一致，连接相邻刚体单元平面之间的弹簧将发生压缩或者拉伸，从而产生了相
互作用力。以第j个刚体单元的β平面和第j+1个刚体单元的α平面的相互作用为例，使用上
标i表示惯性坐标系。在惯性坐标系中，Nj和Mj+1分别由向量和确定。由拉伸弹簧产生
的作用在第j个刚体单元上的相互作用力和作用在第j+1个刚体单元上的
相互作用力表示为：

$\begin{array}{c}{F}_{{\mathrm{\β}}_j}^i=k_{{\mathrm{\β}}_j,a_{j+1}}^s\left(r_{M_{j+1}}^i-r_{N_j}^i\right)=\frac{E\·\min \left(A_j,A_{j+1}\right)}{\left(\frac{l_j+l_{j+1}}{2}\right)}\left(r_{M_{j+1}}^i-r_{N_j}^i\right)\\ F_{{\mathrm{\α}}_{j+1}}^i=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^s\left(r_{N_j}^i-r_{M_{j+1}}^i\right)=\frac{E\·\min \left(A_j,A_{j+1}\right)}{\left(\frac{l_j+l_{j+1}}{2}\right)}\left(r_{N_j}^i-r_{M_{j+1}}^i\right)\end{array}---\left(1\right)$

式中：——第j个刚体单元β面与第j+1个刚体单元α面之间的拉伸弹簧的刚
度/N.m-1；E——转子的弹性模量/Pa；lj和lj+1——第j个和第j+1个刚体单元的长度/m；Aj和
Aj+1——第j个和第j+1个刚体单元端面的面积/m2；函数min(a,b)计算a和b的最小值。

除了拉伸弹簧的相互作用力，同时还有由三个方向上的扭转弹簧产生的作用在第
j个刚体单元的作用力矩和作用在第j+1个刚体单元的作用力矩
相邻刚体单元的姿态角为(ηj,ξj,λj)和(ηj+1,ξj+1,λj+1)，相互作用
力矩和由下式确定：

$\left\{\begin{array}{c}{M}_j^{\mathrm{\β}x}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rx}\frac{{\mathrm{\η}}_{j+1}-{\mathrm{\η}}_j}{l_j+l_{j+1}}=2G\·min(I_{pj},I_{pj+1})\frac{{\mathrm{\η}}_{j+1}-{\mathrm{\η}}_j}{l_j+l_{j+1}}\\ M_{j+1}^{\mathrm{\α}x}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rx}\frac{{\mathrm{\η}}_j-{\mathrm{\η}}_{j+1}}{l_{j+1}+l_j}=2G\·min(I_{pj},I_{pj+1})\frac{{\mathrm{\η}}_j-{\mathrm{\η}}_{j+1}}{l_{j+1}+l_j}\end{array}\right.---\left(2\right)$

式中：——第j个刚体单元β面与第j+1个刚体单元α面之间x方向扭转弹簧的
刚度/N.m.rad-1，G——转子的剪切模量/Pa；Ipj和Ipj+1——第j个和第j+1个刚体单元的极惯
性矩/m4；

力矩分量和通过下式确定：

$\left\{\begin{array}{c}{M}_j^{\mathrm{\β}y}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{ry}\frac{{\mathrm{\ξ}}_j-{\mathrm{\ξ}}_{j+1}}{l_j+l_{j+1}}=2E\·\min \left(I_{yj},I_{yj+1}\right)\frac{{\mathrm{\ξ}}_j-{\mathrm{\ξ}}_{j+1}}{l_j+l_{j+1}}\\ M_{j+1}^{\mathrm{\α}y}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{ry}\frac{{\mathrm{\ξ}}_{j+1}-{\mathrm{\ξ}}_j}{l_j+l_{j+1}}=2E\·\min \left(I_{yj},I_{yj+1}\right)\frac{{\mathrm{\ξ}}_{j+1}-{\mathrm{\ξ}}_j}{l_j+l_{j+1}}\end{array}\right.---\left(3\right)$

式中：——第j个刚体单元β面与第j+1个刚体单元α面之间y方向扭转弹簧的
刚度/N.m.rad-1；Iyj和Iyj+1——第j个刚体单元和第j+1个刚体单元相对于y轴的惯性矩/m4；

力矩分量和表示为：

$\left\{\begin{array}{c}{M}_j^{\mathrm{\β}z}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rz}\frac{{\mathrm{\λ}}_j-{\mathrm{\λ}}_{j+1}}{l_j+l_{j+1}}=2E\·min(I_{zj},I_{zj+1})\frac{{\mathrm{\λ}}_j-{\mathrm{\λ}}_{j+1}}{l_j+l_{j+1}}\\ M_{j+1}^{\mathrm{\α}z}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rz}\frac{{\mathrm{\λ}}_{j+1}-{\mathrm{\λ}}_j}{l_j+l_{j+1}}=2E\·min(I_{zj},I_{zj+1})\frac{{\mathrm{\λ}}_{j+1}-{\mathrm{\λ}}_j}{l_j+l_{j+1}}\end{array}\right.---\left(4\right)$

式中：——第j个刚体单元β面与第j+1个刚体单元α面之间z方向扭转弹簧的
刚度/N.m.rad-1；Izj和Izj+1——第j个刚体单元和第j+1个刚体单元相对于z轴的惯性矩/
m4。

进一步的，步骤3)中，将每个刚体单元的运动分解成为质心的平动运动和绕其质
心的转动；和分别是由不平衡质量产生的不平衡力和施加
在第j个刚体单元上的外力，则第j个刚体单元的平动运动方程为：

$\left\{\begin{array}{c}{m}_j{\stackrel{\·\·}{x}}_j=F_{{\mathrm{\α}}_j}^x+F_{{\mathrm{\β}}_j}^x+F_{cj}^x+F_{ej}^x\\ m_j{\stackrel{\·\·}{y}}_j=F_{{\mathrm{\α}}_j}^y+F_{{\mathrm{\β}}_j}^y+F_{cj}^y+F_{ej}^y\\ m_j{\stackrel{\·\·}{z}}_j=F_{{\mathrm{\α}}_j}^z+F_{{\mathrm{\β}}_j}^z+F_{cj}^z+F_{ej}^z+G_j\end{array}\right.---\left(5\right)$

式中：mj——第j个刚体单元的质量/kg；(xj,yj,zj)——质心Orj在惯性坐标系中的
坐标分量；Gj——刚体单元的重力/N；不平衡力表示为：

$F_{cj}^i=T_{rj,i}\left\{\begin{array}{c}0\\ m_{uj}{r}_{uj}{\stackrel{\·}{\mathrm{\η}}}_j^2{\mathrm{cos\α}}_j\\ m_{uj}{r}_{uj}{\stackrel{\·}{\mathrm{\η}}}_j^2{\mathrm{sin\α}}_j\end{array}\right\}---\left(6\right)$

式中：mujruj——不平衡质量和其半径的乘积/kg.m；αj——不平衡质量的初始方位
角/rad；Trj,i——从刚体单元定体坐标系到惯性坐标系的变换矩阵；

第j个刚体单元的角速度为ωj(ωjx,ωjy,ωjz)，则刚体单元的旋转运动方程表示
为：

$\left\{\begin{array}{c}{I}_{jx}{\stackrel{\·}{\mathrm{\ω}}}_{jx}-(I_{jy}-I_{jz}){\mathrm{\ω}}_{jy}{\mathrm{\ω}}_{jz}=M_j^{\mathrm{\α}x}+M_j^{\mathrm{\β}x}+M_{Fj}^{\mathrm{\α}x}+M_{Fj}^{\mathrm{\β}x}+M_{ej}^x\\ I_{jy}{\stackrel{\·}{\mathrm{\ω}}}_{jy}-(I_{jz}-I_{jx}){\mathrm{\ω}}_{jz}{\mathrm{\ω}}_{jx}=M_j^{ay}+M_j^{\mathrm{\β}y}+M_{Fj}^{\mathrm{\α}y}+M_{Fj}^{\mathrm{\β}y}+M_{ej}^y\\ I_{jz}{\stackrel{\·}{\mathrm{\ω}}}_{jz}-(I_{jx}-I_{jy}){\mathrm{\ω}}_{jx}{\mathrm{\ω}}_{jy}=M_j^{az}+M_j^{\mathrm{\β}z}+M_{Fj}^{\mathrm{\α}z}+M_{Fj}^{\mathrm{\β}z}+M_{ej}^z\end{array}\right.---\left(7\right)$

式中施加在第j个刚体单元上的外部力矩和
是由和产生的力矩，表示为：

式中：Ti,rj——从惯性坐标系到刚体单元定体坐标系的变换矩阵。

所述步骤4)中，使用变步长四阶龙格-库塔-费尔伯格方法进行积分，可得到下一
时刻各刚体单元的速度和位移。重复此步骤可求得每个刚体单元的振动响应，对应于转子
不同位置的振动响应。

本发明的进一步改进在于，转子模型的求解步骤为：首先根据转子的几何参数、材
料属性、运行状态计算各转子单元初始速度和位移，然后根据转子相邻刚体单元之间相互
作用模型计算作用在每个刚体单元上的合力和合力矩，并建立每个刚体单元的动力学微分
方程；使用变步长的四阶龙格-库塔-费尔伯格方法求解各个刚体单元的动力学方程来获得
刚体单元的振动响应。

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

1、提出了一个新的转子系统动力学建模方法(DM)，将转子分隔成离散的刚体单
元，通过求解各个刚体单元的动力学方程来获得刚体单元的动力学响应；

2、利用本发明的方法建立的转子动力学模型，不需要做任何假设，就能够与Gupta
轴承动力学模型实现耦合，建立整个轴承-转子系统的完全动力学模型。

本发明通过将转子分隔成由弹簧相连接的离散的刚体单元，根据转子相邻刚体单
元之间相互作用模型计算作用在每个刚体单元上的合力和合力矩，并建立每个刚体单元的
动力学微分方程，使用变步长的四阶龙格-库塔-费尔伯格方法求解各个刚体单元的动力学
方程来获得刚体单元的振动响应。利用该方法建立的转子系统动力学模型，不需要任何假
设，就能够与Gupta轴承动力学模型实现耦合，建立整个滚动轴承-转子系统的完全动力学
模型。

附图说明

图1是本发明中转子的任意两个相连接刚体单元之间的相互作用：坐标系分别为
惯性坐标系i(xi,yi,zi)和第j个刚体单元的定体坐标系rj(xrj,yrj,zrj)；

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

图3是本发明的转子模型，其中图3(a)为转子系统动力学模型，图3(b)为转子系统
有限元模型；

图4是本发明中节点4在y方向的时域响应图；

图5是本发明中节点4在z方向的时域响应图；

图6是本发明中转子各节点的位移图；

图7是本发明中转子转速为10000r/min时不同节点的轨迹。

具体实施方式

请参阅图1及图2所示，本发明一种用于转子系统动力学建模的刚体单元法，包括
以下步骤：

1)收集转子的几何参数、材料属性和运行状态，为转子系统动力学建模提供数据
支持；

2)将转子划分为若干个离散的刚体单元，任意相邻的两个刚体单元之间通过四个
假想的弹簧相连接，即一个约束平移运动的拉伸弹簧和三个约束旋转运动的扭转弹簧；

3)根据转子各个相邻刚体单元之间相互作用模型，计算作用在每个转子刚体单元
上的合力和合力矩，并建立每个刚体单元的动力学微分方程；

4)使用变步长的四阶龙格-库塔-费尔伯格方法求解各个刚体单元的动力学方程
来获得每个刚体单元的振动响应。

步骤1)中，转子的材料特性参数包括转子长度、直径，以及弹性模量、泊松比和密
度。

步骤2)中，在将转子划分为若干个离散的刚体单元时，为了使建立的转子动力学
模型足够合理，并方便计算作用在刚体单元上的合力和合力矩，需要遵循以下标准：

(1)划分的每一个刚体单元的内径或外径必须一致，同一个刚体单元不能具有不
同的内径或外径；

(2)为了保证各个节点振幅相对误差在合理范围内，建议划分后的每一个刚体单
元的长径比不超过5；

(3)每个轴承在转子上的支承位置的中心必须与它所对应的刚体单元的节点对
齐，而不能与两个相邻刚体单元之间的面对齐。

步骤4)中，使用变步长四阶龙格-库塔-费尔伯格方法进行积分，可得到下一时刻
各刚体单元的速度和位移。重复此步骤可求得每个刚体单元的振动响应。

下面结合一个实例对本发明作进一步详细说明，但本实例并不用于限制本发明。
其步骤为：

1)收集转子的几何参数、材料属性和运行状态，为转子系统动力学建模提供数据
支持，转子的参数如表1所示。

表1转子参数

2)将转子等分成9个刚体单元(其中第2个和第8个刚体单元通过弹簧和阻尼支
承)，任意相邻的两个刚体单元之间通过四个假想的弹簧相连接，即一个约束平移运动的拉
伸弹簧和三个约束旋转运动的扭转弹簧。

3)根据相邻刚体单元之间的相互作用模型，对转子各刚体单元之间的相互作用进
行分析，计算作用在每个转子单元上的力和力矩，建立每个刚体单元的运动微分方程，转子
动力学模型(DM)如图3(a)所示。具体计算方法如下：

选取任意两个相邻刚体单元(第j个和第j+1个)，定义惯性坐标系i(xi,yi,zi)和第
j个刚体单元的定体坐标系rj(xrj,yrj,zrj)，可通过刚体单元的姿态角确定。每个刚体单元
的左端面被命名为α平面，右端面被命名为β平面。Mj和Nj分别是第j个刚体单元α平面和β平
面的中心。同样，Mj+1和Nj+1分别是第j+1个刚体单元α平面和β平面的中心。使用上标i表示惯
性坐标系，在惯性坐标系中，Nj和Mj+1分别由向量和确定。当转子受到外力时，由拉伸
弹簧产生的作用在第j个刚体单元上的相互作用力和作用在第j+1个刚体
单元上的相互作用力可根据下式来计算

$\left\{\begin{array}{c}{F}_{{\mathrm{\β}}_j}^i=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^s\left(r_{M_{j+1}}^i-r_{N_j}^i\right)=\frac{E\·\min \left(A_j,A_{j+1}\right)}{\left(\frac{l_j+l_{j+1}}{2}\right)}\left(r_{M_{j+1}}^i-r_{N_j}^i\right)\\ F_{{\mathrm{\α}}_{j+1}}^i=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^s\left(r_{N_j}^i-r_{M_{j+1}}^i\right)=\frac{E\·\min \left(A_j,A_{j+1}\right)}{\left(\frac{l_j+l_{j+1}}{2}\right)}\left(r_{N_j}^i-r_{M_{j+1}}^i\right)\end{array}\right.---\left(1\right)$

式中：——第j个刚体单元β面与第j+1个刚体单元α面之间的拉伸弹簧的刚
度/N.m-1；E——转子的弹性模量/Pa；lj和lj+1——第j个和第j+1个刚体单元的长度/m；Aj和
Aj+1——第j个和第j+1个刚体单元端面的面积/m2。函数min(a,b)计算a和b的最小值。除了
拉伸弹簧的相互作用力，同时还有由三个方向上的扭转弹簧产生的作用在第j个刚体
单元的作用力矩和作用在第j+1个刚体单元的作用力矩
相邻刚体单元的姿态角为(ηj,ξj,λj)和(ηj+1,ξj+1,λj+1)，则相互作
用力矩和可由下式确定

$\left\{\begin{array}{c}{M}_j^{\mathrm{\β}x}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rx}\frac{{\mathrm{\η}}_{j+1}-{\mathrm{\η}}_j}{l_j+l_{j+1}}=2G\·min(I_{pj},I_{pj+1})\frac{{\mathrm{\η}}_{j+1}-{\mathrm{\η}}_j}{l_j+l_{j+1}}\\ M_{j+1}^{\mathrm{\α}x}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rx}\frac{{\mathrm{\η}}_j-{\mathrm{\η}}_{j+1}}{l_{j+1}+l_j}=2G\·min(I_{pj},I_{pj+1})\frac{{\mathrm{\η}}_j-{\mathrm{\η}}_{j+1}}{l_{j+1}+l_j}\end{array}\right.---\left(2\right)$

式中：——第j个刚体单元β面与第j+1个刚体单元α面之间x方向扭转弹簧的
刚度/N.m.rad-1，G——转子的剪切模量/Pa；Ipj和Ipj+1——第j个和第j+1个刚体单元的极惯
性矩/m4。

力矩分量和可以通过下式确定

$\left\{\begin{array}{c}{M}_j^{\mathrm{\β}y}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{ry}\frac{{\mathrm{\ξ}}_j-{\mathrm{\ξ}}_{j+1}}{l_j+l_{j+1}}=2E\·\min \left(I_{yj},I_{yj+1}\right)\frac{{\mathrm{\ξ}}_j-{\mathrm{\ξ}}_{j+1}}{l_j+l_{j+1}}\\ M_{j+1}^{\mathrm{\α}y}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{ry}\frac{{\mathrm{\ξ}}_{j+1}-{\mathrm{\ξ}}_j}{l_j+l_{j+1}}=2E\·\min \left(I_{yj},I_{yj+1}\right)\frac{{\mathrm{\ξ}}_{j+1}-{\mathrm{\ξ}}_j}{l_j+l_{j+1}}\end{array}\right.---\left(3\right)$

式中：——第j个刚体单元β面与第j+1个刚体单元α面之间y方向扭转弹簧的
刚度/N.m.rad-1；Iyj和Iyj+1——第j个刚体单元和第j+1个刚体单元相对于y轴的惯性矩/m4。

力矩分量和可表示为

$\left\{\begin{array}{c}{M}_j^{\mathrm{\β}z}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rz}\frac{{\mathrm{\λ}}_j-{\mathrm{\λ}}_{j+1}}{l_j+l_{j+1}}=2E\·min(I_{zj},I_{zj+1})\frac{{\mathrm{\λ}}_j-{\mathrm{\λ}}_{j+1}}{l_j+l_{j+1}}\\ M_{j+1}^{\mathrm{\α}z}=k_{{\mathrm{\β}}_j,{\mathrm{\α}}_{j+1}}^{rz}\frac{{\mathrm{\λ}}_{j+1}-{\mathrm{\λ}}_j}{l_j+l_{j+1}}=2E\·min(I_{zj},I_{zj+1})\frac{{\mathrm{\λ}}_{j+1}-{\mathrm{\λ}}_j}{l_j+l_{j+1}}\end{array}\right.---\left(4\right)$

式中：——第j个刚体单元β面与第j+1个刚体单元α面之间z方向扭转弹簧的
刚度/N.m.rad-1；Izj和Izj+1——第j个刚体单元和第j+1个刚体单元相对于z轴的惯性矩/m4。

每个刚体单元的运动可以分解成为质心的平动运动和绕其质心的转动。假设
和分别是由不平衡质量产生的不平衡力和施加在第j个刚体
单元上的外力，则第j个刚体单元的平动运动方程为

$\left\{\begin{array}{c}{m}_j{\stackrel{\·\·}{x}}_j=F_{{\mathrm{\α}}_j}^x+F_{{\mathrm{\β}}_j}^x+F_{cj}^x+F_{ej}^x\\ m_j{\stackrel{\·\·}{y}}_j=F_{{\mathrm{\α}}_j}^y+F_{{\mathrm{\β}}_j}^y+F_{cj}^y+F_{ej}^y\\ m_j{\stackrel{\·\·}{z}}_j=F_{{\mathrm{\α}}_j}^z+F_{{\mathrm{\β}}_j}^z+F_{cj}^z+F_{ej}^z+G_j\end{array}\right.---\left(5\right)$

式中：mj——第j个刚体单元的质量/kg；(xj,yj,zj)——质心Orj在惯性坐标系中的
坐标分量；Gj——刚体单元的重力/N。不平衡力可以表示为

$F_{cj}^i=T_{rj,i}\left\{\begin{array}{c}0\\ m_{uj}{r}_{uj}{\stackrel{\·}{\mathrm{\η}}}_j^2{\mathrm{cos\α}}_j\\ m_{uj}{r}_{uj}{\stackrel{\·}{\mathrm{\η}}}_j^2{\mathrm{sin\α}}_j\end{array}\right\}---\left(6\right)$

式中：mujruj——不平衡质量和其半径的乘积/kg.m；αj——不平衡质量的初始方位
角/rad；Trj,i——从刚体单元定体坐标系到惯性坐标系的变换矩阵。

假设第j个刚体单元的角速度为ωj(ωjx,ωjy,ωjz)，则刚体单元的旋转运动方程
可以表示为

$\left\{\begin{array}{c}{I}_{jx}{\stackrel{\·}{\mathrm{\ω}}}_{jx}-(I_{jy}-I_{jz}){\mathrm{\ω}}_{jy}{\mathrm{\ω}}_{jz}=M_j^{\mathrm{\α}x}+M_j^{\mathrm{\β}x}+M_{Fj}^{\mathrm{\α}x}+M_{Fj}^{\mathrm{\β}x}+M_{ej}^x\\ I_{jy}{\stackrel{\·}{\mathrm{\ω}}}_{jy}-(I_{jz}-I_{jx}){\mathrm{\ω}}_{jz}{\mathrm{\ω}}_{jx}=M_j^{ay}+M_j^{\mathrm{\β}y}+M_{Fj}^{\mathrm{\α}y}+M_{Fj}^{\mathrm{\β}y}+M_{ej}^y\\ I_{jz}{\stackrel{\·}{\mathrm{\ω}}}_{jz}-(I_{jx}-I_{jy}){\mathrm{\ω}}_{jx}{\mathrm{\ω}}_{jy}=M_j^{az}+M_j^{\mathrm{\β}z}+M_{Fj}^{\mathrm{\α}z}+M_{Fj}^{\mathrm{\β}z}+M_{ej}^z\end{array}\right.---\left(7\right)$

式中施加在第j个刚体单元上的外部力矩和
是由和产生的力矩，可以表示为

式中：Ti,rj——从惯性坐标系到刚体单元定体坐标系的变换矩阵。

4)使用变步长四阶龙格-库塔-费尔伯格方法求解各个刚体单元的动力学方程，最
大允许截断误差被设定为1.0E-5，可以得到任意一个或所关注刚体单元的振动响应。

本实例取转子动力学模型的节点4为研究对象，当转子转速为1000r/min时，在节
点4的z方向施加一个大小为100N的正弦激振力，通过仿真分析，可以得到它的y方向和z方
向的时域响应，分别如图4和图5所示，可以看出节点4在y方向和z方向的时域响应曲线均为
正弦曲线。重复以上方法，就可以求解任意一个或所关注刚体单元的振动响应。

另外，图3(b)给出了包含10个单元的转子的有限元模型(FEM)，第1个和第9个单元
的长度均为20mm，其他单元的长度都为40mm。转子的参数如表1示所列。有限元模型使用
Newmarkβ方法求解，计算时间步长设定为2.0E-4s。

1)假设当转子转速为0时，在节点4的z方向施加一个径向力，图6为不同大小的径
向力作用下转子各节点在z方向的位移。从图6可以看出当转子受到静载荷时，FEM的计算结
果和DM的计算结果匹配较好。

2)假设在节点4处存在不平衡质量，不平衡质量和其半径的乘积是0.001kg.m。当
转子旋转时，转子在不平衡力作用下会产生振动响应。图7是当转速为10000r/min时各个节
点的轨迹(图中：FEM——有限元模型计算的结果，DM——动力学模型计算的结果)，由于支
承弹簧的限制，节点2和节点8的振幅要小于其他节点。