技术特征：
1.一种滑动轴承转子系统非线性动力学特性快速计算方法，其特征在于，包括以下步骤：1)根据给定的转子参数、圆盘参数、滑动轴承参数、材料参数，建立滑动轴承转子系统euler
‑
bernoulli梁有限元模型和运动微分方程，滑动轴承转子系统euler
‑
bernoulli梁有限元模型的每个节点四个自由度；2)根据步骤1)中有限元模型和运动微分方程，计算总体质量矩阵m、总体刚度矩阵k、总体陀螺矩阵g、总体重力矢量fg、总体不平衡力矢量f1和f2；3)根据给定的滑动轴承参数计算轴承的油膜厚度h(θ)，令h(θ)等于零求解轴承的轴心轨迹可行域s；4)采用delaunay三角剖分算法对可行域s进行网格划分，且靠近边界处，网格细化；5)利用有限差分法求解reynolds方程，计算滑动轴承的油膜压力分布p；6)利用步骤5)中得到的油膜压力分布p，计算油膜承载力f
x
、f
y
；7)利用步骤6)中得到的油膜承载力，采用小扰动法计算可行域s内每个节点处滑动轴承的刚度系数阻尼系数i为可行域s内节点的编号；8)采用二维曲面散乱插值函数griddata分别将步骤7)中计算得到刚度、阻尼系数拟合为关于可行域s内坐标的函数，进而得到轴心轨迹可行域s内任意位置处滑动轴承的刚度系数k
xx
、k
xy
、k
yx
、k
yy
，阻尼系数c
xx
、c
xy
、c
yx
、c
yy
；9)根据偏心参数确定初始时刻轴承的偏心位置(x
t0
,y
t0
)和速度根据给定工况确定迭代次数以及迭代终止时间t
max
，并令t＝t0；10)根据t时刻轴承的偏心位置确定所在的delaunay三角形区域，并利用voronoi图确定delaunay三角形区域内距离偏心位置最近的节点q；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
；12)调用步骤2)中计算得到的m、k、g、fg、f1、f2和步骤11)中计算得到的非线性油膜力矢量f
oil
，并将其代入步骤1)所建立的系统运动微分方程中，求解给定工况下t时刻系统的瞬态动力学响应；13)假设时间间隔为δt，根据初始偏心位置、外部激励载荷、步骤5)中油膜压力分布，计算t＝t+δt时刻轴承的偏心位置(x
t+δt
,y
t+δt
)和速度14)重复步骤10)至13)，并判断t是否大于等于t
max
，若t﹤t
max
，则继续迭代；若t≧t
max
，则迭代结束，得到轴承的轴心轨迹以及各个时刻系统的瞬态动力学响应，进而用于对系统的非线性动力学特性进行分析。2.根据权利要求1所述的一种滑动轴承转子系统非线性动力学特性快速计算方法，其特征在于，步骤1)中，运动微分方程为：式中：m是总体质量矩阵，g是总体陀螺矩阵，k是总体刚度矩阵，fg是总体重力矢量；f1和f2分别为系统沿x、y方向的总体不平衡力矢量，f
oil
是滑动轴承非线性油膜力矢量，q(t)是
广义节点位移，ω是转子转速。3.根据权利要求2所述的一种滑动轴承转子系统非线性动力学特性快速计算方法，其特征在于，步骤2)中，总体质量矩阵m、总体刚度矩阵k、总体陀螺矩阵g、总体重力矢量fg、总体不平衡力矢量f1和f2分别由各单元矩阵组装得到，其中单元质量矩阵m
e
、单元刚度矩阵k
e
、单元陀螺矩阵g
e
可根据转子动力学理论直接获得，单元力矢量fg
e
、f
1e
和f
2e
表达式如下所示：单元重力矢量fg
e
：式中：m
i
是第i个梁单元处转子轴段的质量；单元不平衡力矢量f
1e
和f
2e
：式中：m
e
d是不平衡质量，是不平衡相位角，ω是转子转速，f
1e
表示x方向单元不平衡力矢量，f
2e
表示y方向单元不平衡力矢量；组装得到的总体力矢量fg、总体不平衡力矢量f1和f2表达式如下所示：总体重力矢量fg：f
g
＝[0,m1g,0,0,
……
0,m
i
g,0,0,
……
,0,m
n
‑1g,0,0]
t
ꢀꢀ
(4)式中：m
i
是第i个梁单元处转子轴段的质量，n是有限元模型中的节点总数；总体不平衡力矢量f1和f2：式中：m
e
d是不平衡质量，是不平衡相位角，ω是转子转速，j是不平衡质量位置处的节点编号，f1表示x方向不平衡力矢量，f2表示y方向不平衡力矢量。4.根据权利要求3所述的一种滑动轴承转子系统非线性动力学特性快速计算方法，其特征在于，步骤3)中，对于固定瓦轴承，其油膜厚度表达式如式(6)所示；对于可倾瓦轴承，其油膜厚度表达式如式(7)所示：h(θ)＝c
p
‑
x
j
cos(θ)
‑
y
j
sin(θ)
‑
c
p
mcos(θ
‑
θ
p
)
ꢀꢀꢀ
(6)h(θ)＝c
p
‑
x
j
cos(θ)
‑
y
j
sin(θ)
‑
c
p
mcos(θ
‑
θ
p
)
‑
(r+t)δsin(θ
‑
θ
p
)
ꢀꢀꢀ
(7)式中：θ为周向位置角，c
p
为滑动轴承的半径间隙，x
j
是轴颈周向位移，y
j
是轴颈轴向位移，r为轴承半径，m为预负荷系数，t为预偏心距。5.根据权利要求4所述的一种滑动轴承转子系统非线性动力学特性快速计算方法，其特征在于，步骤5)中reynolds方程为：
式中：g
x
是周向紊流修正因子，g
y
是轴向紊流修正因子，η为润滑油粘度，ρ为润滑油密度，p为油膜压力，x、y分别为周向和轴向坐标；边界条件为：式中：p
s
为供油压力，b为轴承宽度。6.根据权利要求5所述的一种滑动轴承转子系统非线性动力学特性快速计算方法，其特征在于，步骤6)中，x、y方向滑动轴承油膜承载力计算公式为；式中：θ为周向位置角，θ1油膜起始角，θ2油膜破裂角。7.根据权利要求6所述的一种滑动轴承转子系统非线性动力学特性快速计算方法，其特征在于，步骤7)中，假设轴承轴心位于可行域内第i个节点，轴颈在x、y方向产生微小扰动，此时轴颈轴心位置从(x，y)移动到(x0，y0)，则f
xi
、f
yi
、f
x0i
、f
y0i
之间的表达式为：式中：f
xi
、f
yi
是轴颈轴心位置位于(x，y)处的油膜承载力，f
x0i
、f
y0i
是轴颈轴心位置位于(x0，y0)处的油膜承载力；可行域内第i个节点处的刚度阻尼系数为：式中：分别表示第i个节点处x、y方向直接刚度，表示第i个节点处交叉刚度；分别表示第i个节点处x、y方向直接阻尼，表示第i个节点处交叉阻尼。8.根据权利要求7所述的一种滑动轴承转子系统非线性动力学特性快速计算方法，其特征在于，步骤8)中，采用二维曲面散乱插值函数griddata拟合得到的函数如下所示：
式中：(x，y)分别是可行域s内任意一点的坐标，且拟合函数都满足步骤7)中节点处的数据点。9.根据权利要求8所述的一种滑动轴承转子系统非线性动力学特性快速计算方法，其特征在于，步骤11)中，t时刻轴承的油膜力f
oilx
、f
oily
为：式中：(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方向的距离；组装得到的非线性油膜力矢量f
oil
表达式为：式中：n为有限元模型中的节点总数，p、q分别为有限元模型中左、右滑动轴承位置处的节点编号。10.根据权利要求9所述的一种滑动轴承转子系统非线性动力学特性快速计算方法，其特征在于，步骤13)中，t时刻和t+δt时刻轴承偏心位置表达式如下所示：式中：(x
t
，y
t
)是t时刻轴承的偏心位置，(x
t+δt
，y
t+δt
)是t+δt时刻轴承的偏心位置，δt是时间间隔，分别是t时刻偏心位置处x、y方向的速度，和可根据压力分布p等于外部激励载荷求得。
技术总结
本发明公开了一种滑动轴承转子系统非线性动力学特性快速计算方法，该方法解决了滑动轴承转子系统在爆炸、冲击、摇摆、急加/减速等特殊工况下，传统方法计算轴承非线性油膜力效率低、适用性范围局限性的问题；本发明将轴心轨迹可行域内节点处的刚度阻尼系数拟合为可行域内任意位置处的函数，在后续求解转子非线性动力学时，可快速得到给定工况下t时刻系统的非线性油膜力矢量，避免了重复迭代、调用的繁琐过程，使得非线性动力学计算过程效率更高；本发明提出的计算方法包含非线性油膜力计算和转子非线性动力学计算两大部分，详细的考虑了系统的各项参数指标，为滑动轴承转子系统非线性动力学特性研究提供了一种全面、规范的快速计算方法。快速计算方法。快速计算方法。
技术研发人员：裴世源 张庆山 洪军
受保护的技术使用者：西安交通大学
技术研发日：2021.07.07
技术公布日：2021/9/24