弹流润滑条件下计算滚动轴承载荷和压力的边界元法的制作方法

本发明属于轧机多列滚动轴承数值分析领域，具体涉及一种弹流润滑条件下计算滚动轴承载荷和压力的边界元法。

背景技术

轧机多列滚动轴承作为轧机的关键部件，在工作中承受着较大的径向负荷、轴向负荷和热负荷，致使多列滚动轴承极易产生偏载，使轴承的寿命急剧下降，并频繁发生异常烧损和大面积疲劳剥落等事故，严重制约了生产效率的提高。其润滑性能好坏对轧机滚动轴承的寿命、振动等性能有着很大的影响，在润滑良好时，轴承的载荷及其分布成为影响其运行行为的主要因素。因此对轧机多列滚动轴承考虑弹流润滑情况下的载荷分布进行研究，研究具有重大意义。

对轧机多列滚动轴承载荷分布的研究中，常采用的方法是有限元法，但基本上都没有考虑滚子弹流润滑对轴承载荷分布的影响，由于有限元法划分单元较多，在高度非线性问题下，有限元计算结果并不理想；

三维有摩擦弹性接触边界元法部分，参见【杨霞，轧机四列圆锥滚子轴承的边界元法理论与实验研究[d]，燕山大学，2011】。

技术实现要素：

鉴于现有技术的上述缺陷，本发明所要解决的技术问题是提供一种弹流润滑条件下计算滚动轴承载荷和压力的边界元法，以满足计算要求并解决现有方法的不足。

为实现上述目的，本发明提供了一种弹流润滑条件下计算滚动轴承载荷和压力的边界元法，包括以下步骤，

步骤1.在三维有摩擦弹性接触边界元法基础上引入轴承边界单元、接触宽度赫兹修正理论和板单元，建立轴承边界元法；

步骤2.根据有限长线接触弹流润滑理论，推导摩擦系数方程，并对基本方程进行量纲一化；

步骤3.用fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序；

步骤4.将轴承边界元法和有限长线接触弹流润滑理论耦合，建立滚动轴承弹流润滑条件下的轴承边界元法；

步骤5.编写滚动轴承弹流润滑条件下的轴承边界元法的fortran计算程序；

步骤6.建立轧辊、轴承内圈、轴承外圈和轴承座的表面单元离散模型，并进行数据的前处理；

步骤7.将前处理好的轧辊、轴承内圈、轴承外圈和轴承座的节点坐标和单元组成信息代入编制的考动轴承弹流润滑条件下的轴承边界元法的fortran计算程序进行计算；

步骤8.对程序计算的数据进行后处理，得到滚动轴承的载荷和压力分布数据。

进一步，所述步骤1中的三维有摩擦弹性接触边界元法为，

(1)增量形式的边界积分方程

考虑两个线弹性体相互接触，a物体的边界记为γa，b物体的边界记为γb，接触区边界记为γc，位移已知边界记为γu，面力已知边界记为γt；

采用等增量加载方法，假设总的增量步为n步，当第m步增量加载时，用增量表示的边界积分方程为：

其中，k表示a，b两个物体，x为源点，y为场点；cij(x)是与边界点x处几何形状有关的函数。如果边界点x处是光滑的，也就是说x点处的外法线矢量是连续的，那么cij(x)＝δij/2。分别表示k物体在第m步增量加载时j方向的位移和面力；积分核函和分别是弹性问题位移和面力的基本解；

下面是基本解函数的表达式：

式中，i,j,k,l＝1，2，3，ri＝yi-xi，其中δij为kronecker函数：

总的量用增量的和来计算，因此第m步增量加载后总的位移和面力分别为：

其中，分别为第m步增量加载后j方向总的位移和面力；分别为第m步增量加载时j方向的位移和面力；

(2)离散形式的边界积分方程

由于工程问题中通常会遇到几何形状及边界条件都比较复杂的问题，因此，需要采用离散技术，进行数值求解。本发明采用4节点等参单元对边界进行离散。

对于四边形4节点线性单元，其插值函数为：

n^β(ξ1,ξ2)＝(1+rβξ1)(1+sβξ2)/4β＝1,2,3,4；(7)

其中，rβ为第β个节点ξ1方向的局部坐标分量；sβ为第β个节点ξ2方向的局部坐标分量。

单元内部任意一点的整体坐标、位移和面力分别由单元节点上的坐标、位移和面力来描述，设分别表示第ne单元中的β节点的坐标、位移和面力分量，则：

其中，q为单元的节点个数；n^β为形函数，即插值函数。

在边界元法中，将k物体的整个边界划分为n^k个单元，然后将式(8)代入到增量边界积分方程式(1)中，得到离散形式的增量边界积分方程：

