一种定位预紧下角接触球轴承非线性刚度的计算系统及方法与流程

本发明涉及机床主轴设计过程中轴承的动刚度计算领域，具体为一种定位预紧下角接触球轴承非线性刚度的计算系统及方法。

背景技术：

随着高速，高精密，高刚度机床的发展，主轴作为高档数控机床的核心部件，其动静态特性直接影响机床的加工精度，而轴承又是主轴的关键部件之一，其动静态特性对主轴的性能有很大影响。在设计一款主轴时，轴承的选择与设计计算是必要的步骤之一，其中包括轴承的刚度的设计计算。轴承的刚度特性与轴承的预紧方式有很大关系，预紧方式分为定位预紧和定压预紧，不同预紧方式下轴承的动态特性有很大差别。对于定压预紧方式下轴承刚度及动态特性的计算已经有很多学者进行研究，在此不做过多叙述；定位预紧方式如图1所示，二个相互配对的角接触球轴承在出厂时就确定了内外圈之间的宽度差，即压入量，将轴承通过轴承座1组装在轴2上之后用锁紧螺母3拧紧消除间隙δx即可使得二个轴承均处于预紧状态，这是定位预紧区别于定压预紧的关键，定压预紧是用弹簧对轴承进行压紧，弹簧的力在轴承运转过程中保持不变，而定位预紧下，轴承在运转过程中，若不考虑热变形的影响，轴承初始压入量δx保持不变，即轴承运行过程中，内外圈的相对位置保持不变。目前对于计算轴承刚度矩阵的计算大多是基于定压预紧进行计算，也有把定位预紧当作定压预紧来计算的，此时在高速下刚度数值会产生很大差别；由于轴承拟静力学模型涉及的方程是高度的非线性方程，而且所要求解的未知量数目较多，求解时经常遇到不收敛，而且求解速度慢的情况，因此有一些学者提出一些刚度计算的经验公式，但是这些经验公式只能对刚度的大小进行粗略估计，而且高速下计算精度和刚度的真实值有很大区别，并且只能粗略估计轴承的主刚度，不能对轴承交叉刚度进行计算，交叉刚度在主轴设计中有非常重要的作用；总之对于定位预紧方式下计算轴承刚度矩阵的计算目前没有提出一个完整的清晰的计算流程。

技术实现要素：

针对现有技术中存在的问题，本发明提供一种定位预紧下角接触球轴承非线性刚度的计算系统及方法，求解速度快、求解精度高和易于收敛，不仅能够得到轴承的主刚度，而且能够求出轴承的刚度矩阵，为电主轴的设计过程中轴承的预紧方式及刚度计算提供了一些基础的理论指导。

本发明是通过以下技术方案来实现：

一种定位预紧下角接触球轴承非线性刚度的计算方法，包括如下步骤，

步骤1，建立滚动轴承的五自由度拟静力学模型；得到滚动体受力平衡方程、定位预紧下轴承变形的几何相容方程和轴承总体受力平衡方程；

步骤2，确定待计算角接触球轴承的结构参数、材料参数、轴承的初始压入量δx、所受的外载荷以及转速参数；

步骤3，根据步骤2中确定的参数，通过求解滚动轴承的静力学模型，得到根据拟静力学模型待计算角接触球轴承整体变形量的迭代初值d⁰＝{δy,δz,θy,θz}⁰和每个滚动体的迭代初值j＝1,2...n，n为滚动体个数；其中，x1,x2,δi,δo分别为滚动体中心与外沟道曲率中心之间的轴向距离、径向距离，滚动体与内、外圈之间的接触变形量，；δy,δz,θy,θz分别为角接触球轴承受力后y轴，z轴变形量以及绕y轴和z轴转角变形；

步骤4，根据步骤3中得到的迭代初值，求解步骤1中得到的滚动体受力平衡方程和几何相容方程，得到满足求解精度的xj＝{x1,x2,δi,δo}j；

