综合拟静力学与动力学的轴承摩擦动力学响应分析方法

技术领域

本发明属于机械工程技术摩擦动力学领域，具体涉及综合拟静力学与动力学的轴承摩擦动力学响应分析方法。

背景技术

固体润滑球轴承广泛应用于航天航空等高端装备，是设备转子系统的关键基础零部件。特别是在液体火箭涡轮泵中，面临高轴向力、高径向力、高转速和固体润滑的极端工况，固体润滑轴承摩擦打滑、自旋和振动加速度等摩擦学和振动特性均严重影响发动机转子系统的稳定运行和工作寿命。作为极端工况下旋转设备安全、稳定和长寿命运行的基础，固体润滑球轴承摩擦动力学特性的准确分析是系统摩擦学、结构学研究和设计的重要前提。现有的滚动轴承分析技术主要集中于油润滑轴承的建模，难以应用于固体润滑轴承；而少有的固体润滑球轴承分析模型也仅局限于静态载荷分析，无法满足工程实际中对变工况/变转速等时变条件下的摩擦动力学分析需求，为关键零部件中轴承的设计和分析带来了一定的挑战。

现有的针对固体润滑球轴承的分析模型仅能考虑轴向载荷，对模型进行了过多的简化，分析能力有限；另外目前的技术并未涉及固体润滑球轴承的动力学建模与分析，主要原因是模型较为复杂，建模和求解均具有一定难度，因此无法满足实际工程应用中对固体润滑球轴承在多载荷、变载荷条件下的摩擦学及动力学分析。

发明内容

本发明的目的是针对目前的滚动轴承分析技术无法满足于工程中对固体润滑球轴承的摩擦动力学分析的问题，提出一种综合拟静力学与动力学的轴承摩擦动力学响应分析方法，该方法首先根据球轴承的安装方向和转动方向建立参考坐标系，在此基础上建立具有一致位移变量的拟静力学和动力学方程；其次根据接触角与接触参数的关系实现固体润滑角接触球轴承摩擦力的快速求解，进而计算得到拟静力学和动力学方程中力和力矩；最后求解拟静力学平衡方程，根据拟静力学方程与动力学方程中位移变量的一致性，将稳态解映射为动力学模型的初值并求解动力学微分方程，获取轴承的摩擦动力学响应。

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

一种综合拟静力学与动力学的轴承摩擦动力学响应分析方法，包括以下步骤：

步骤一：首先根据固体润滑角接触球轴承的安装方向和转动方向建立参考坐标系，根据参考坐标系建立固体润滑球轴承具有一致位移变量的拟静力学平衡方程和动力学微分方程；

步骤二：根据接触角与接触参数的关系，计算得到拟静力学平衡方程中的摩擦力和摩擦力矩；

步骤三：根据摩擦力和摩擦力矩，求解步骤一建立的拟静力学平衡方程，得到稳态分析结果，然后将稳态分析结果作为动力学微分方程的初值，采用四阶龙格库塔数值迭代算法求解动力学微分方程，获取轴承的摩擦动力学响应。

进一步的，根据固体润滑角接触球轴承的安装方向和转动方向，建立参考坐标系，包括：

1)惯性坐标系O-XYZ，将Z轴的正方向设置为角接触轴承的正面方向；

2)钢球随动坐标系o

3)钢球与内、外滚道间的局部接触坐标系o

4)保持架兜孔坐标系o

进一步的，拟静力学平衡方程包括钢球的拟静力学平衡方程与内圈的平衡方程，钢球的拟静力学平衡方程为：

F

式中，F

式中，

当球轴承受到纯轴向力F

F

式中，N

当内圈承受5自由度载荷时，内圈的平衡方程为：

T

其中，T

进一步的，将随动参考系施加的力转换为惯性参考系的矩阵如下：

进一步的，球在随动坐标系中作用在内滚道上的力如下：

进一步的，作用在内圈上的5自由度载荷如下：