其中，n^k是k物体离散单元总数；分别为k物体在第m步增量加载中第n单元上的β节点的位移和面力；n^β(ξ1,ξ2)为单元上的插值函数；g(ξ1,ξ2)为单元的雅可比变换系数。

(3)建立局部坐标系

本发明在求解弹性接触问题的边界积分方程时采用凝缩解法，离散模型的接触模式采用点对点接触模式。由于未知量数目多于方程数目，为了补充方程，从而可以求得未知量，故需在接触区γc上建立局部坐标系(ξ1,ξ2,ξ3)。

对于接触区为曲面的问题，在接触区上，a和b物体相对应的节点的曲率不同，为了保证计算结果的精度，只在a物体的接触区上建立局部坐标系，并且保证局部坐标系的ξ3方向为a物体接触边界的外法线方向，ξ1和ξ2方向只需要满足右手定则即可，

局部坐标系(ξ1,ξ2,ξ3)相对于整体坐标系(x1,x2,x3)下的方向余弦为αlj，

在接触区γc上，节点的位移和面力在两坐标系下的关系如下式：

其中，分别为k物体i节点j方向的位移和面力；分别为k物体i节点在局部坐标系的ξl方向的位移和面力；为a物体i节点在局部坐标系下ξl与整体坐标系下xj的方向余弦。

在接触区建立局部坐标系以后，增量边界积分方程式(9)可表示为：

其中，αlj为局部坐标系ξl在整体坐标系xj上的方向余弦；n^kc、n^kf分别为k物体的接触区单元数和非接触区单元数；分别为k物体第m步增量加载时第n单元的β节点在局部坐标系下的ξl方向的增量位移。

对于给定的两接触物体a和b，令ya和yb为相互接触的节点对，初始间隙为δ0，第m步增量加载状态的起始间隙为δm，则：

其中，为k物体第i步增量加载时在局部坐标系下的ξ3方向上的位移。

(4)接触节点状态

更进一步，所述接触区γc上的接触节点对ya和yb在第m步增量加载时，可能处于分离状态、粘着状态和滑移状态三种状态之一；

当处于分离状态时，ya和yb暂时没有接触，但随着外载荷的增大或其它状态的改变，可能会进入接触状态，处于分离状态的节点满足关系式：

当处于粘着状态时，ya和yb已经接触，但没有发生相对滑动，处于粘着状态的节点满足关系式：

当处于滑移状态时，ya和yb已经接触，且发生了相对滑动，处于滑移状态的节点满足关系式：

当摩擦系数μ一定时，如果接触区上的某一点从第m-1步增量的粘着状态变为第m步增量的滑移状态，且ya点相对于yb点的滑移方向与ξ1轴的夹角为φ^k，得到滑移角之间的关系；

(5)耦合的增量边界积分方程

将上述接触状态关系式(13)～(15)和滑移角关系式(16)代入到增量边界积分方程(1)中，分别得到耦合的物体a的边界积分方程：

同理，将滑移角的关系式代入增量积分方程式(1)中，得到物体b的简化的增量边界积分方程：

进一步，所述步骤1中的轴承边界单元理论是假设轴承内圈有n个滚动体，则沿圆周方向轴承内圈被划分为2n个单元；同理，对于轴承外圈按一样的模式进行划分；在轴承内圈的离散模型中任取一个轴承边界单元i，此单元一定处于接触状态或者非接触状态；若此单元处于接触状态，则其上的接触面力等于滚动体上的接触载荷，因而在该单元上具有位移连续但面力不连续的特点；

此时，一个轴承边界单元分为两个子单元，子单元上作用有连续面力，上面力为零；假定子单元上法向面力沿宽度方向呈抛物线分布，在长度方向上呈线性分布；

子单元的面力增量形式表示如下：

其中，形函数表达式如下，