步骤5，根据步骤1中的轴承总体受力平衡方程得到轴承刚度矩阵k的具体解析表达式；并且根据步骤4中求出的xj＝{x1,x2,δi,δo}j和初始条件d⁰＝{δy,δz,θy,θz}⁰得到该次迭代的轴承刚度矩阵的具体数值kn；

步骤6，根据轴承总体受力平衡方程得到迭代初值为d⁰＝{δy,δz,θy,θz}⁰时轴承内圈所受的全局载荷向量fn，并将fn和给定的外载荷向量f＝{fy,fz,my,mz}对比，并计算误差值ε2＝||fn-f||；

若满足精度，则执行步骤7，输出结果；

若不满足给定的精度，用步骤5中得到的轴承刚度矩阵的第2行到第5行以及第2列到第5列，形成的4×4的刚度修正矩阵对迭代初值进行如下修正，然后返回步骤4；

步骤7，根据求解所得到的d＝{δy,δz,θy,θz}，xj＝{x1,x2,δi,δo}j，以及初始压入量δx即可求得轴承不同转速，不同载荷下的刚度矩阵；从而得到定位预紧下角接触球轴承非线性刚度。

优选的，步骤1中，滚动体受力平衡方程包括如下的滚动体水平方向平衡方程和竖直方向平衡方程；

式中，qij,qoj分别为第j个滚动体与内、外圈之间的法向接触载荷，αij,αoj为第j个滚动体的内、外接触角，λij,λoj为第j个滚动体的内、外圈载荷分配系数，mgj,fcj分别为第j个滚动体的陀螺力矩和离心力。

优选的，步骤1中，定位预紧下轴承变形的几何相容方程如下，

式中，a1j和a2j为第j个滚动体内、外沟道曲率中心之间的轴向距离和径向距离，x1j,x2j为第j个滚动体中心与外沟道曲率中心之间的轴向距离和径向距离，δij,δoj分别为第j个滚动体与内、外圈之间的接触变形量，fi,fo分别为内、外圈沟曲率半径系数，d为滚动体直径；

其中，

式中，bd为轴承为受力前内、外沟曲率中心之间的距离，α⁰为初始接触角，δx为初始压入量，δy,δz,θy,θz分别为轴承受力后y轴，z轴变形量以及绕y轴和z轴转角变形，ri为内沟道曲率中心轨迹半径，ψj为第j个滚动体的位置角。

优选的，步骤1中，轴承总体受力平衡方程如下，

式中：n为滚动体个数，qij为第j个滚动体与内圈之间的法向接触载荷，αij为第j个滚动体的内接触角，λij为第j个滚动体的内圈载荷分配系数，mgj为第j个滚动体的陀螺力矩，d为滚动体直径，fi为内圈沟曲率半径系数，ψj为第j个滚动体的位置角。

进一步，步骤5中，得到轴承刚度矩阵k的具体解析表达式如下，

优选的，步骤5中求解轴承刚度矩阵时，运用全解析法进行计算。

一种定位预紧下角接触球轴承非线性刚度的计算系统，包括处理器和存储器，以及存储在存储器中通过处理器执行的程序；所述的程序如本发明的方法所述，用于计算定位预紧下角接触球轴承非线性刚度。

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

本发明提出一种完整的定位预紧下轴承刚度矩阵的计算系统及方法，区别与定压预紧下轴承刚度矩阵的计算方法，也区别与利用经验公式计算轴承刚度的方法，通过拟静力学模型、受力平衡方程、定位预紧下轴承变形的几何相容方程和轴承总体受力平衡方程的建立及求解，清晰的分辨出所要求的未知量在轴承平衡方程中的分布情况，在求解过程中巧妙地运用一个4×4的刚度修正矩阵对方程进行求解。从而能够针对实际情况中所有的定位预紧情况进行分析计算，只需给出轴承的相关结构尺寸，材料参数，转速及载荷，就可以对定位预紧下轴承的刚度矩阵进行计算，该方法特指定位预紧下轴承刚度矩阵的求解，计算流程清晰完整，易于编程，给主轴设计过程中定位预紧情况下计算轴承刚度矩阵带来很大的方便，在实际设计中，可以用该程序很方便的计算出轴承刚度是否满足要求。

