一种滑动轴承转子系统非线性动力学特性快速计算方法

1.本发明属于转子动力学领域，具体涉及一种滑动轴承转子系统非线性动力学特性快速计算方法。


背景技术：

2.滑动轴承转子系统是旋转机械设备非常重要的组成部分，广泛应用于航空航天、工程机械、国防军工等领域。随着科技的不断进步和市场的特殊要求，旋转机械设备的工作环境逐渐趋向于高速、高温、重载，使得滑动轴承转子系统经常工作在爆炸、冲击、摇摆、急加/减速等特殊工况下，此时作用在轴承上的载荷随时间不断变化，每一时刻转子轴心的平衡位置、轴承油膜厚度也都发生变化，轴承会产生非线性油膜力，导致滑动轴承转子系统表现出非常显著的非线性特征。相关研究发现，非线性特征会对旋转机械设备的可靠性、安全性产生很大的影响，所以准确、有效、快速地分析滑动轴承转子系统的非线性动力学行为显得越来越重要。
3.在滑动轴承转子系统中，轴承的非线性油膜力是系统表现出非线性特征的主要因素之一，所以获得系统非线性特征的关键在于轴承非线性油膜力的计算。然而传统的刚度、阻尼计算方法，仅适用于轴承轴心轨迹在小范围内变化的情况，且计算时间过长，计算效率低下。而迁移速率法和阻抗法，仅适用于圆对称滑动轴承，对于椭圆、三油叶、阶梯面等非圆对称滑动轴承，具有一定的局限性。因此为了求解滑动轴承转子系统在特殊工况下的非线性动力学特性，解决传统方法计算效率低下、适用范围局限性的问题，急需一种快速求解滑动轴承转子系统非线性动力学特性的计算方法。


技术实现要素：