F＝(F

其中，

式中，M

进一步的，内圈的平移运动的微分方程为：

式中，

进一步的，保持架的运动微分方程：

式中，I

进一步的，接触参数包括钢球的内和外接触角，接触曲率和以及接触曲率差，其中，钢球的内和外接触角为：

式中，α

接触曲率和与接触曲率差分别表示为：

式中，∑ρ

进一步的，钢球相对于内/外滚道的滑动率和自旋率为：

式中，ξ

根据钢球相对于外滚道的纵向滑动率、横向滑动率和自旋率，通过FASTSIM算法得到钢球受到的摩擦力和摩擦力矩。

与现有技术相比，本发明具有的有益效果如下：

本发明中将相同的位移变量同时用于拟静力学和动力学微分方程的建立，一方面便于技术人员在实际工程应用中根据需求选择模型，另一方面也使得拟静力学的解能够轻易地转换为动力学方程的初值，从而提高动力学模型的求解效率；用一致位移变量来表示滚动体的接触角，根据接触角与接触参数之间的对应关系，预先计算出不同接触角对应的接触参数，在求解方程的过程中可直接利用插值快速计算接触力和摩擦力，能够节约大量的计算时间；最后实现对于轴承摩擦动力学响应的分析。本发明可有效、准确地分析固体润滑球轴承在稳定或时变工况下的摩擦学和动力学特性，为固体润滑球轴承的摩擦学分析设计提供指导。

进一步的，固体润滑球轴承中的钢球受到的摩擦力主要为滚动摩擦，与现有基于油润滑的动力学模型不同，本发明由于引入了滚动接触摩擦力的FASTSIM算法，计算得到钢球受到的摩擦力和摩擦力矩，所提出的动力学模型对于固体润滑球轴承而言，具有更高的分析精度。

进一步的，所提出的拟静力学模型和动力学模型能够同时考虑轴向力、径向力和倾覆力矩的加载情况，分析范围更加广泛，而所提出的拟静力学模型载荷迭代求解框架也保证了模型在复杂工况下也能够有效收敛。

附图说明

图1为建立平衡方程所使用的参考坐标系；其中，(a)为惯性坐标系，(b)为钢球随动坐标系，(c)为保持架接触坐标系；

图2为内圈与钢球的受力示意图；其中，(a)为内圈受力示意图，(b)为钢球受力示意图；

图3为接触区域示意图；其中(a)，为接触区域说明，(b)为采用滚动接触FASTSIM算法时对接触区域进行的有限元网格划分；

图4为拟静力学模型的综合求解策略流程图；

图5为拟静力学分析的接触力稳态结果及与Jones-Harris经典方法对比；其中，(a)为钢球-内滚道接触力(N)，(b)为钢球-外滚道接触力(N)。

图6为动力学分析得到的内圈X和Y方向的振动加速度；其中，(a)为内圈X方向振动加速度，(b)为内圈Y方向振动加速度。

图7为动力学分析得到的钢球角速度动力学响应以及动态摩擦力矩计算结果。其中，(a)为钢球角速度动态响应，(b)为动态摩擦力矩计算结果。

具体实施方式

为了更清楚地描述本发明的具体实施过程，下面结合附图对实施流程进行详细的描述。

本发明的一种综合拟静力学与动力学的轴承摩擦动力学响应分析方法，包括以下步骤：

步骤一：首先根据球轴承的安装方向和转动方向建立参考坐标系，在此基础上建立固体润滑球轴承具有一致位移变量的拟静力学平衡方程和动力学微分方程。

根据固体润滑角接触球轴承的安装方向和转动方向，建立如图1中(a)、(b)和(c)所示四种参考坐标系，包括：

1)惯性坐标系O-XYZ，将Z轴的正方向设置为角接触轴承的正面方向；

2)钢球随动坐标系o

3)钢球与内、外滚道间的局部接触坐标系o

4)保持架兜孔坐标系o

在惯性坐标系和随动坐标系的基础上，定义钢球和内圈的位移变量，作为拟静力学和动力学模型的统一变量。以惯性坐标系为柱坐标基准，定义钢球沿y