附图说明

图1为现有技术中定位预紧的结构示意图。

图中：轴承座1，轴2，锁紧螺母3。

图2为本发明实例中所述滚动体载荷分析图。

图3为本发明实例中所述载荷作用前后球中心和沟道曲率中心位置图。

图4为本发明实例中所述整体求解流程图。

图5为本发明实例中所述轴向刚度随转速变化图。

图6为本发明实例中所述径向刚度随径向载荷变化图。

具体实施方式

下面结合具体的实施例对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。

本发明一种定位预紧下角接触球轴承非线性刚度的计算方法，特指在定位预紧方式下求解轴承五自由度拟静力学模型，其中已知条件为轴承初始压入量δx以及外载荷向量f＝{fy,fz,my,mz}，所要求解的未知量为轴向力fx、轴承整体变形量d＝{δy,δz,θy,θz}和每个滚动体的变形量xj＝{x1,x2,δi,δo}j，j＝1,2...n，n为滚动体个数。其中：下标x，y，z分别表示与坐标轴x轴，y轴，z轴有关的量；下标i表示和内圈有关的变量，下标o表示和外圈有关的变量；下标1，2是为了编程方便而设置的变量区别符号，my是绕y轴的弯矩，mz是绕z轴的弯矩。

通过对除轴向平衡方程以外的总体平衡方程、滚动体受力平衡方程和几何相容方程进行迭代求解，得到每次迭代的所有滚动体的求解变量，及载荷误差ε2。

再运用一个4×4的刚度修正矩阵来对轴承的迭代初值进行修正，使得整体求解过程易于收敛。

最后运用全解析法对轴承刚度矩阵进行计算，最终得到定位预紧下角接触球轴承非线性刚度。

如图4所示，具体的方法流程如下所述。

求解轴承刚度矩阵时，首先需要对轴承进行受力分析及建模，受力分析包括对每个滚动体受力分析和对整个轴承进行受力分析，通过受力分析得出滚动体和整个轴承的受力平衡方程。

步骤1,建立滚动轴承的五自由度拟静力学模型；

1.1滚动体受力平衡方程

如图2所示滚动体受力分析图，滚动体水平方向和竖直方向上平衡方程为：

式中，qij,qoj分别为第j个滚动体与内、外圈之间的法向接触载荷，αij,αoj为钢球的内、外接触角，λij,λoj为内、外圈载荷分配系数，mgj,fcj分别为第j个滚动体的陀螺力矩和离心力。

1.2定位预紧下轴承变形的几何相容方程

如图3所示，由勾股定理可得：

式中，a1j和a2j为第j个滚动体内、外沟道曲率中心之间的轴向距离和径向距离，x1j,x2j为第j个滚动体中心与外沟道曲率中心之间的轴向距离和径向距离，δij,δoj分别为第j个滚动体与内、外圈之间的接触变形量，fi,fo分别为内外圈沟曲率半径系数，d为滚动体直径。

其中，

式中，bd为轴承为受力前内、外沟曲率中心之间的距离，α⁰为初始接触角，δx为初始压入量，δy,δz,θy,θz分别为轴承受力后y轴，z轴变形量以及绕y轴和z轴转角变形，ri为内沟道曲率中心轨迹半径，ψj为第j个滚动体的位置角。

1.3轴承总体受力平衡方程；

一般条件下可认为轴承所受外力是通过轴传递到内圈，内圈通过滚动体传递到外圈，在传递到轴承座上面，因此可以把所有滚动体对内圈的力相加就可以得到内圈的平衡方程。

式中：n为滚动体个数。

步骤2，给定一款轴承的具体的结构参数、材料参数、轴承的初始压入量δx、所受的外载荷f＝{fy,fz,my,mz}以及转速参数。