其中，为i单元为轴承边界单元i时，第m步增量加载时子单元在法线方向上的面力；为i单元为轴承边界单元ii时，第m步增量加载时子单元在法线方向上的面力；为i单元β节点在第m步增量加载时法线方向上的力；分别为轴承边界单元i和轴承边界单元ii的形函数；是轴承边界单元的子单元的内侧边在ξ2方向的坐标值(为接触区的半宽度；所述划分为至少两个微单元，并使微单元的长宽比小于3，且保证奇异点所在的微单元的长宽度比等于1；

所述微单元进行坐标变换如下：

所述微单元区域为ξa≤ξ1≤ξb，ηa≤ξ2≤ηb，采用高斯积分时，应在区域-1≤ξ≤1，-1≤η≤1下进行；

进一步，所述步骤1中的接触宽度赫兹修正理论是，滚动轴承系统的高度非线性，在计算过程中需要初步假设滚动体与轴承外圈的接触宽度，然后再用赫兹接触理论对该接触宽度进行修正；对于轴承边界单元i，经过计算得到各个接触节点的压力值为t^i,β，β表示i单元的节点编号。若此单元为轴承边界单元i，设子单元上的合力为tiⁱ，需要逐个在微单元上进行积分求和，

计算tiⁱ的离散形式如下，

若为轴承边界单元ii，则子单元上的合力为tiⁱⁱ,计算公式如下，

其中，n为微单元个数；l为高斯积分点数；γ1、γ2为第m个微单元对应的坐标转换系数；为第r个积分点对应的权系数；ξ1、ξ2为高斯积分点坐标；分别为轴承边界单元i、ii的形函数；g为雅可比值；

对于圆锥或圆柱滚子轴承，利用赫兹接触理论，滚动体与内、外圈接触的接触半宽度按下式计算，

其中，

其中，bi，bo分别为滚动体与轴承内、外圈的接触半宽度；ri为内圈外径；ro为外圈内径；r为滚动体半径；w为单位长度上的载荷；l为滚动体有效长度。

进一步，所述步骤1中的板单元理论是，滚动体在滚动轴承中对接触状态的影响只关系到径向位移，而其余方向位移对接触状态均无影响；此时，两个圆柱体接触，弹性趋变形量为：

而对于圆柱或圆锥滚子轴承，滚动体的本构方程为：

其中，r为滚动体半径；ri为内圈半径；ro为外圈半径；l为是滚动体长度；e、ei和eo分别为滚动体、内圈和外圈的弹性模量；υ、υi和υo分别为滚动体、内圈和外圈的泊松比；分别为点和yi^r在局部坐标系下ξ3方向的位移；分别为点yi^a和在局部坐标系下ξ3方向的力，bi和bo值如式(24)所示；

此时，考虑滚动体弹性变形、润滑油膜厚度和粗糙度幅值后，第m步增量加载时的间隙为，

其中，δ0为初始间隙，ham-1为第m-1步增量加载时的平均油膜厚度值，其初始值为0；ra为粗糙度幅值；

进一步，所述步骤1中的轴承边界元法，将轴承边界单元、接触宽度赫兹修正理论和板单元引入到三维有摩擦弹性接触边界元法中，并假定板单元固接于轴承内圈上，离散边界以后，得到轴承边界元法离散的边界积分方程如下：

其中，k为第k个物体；n^f为非接触单元数；n^bi为轴承边界单元i的数量；n^bii为轴承边界单元ιι的数量；n^c为接触单元总数。

进一步，所述步骤2中，根据有限长线接触弹流润滑理论，其基本方程如下：

(1)reynolds方程：

(2)reynolds方程的边界条件：

(3)油膜厚度方程：

式中，h0为刚体油膜厚度，r为当量曲率半径，-”为滚动体与内圈，“+”代表滚动体与外圈；

(4)弹性变形方程：

其中，e′为综合弹性模量，e1和e2分别代表两接触曲面的弹性模量，υ1和υ2代表两接触材料的泊松比；

(5)粗糙度方程：

式中，ra表示粗糙度波峰幅值，l′表示粗糙度波长，β′表示粗糙度纹理方向；

(6)粘压方程：

η＝η0exp{(lnη0+9.67)[-1+(1+p0p)^z]}；(35)

式中，α为粘压系数，z取0.68，p0为压力系数，取5.1×10^-9；

(7)密压方程：

(8)载荷方程：

其中w的表达式如式(25)所示；

(9)摩擦系数方程：

采用ree-eyring型非牛顿流体来对圆柱滚子轴承接触区域的摩擦力进行计算，进而得到摩擦系数；

ree-eyring模型本构方程是：

式中，τ0为润滑油的极限剪切力，ha为滚子弹流接触区域的平均油膜厚度，η为接触区的粘度参数，u和v分别表示滚动方向的吸卷速度和润滑油的端泄速度；

所以，可以推导得出摩擦系数的计算公式为：

其中w的表达式如式(25)所示；

进一步，所述步骤2中，将有限长线接触弹流润滑理论的一系列方程进行量纲一化，

p为无量纲油膜压力，ph为赫兹接触压力，b1为赫兹修正后的滚动体接触宽度。

(1)无量纲reynolds方程：

(2)无量纲reynolds方程的边界条件：

(3)无量纲油膜厚度：

(4)无量纲弹性变形方程：

(5)无量纲粗糙度方程：

式中，为无量纲粗糙度幅值，为无量纲粗糙度波长，

(6)无量纲粘压方程：

η^*＝exp{(lnη0+9.67)[-1+(1+p0php)^z]}；(46)

式中，α为粘压系数，z通常取0.68，p0为压力系数，通常取5.1×10^-9。

(7)密压公式：

(8)载荷方程：

式中，为圆柱滚子母线无量纲长度，l为母线长度。

(9)摩擦系数方程：

进一步，所述步骤3中，用fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序，将量纲一化的有限长线接触reynolds方程(41)等式左边中心差分，等式右面向后差分进行差分格式的离散得到的reynolds离散形式的差分方程：

整理后得，

式中，

油膜厚度离散方程可写成：

线接触弹性变形的离散方程：

载荷离散方程：

并根据上述离散方程，用fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序。

更进一步，所述有限长线接触弹流润滑计算程序包括以下步骤，

step1，输入轴承内径、外径、滚动体直径、长度、载荷值和接触宽度值；设定初始油膜厚度h0和初始油膜压力p0；

step2，设定粗糙度幅值、波长和纹理角度；

step3，计算粘度、密度、弹性变形量v和膜厚h；

step4,采用有限差分法求解reynolds方程计算出新的油膜厚度h和油膜压力p；

step5，对比迭代前后压力相对误差，判断是否收敛；

step6，对比迭代前后荷载相对误差，判断是否收敛；

step7，计算量纲化后的油膜厚度和压力值以及弹流接触区域的平均有膜厚度；

step8，计算摩擦系数并输出计算结果。

进一步，当所述step5为不收敛时，修正初始油膜压力p0，并重新进行step3—step8；

当所述step6为不收敛时，修正初始油膜厚度h0，并重新进行step3—step8。

进一步，所述步骤4中，将轴承边界元法和有限长线接触弹流润滑理论耦合，建立滚动轴承弹流润滑条件下的轴承边界元法，其具体的方法为，

将弹流润滑理论计算出的油膜厚度(32)代入间隙值的计算公式(28)，把摩擦系数方程(39)代入边界积分方程(17)和(18)，并假定板单元固接于轴承内圈上，轴承内圈与轧辊固定，轴承外圈与轴承座固定，轴承系统简化为两个接触物体，耦合的离散边界积分方程如下：

其中，间隙δm表达式如式(28)所示；摩擦系数公式(39)中载荷w由式(25)、(22)和(23)确定。k为第k个物体，k＝1,2；n^f为非接触单元数；n^bι为k物体轴承边界单元i的数量；n^bii为k物体轴承边界单元ii的数量；为预接触区单元数量，为黏着区单元数量，为滑移区单元数量，为处于黏着状态的轴承边界单元i的数量，为处于滑移状态的轴承边界单元i的数量，为处于黏着状态的轴承边界单元ii的数量，为处于滑移状态的轴承边界单元ii的数量。

进一步，所述步骤5中，所述编写滚动轴承弹流润滑条件下的轴承边界元法的fortran计算程序包括以下步骤，

step1,读入轧辊和轴承座的节点坐标、单元组成信息及边界条件；

step2,假设接触半宽b0,初始摩擦系数μ0和初始间隙δ0；

step3,计算所有物体的积分系数；

step4,对系数矩阵方程进行高斯消去，得到接触区矩阵方程；

step5,对接触区系数矩阵耦合得到耦合区的系数矩阵；

step6,用边界元法计算接触区的面力和位移；

step7,对接触区进行接触状态和摩擦状态判断；

step8,满足接触状态收敛准则和摩擦状态收敛准则；

step9，修正滚动体的接触宽度b1，并计算所有滚动体上的载荷w；

step10,判断是否满足|b1-b0|≤ε；ε为设定的精度值；

step11,判断是否为第一次迭代计算。

进一步，在所述step10中，

当满足时，执行step11；

不满足时，调用弹流润滑程序，利用接触宽度b1和载荷w计算油膜厚度平均值ha和摩擦系数μ；并假设接触半宽b1,摩擦系数μ，间隙δ0+ha-ra，重新进行step4-step10。

进一步，在所述step11中，

当满足时，调用弹流润滑程序，利用接触宽度b1和载荷w计算油膜厚度平均值ha和摩擦系数μ；并假设接触半宽b1,摩擦系数μ，间隙δ0+ha-ra，重新进行step4-step10；

不满足时，将接触区未知量代入系数矩阵方程计算非接触区的未知量。

所述步骤6中，建立轧辊和轴承座的表面单元离散模型，并进行数据的前处理；利用marc有限元软件，按照物体的实际尺寸建立表面单元模型并导出各物体的节点信息和单元组成。由于滚动轴承弹流润滑条件下的轴承边界元法计算程序采用凝缩解法，要求输入数据所有物体的接触节点和单元排在后面，且一一对应。需要对marc导出的数据进行前处理。

进一步，步骤7.将前处理好的轧辊和轴承座的节点坐标和单元组成信息代入编制的滚动轴承弹流润滑条件下的轴承边界元法的fortran计算程序进行计算；计算结束后输出数据格式的结果。

进一步，步骤8.对程序计算的数据进行后处理，得到滚动轴承的载荷和压力分布数据。

本发明的有益效果是：建立了采用有限长线接触弹流润滑理论，推导出摩擦系数公式，将油膜厚度作为间隙值，摩擦系数公式与边界积分方程耦合，用板单元模拟滚动体，轴承边界单元实现轴承面力不连续现象，用赫兹接触理论修正滚动体与轴承滚道的接触半宽，编制弹流润滑条件下的轴承边界元法程序；与现有轴承载荷数值分析方法相比，本发明方法由于将弹流润滑理论与边界元法相结合，能够分析多列滚动轴承、轧辊和轴承座全模型下的轴承载荷和压力分布；能够通过改变粗糙度影响因素，分析不同润滑条件下轴承的载荷和压力的分布；只需对轴承、轧辊和轴承座进行表面单元的划分，单元划分少，计算精度高。

附图说明

图1为两物体接触模型；

图2为4节点线性单元；

图3为局部坐标系示意图；

图4为滑移方向示意图；

图5为轴承内圈离散示意图

图6为轴承边界子单元；

图7为轴承边界微单元；

图8为板单元示意图；

图9有限长线接触模型；

图10线接触弹流润滑问题计算流程图；

图11为考虑轴承弹流润滑的边界元法计算流程图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明：一种弹流润滑条件下计算滚动轴承载荷和压力的边界元法，包括以下步骤，

步骤1.在三维有摩擦弹性接触边界元法基础上引入轴承边界单元、接触宽度赫兹修正理论和板单元，建立轴承边界元法；

步骤2.根据有限长线接触弹流润滑理论，推导摩擦系数方程，并对基本方程进行量纲一化；

步骤3.用fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序；

步骤4.将轴承边界元法和有限长线接触弹流润滑理论耦合，建立滚动轴承弹流润滑条件下的轴承边界元法；

步骤5.编写滚动轴承弹流润滑条件下的轴承边界元法的fortran计算程序；

步骤6.建立轧辊、轴承内圈、轴承外圈和轴承座的表面单元离散模型，并进行数据的前处理；

步骤7.将前处理好的轧辊、轴承内圈、轴承外圈和轴承座的节点坐标和单元组成信息代入编制的考动轴承弹流润滑条件下的轴承边界元法的fortran计算程序进行计算；

步骤8.对程序计算的数据进行后处理，得到滚动轴承的载荷和压力分布数据。

进一步，所述步骤1中的三维有摩擦弹性接触边界元法为，

(1)增量形式的边界积分方程

考虑两个线弹性体相互接触，如图1所示，物体的边界记为γa，b物体的边界记为γb，接触区边界记为γc，位移已知边界记为γu，面力已知边界记为γt；

采用等增量加载方法，假设总的增量步为n步，当第m步增量加载时，用增量表示的边界积分方程为：

其中，k表示a，b两个物体，x为源点，y为场点；cij(x)是与边界点x处几何形状有关的函数。如果边界点x处是光滑的，也就是说x点处的外法线矢量是连续的，那么cij(x)＝δij/2。分别表示k物体在第m步增量加载时j方向的位移和面力；积分核函和分别是弹性问题位移和面力的基本解；

下面是基本解函数的表达式：

式中，i,j,k,l＝1，2，3，ri＝yi-xi，其中δij为kronecker函数：

总的量用增量的和来计算，因此第m步增量加载后总的位移和面力分别为：

其中，分别为第m步增量加载后j方向总的位移和面力；分别为第m步增量加载时j方向的位移和面力；

(2)离散形式的边界积分方程

由于工程问题中通常会遇到几何形状及边界条件都比较复杂的问题，因此，需要采用离散技术，进行数值求解。本发明采用4节点等参单元对边界进行离散，离散单元模型如下图2所示。

对于四边形4节点线性单元，其插值函数为：

n^β(ξ1,ξ2)＝(1+rβξ1)(1+sβξ2)/4β＝1,2,3,4；(7)

其中，rβ为第β个节点ξ1方向的局部坐标分量；sβ为第β个节点ξ2方向的局部坐标分量。

单元内部任意一点的整体坐标、位移和面力分别由单元节点上的坐标、位移和面力来描述，设分别表示第ne单元中的β节点的坐标、位移和面力分量，则：

其中，q为单元的节点个数；n^β为形函数，即插值函数。

在边界元法中，将k物体的整个边界划分为n^k个单元，然后将式(8)代入到增量边界积分方程式(1)中，得到离散形式的增量边界积分方程：

其中，n^k是k物体离散单元总数；分别为k物体在第m步增量加载中第n单元上的β节点的位移和面力；n^β(ξ1,ξ2)为单元上的插值函数；g(ξ1,ξ2)为单元的雅可比变换系数。

(3)建立局部坐标系

本发明在求解弹性接触问题的边界积分方程时采用凝缩解法，离散模型的接触模式采用点对点接触模式。由于未知量数目多于方程数目，为了补充方程，从而可以求得未知量，故需在接触区γc上建立局部坐标系(ξ1,ξ2,ξ3)。

对于接触区为曲面的问题，在接触区上，a和b物体相对应的节点的曲率不同，为了保证计算结果的精度，只在a物体的接触区上建立局部坐标系，并且保证局部坐标系的ξ3方向为a物体接触边界的外法线方向，ξ1和ξ2方向只需要满足右手定则即可，如图3所示为接触区上建立局部坐标系的示意图。

局部坐标系(ξ1,ξ2,ξ3)相对于整体坐标系(x1,x2,x3)下的方向余弦为αlj，

在接触区γc上，节点的位移和面力在两坐标系下的关系如下式：

其中，分别为k物体i节点j方向的位移和面力；分别为k物体i节点在局部坐标系的ξl方向的位移和面力；为a物体i节点在局部坐标系下ξl与整体坐标系下xj的方向余弦。

在接触区建立局部坐标系以后，增量边界积分方程式(9)可表示为：

其中，αlj为局部坐标系ξl在整体坐标系xj上的方向余弦；n^kc、n^kf分别为k物体的接触区单元数和非接触区单元数；分别为k物体第m步增量加载时第n单元的β节点在局部坐标系下的ξl方向的增量位移。

对于给定的两接触物体a和b，令ya和yb为相互接触的节点对，初始间隙为δ0，第m步增量加载状态的起始间隙为δm，则：

其中，为k物体第i步增量加载时在局部坐标系下的ξ3方向上的位移。

(4)接触节点状态

更进一步，所述接触区γc上的接触节点对ya和yb在第m步增量加载时，可能处于分离状态、粘着状态和滑移状态三种状态之一；

当处于分离状态时，ya和yb暂时没有接触，但随着外载荷的增大或其它状态的改变，可能会进入接触状态，处于分离状态的节点满足关系式：

当处于粘着状态时，ya和yb已经接触，但没有发生相对滑动，处于粘着状态的节点满足关系式：

当处于滑移状态时，ya和yb已经接触，且发生了相对滑动，处于滑移状态的节点满足关系式：

当摩擦系数μ一定时，如果接触区上的某一点从第m-1步增量的粘着状态变为第m步增量的滑移状态，且ya点相对于yb点的滑移方向与ξ1轴的夹角为φ^k，如图4所示，得到滑移角之间的关系；

(5)耦合的增量边界积分方程

将上述接触状态关系式(13)～(15)和滑移角关系式(16)代入到增量边界积分方程(1)中，分别得到耦合的物体a的边界积分方程：

同理，将滑移角的关系式代入增量积分方程式(1)中，得到物体b的简化的增量边界积分方程：

进一步，所述步骤1中的轴承边界单元理论是以轴承内圈为例，其离散模型及滚动体的位置关系如图5所示。假设轴承内圈有n个滚动体，则沿圆周方向轴承内圈被划分为2n个单元；同理，对于轴承外圈按一样的模式进行划分；在轴承内圈的离散模型中任取一个轴承边界单元i，此单元一定处于接触状态或者非接触状态；若此单元处于接触状态，则其上的接触面力等于滚动体上的接触载荷，因而在该单元上具有位移连续但面力不连续的特点；

此时，一个轴承边界单元分为两个子单元，子单元上作用有连续面力，上面力为零；假定子单元上法向面力沿宽度方向呈抛物线分布，在长度方向上呈线性分布,如图6所示；

子单元的面力增量形式表示如下：

其中，形函数表达式如下，

其中，为i单元为轴承边界单元i时，第m步增量加载时子单元在法线方向上的面力；为i单元为轴承边界单元ii时，第m步增量加载时子单元在法线方向上的面力；为i单元β节点在第m步增量加载时法线方向上的力；分别为轴承边界单元i和轴承边界单元ii的形函数；是轴承边界单元的子单元的内侧边在ξ2方向的坐标值(为接触区的半宽度；所述划分为至少两个微单元，并使微单元的长宽比小于3，且保证奇异点所在的微单元的长宽度比等于1，如图7所示；

所述微单元进行坐标变换如下：

所述微单元区域为ξa≤ξ1≤ξb，ηa≤ξ2≤ηb，采用高斯积分时，应在区域-1≤ξ≤1，-1≤η≤1下进行；

进一步，所述步骤1中的接触宽度赫兹修正理论是，滚动轴承系统的高度非线性，在计算过程中需要初步假设滚动体与轴承外圈的接触宽度，然后再用赫兹接触理论对该接触宽度进行修正；对于轴承边界单元i，经过计算得到各个接触节点的压力值为t^i,β，β表示i单元的节点编号。若此单元为轴承边界单元i，设子单元上的合力为tiⁱ，需要逐个在微单元上进行积分求和，

计算tiⁱ的离散形式如下，

若为轴承边界单元ii，则子单元上的合力为ti^ιι,计算公式如下，

其中，n为微单元个数；l为高斯积分点数；γ1、γ2为第m个微单元对应的坐标转换系数；为第r个积分点对应的权系数；ξ1、ξ2为高斯积分点坐标；分别为轴承边界单元i、ii的形函数；g为雅可比值；

对于圆锥或圆柱滚子轴承，利用赫兹接触理论，滚动体与内、外圈接触的接触半宽度按下式计算，

其中，

其中，bi，bo分别为滚动体与轴承内、外圈的接触半宽度；ri为内圈外径；ro为外圈内径；r为滚动体半径；w为单位长度上的载荷；l为滚动体有效长度。

进一步，所述步骤1中的板单元理论是，滚动体在滚动轴承中对接触状态的影响只关系到径向位移，而其余方向位移对接触状态均无影响；图8为板单元示意图，此时，两个圆柱体接触，弹性趋变形量为：

而对于圆柱或圆锥滚子轴承，滚动体的本构方程为：

其中，r为滚动体半径；ri为内圈半径；ro为外圈半径；l为是滚动体长度；e、ei和eo分别为滚动体、内圈和外圈的弹性模量；υ、υi和υo分别为滚动体、内圈和外圈的泊松比；分别为点和yi^r在局部坐标系下ξ3方向的位移；分别为点yi^a和在局部坐标系下ξ3方向的力，bi和bo值如式(24)所示；

此时，考虑滚动体弹性变形、润滑油膜厚度和粗糙度幅值后，第m步增量加载时的间隙为，

其中，δ0为初始间隙，ham-1为第m-1步增量加载时的平均油膜厚度值，其初始值为0；ra为粗糙度幅值；

进一步，所述步骤1中的轴承边界元法，将轴承边界单元、接触宽度赫兹修正理论和板单元引入到三维有摩擦弹性接触边界元法中，并假定板单元固接于轴承内圈上，离散边界以后，得到轴承边界元法离散的边界积分方程如下：

其中，k为第k个物体；n^f为非接触单元数；n^bi为轴承边界单元i的数量；n^bii为轴承边界单元ii的数量；n^c为接触单元总数。

进一步，所述步骤2中，根据有限长线接触弹流润滑理论，其基本方程如下：

(1)reynolds方程：

(2)reynolds方程的边界条件：

(3)油膜厚度方程：

式中，h0为刚体油膜厚度，r为当量曲率半径，-”为滚动体与内圈，“+”代表滚动体与外圈；

(4)弹性变形方程：

其中，e′为综合弹性模量，e1和e2分别代表两接触曲面的弹性模量，υ1和υ2代表两接触材料的泊松比；

(5)粗糙度方程：

式中，ra表示粗糙度波峰幅值，l′表示粗糙度波长，β′表示粗糙度纹理方向；

(6)粘压方程：

η＝η0exp{(lnη0+9.67)[-1+(1+p0p)^z]}；(35)

式中，α为粘压系数，z取0.68，p0为压力系数，取5.1×10^-9；

(7)密压方程：

(8)载荷方程：

其中w的表达式如式(25)所示；

(9)摩擦系数方程：

采用ree-eyring型非牛顿流体来对圆柱滚子轴承接触区域的摩擦力进行计算，进而得到摩擦系数；

ree-eyring模型本构方程是：

式中，τ0为润滑油的极限剪切力，ha为滚子弹流接触区域的平均油膜厚度，η为接触区的粘度参数，u和v分别表示滚动方向的吸卷速度和润滑油的端泄速度；

所以，可以推导得出摩擦系数的计算公式为：

其中w的表达式如式(25)所示；

进一步，所述步骤2中，将有限长线接触弹流润滑理论的一系列方程进行量纲一化，其有限长线接触模型如附图9所示，

p为无量纲油膜压力，ph为赫兹接触压力，b1为赫兹修正后的滚动体接触宽度。

(1)无量纲reynolds方程：

(2)无量纲reynolds方程的边界条件：

(3)无量纲油膜厚度：

(4)无量纲弹性变形方程：

(5)无量纲粗糙度方程：

式中，为无量纲粗糙度幅值，为无量纲粗糙度波长，

(6)无量纲粘压方程：

η^*＝exp{(lnη0+9.67)[-1+(1+p0php)^z]}；(46)

式中，α为粘压系数，z通常取0.68，p0为压力系数，通常取5.1×10^-9。

(7)密压公式：

(8)载荷方程：

式中，为圆柱滚子母线无量纲长度，l为母线长度。

(9)摩擦系数方程：

进一步，所述步骤3中，用fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序，将量纲一化的有限长线接触reynolds方程(41)等式左边中心差分，等式右面向后差分进行差分格式的离散得到的reynolds离散形式的差分方程：

整理后得，

式中，

油膜厚度离散方程可写成：

线接触弹性变形的离散方程：

载荷离散方程：

并根据上述离散方程，用fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序。

更进一步，所述有限长线接触弹流润滑计算程序包括以下步骤，如附图10所示，

step1，输入轴承内径、外径、滚动体直径、长度、载荷值和接触宽度值；设定初始油膜厚度h0和初始油膜压力p0；

step2，设定粗糙度幅值、波长和纹理角度；

step3，计算粘度、密度、弹性变形量v和膜厚h；

step4,采用有限差分法求解reynolds方程计算出新的油膜厚度h和油膜压力p；

step5，对比迭代前后压力相对误差，判断是否收敛；

step6，对比迭代前后荷载相对误差，判断是否收敛；

step7，计算量纲化后的油膜厚度和压力值以及弹流接触区域的平均有膜厚度；

step8，计算摩擦系数并输出计算结果。

进一步，当所述step5为不收敛时，修正初始油膜压力p0，并重新进行step3—step8；

当所述step6为不收敛时，修正初始油膜厚度h0，并重新进行step3—step8。

进一步，所述步骤4中，将轴承边界元法和有限长线接触弹流润滑理论耦合，建立滚动轴承弹流润滑条件下的轴承边界元法，其具体的方法为，

将弹流润滑理论计算出的油膜厚度(32)代入间隙值的计算公式(28)，把摩擦系数方程(39)代入边界积分方程(17)和(18)，并假定板单元固接于轴承内圈上，轴承内圈与轧辊固定，轴承外圈与轴承座固定，轴承系统简化为两个接触物体，耦合的离散边界积分方程如下：

其中，间隙δm表达式如式(28)所示；摩擦系数公式(39)中载荷w由式(25)、(22)和(23)确定。k为第k个物体，k＝1,2；n^f为非接触单元数；n^bi为k物体轴承边界单元i的数量；n^bii为k物体轴承边界单元ii的数量；为预接触区单元数量，为黏着区单元数量，为滑移区单元数量，为处于黏着状态的轴承边界单元i的数量，为处于滑移状态的轴承边界单元i的数量，为处于黏着状态的轴承边界单元ii的数量，为处于滑移状态的轴承边界单元ii的数量。

进一步，所述步骤5中，所述编写滚动轴承弹流润滑条件下的轴承边界元法的fortran计算程序包括以下步骤，如附图11所示，

step1,读入轧辊和轴承座的节点坐标、单元组成信息及边界条件；

step2,假设接触半宽b0,初始摩擦系数μ0和初始间隙δ0；

step3,计算所有物体的积分系数；

step4,对系数矩阵方程进行高斯消去，得到接触区矩阵方程；

step5,对接触区系数矩阵耦合得到耦合区的系数矩阵；

step6,用边界元法计算接触区的面力和位移；

step7,对接触区进行接触状态和摩擦状态判断；

step8,满足接触状态收敛准则和摩擦状态收敛准则；

step9，修正滚动体的接触宽度b1，并计算所有滚动体上的载荷w；

step10,判断是否满足|b1-b0|≤ε；ε为设定的精度值；

step11,判断是否为第一次迭代计算。

进一步，在所述step10中，

当满足时，执行step11；

不满足时，调用弹流润滑程序，利用接触宽度b1和载荷w计算油膜厚度平均值ha和摩擦系数μ；并假设接触半宽b1,摩擦系数μ，间隙δ0+ha-ra，重新进行step4-step10。

进一步，在所述step11中，

当满足时，调用弹流润滑程序，利用接触宽度b1和载荷w计算油膜厚度平均值ha和摩擦系数μ；并假设接触半宽b1,摩擦系数μ，间隙δ0+ha-ra，重新进行step4-step10；

不满足时，将接触区未知量代入系数矩阵方程计算非接触区的未知量。

所述步骤6中，建立轧辊和轴承座的表面单元离散模型，并进行数据的前处理；利用marc有限元软件，按照物体的实际尺寸建立表面单元模型并导出各物体的节点信息和单元组成。由于滚动轴承弹流润滑条件下的轴承边界元法计算程序采用凝缩解法，要求输入数据所有物体的接触节点和单元排在后面，且一一对应。需要对marc导出的数据进行前处理。

进一步，步骤7.将前处理好的轧辊和轴承座的节点坐标和单元组成信息代入编制的滚动轴承弹流润滑条件下的轴承边界元法的fortran计算程序进行计算；计算结束后输出数据格式的结果。

进一步，步骤8.对程序计算的数据进行后处理，得到滚动轴承的载荷和压力分布数据。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。