4.本发明的目的在于针对滑动轴承转子系统工作在爆炸、冲击、摇摆、急加/减速等特殊工况下，系统非线性油膜力计算时间过长的问题，提供了一种滑动轴承转子系统非线性动力学特性快速计算方法，该方法包含系统质量矩阵m、刚度矩阵k、陀螺矩阵g、重力矢量fg、不平衡力矢量f1和f2、轴心轨迹可行域s、轴承刚度阻尼系数、非线性油膜力矢量f
oil
的计算。已知条件包括转子参数、圆盘参数、滑动轴承参数、材料参数和润滑油参数。
5.本发明采用如下技术方案来实现的：
6.一种滑动轴承转子系统非线性动力学特性快速计算方法，包括以下步骤：
7.1)根据给定的转子参数、圆盘参数、滑动轴承参数、材料参数，建立滑动轴承转子系统euler
‑
bernoulli梁有限元模型和运动微分方程，滑动轴承转子系统euler
‑
bernoulli梁有限元模型的每个节点四个自由度；
8.2)根据步骤1)中有限元模型和运动微分方程，计算总体质量矩阵m、总体刚度矩阵k、总体陀螺矩阵g、总体重力矢量fg、总体不平衡力矢量f1和f2；
9.3)根据给定的滑动轴承参数计算轴承的油膜厚度h(θ)，令h(θ)等于零求解轴承的轴心轨迹可行域s；
10.4)采用delaunay三角剖分算法对可行域s进行网格划分，且靠近边界处，网格细化；
11.5)利用有限差分法求解reynolds方程，计算滑动轴承的油膜压力分布p；
12.6)利用步骤5)中得到的油膜压力分布p，计算油膜承载力f
x
、f
y
；
13.7)利用步骤6)中得到的油膜承载力，采用小扰动法计算可行域s内每个节点处滑动轴承的刚度系数阻尼系数i为可行域s内节点的编号；
14.8)采用二维曲面散乱插值函数griddata分别将步骤7)中计算得到刚度、阻尼系数拟合为关于可行域s内坐标的函数，进而得到轴心轨迹可行域s内任意位置处滑动轴承的刚度系数k
xx
、k
xy
、k
yx
、k
yy
，阻尼系数c
xx
、c
xy
、c
yx
、c
yy
；
15.9)根据偏心参数确定初始时刻轴承的偏心位置(x
t0
,y
t0
)和速度根据给定工况确定迭代次数以及迭代终止时间t
max
，并令t＝t0；
16.10)根据t时刻轴承的偏心位置确定所在的delaunay三角形区域，并利用voronoi图确定delaunay三角形区域内距离偏心位置最近的节点q；
17.11)调用步骤8)中拟合得到的刚度系数k
xx
、k
xy
、k
yx
、k
yy
和阻尼系数c
xx
、c
xy
、c
yx
、c
yy
，步骤10)中t时刻轴承的偏心位置和节点q，求解t时刻轴承的油膜力f
oilx
、f
oily
，并组装得到非线性油膜力矢量f
oil
；
18.12)调用步骤2)中计算得到的m、k、g、fg、f1、f2和步骤11)中计算得到的非线性油膜力矢量f
oil
，并将其代入步骤1)所建立的系统运动微分方程中，求解给定工况下t时刻系统的瞬态动力学响应；
19.13)假设时间间隔为δt，根据初始偏心位置、外部激励载荷、步骤5)中油膜压力分布，计算t＝t+δt时刻轴承的偏心位置(x
t+δt
,y
t+δt
)和速度
20.14)重复步骤10)至13)，并判断t是否大于等于t
max
，若t﹤t
max
，则继续迭代；若t≧t
max
，则迭代结束，得到轴承的轴心轨迹以及各个时刻系统的瞬态动力学响应，进而用于对系统的非线性动力学特性进行分析。
21.本发明进一步的改进在于，步骤1)中，运动微分方程为：
[0022][0023]
式中：m是总体质量矩阵，g是总体陀螺矩阵，k是总体刚度矩阵，fg是总体重力矢量；f1和f2分别为系统沿x、y方向的总体不平衡力矢量，f
oil
是滑动轴承非线性油膜力矢量，q(t)是广义节点位移，ω是转子转速。
[0024]
本发明进一步的改进在于，步骤2)中，总体质量矩阵m、总体刚度矩阵k、总体陀螺矩阵g、总体重力矢量fg、总体不平衡力矢量f1和f2分别由各单元矩阵组装得到，其中单元质量矩阵m
e
、单元刚度矩阵k
e
、单元陀螺矩阵g
e
可根据转子动力学理论直接获得，单元力矢量fg
e
、f
1e
和f
2e
表达式如下所示：
[0025]
单元重力矢量fg
e
：
[0026]
f
ge
＝[0,m
i
g,0,0]
t
ꢀꢀ
(2)
[0027]
式中：m
i
是第i个梁单元处转子轴段的质量；
[0028]
单元不平衡力矢量f
1e
和f
2e
：
[0029][0030]
式中：m
e
d是不平衡质量，是不平衡相位角，ω是转子转速，f
1e
表示x方向单元不平衡力矢量，f
2e
表示y方向单元不平衡力矢量；
[0031]
组装得到的总体力矢量fg、总体不平衡力矢量f1和f2表达式如下所示：
[0032]
总体重力矢量fg：
[0033]
f
g
＝[0,m1g,0,0,
……
0,m
i
g,0,0,
……
,0,m
n
‑1g,0,0]
t
ꢀꢀ
(4)
[0034]
式中：m
i
是第i个梁单元处转子轴段的质量，n是有限元模型中的节点总数；
[0035]
总体不平衡力矢量f1和f2：
[0036][0037]
式中：m
e
d是不平衡质量，是不平衡相位角，ω是转子转速，j是不平衡质量位置处的节点编号，f1表示x方向不平衡力矢量，f2表示y方向不平衡力矢量。
[0038]
本发明进一步的改进在于，步骤3)中，对于固定瓦轴承，其油膜厚度表达式如式(6)所示；对于可倾瓦轴承，其油膜厚度表达式如式(7)所示：
[0039]
h(θ)＝c
p
‑
x
j
cos(θ)
‑
y
j
sin(θ)
‑
c
p
mcos(θ
‑
θ
p
)
ꢀꢀ
(6)
[0040]
h(θ)＝c
p
‑
x
j
cos(θ)
‑
y
j
sin(θ)
‑
c
p
mcos(θ
‑
θ
p
)
‑
(r+t)δsin(θ
‑
θ
p
)
ꢀꢀ
(7)
[0041]
式中：θ为周向位置角，c
p
为滑动轴承的半径间隙，x
j
是轴颈周向位移，y
j
是轴颈轴向位移，r为轴承半径，m为预负荷系数，t为预偏心距。
[0042]
本发明进一步的改进在于，步骤5)中reynolds方程为：
[0043][0044]
式中：g
x
是周向紊流修正因子，g
y
是轴向紊流修正因子，η为润滑油粘度，ρ为润滑油密度，p为油膜压力，x、y分别为周向和轴向坐标；
[0045]
边界条件为：
[0046][0047]
式中：p
s
为供油压力，b为轴承宽度。
[0048]
本发明进一步的改进在于，步骤6)中，x、y方向滑动轴承油膜承载力计算公式为；
[0049]
[0050]
式中：θ为周向位置角，θ1油膜起始角，θ2油膜破裂角。
[0051]
本发明进一步的改进在于，步骤7)中，假设轴承轴心位于可行域内第i个节点，轴颈在x、y方向产生微小扰动，此时轴颈轴心位置从(x，y)移动到(x0，y0)，则)，则之间的表达式为：
[0052][0053]
式中：是轴颈轴心位置位于(x，y)处的油膜承载力，是轴颈轴心位置位于(x0，y0)处的油膜承载力；
[0054]
可行域内第i个节点处的刚度阻尼系数为：
[0055][0056]
式中：分别表示第i个节点处x、y方向直接刚度，表示第i个节点处交叉刚度；分别表示第i个节点处x、y方向直接阻尼，表示第i个节点处交叉阻尼。
[0057]
本发明进一步的改进在于，步骤8)中，采用二维曲面散乱插值函数griddata拟合得到的函数如下所示：
[0058][0059]
式中：(x，y)分别是可行域s内任意一点的坐标，且拟合函数都满足步骤7)中节点处的数据点。
[0060]
本发明进一步的改进在于，步骤11)中，t时刻轴承的油膜力f
oilx
、f
oily
为：
[0061][0062]
式中：(x
t
，y
t
)是t时刻轴承的偏心位置，分别是t时刻偏心位置处x、y方向的速度，k
xx
、k
xy
、k
yx
、k
yy
是(x
t
，y
t
)位置处的刚度系数，c
xx
、c
xy
、c
yx
、c
yy
是(x
t
，y
t
)位置处的阻尼系数，f
qx
、f
qy
是节点q处的油膜承载力，δx
t
、δy
t
分别是偏心位置与节点q在x、y方向的距离；
[0063]
组装得到的非线性油膜力矢量f
oil
表达式为：
[0064][0065]
式中：n为有限元模型中的节点总数，p、q分别为有限元模型中左、右滑动轴承位置处的节点编号。
[0066]
本发明进一步的改进在于，步骤13)中，t时刻和t+δt时刻轴承偏心位置表达式如下所示：
[0067][0068]
式中：(x
t
，y
t
)是t时刻轴承的偏心位置，(x
t+δt
，y
t+δt
)是t+δt时刻轴承的偏心位置，δt是时间间隔，分别是t时刻偏心位置处x、y方向的速度，和可根据压力分布p等于外部激励载荷求得。
[0069]
本发明至少具有如下有益的技术效果：
[0070]
1)本发明提供的一种滑动轴承转子系统非线性动力学特性快速计算方法，解决了传统计算方法求解爆炸、冲击、摇摆、急加/减速等特殊工况下轴承非线性油膜力计算效率低下、适用范围局限性的问题。
[0071]
2)本发明将轴心轨迹可行域内节点处的刚度、阻尼系数拟合为可行域内任意位置处的二维曲面函数，在后续求解转子非线性动力学时，结合上述拟合结果、偏心位置、速度以及voronoi图可快速得到给定工况下t时刻系统的非线性油膜力矢量，避免了重复迭代、调用的繁琐过程，使得非线性动力学计算过程更简单、效率更高。
[0072]
3)本发明提出的计算方法包含非线性油膜力计算和转子非线性动力学计算两大部分，具体包括系统质量矩阵m、刚度矩阵k、陀螺矩阵g、重力矢量fg、不平衡力矢量f1和f2、轴心轨迹可行域s、轴承刚度阻尼系数、非线性油膜力矢量f
oil
的计算，详细的考虑了系统的各项参数指标，为滑动轴承转子系统非线性动力学特性研究提供了一种全面、规范的快速计算方法。
附图说明
[0073]
图1为滑动轴承转子系统非线性动力学特性快速计算方法流程图。
[0074]
图2为滑动轴承转子系统euler
‑
bernoulli梁有限元模型。
[0075]
图3为矩阵组装示意图。
[0076]
图4为轴心轨迹可行域网格划分示意图。
[0077]
图5为轴心轨迹示意图。
具体实施方式
[0078]
下面结合附图对本发明做进一步的详细描述：
[0079]
如图1所示，其为快速计算方法的流程图，本发明提供的一种滑动轴承转子系统非线性动力学特性快速计算方法，包括以下步骤：
[0080]
1)根据给定的转子参数、圆盘参数、材料参数，建立euler
‑
bernoulli梁有限元模
型(每个节点4个自由度)和运动微分方程。euler
‑
bernoulli梁有限元模型如图2所示，运动微分方程如下所示：
[0081][0082]
式中：m是总体质量矩阵，g是总体陀螺矩阵，k是总体刚度矩阵，fg是总体重力矢量；f1和f2分别为系统沿x、y方向的总体不平衡力矢量，f
oil
是滑动轴承非线性油膜力矢量，q(t)是广义节点位移，ω是转子转速；
[0083]
2)根据步骤1)中有限元模型和运动微分方程，计算总体质量矩阵m、总体刚度矩阵k、总体陀螺矩阵g、总体重力矢量fg、总体不平衡力矢量f1和f2。各个总体矩阵分别由各单元矩阵组装得到，其中单元质量矩阵m
e
、单元刚度矩阵k
e
、单元陀螺矩阵g
e
可根据转子动力学理论直接获得，单元力矢量fg
e
、f
1e
和f
2e
表达式如下所示：
[0084]
单元重力矢量fg
e
：
[0085]
f
ge
＝[0,m
i
g,0,0]
t
ꢀꢀ
(2)
[0086]
式中：m
i
是第i个梁单元处转子轴段的质量。
[0087]
单元不平衡力矢量f
1e
和f
2e
：
[0088][0089]
式中：m
e
d是不平衡质量，是不平衡相位角，ω是转子转速，f
1e
表示x方向单元不平衡力矢量，f
2e
表示y方向单元不平衡力矢量。
[0090]
总体质量矩阵m、总体刚度矩阵k、总体陀螺矩阵g组装方式如图3所示，总体重力矢量fg、总体不平衡力矢量f1和f2组装方式如下所示：
[0091]
总体重力矢量fg：
[0092]
f
g
＝[0,m1g,0,0,
……
0,m
i
g,0,0,
……
,0,m
n
‑1g,0,0]
t
ꢀꢀ
(4)
[0093]
式中：m
i
是第i个梁单元处转子轴段的质量，n是有限元模型中的节点总数。
[0094]
总体不平衡力矢量f1和f2：
[0095][0096]
式中：m
e
d是不平衡质量，是不平衡相位角，ω是转子转速，j是不平衡质量位置处的节点编号，f1表示x方向不平衡力矢量，f2表示y方向不平衡力矢量。
[0097]
3)根据给定的滑动轴承参数计算轴承的油膜厚度h(θ)，令h(θ)等于零求解轴承的轴心轨迹可行域s。对于固定瓦轴承，其油膜厚度表达式如式(6)所示；对于可倾瓦轴承，其油膜厚度表达式如式(7)所示；
[0098]
h(θ)＝c
p
‑
x
j
cos(θ)
‑
y
j
sin(θ)
‑
c
p
mcos(θ
‑
θ
p
)
ꢀꢀ
(6)
[0099]
h(θ)＝c
p
‑
x
j
cos(θ)
‑
y
j
sin(θ)
‑
c
p
mcos(θ
‑
θ
p
)
‑
(r+t)δsin(θ
‑
θ
p
)
ꢀꢀ
(7)
[0100]
式中：θ为周向位置角，c
p
为滑动轴承的半径间隙，x
j
是轴颈周向位移，y
j
是轴颈轴
向位移，r为轴承半径，m为预负荷系数，t为预偏心距。
[0101]
4)采用delaunay三角剖分算法对可行域进行网格划分。为了提高后续计算结果的准确性，靠近边界处，网格需细化，如图4所示；
[0102]
5)利用有限差分法求解reynolds方程，计算滑动轴承的油膜压力分布p，reynolds方程表达式如下所示：
[0103][0104]
式中：g
x
是周向紊流修正因子，g
y
是轴向紊流修正因子，η为润滑油粘度，ρ为润滑油密度，p为油膜压力，x、y分别为周向和轴向坐标。
[0105]
边界条件为：
[0106][0107]
式中：p
s
为供油压力，b为轴承宽度。
[0108]
6)利用步骤5)中得到的油膜压力分布p，计算油膜承载力f
x
、f
y
，具体计算公式如下所示：
[0109][0110]
式中：θ为周向位置角，θ1油膜起始角，θ2油膜破裂角
[0111]
7)利用步骤6)中得到的油膜承载力，采用小扰动法计算可行域内每个节点处滑动轴承的刚度系数阻尼系数i为可行域s内节点的编号。假设轴承轴心位于可行域内第i个节点，轴颈在x、y方向产生微小扰动，此时轴颈轴心位置从(x，y)移动到(x0，y0)，则之间的表达式为：
[0112][0113]
式中：是轴颈轴心位置位于(x，y)处的油膜承载力，是轴颈轴心位置位于(x0，y0)处的油膜承载力。
[0114]
则可行域内第i个节点处的刚度阻尼系数为：
[0115]
[0116]
式中：分别表示第i个节点处x、y方向直接刚度，表示第i个节点处交叉刚度；分别表示第i个节点处x、y方向直接阻尼，表示第i个节点处交叉阻尼。
[0117]
8)采用二维曲面散乱插值函数griddata分别将步骤7)中计算得到节点处的刚度、阻尼系数拟合为关于可行域内任意坐标的函数，进而得到轴心轨迹可行域内任意位置处滑动轴承的刚度系数k
xx
、k
xy
、k
yx
、k
yy
，阻尼系数c
xx
、c
xy
、c
yx
、c
yy
。上述拟合得到的函数如下所示：
[0118][0119]
式中：(x，y)分别是可行域s内任意一点的坐标，且拟合函数都满足步骤7)中节点处的数据点。
[0120]
9)根据偏心参数确定初始时刻轴承的偏心位置(x
t0
,y
t0
)和速度根据给定工况确定迭代次数以及迭代终止时间t
max
，并令t＝t0；
[0121]
10)根据t时刻轴承的偏心位置确定所在的delaunay三角形区域，并利用voronoi图确定delaunay三角形区域内距离偏心位置最近的节点q。
[0122]
11)调用步骤8)中拟合得到的刚度系数k
xx
、k
xy
、k
yx
、k
yy
和阻尼系数c
xx
、c
xy
、c
yx
、c
yy
，步骤10)中t时刻轴承的偏心位置和节点q，求解t时刻轴承的油膜力f
oilx
、f
oily
，并组装得到非线性油膜力矢量f
oil
。t时刻轴承的油膜力f
oilx
、f
oily
计算公式为：
[0123][0124]
式中：(x
t
，y
t
)是t时刻轴承的偏心位置，分别是t时刻偏心位置处x、y方向的速度，k
xx
、k
xy
、k
yx
、k
yy
是(x
t
，y
t
)位置处的刚度系数，c
xx
、c
xy
、c
yx
、c
yy
是(x
t
，y
t
)位置处的阻尼系数，f
qx
、f
qy
节点q处的油膜承载力，δx
t
、δy
t
分别是偏心位置与节点q在x、y方向的距离。
[0125]
组装得到的非线性油膜力矢量f
oil
表达式为：
[0126][0127]
式中：n为有限元模型中的节点总数，p、q分别为有限元模型中左、右滑动轴承位置处的节点编号。
[0128]
12)调用步骤2)中计算得到的m、k、g、fg、f1、f2、步骤11)中计算得到的非线性油膜力矢量f
oil
，并将其代入步骤1)所建立的系统运动微分方程中，求解给定工况下t时刻系统的瞬态动力学响应；
[0129]
13)假设时间间隔为δt，根据初始偏心位置、外部激励载荷、步骤3)中油膜压力分布，计算t＝t+δt时刻轴承的偏心位置(x
t+δt
,y
t+δt
)和速度t时刻和t+δt时刻
轴承偏心位置表达式如下所示：
[0130][0131]
式中：(x
t
，y
t
)是t时刻轴承的偏心位置，(x
t+δt
，y
t+δt
)是t+δt时刻轴承的偏心位置，δt是时间间隔，分别是t时刻偏心位置处x、y方向的速度，和可根据压力分布p等于外部激励载荷求得。
[0132]
14)重复步骤10)至13)，并判断t是否大于等于t
max
。若t﹤t
max
，则继续迭代；若t≧t
max
，则迭代结束。至此便可得到如图5所示的轴心轨迹以及各个时刻系统的瞬态动力学响应，进而就可对系统的非线性动力学特性进行分析。
[0133]
综上所述，本发明提供的一种滑动轴承转子系统非线性动力学特性快速计算方法，通过delaunay三角剖分算法对轴承轴心轨迹可行域区域剖分，采用二维曲面散乱插值函数griddata将轴心轨迹可行域内节点处的刚度、阻尼系数拟合为可行域内任意坐标的函数，在后续求解转子非线性动力学时，结合上述拟合结果、偏心位置、速度以及voronoi图可快速得到给定工况下t时刻系统的非线性油膜力矢量，避免了重复迭代、调用的繁琐过程，解决了传统方法计算效率低下、适用性差等缺点，为滑动轴承转子系统在爆炸、冲击、摇摆、急加/减速等特殊工况下非线性动力学特性的快速计算提供了一种行之有效的方法。