步骤3，根据步骤2给出的条件，先求解轴承的静力学模型，计算出拟静力学模型计算的轴承整体变形量的迭代初值d⁰＝{δy,δz,θy,θz}⁰和每个滚动体的迭代初值j＝1,2...n，n为滚动体个数。

步骤4，根据步骤3给出的迭代初值，求解滚动体受力平衡方程(1)和几何相容方程(2)，求解出在初始条件为d⁰＝{δy,δz,θy,θz}⁰时满足求解精度的xj＝{x1,x2,δi,δo}j，j＝1,2...n，n为滚动体个数。

步骤5，为了对轴承的迭代初值d⁰＝{δy,δz,θy,θz}⁰进行修正，需要根据方程(4)求解轴承刚度矩阵的具体解析表达式，

并且根据步骤4中求出的xj＝{x1,x2,δi,δo}j和初始条件d⁰＝{δy,δz,θy,θz}⁰来求得该次迭代的轴承刚度矩阵的具体数值kn。

步骤6，根据方程(4)的后四个方程可以求得在迭代初值为d⁰＝{δy,δz,θy,θz}⁰时轴承内圈所受的全局载荷向量fn，fn和给定的外载荷向量f＝{fy,fz,my,mz}对比，并计算误差值ε2＝||fn-f||，若满足给定的求解精度tol2，则执行步骤7，输出结果；若不满足给定的求解精度tol2，用步骤5求出的轴承刚度矩阵的第2行到第5行以及第2列到第5列，形成的4×4的刚度修正矩阵对迭代初值进行修正，然后返回步骤4，直到ε2小于给定的精度。

步骤7，根据求解所得到的d＝{δy,δz,θy,θz}，xj＝{x1,x2,δi,δo}j，j＝1,2...n，n为滚动体个数，以及初始压入量δx即可求得轴承不同转速，不同载荷下的刚度矩阵。

总的来说，轴承在定位预紧时，已知的是轴承的压入量δx，和轴承所受载荷f＝{fy,fz,my,mz}，所要求解的变量包括轴向力fx，轴承的整体变量d＝{δy,δz,θy,θz}，以及每个滚动体的未知量xj＝{x1,x2,δi,δo}j，j＝1,2...n，n为滚动体个数。此时由于轴向载荷未知量fx只在整体方程(4)中的轴向平衡方程中含有，滚动体受力平衡方程(1)和几何相容方程(2)以及整体方程(4)中的后4个方程正好含有的未知量个数为4n+4个，不包含轴向载荷求解变量fx，方程总数也是4n+4个；因此只需先求解滚动体的4n个未知量x＝{x1,x2,δi,δo}j，j＝1,2...n，n为滚动体个数，以及整个轴承的4个未知量d＝{δy,δz,θy,θz}，由于轴承的压入量已知，因此在迭代过程中不需要对轴向位移进行求解，因此取合理的初始迭代向量d⁰＝{δy,δz,θy,θz}⁰，不满足收敛条件时，用轴承刚度矩阵的第2行到第5行以及第2列到第5列，形成的4×4的刚度修正矩阵进行修正即可，满足收敛条件时，把所有的求解变量带入5×5的刚度矩阵即可求出轴承的刚度矩阵。

上述计算过程可以通过应用matlab编程实现。为了更具体证明本方法的有效性，本发明提供一个具体的计算实例。以nsk的一款角接触球轴承为例，具体参数如下表1所示，

表1轴承参数

分别计算了该轴承在初始压入量为10μm,17μm,25μm，转速从0rpm到20000rpm，间隔为1000rpm，计算轴承的轴向刚度，如图5所示，可以看到随转速的增加，刚度有衰减线性，符合一般的规律；

另外计算了初始压入量在17μm，转速为5000rpm时，径向力分别为0n，100n，200n，300n，400n，500n，600n，轴承的径向刚度，如图6所示，可以看到随着径向载荷的增加，径向刚度也是增加的，也是符合一般规律，则证明该方法是有效的。