钢球通过接触力、摩擦力、离心力、摩擦力矩和陀螺力矩保持受力平衡状态，而这些力和力矩均可由上述定义的变量计算得到。根据如图2中(a)和(b)所示的受力示意图，钢球的拟静力学平衡方程表示为：

F

式中，F

式中，

当球轴承受到纯轴向力F

F

式中，N

T

其中，T

F＝(F

其中，

以上完成了轴承的拟静力学建模，接下来进行动力学方程的建立。

动力学方程中，轴承外圈固定，将内圈简化为2个自由度，即X方向和Y方向的平移运动，保持架简化为具有1个转动自由度，如图1中(c)所示。由于保持架兜孔具有一定的间隙，因此仅当钢球与兜孔壁接触时才产生弹性力。假定某个钢球的旋转位移为θ

式中，

式中，N

式中，M

式中，

式中，I

步骤二：根据接触角与接触参数的关系，实现固体润滑角接触球轴承摩擦力的快速求解，进而计算得到拟静力学和动力学方程中力和力矩。

基于步骤一中提出的钢球与内圈的11个变量，推导接触角，并根据接触角计算接触参数。由运动位移关系将钢球的内、外接触角表达为：

式中，α

接触角与接触斑点的长半轴a、短半轴b、接触刚度κ、椭圆率k紧密相关，计算公式分别为：

a＝[6k

b＝[6εRQ/(πkE')]

式中，ε为椭圆率k的第二类椭圆积分，Q＝Q

式中，ρ

对于椭圆率k，可通过以下公式进行迭代计算：

式中，n和(n+1)表述迭代步数，F

式中，∑ρ

根据式(34)～式(42)可知，a、b、κ均与椭圆率k有关，而椭圆率又与接触角相关，当接触角变化，就必须用迭代的方法求解椭圆率，为了避免在解方程过程中反复求解椭圆率而浪费时间，可以事先计算出接触角在0～90°范围内对应的椭圆率，然后在解方程组时只需要根据不同的接触角插值便可快速获得椭圆率，进而计算a、b、κ的值。

滑动率和自旋率表示钢球相对于内、外滚道的相对滑动和自旋速率，是计算固体润滑球轴承中钢球摩擦力的基础条件。基于坐标系规定的方向，钢球相对于内/外滚道的滑动率和自旋率为：

式中，ξ

利用公式(43)-(48)计算出来的6个滑动率和自旋率，结合轴承的材料参数，可利用滚动接触摩擦力算法FASTSIM得到钢球受到的摩擦力和摩擦力矩，具体如下：

算法假定滚动接触区域前沿的滑动为零，即前沿处的摩擦力为零，然后仅需从前沿向后沿逐渐递推求出各质点的摩擦力，最终求和得到整个接触区域的总摩擦力。接触区域为椭圆形，同时包括滑动区和粘着区，如图3中(a)所示。为了便于计算，对接触区域进行有限元网格划分并对节点进行编号，如图3中(b)所示。为了提高计算程序的精度和通用性，首先对一些参数进行无量纲化处理：

P

U

U

上式中，X和Y为接触区域内质点的无量纲坐标，坐标原点在椭圆的正中心，x和y为质点的有量纲坐标，P

式中，C

式中，

式中，T

式中，P

通过上述FASTSIM算法，计算出钢球与内、外滚道之间的摩擦力和摩擦力矩，最终得到钢球在如式(21)中所示的6个摩擦力和摩擦力矩。摩擦力还会对钢球造成额外的力矩，对应于纵向摩擦力T

式中，M

式中，κ

为了使计算更为精确，还应考虑钢球的惯性效应，钢球的离心力F

式中，m

式中，I

步骤三：首先利用载荷迭代的求解框架求解步骤二建立的拟静力学平衡方程，得到稳态分析结果。然后将拟静力学平衡方程的解(即稳态分析结果)作为动力学微分方程的初值，采用四阶龙格库塔数值迭代算法求解微分方程，获取轴承的摩擦动力学响应。具体过程为：

首先求解拟静力学方程。当轴承仅受轴向载荷时，拟静力学方程应选用式(20)和式(22)，一共7个非线性平衡方程，由经典Jones-Harris模型给定方程初值，采用信赖域算法求解。

当轴承除了受到轴向力作用，还有径向力等其他载荷加载时，各钢球受载不同，每个钢球均需要建立各自的式(20)所示的平衡方程组，加上内圈的5自由度运动，总共6N

针对方程组数量较大，难以一次性求解的问题，采用模块化求解方法。即先逐个求解每个钢球的平衡方程，一次仅需解6个方程，当所有钢球求解完毕，最后求解内圈的平衡方程，一次仅需解5个方程。在每个载荷步内反复迭代，直到内圈和所有钢球均收敛，则停止本载荷步的计算。上述求解策略的详细流程如图4所示：①首先，判断负载类别是纯轴向载荷还是联合载荷。对于纯轴向载荷，可通过处理七个平衡方程直接求解。此外，方程组的初始值最初可以从经典JH模型中获得，该模型不考虑摩擦效应。②其次，对联合负载采用渐进求解策略。由于组合荷载对内圈位移有很大影响，收敛对初始值非常敏感，因此径向荷载和力矩荷载必须从零开始逐渐增加，初始值来自纯轴向载荷模型。在荷载的迭代增加中，在先前荷载F下获得的解应迭代用作增加荷载F'的初始值，直到增加的荷载等于给定荷载F

经过拟静力学方程组的求解，最终可以得到：每个钢球的位移y

式中，x

式中，ω

式中，N

将拟静力学得到的变量作为初值，带入到式(29)～式(31)，采用四阶定步长龙格库塔积分方法求解动力学微分方程，最终可以得到各变量随时间变化的动态结果。相较于拟静力学模型，动力学模型可以计算出变载荷(如变转速、变轴向力等)条件下的分析结果，能够用于轴承装备的加速、平稳、减速的全过程分析，因此适用范围更广。将上述得到的变量时间序列利用式(56)～式(62)可以最终计算出轴承状态角、摩擦力、接触力、振动加速度等动态响应。

利用上述建立的摩擦动力学分析模型和方法，可实现工程实践中不同类型载荷(轴向力、径向力、倾覆力矩、稳态、动态)下固体润滑球轴承的摩擦学和动力学分析，满足真实工况下轴承系统在摩擦学和结构力学等方面的分析、设计和优化需求。

实施例1

本实施例利用如表1所示的B7218角接触球轴承参数进行计算，验证了本发明的有效性。显然，本发明中所开展的实施例并不是全部的实施例，而仅作为代表性的一个实施例，本领域技术人员在本发明基础上所得到的实施例，除非取得创造性成果，其他均属于本发明保护的范围。在如表1所示的联合载荷作用下，轴承拟静力学选用式(20)和式(23)建模，总共16个钢球，每个钢球有6个平衡方程，内圈有5个平衡方程，系统共101个平衡方程。根据图4所示的求解流程，采用信赖域方法求解，最终得到的变量如表2所示，其中仅列出了1号、5号和9号钢球的变量。根据拟静力学结果，结合式(56)-(62)可以得到轴承的稳态分析结果。其中，将钢球公转一周过程中与内、外滚道的接触力计算结果与经典Jones-Harris拟静力学计算结果的对比如图5中(a)和(b)所示，从图中可以发现最大接触载荷出现在X正轴位置，与径向力加载方向相对应，而且本发明模型的计算结果与JH模型基本一致，证明了本发明模型的有效性。此外，用式(68)计算出的内圈稳态摩擦力矩大小为4.43Nm。

表1实施例1使用的角接触球轴承参数与工况

表2实施例1中轴承拟静力学方程的求解结果

动力学计算：设定轴承的动态转速为：0到0.5秒内圈怠速(转速稳定在1000rpm)；然后在0.5秒到1.4秒时间段内转速从1000rpm均匀增加到10000rpm；最后在1.4秒之后稳定在10000rpm，动力学总计算时长为2.5秒。将如图5所示的拟静力学解作为动力学方程的初值，采用四阶定步长龙格库塔方法求解微分方程，步长为10