CN108984933B - 弹流润滑条件下计算滚动轴承载荷和压力的边界元法 - Google Patents

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

Info

Publication number
CN108984933B
CN108984933B CN201810829248.7A CN201810829248A CN108984933B CN 108984933 B CN108984933 B CN 108984933B CN 201810829248 A CN201810829248 A CN 201810829248A CN 108984933 B CN108984933 B CN 108984933B
Authority
CN
China
Prior art keywords
bearing
contact
boundary
rolling
element method
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810829248.7A
Other languages
English (en)
Other versions
CN108984933A (zh
Inventor
杨霞
李亚超
桂海莲
刘光明
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Taiyuan University of Science and Technology
Original Assignee
Taiyuan University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Taiyuan University of Science and Technology filed Critical Taiyuan University of Science and Technology
Priority to CN201810829248.7A priority Critical patent/CN108984933B/zh
Publication of CN108984933A publication Critical patent/CN108984933A/zh
Application granted granted Critical
Publication of CN108984933B publication Critical patent/CN108984933B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Rolling Contact Bearings (AREA)

Abstract

本发明属于轧机多列滚动轴承数值分析领域，具体涉及一种弹流润滑条件下计算滚动轴承载荷和压力的边界元法，包括在三维有摩擦弹性接触边界元法基础上引入轴承边界单元、接触宽度赫兹修正理论和板单元，建立轴承边界元法；根据有限长线接触弹流润滑理论，推导摩擦系数方程；将轴承边界元法和有限长线接触弹流润滑理论耦合，建立弹流润滑条件下的轴承边界元法；编写滚动轴承弹流润滑条件下的轴承边界元法的Fortran计算程序等步骤；该方法能够分析多列滚动轴承、轧辊和轴承座全模型下的载荷和压力分布；能够通过改变粗糙度影响因素，分析不同润滑条件下轴承的载荷和压力的分布；只需对轴承、轧辊和轴承座进行表面单元的划分，单元划分少，计算精度高。

Description

弹流润滑条件下计算滚动轴承载荷和压力的边界元法
技术领域
本发明属于轧机多列滚动轴承数值分析领域，具体涉及一种弹流润滑条件下计算滚动轴承载荷和压力的边界元法。
背景技术
轧机多列滚动轴承作为轧机的关键部件，在工作中承受着较大的径向负荷、轴向负荷和热负荷，致使多列滚动轴承极易产生偏载，使轴承的寿命急剧下降，并频繁发生异常烧损和大面积疲劳剥落等事故，严重制约了生产效率的提高。其润滑性能好坏对轧机滚动轴承的寿命、振动等性能有着很大的影响，在润滑良好时，轴承的载荷及其分布成为影响其运行行为的主要因素。因此对轧机多列滚动轴承考虑弹流润滑情况下的载荷分布进行研究，研究具有重大意义。
对轧机多列滚动轴承载荷分布的研究中，常采用的方法是有限元法，但基本上都没有考虑滚子弹流润滑对轴承载荷分布的影响，由于有限元法划分单元较多，在高度非线性问题下，有限元计算结果并不理想；
三维有摩擦弹性接触边界元法部分，参见【杨霞，轧机四列圆锥滚子轴承的边界元法理论与实验研究[D]，燕山大学，2011】。
发明内容
鉴于现有技术的上述缺陷，本发明所要解决的技术问题是提供一种弹流润滑条件下计算滚动轴承载荷和压力的边界元法，以满足计算要求并解决现有方法的不足。
为实现上述目的，本发明提供了一种弹流润滑条件下计算滚动轴承载荷和压力的边界元法，包括以下步骤，
步骤1.在三维有摩擦弹性接触边界元法基础上引入轴承边界单元、接触宽度赫兹修正理论和板单元，建立轴承边界元法；
步骤2.根据有限长线接触弹流润滑理论，推导摩擦系数方程，并对基本方程进行量纲一化；
步骤3.用Fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序；
步骤4.将轴承边界元法和有限长线接触弹流润滑理论耦合，建立滚动轴承弹流润滑条件下的轴承边界元法；
步骤5.编写滚动轴承弹流润滑条件下的轴承边界元法的Fortran计算程序；
步骤6.建立轧辊、轴承内圈、轴承外圈和轴承座的表面单元离散模型，并进行数据的前处理；
步骤7.将前处理好的轧辊、轴承内圈、轴承外圈和轴承座的节点坐标和单元组成信息代入编制的考动轴承弹流润滑条件下的轴承边界元法的Fortran计算程序进行计算；
步骤8.对程序计算的数据进行后处理，得到滚动轴承的载荷和压力分布数据。
进一步，所述步骤1中的三维有摩擦弹性接触边界元法为，
(1)增量形式的边界积分方程
考虑两个线弹性体相互接触，A物体的边界记为ΓA，B物体的边界记为ΓB，接触区边界记为ΓC，位移已知边界记为Γu，面力已知边界记为Γt；
采用等增量加载方法，假设总的增量步为n步，当第m步增量加载时，用增量表示的边界积分方程为：
Figure GDA0003586938950000021
其中，k表示A，B两个物体，X为源点，Y为场点；Cij(X)是与边界点X处几何形状有关的函数。如果边界点X处是光滑的，也就是说X点处的外法线矢量是连续的，那么
Figure GDA0003586938950000022
分别表示k物体在第m步增量加载时j方向的位移和面力；积分核函
Figure GDA0003586938950000023
和
Figure GDA0003586938950000024
分别是弹性问题位移和面力的基本解；
下面是基本解函数的表达式：
Figure GDA0003586938950000025
Figure GDA0003586938950000026
式中，i,j,k,l＝1，2，3，ri＝yi-xi，
Figure GDA0003586938950000027
其中δij为Kronecker函数：
Figure GDA0003586938950000028
总的量用增量的和来计算，因此第m步增量加载后总的位移和面力分别为：
Figure GDA0003586938950000029
Figure GDA00035869389500000210
其中，
Figure GDA00035869389500000211
分别为第m步增量加载后j方向总的位移和面力；
Figure GDA00035869389500000212
分别为第m步增量加载时j方向的位移和面力；
(2)离散形式的边界积分方程
由于工程问题中通常会遇到几何形状及边界条件都比较复杂的问题，因此，需要采用离散技术，进行数值求解。本发明采用4节点等参单元对边界进行离散。
对于四边形4节点线性单元，其插值函数为：
Nβ(ξ1,ξ2)＝(1+rβξ1)(1+sβξ2)/4β＝1,2,3,4； (7)
其中，rβ为第β个节点ξ1方向的局部坐标分量；sβ为第β个节点ξ2方向的局部坐标分量。
单元内部任意一点的整体坐标、位移和面力分别由单元节点上的坐标、位移和面力来描述，设
Figure GDA00035869389500000213
Figure GDA0003586938950000031
分别表示第ne单元中的β节点的坐标、位移和面力分量，则：
Figure GDA0003586938950000032
其中，q为单元的节点个数；Nβ为形函数，即插值函数。
在边界元法中，将k物体的整个边界划分为Nk个单元，然后将式(8)代入到增量边界积分方程式(1)中，得到离散形式的增量边界积分方程：
Figure GDA0003586938950000033
其中，Nk是k物体离散单元总数；
Figure GDA0003586938950000034
分别为k物体在第m步增量加载中第n单元上的β节点的位移和面力；Nβ(ξ1,ξ2)为单元上的插值函数；G(ξ1,ξ2)为单元的雅可比变换系数。
(3)建立局部坐标系
本发明在求解弹性接触问题的边界积分方程时采用凝缩解法，离散模型的接触模式采用点对点接触模式。由于未知量数目多于方程数目，为了补充方程，从而可以求得未知量，故需在接触区ΓC上建立局部坐标系(ξ1,ξ2,ξ3)。
对于接触区为曲面的问题，在接触区上，A和B物体相对应的节点的曲率不同，为了保证计算结果的精度，只在A物体的接触区
Figure GDA00035869389500000310
上建立局部坐标系，并且保证局部坐标系的ξ3方向为A物体接触边界的外法线方向，ξ1和ξ2方向只需要满足右手定则即可，
局部坐标系(ξ1,ξ2,ξ3)相对于整体坐标系(x1,x2,x3)下的方向余弦为αlj，
在接触区ΓC上，节点的位移和面力在两坐标系下的关系如下式：
Figure GDA0003586938950000035
其中，
Figure GDA0003586938950000036
分别为k物体i节点j方向的位移和面力；
Figure GDA0003586938950000037
分别为k物体i节点在局部坐标系的ξl方向的位移和面力；
Figure GDA0003586938950000038
为A物体i节点在局部坐标系下ξl与整体坐标系下xj的方向余弦。
在接触区建立局部坐标系以后，增量边界积分方程式(9)可表示为：
Figure GDA0003586938950000039
其中，αlj为局部坐标系ξl在整体坐标系xj上的方向余弦；Nkc、Nkf分别为k物体的接触区单元数和非接触区单元数；
Figure GDA0003586938950000041
分别为k物体第m步增量加载时第n单元的β节点在局部坐标系下的ξl方向的增量位移。
对于给定的两接触物体A和B，令YA和YB为相互接触的节点对，初始间隙为δ0，第m步增量加载状态的起始间隙为δm，则：
Figure GDA0003586938950000042
其中，
Figure GDA0003586938950000043
为k物体第i步增量加载时在局部坐标系下的ξ3方向上的位移。
(4)接触节点状态
更进一步，所述接触区ΓC上的接触节点对YA和YB在第m步增量加载时，可能处于分离状态、粘着状态和滑移状态三种状态之一；
当处于分离状态时，YA和YB暂时没有接触，但随着外载荷的增大或其它状态的改变，可能会进入接触状态，处于分离状态的节点满足关系式：
Figure GDA0003586938950000044
当处于粘着状态时，YA和YB已经接触，但没有发生相对滑动，处于粘着状态的节点满足关系式：
Figure GDA0003586938950000045
当处于滑移状态时，YA和YB已经接触，且发生了相对滑动，处于滑移状态的节点满足关系式：
Figure GDA0003586938950000046
当摩擦系数μ一定时，如果接触区上的某一点从第m-1步增量的粘着状态变为第m步增量的滑移状态，且YA点相对于YB点的滑移方向与ξ1轴的夹角为φk，得到滑移角之间的关系；
Figure GDA0003586938950000047
(5)耦合的增量边界积分方程
将上述接触状态关系式(13)～(15)和滑移角关系式(16)代入到增量边界积分方程(1)中，分别得到耦合的物体A的边界积分方程：
Figure GDA0003586938950000051
同理，将滑移角的关系式代入增量积分方程式(1)中，得到物体B的简化的增量边界积分方程：
Figure GDA0003586938950000052
进一步，所述步骤1中的轴承边界单元理论是假设轴承内圈有n个滚动体，则沿圆周方向轴承内圈被划分为2n个单元；同理，对于轴承外圈按一样的模式进行划分；在轴承内圈的离散模型中任取一个轴承边界单元i，此单元一定处于接触状态或者非接触状态；若此单元处于接触状态，则其上的接触面力等于滚动体上的接触载荷，因而在该单元上具有位移连续但面力不连续的特点；
此时，一个轴承边界单元分为两个子单元，子单元
Figure GDA0003586938950000053
上作用有连续面力，
Figure GDA0003586938950000054
上面力为零；假定
Figure GDA0003586938950000055
子单元上法向面力沿宽度方向呈抛物线分布，在长度方向上呈线性分布；
子单元
Figure GDA0003586938950000056
的面力增量形式表示如下：
Figure GDA0003586938950000057
其中，形函数表达式如下，
Figure GDA0003586938950000058
其中，
Figure GDA0003586938950000059
为i单元为轴承边界单元Ι时，第m步增量加载时子单元
Figure GDA00035869389500000510
在法线方向上的面力；
Figure GDA00035869389500000511
为i单元为轴承边界单元ΙΙ时，第m步增量加载时子单元
Figure GDA00035869389500000512
在法线方向上的面力；
Figure GDA00035869389500000513
为i单元β节点在第m步增量加载时法线方向上的力；
Figure GDA00035869389500000514
分别为轴承边界单元Ι和轴承边界单元ΙΙ的形函数；
Figure GDA00035869389500000515
是轴承边界单元的子单元
Figure GDA00035869389500000516
的内侧边在ξ2方向的坐标值
Figure GDA00035869389500000517
为接触区的半宽度；所述
Figure GDA00035869389500000518
划分为至少两个微单元，并使微单元的长宽比小于3，且保证奇异点所在的微单元的长宽度比等于1；
所述微单元进行坐标变换如下：
Figure GDA0003586938950000061
所述微单元区域为ξa≤ξ1≤ξb，ηa≤ξ2≤ηb，采用高斯积分时，应在区域-1≤ξ≤1，-1≤η≤1下进行；
进一步，所述步骤1中的接触宽度赫兹修正理论是，滚动轴承系统的高度非线性，在计算过程中需要初步假设滚动体与轴承外圈的接触宽度，然后再用赫兹接触理论对该接触宽度进行修正；对于轴承边界单元i，经过计算得到各个接触节点的压力值为ti,β，β表示i单元的节点编号。若此单元为轴承边界单元Ι，设子单元
Figure GDA0003586938950000062
上的合力为
Figure GDA0003586938950000063
需要逐个在微单元上进行积分求和，
计算
Figure GDA0003586938950000064
的离散形式如下，
Figure GDA0003586938950000065
若为轴承边界单元ΙΙ，则子单元
Figure GDA0003586938950000066
上的合力为
Figure GDA0003586938950000067
计算公式如下，
Figure GDA0003586938950000068
其中，n为微单元个数；l为高斯积分点数；γ1、γ2为第m个微单元对应的坐标转换系数；
Figure GDA0003586938950000069
为第r个积分点对应的权系数；ξ1、ξ2为高斯积分点坐标；
Figure GDA00035869389500000610
分别为轴承边界单元Ι、ΙΙ的形函数；G为雅可比值；
对于圆锥或圆柱滚子轴承，利用赫兹接触理论，滚动体与内、外圈接触的接触半宽度按下式计算，
Figure GDA00035869389500000611
其中，
Figure GDA00035869389500000612
其中，bi，bo分别为滚动体与轴承内、外圈的接触半宽度；Ri为内圈外径；Ro为外圈内径；R为滚动体半径；w为单位长度上的载荷；D'为滚动体有效长度。
进一步，所述步骤1中的板单元理论是，滚动体在滚动轴承中对接触状态的影响只关系到径向位移，而其余方向位移对接触状态均无影响；此时，两个圆柱体接触，弹性趋变形量为：
Figure GDA00035869389500000613
而对于圆柱或圆锥滚子轴承，滚动体的本构方程为：
Figure GDA0003586938950000071
其中，R为滚动体半径；Ri为内圈半径；Ro为外圈半径；D为是滚动体长度；E、Ei和Eo分别为滚动体、内圈和外圈的弹性模量；υ、υi和υo分别为滚动体、内圈和外圈的泊松比；
Figure GDA0003586938950000072
分别为点
Figure GDA0003586938950000073
和
Figure GDA0003586938950000074
在局部坐标系下ξ3方向的位移；
Figure GDA0003586938950000075
分别为点
Figure GDA0003586938950000076
和
Figure GDA0003586938950000077
在局部坐标系下ξ3方向的力，bi和bo值如式(24)所示；
此时，考虑滚动体弹性变形、润滑油膜厚度和粗糙度幅值后，第m步增量加载时的间隙为，
Figure GDA0003586938950000078
其中，δ0为初始间隙，ham-1为第m-1步增量加载时的平均油膜厚度值，其初始值为0；Ra为粗糙度幅值；
进一步，所述步骤1中的轴承边界元法，将轴承边界单元、接触宽度赫兹修正理论和板单元引入到三维有摩擦弹性接触边界元法中，并假定板单元固接于轴承内圈上，离散边界以后，得到轴承边界元法离散的边界积分方程如下：
Figure GDA0003586938950000079
其中，k为第k个物体；Nf为非接触单元数；NbΙ为轴承边界单元Ι的数量；NbΙΙ为轴承边界单元ΙΙ的数量；Nc为接触单元总数。
进一步，所述步骤2中，根据有限长线接触弹流润滑理论，其基本方程如下：
(1)Reynolds方程：
Figure GDA0003586938950000081
(2)Reynolds方程的边界条件：
Figure GDA0003586938950000082
(3)油膜厚度方程：
Figure GDA0003586938950000083
式中，h0为刚体油膜厚度，R为当量曲率半径，
Figure GDA0003586938950000084
“-”为滚动体与内圈，“+”代表滚动体与外圈；
(4)弹性变形方程：
Figure GDA0003586938950000085
其中，E′为综合弹性模量，
Figure GDA0003586938950000086
E1和E2分别代表两接触曲面的弹性模量，υ1和υ2代表两接触材料的泊松比；
(5)粗糙度方程：
Figure GDA0003586938950000087
式中，Ra表示粗糙度波峰幅值，l′表示粗糙度波长，β′表示粗糙度纹理方向；
(6)粘压方程：
η＝η0exp{(lnη0+9.67)[-1+(1+p0p)z]}； (35)
式中，
Figure GDA0003586938950000088
α为粘压系数，z取0.68，p0为压力系数，取5.1×10-9；
(7)密压方程：
Figure GDA0003586938950000089
(8)载荷方程：
Figure GDA00035869389500000810
其中w的表达式如式(25)所示；
(9)摩擦系数方程：
采用Ree-Eyring型非牛顿流体来对圆柱滚子轴承接触区域的摩擦力进行计算，进而得到摩擦系数；
Ree-Eyring模型本构方程是：
Figure GDA0003586938950000091
式中，τ0为润滑油的极限剪切力，ha为滚子弹流接触区域的平均油膜厚度，η为接触区的粘度参数，u和v分别表示滚动方向的吸卷速度和润滑油的端泄速度；
所以，可以推导得出摩擦系数的计算公式为：
Figure GDA0003586938950000092
其中w的表达式如式(25)所示；
进一步，所述步骤2中，将有限长线接触弹流润滑理论的一系列方程进行量纲一化，
Figure GDA0003586938950000093
P为无量纲油膜压力，pH为赫兹接触压力，b1为赫兹修正后的滚动体接触宽度。
(1)无量纲Reynolds方程：
Figure GDA0003586938950000094
(2)无量纲Reynolds方程的边界条件：
Figure GDA0003586938950000095
(3)无量纲油膜厚度：
Figure GDA0003586938950000096
(4)无量纲弹性变形方程：
Figure GDA0003586938950000097
(5)无量纲粗糙度方程：
Figure GDA0003586938950000098
式中，
Figure GDA0003586938950000099
为无量纲粗糙度幅值，
Figure GDA00035869389500000910
为无量纲粗糙度波长，
Figure GDA00035869389500000911
(6)无量纲粘压方程：
Figure GDA0003586938950000101
式中，
Figure GDA0003586938950000102
α为粘压系数，z通常取0.68，p0为压力系数，通常取5.1×10-9。
(7)密压公式：
Figure GDA0003586938950000103
(8)载荷方程：
Figure GDA0003586938950000104
式中，
Figure GDA0003586938950000105
为圆柱滚子母线无量纲长度，
Figure GDA0003586938950000106
L为母线长度。
(9)摩擦系数方程：
Figure GDA0003586938950000107
进一步，所述步骤3中，用Fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序，将量纲一化的有限长线接触Reynolds方程(41)等式左边中心差分，等式右面向后差分进行差分格式的离散得到的Reynolds离散形式的差分方程：
Figure GDA0003586938950000108
整理后得，
Figure GDA0003586938950000109
式中，
Figure GDA00035869389500001010
油膜厚度离散方程可写成：
Figure GDA00035869389500001011
线接触弹性变形的离散方程：
Figure GDA00035869389500001012
载荷离散方程：
Figure GDA0003586938950000111
并根据上述离散方程，用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)，并假定板单元固接于轴承内圈上，轴承内圈与轧辊固定，轴承外圈与轴承座固定，轴承系统简化为两个接触物体，耦合的离散边界积分方程如下：
Figure GDA0003586938950000121
其中，间隙δm表达式如式(28)所示；摩擦系数公式(39)中载荷w由式(25)、(22)和(23)确定。k为第k个物体，k＝1,2；Nf为非接触单元数；NbΙ为k物体轴承边界单元Ι的数量；Nb ΙΙ为k物体轴承边界单元ΙΙ的数量；
Figure GDA0003586938950000122
为预接触区单元数量，
Figure GDA0003586938950000123
为黏着区单元数量，
Figure GDA0003586938950000124
为滑移区单元数量，
Figure GDA0003586938950000125
为处于黏着状态的轴承边界单元Ι的数量，
Figure GDA0003586938950000126
为处于滑移状态的轴承边界单元Ι的数量，
Figure GDA0003586938950000127
为处于黏着状态的轴承边界单元ΙΙ的数量，
Figure GDA0003586938950000128
为处于滑移状态的轴承边界单元ΙΙ的数量。
进一步，所述步骤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步增量加载时，用增量表示的边界积分方程为：
Figure GDA0003586938950000141
其中，k表示A，B两个物体，X为源点，Y为场点；Cij(X)是与边界点X处几何形状有关的函数。如果边界点X处是光滑的，也就是说X点处的外法线矢量是连续的，那么Cij(X)＝δij/2。
Figure GDA0003586938950000151
分别表示k物体在第m步增量加载时j方向的位移和面力；积分核函
Figure GDA0003586938950000152
和
Figure GDA0003586938950000153
分别是弹性问题位移和面力的基本解；
下面是基本解函数的表达式：
Figure GDA0003586938950000154
Figure GDA0003586938950000155
式中，i,j,k,l＝1，2，3，ri＝yi-xi，
Figure GDA0003586938950000156
其中δij为Kronecker函数：
Figure GDA0003586938950000157
总的量用增量的和来计算，因此第m步增量加载后总的位移和面力分别为：
Figure GDA0003586938950000158
Figure GDA0003586938950000159
其中，
Figure GDA00035869389500001510
分别为第m步增量加载后j方向总的位移和面力；
Figure GDA00035869389500001511
分别为第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方向的局部坐标分量。
单元内部任意一点的整体坐标、位移和面力分别由单元节点上的坐标、位移和面力来描述，设
Figure GDA00035869389500001512
Figure GDA00035869389500001513
分别表示第ne单元中的β节点的坐标、位移和面力分量，则：
Figure GDA00035869389500001514
其中，q为单元的节点个数；Nβ为形函数，即插值函数。
在边界元法中，将k物体的整个边界划分为Nk个单元，然后将式(8)代入到增量边界积分方程式(1)中，得到离散形式的增量边界积分方程：
Figure GDA00035869389500001515
其中，Nk是k物体离散单元总数；
Figure GDA0003586938950000161
分别为k物体在第m步增量加载中第n单元上的β节点的位移和面力；Nβ(ξ1,ξ2)为单元上的插值函数；G(ξ1,ξ2)为单元的雅可比变换系数。
(3)建立局部坐标系
本发明在求解弹性接触问题的边界积分方程时采用凝缩解法，离散模型的接触模式采用点对点接触模式。由于未知量数目多于方程数目，为了补充方程，从而可以求得未知量，故需在接触区ΓC上建立局部坐标系(ξ1,ξ2,ξ3)。
对于接触区为曲面的问题，在接触区上，A和B物体相对应的节点的曲率不同，为了保证计算结果的精度，只在A物体的接触区
Figure GDA00035869389500001610
上建立局部坐标系，并且保证局部坐标系的ξ3方向为A物体接触边界的外法线方向，ξ1和ξ2方向只需要满足右手定则即可，如图3所示为接触区上建立局部坐标系的示意图。局部坐标系(ξ1,ξ2,ξ3)相对于整体坐标系(x1,x2,x3)下的方向余弦为αlj，
在接触区ΓC上，节点的位移和面力在两坐标系下的关系如下式：
Figure GDA0003586938950000162
其中，
Figure GDA0003586938950000163
分别为k物体i节点j方向的位移和面力；
Figure GDA0003586938950000164
分别为k物体i节点在局部坐标系的ξl方向的位移和面力；
Figure GDA0003586938950000165
为A物体i节点在局部坐标系下ξl与整体坐标系下xj的方向余弦。
在接触区建立局部坐标系以后，增量边界积分方程式(9)可表示为：
Figure GDA0003586938950000166
其中，αlj为局部坐标系ξl在整体坐标系xj上的方向余弦；Nkc、Nkf分别为k物体的接触区单元数和非接触区单元数；
Figure GDA0003586938950000167
分别为k物体第m步增量加载时第n单元的β节点在局部坐标系下的ξl方向的增量位移。
对于给定的两接触物体A和B，令YA和YB为相互接触的节点对，初始间隙为δ0，第m步增量加载状态的起始间隙为δm，则：
Figure GDA0003586938950000168
其中，
Figure GDA0003586938950000169
为k物体第i步增量加载时在局部坐标系下的ξ3方向上的位移。
(4)接触节点状态
更进一步，所述接触区ΓC上的接触节点对YA和YB在第m步增量加载时，可能处于分离状态、粘着状态和滑移状态三种状态之一；
当处于分离状态时，YA和YB暂时没有接触，但随着外载荷的增大或其它状态的改变，可能会进入接触状态，处于分离状态的节点满足关系式：
Figure GDA0003586938950000171
当处于粘着状态时，YA和YB已经接触，但没有发生相对滑动，处于粘着状态的节点满足关系式：
Figure GDA0003586938950000172
当处于滑移状态时，YA和YB已经接触，且发生了相对滑动，处于滑移状态的节点满足关系式：
Figure GDA0003586938950000173
当摩擦系数μ一定时，如果接触区上的某一点从第m-1步增量的粘着状态变为第m步增量的滑移状态，且YA点相对于YB点的滑移方向与ξ1轴的夹角为φk，如图4所示，得到滑移角之间的关系；
Figure GDA0003586938950000174
(5)耦合的增量边界积分方程
将上述接触状态关系式(13)～(15)和滑移角关系式(16)代入到增量边界积分方程(1)中，分别得到耦合的物体A的边界积分方程：
Figure GDA0003586938950000175
同理，将滑移角的关系式代入增量积分方程式(1)中，得到物体B的简化的增量边界积分方程：
Figure GDA0003586938950000181
进一步，所述步骤1中的轴承边界单元理论是以轴承内圈为例，其离散模型及滚动体的位置关系如图5所示。假设轴承内圈有n个滚动体，则沿圆周方向轴承内圈被划分为2n个单元；同理，对于轴承外圈按一样的模式进行划分；在轴承内圈的离散模型中任取一个轴承边界单元i，此单元一定处于接触状态或者非接触状态；若此单元处于接触状态，则其上的接触面力等于滚动体上的接触载荷，因而在该单元上具有位移连续但面力不连续的特点；
此时，一个轴承边界单元分为两个子单元，子单元
Figure GDA0003586938950000182
上作用有连续面力，
Figure GDA0003586938950000183
上面力为零；假定
Figure GDA0003586938950000184
子单元上法向面力沿宽度方向呈抛物线分布，在长度方向上呈线性分布,如图6所示；
子单元
Figure GDA0003586938950000185
的面力增量形式表示如下：
Figure GDA0003586938950000186
其中，形函数表达式如下，
Figure GDA0003586938950000187
其中，
Figure GDA0003586938950000188
为i单元为轴承边界单元Ι时，第m步增量加载时子单元
Figure GDA00035869389500001818
在法线方向上的面力；
Figure GDA0003586938950000189
为i单元为轴承边界单元ΙΙ时，第m步增量加载时子单元
Figure GDA00035869389500001810
在法线方向上的面力；
Figure GDA00035869389500001811
为i单元β节点在第m步增量加载时法线方向上的力；
Figure GDA00035869389500001812
分别为轴承边界单元Ι和轴承边界单元ΙΙ的形函数；
Figure GDA00035869389500001813
是轴承边界单元的子单元
Figure GDA00035869389500001814
的内侧边在ξ2方向的坐标值
Figure GDA00035869389500001815
为接触区的半宽度；所述
Figure GDA00035869389500001816
划分为至少两个微单元，并使微单元的长宽比小于3，且保证奇异点所在的微单元的长宽度比等于1，如图7所示；
所述微单元进行坐标变换如下：
Figure GDA00035869389500001817
所述微单元区域为ξa≤ξ1≤ξb，ηa≤ξ2≤ηb，采用高斯积分时，应在区域-1≤ξ≤1，-1≤η≤1下进行；
进一步，所述步骤1中的接触宽度赫兹修正理论是，滚动轴承系统的高度非线性，在计算过程中需要初步假设滚动体与轴承外圈的接触宽度，然后再用赫兹接触理论对该接触宽度进行修正；对于轴承边界单元i，经过计算得到各个接触节点的压力值为ti,β，β表示i单元的节点编号。若此单元为轴承边界单元Ι，设子单元
Figure GDA0003586938950000191
上的合力为
Figure GDA0003586938950000192
需要逐个在微单元上进行积分求和，
计算
Figure GDA0003586938950000193
的离散形式如下，
Figure GDA0003586938950000194
若为轴承边界单元ΙΙ，则子单元
Figure GDA0003586938950000195
上的合力为
Figure GDA0003586938950000196
计算公式如下，
Figure GDA0003586938950000197
其中，n为微单元个数；l为高斯积分点数；γ1、γ2为第m个微单元对应的坐标转换系数；
Figure GDA0003586938950000198
为第r个积分点对应的权系数；ξ1、ξ2为高斯积分点坐标；
Figure GDA0003586938950000199
分别为轴承边界单元Ι、ΙΙ的形函数；G为雅可比值；
对于圆锥或圆柱滚子轴承，利用赫兹接触理论，滚动体与内、外圈接触的接触半宽度按下式计算，
Figure GDA00035869389500001910
其中，
Figure GDA00035869389500001911
其中，bi，bo分别为滚动体与轴承内、外圈的接触半宽度；Ri为内圈外径；Ro为外圈内径；R为滚动体半径；w为单位长度上的载荷；D'为滚动体有效长度。
进一步，所述步骤1中的板单元理论是，滚动体在滚动轴承中对接触状态的影响只关系到径向位移，而其余方向位移对接触状态均无影响；图8为板单元示意图，此时，两个圆柱体接触，弹性趋变形量为：
Figure GDA00035869389500001912
而对于圆柱或圆锥滚子轴承，滚动体的本构方程为：
Figure GDA00035869389500001913
其中，R为滚动体半径；Ri为内圈半径；Ro为外圈半径；D为是滚动体长度；E、Ei和Eo分别为滚动体、内圈和外圈的弹性模量；υ、υi和υo分别为滚动体、内圈和外圈的泊松比；
Figure GDA00035869389500001914
分别为点
Figure GDA0003586938950000201
和Yi R在局部坐标系下ξ3方向的位移；
Figure GDA0003586938950000202
分别为点Yi A和
Figure GDA0003586938950000203
在局部坐标系下ξ3方向的力，bi和bo值如式(24)所示；
此时，考虑滚动体弹性变形、润滑油膜厚度和粗糙度幅值后，第m步增量加载时的间隙为，
Figure GDA0003586938950000204
其中，δ0为初始间隙，ham-1为第m-1步增量加载时的平均油膜厚度值，其初始值为0；Ra为粗糙度幅值；
进一步，所述步骤1中的轴承边界元法，将轴承边界单元、接触宽度赫兹修正理论和板单元引入到三维有摩擦弹性接触边界元法中，并假定板单元固接于轴承内圈上，离散边界以后，得到轴承边界元法离散的边界积分方程如下：
Figure GDA0003586938950000205
其中，k为第k个物体；Nf为非接触单元数；NbΙ为轴承边界单元Ι的数量；NbΙΙ为轴承边界单元ΙΙ的数量；Nc为接触单元总数。
进一步，所述步骤2中，根据有限长线接触弹流润滑理论，其基本方程如下：
(1)Reynolds方程：
Figure GDA0003586938950000206
(2)Reynolds方程的边界条件：
Figure GDA0003586938950000207
(3)油膜厚度方程：
Figure GDA0003586938950000208
式中，h0为刚体油膜厚度，R为当量曲率半径，
Figure GDA0003586938950000211
“-”为滚动体与内圈，“+”代表滚动体与外圈；
(4)弹性变形方程：
Figure GDA0003586938950000212
其中，E′为综合弹性模量，
Figure GDA0003586938950000213
E1和E2分别代表两接触曲面的弹性模量，υ1和υ2代表两接触材料的泊松比；
(5)粗糙度方程：
Figure GDA0003586938950000214
式中，Ra表示粗糙度波峰幅值，l′表示粗糙度波长，β′表示粗糙度纹理方向；
(6)粘压方程：
η＝η0exp{(lnη0+9.67)[-1+(1+p0p)z]}； (35)
式中，
Figure GDA0003586938950000215
α为粘压系数，z取0.68，p0为压力系数，取5.1×10-9；
(7)密压方程：
Figure GDA0003586938950000216
(8)载荷方程：
Figure GDA0003586938950000217
其中w的表达式如式(25)所示；
(9)摩擦系数方程：
采用Ree-Eyring型非牛顿流体来对圆柱滚子轴承接触区域的摩擦力进行计算，进而得到摩擦系数；
Ree-Eyring模型本构方程是：
Figure GDA0003586938950000218
式中，τ0为润滑油的极限剪切力，ha为滚子弹流接触区域的平均油膜厚度，η为接触区的粘度参数，u和v分别表示滚动方向的吸卷速度和润滑油的端泄速度；
所以，可以推导得出摩擦系数的计算公式为：
Figure GDA0003586938950000221
其中w的表达式如式(25)所示；
进一步，所述步骤2中，将有限长线接触弹流润滑理论的一系列方程进行量纲一化，其有限长线接触模型如附图9所示，
Figure GDA0003586938950000222
P为无量纲油膜压力，pH为赫兹接触压力，b1为赫兹修正后的滚动体接触宽度。
(1)无量纲Reynolds方程：
Figure GDA0003586938950000223
(2)无量纲Reynolds方程的边界条件：
Figure GDA0003586938950000224
(3)无量纲油膜厚度：
Figure GDA0003586938950000225
(4)无量纲弹性变形方程：
Figure GDA0003586938950000226
(5)无量纲粗糙度方程：
Figure GDA0003586938950000227
式中，
Figure GDA0003586938950000228
为无量纲粗糙度幅值，
Figure GDA0003586938950000229
为无量纲粗糙度波长，
Figure GDA00035869389500002210
(6)无量纲粘压方程：
η*＝exp{(lnη0+9.67)[-1+(1+p0pHP)z]}； (46)
式中，
Figure GDA00035869389500002211
α为粘压系数，z通常取0.68，p0为压力系数，通常取5.1×10-9。
(7)密压公式：
Figure GDA0003586938950000231
(8)载荷方程：
Figure GDA0003586938950000232
式中，
Figure GDA0003586938950000233
为圆柱滚子母线无量纲长度，
Figure GDA0003586938950000234
L为母线长度。
(9)摩擦系数方程：
Figure GDA0003586938950000235
进一步，所述步骤3中，用Fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序，将量纲一化的有限长线接触Reynolds方程(41)等式左边中心差分，等式右面向后差分进行差分格式的离散得到的Reynolds离散形式的差分方程：
Figure GDA0003586938950000236
整理后得，
Figure GDA0003586938950000237
式中，
Figure GDA0003586938950000238
油膜厚度离散方程可写成：
Figure GDA0003586938950000239
线接触弹性变形的离散方程：
Figure GDA00035869389500002310
载荷离散方程：
Figure GDA00035869389500002311
并根据上述离散方程，用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)，并假定板单元固接于轴承内圈上，轴承内圈与轧辊固定，轴承外圈与轴承座固定，轴承系统简化为两个接触物体，耦合的离散边界积分方程如下：
Figure GDA0003586938950000251
其中，间隙δm表达式如式(28)所示；摩擦系数公式(39)中载荷w由式(25)、(22)和(23)确定。k为第k个物体，k＝1,2；Nf为非接触单元数；NbΙ为k物体轴承边界单元Ι的数量；Nb ΙΙ为k物体轴承边界单元ΙΙ的数量；
Figure GDA0003586938950000252
为预接触区单元数量，
Figure GDA0003586938950000253
为黏着区单元数量，
Figure GDA0003586938950000254
为滑移区单元数量，
Figure GDA0003586938950000255
为处于黏着状态的轴承边界单元Ι的数量，
Figure GDA0003586938950000256
为处于滑移状态的轴承边界单元Ι的数量，
Figure GDA0003586938950000257
为处于黏着状态的轴承边界单元ΙΙ的数量，
Figure GDA0003586938950000258
为处于滑移状态的轴承边界单元ΙΙ的数量。
进一步，所述步骤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.对程序计算的数据进行后处理，得到滚动轴承的载荷和压力分布数据。
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (10)

1.弹流润滑条件下计算滚动轴承载荷和压力的边界元法，包括以下步骤，
步骤1.在三维有摩擦弹性接触边界元法基础上引入轴承边界单元、接触宽度赫兹修正理论和板单元，建立轴承边界元法；
步骤2.根据有限长线接触弹流润滑理论，推导摩擦系数方程，并对基本方程进行量纲一化；
步骤3.用Fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序；
步骤4.将轴承边界元法和有限长线接触弹流润滑理论耦合，建立滚动轴承弹流润滑条件下的轴承边界元法；
步骤5.编写滚动轴承弹流润滑条件下的轴承边界元法的Fortran计算程序；
步骤6.建立轧辊和轴承座的表面单元离散模型，并进行数据的前处理；
步骤7.将前处理好的轧辊和轴承座的节点坐标和单元组成信息代入编制的考动轴承弹流润滑条件下的轴承边界元法的Fortran计算程序进行计算；
步骤8.对程序计算的数据进行后处理，得到滚动轴承的载荷和压力分布数据。
2.根据权利要求1所述的弹流润滑条件下计算滚动轴承载荷和压力的边界元法，其特征在于：所述步骤1中的接触宽度赫兹修正理论是，滚动轴承系统的高度非线性，在计算过程中需要初步假设滚动体与轴承外圈的接触宽度，然后再用赫兹接触理论对该接触宽度进行修正；对于轴承边界单元i，经过计算得到各个接触节点的压力值为ti,β，β表示i单元的节点编号；若此单元为轴承边界单元Ι，设子单元
Figure FDA0003586938940000011
上的合力为
Figure FDA0003586938940000012
需要逐个在微单元上进行积分求和，
计算
Figure FDA0003586938940000013
的离散形式如下，
Figure FDA0003586938940000014
若为轴承边界单元ΙΙ，则子单元
Figure FDA0003586938940000015
上的合力为Ti ΙΙ,计算公式如下，
Figure FDA0003586938940000016
其中，n为微单元个数；l为高斯积分点数；γ1、γ2为第m个微单元对应的坐标转换系数；
Figure FDA0003586938940000017
为第r个积分点对应的权系数；ξ1、ξ2为高斯积分点坐标；
Figure FDA0003586938940000018
分别为轴承边界单元Ι、ΙΙ的形函数；G为雅可比值；
对于圆锥或圆柱滚子轴承，利用赫兹接触理论，滚动体与内、外圈接触的接触宽度按下式计算，
Figure FDA0003586938940000019
其中，
Figure FDA0003586938940000021
其中，bi，bo分别为滚动体与轴承内、外圈的接触半宽度；Ei和Eo分别为内圈和外圈的弹性模量；Ri为内圈外径；Ro为外圈内径；R为滚动体半径；w为单位长度上的载荷；D'为滚动体有效长度。
3.根据权利要求1所述的弹流润滑条件下计算滚动轴承载荷和压力的边界元法，其特征在于：所述步骤1中的板单元理论是，滚动体在滚动轴承中对接触状态的影响只关系到径向位移，而其余方向位移对接触状态均无影响；此时，两个圆柱体接触，弹性趋变形量为：
Figure FDA0003586938940000022
其中，r1和r2分别代表两个相接触的圆柱体的半径；
而对于圆柱或圆锥滚子轴承，滚动体的本构方程为：
Figure FDA0003586938940000023
其中，R为滚动体半径；Ri为内圈半径；Ro为外圈半径；D为是滚动体长度；E、Ei和Eo分别为滚动体、内圈和外圈的弹性模量；υ、υi和υo分别为滚动体、内圈和外圈的泊松比；
Figure FDA0003586938940000024
分别为点
Figure FDA0003586938940000025
和Yi R在局部坐标系下ξ3方向的位移；
Figure FDA0003586938940000026
分别为点Yi A和Yo B在局部坐标系下ξ3方向的力，bi和bo值如式(24)所示；
此时，考虑滚动体弹性变形、润滑油膜厚度和粗糙度幅值后，第m步增量加载时的间隙为，
Figure FDA0003586938940000027
其中，δ0为初始间隙，ham-1为第m-1步增量加载时的平均油膜厚度值，其初始值为0；Ra为粗糙度幅值。
4.根据权利要求1所述的弹流润滑条件下计算滚动轴承载荷和压力的边界元法，其特征在于：所述步骤2中，根据有限长线接触弹流润滑理论，推导摩擦系数方程，摩擦系数方程如下：
采用Ree-Eyring型非牛顿流体来对圆柱滚子轴承接触区域的摩擦力进行计算，进而得到摩擦系数；
Ree-Eyring模型本构方程是：
Figure FDA0003586938940000028
式中，τ0为润滑油的极限剪切力，τx和τy分别为x、y两个方向的剪应力，ha为滚子弹流接触区域的平均油膜厚度，η为接触区的粘度参数，u和v分别表示滚动方向的吸卷速度和润滑油的端泄速度；
所以，可以推导得出摩擦系数的计算公式为：
Figure FDA0003586938940000031
式中，w为滚动体上所受的单位载荷。
5.根据权利要求1所述的弹流润滑条件下计算滚动轴承载荷和压力的边界元法，其特征在于：所述步骤2中，对限长线接触弹流润滑理论的基本方程进行量纲一化，其中量纲一化后的油膜压力计算公式如下：
Figure FDA0003586938940000032
其中P为无量纲油膜压力，p为油膜压力，pH为赫兹接触压力，w为滚动体上所受的单位载荷，计算公式如式(25)所示，b1为滚动体的接触宽度，如式(24)所示；
对摩擦系数方程进行量纲一化，公式如下：
Figure FDA0003586938940000033
6.根据权利要求1所述的弹流润滑条件下计算滚动轴承载荷和压力的边界元法，其特征在于：所述步骤3中，用Fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序，其中量纲一化的有限长线接触Reynolds方程，等式左边中心差分，等式右面向后差分进行差分格式的离散得到离散后的差分Reynolds方程如下：
Figure FDA0003586938940000034
整理后得，
Figure FDA0003586938940000035
式中，
Figure FDA0003586938940000036
式中：ε为设定的精度值。
7.根据权利要求1所述的弹流润滑条件下计算滚动轴承载荷和压力的边界元法，其特征在于：所述步骤3中，用Fortran语言编写计算滚动轴承弹流润滑的有限差分法计算程序包括以下步骤，
step1，输入轴承内径、外径、滚动体直径、长度、载荷值和接触宽度值；设定初始油膜厚度h0和初始油膜压力p0；
step2，设定粗糙度幅值、波长和纹理角度；
step3，计算粘度、密度、弹性变形量v和膜厚h；
step4,采用有限差分法求解Reynolds方程计算出新的油膜厚度h和油膜压力p；
step5，对比迭代前后压力相对误差，判断是否收敛；
step6，对比迭代前后荷载相对误差，判断是否收敛；
Step7，计算量纲化后的油膜厚度和压力值；
Step8，计算摩擦系数并输出计算结果；
进一步，当所述step5为不收敛时，修正初始油膜压力p0，并重新进行step3—step8；
当所述step6为不收敛时，修正初始油膜厚度h0，并重新进行step3—step8。
8.根据权利要求1所述的弹流润滑条件下计算滚动轴承载荷和压力的边界元法，其特征在于：所述步骤4中，将轴承边界元法和有限长线接触弹流润滑理论耦合，建立滚动轴承弹流润滑条件下的轴承边界元法，其具体的方法为，
将弹流润滑理论计算出的油膜厚度(32)代入间隙值的计算公式(28)，把摩擦系数方程(39)代入边界积分方程(17)和(18)，并假定板单元固接于轴承内圈上，轴承内圈与轧辊固定，轴承外圈与轴承座固定，轴承系统简化为两个接触物体，耦合的离散边界积分方程如下：
Figure FDA0003586938940000051
其中，
Figure FDA0003586938940000052
表示k物体在第m步增量加载时j方向的位移，
Figure FDA0003586938940000053
为k物体在第m步增量加载中第n单元上的β节点的位移，
Figure FDA0003586938940000054
为k物体第m步增量加载时第n单元的β节点在局部坐标系下的ξl方向的增量位移；间隙δm表达式如式(28)所示；摩擦系数公式(39)中载荷w由式(25)、(22)和(23)确定；k为第k个物体，k＝1,2；Nf为非接触单元数；NbΙ为k物体轴承边界单元Ι的数量；NbΙΙ为k物体轴承边界单元ΙΙ的数量；
Figure FDA0003586938940000055
为预接触区单元数量，
Figure FDA0003586938940000056
为黏着区单元数量，
Figure FDA0003586938940000057
为滑移区单元数量，
Figure FDA0003586938940000058
为处于黏着状态的轴承边界单元Ι的数量，
Figure FDA0003586938940000059
为处于滑移状态的轴承边界单元Ι的数量，
Figure FDA00035869389400000510
为处于黏着状态的轴承边界单元ΙΙ的数量，
Figure FDA0003586938940000061
为处于滑移状态的轴承边界单元ΙΙ的数量。
9.根据权利要求1所述的弹流润滑条件下计算滚动轴承载荷和压力的边界元法，其特征在于：所述步骤5中，编写滚动轴承弹流润滑条件下的轴承边界元法的Fortran计算程序包括以下步骤，
Step1,读入轧辊和轴承座的节点坐标、单元组成信息及边界条件；
Step2,假设接触半宽b0,初始摩擦系数μ0和初始间隙δ0；
Step3,计算所有物体的积分系数；
Step4,对系数矩阵方程进行高斯消去，得到接触区矩阵方程；
Step5,对接触区系数矩阵耦合得到耦合区的系数矩阵；
Step6,用边界元法计算接触区的面力和位移；
Step7,对接触区进行接触状态和摩擦状态判断；
Step8,满足接触状态收敛准则和摩擦状态收敛准则；
Step9，修正滚动体的接触宽度b1，并计算所有滚动体上的载荷w；
Step10,判断是否满足|b1-b0|≤ε；ε为设定的精度值；
Step11,判断是否为第一次迭代计算。
10.根据权利要求9所述的弹流润滑条件下计算滚动轴承载荷和压力的边界元法，其特征在于：
所述step10中，
当满足时，执行Step11；
不满足时，调用弹流润滑程序，利用接触宽度b1和载荷w计算油膜厚度平均值ha和摩擦系数μ；并假设接触半宽b1,摩擦系数μ，间隙δ0+ha-Ra，重新进行step4-step10；
在所述step11中，
当满足时，调用弹流润滑程序，利用接触宽度b1和载荷w计算油膜厚度平均值ha和摩擦系数μ；并假设接触半宽b1,摩擦系数μ，间隙δ0+ha-Ra，重新进行step4-step10；
不满足时，将接触区未知量代入系数矩阵方程计算非接触区的未知量。
CN201810829248.7A 2018-07-25 2018-07-25 弹流润滑条件下计算滚动轴承载荷和压力的边界元法 Active CN108984933B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810829248.7A CN108984933B (zh) 2018-07-25 2018-07-25 弹流润滑条件下计算滚动轴承载荷和压力的边界元法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810829248.7A CN108984933B (zh) 2018-07-25 2018-07-25 弹流润滑条件下计算滚动轴承载荷和压力的边界元法

Publications (2)

Publication Number Publication Date
CN108984933A CN108984933A (zh) 2018-12-11
CN108984933B true CN108984933B (zh) 2022-05-20

Family

ID=64550883

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810829248.7A Active CN108984933B (zh) 2018-07-25 2018-07-25 弹流润滑条件下计算滚动轴承载荷和压力的边界元法

Country Status (1)

Country Link
CN (1) CN108984933B (zh)

Families Citing this family (19)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109815548B (zh) * 2018-12-27 2021-01-19 西安交通大学 一种基于Garlerkin思想的流体膜压力计算方法
CN109753723B (zh) * 2019-01-02 2022-10-04 太原理工大学 一种向心滚动轴承疲劳寿命计算方法
CN109829262B (zh) * 2019-04-04 2021-10-19 哈尔滨工程大学 一种转子-轴承系统非线性动力学分析方法
CN110119576A (zh) * 2019-05-15 2019-08-13 重庆大学 一种基于adina二次开发的点接触弹流润滑仿真分析方法
CN110378018B (zh) * 2019-07-18 2023-12-26 上海理工大学 一种液体动静压球轴承的稳态承载能力的计算方法
CN110765617B (zh) * 2019-10-25 2023-07-25 常州市乾憬轴承科技有限公司 基于润滑理论的圆柱滚子轴承滚子对数修形的设计方法
CN111177898B (zh) * 2019-12-16 2024-03-29 重庆大学 基于bp神经网络的滚动轴承-转子系统耦合性能求解方法
CN111209657B (zh) * 2019-12-30 2022-04-22 南京航空航天大学 考虑液体表面张力的固体变形界面计算方法
CN111075920B (zh) * 2020-01-13 2023-02-21 重庆大学 基于fft和润滑影响的rv减速器摆线针轮残余应力求解方法
CN111209686B (zh) * 2020-01-16 2024-04-19 重庆大学 基于复合形法的滚动轴承多体润滑性能求解方法
CN112307571B (zh) * 2020-07-08 2021-06-29 重庆大学 径向推力一体式水润滑轴承及其自适应混合润滑分析方法
CN112131771B (zh) * 2020-09-18 2022-10-11 重庆长安汽车股份有限公司 一种汽车发动机的气门油封机油泄漏量的预测方法
CN112270078B (zh) * 2020-10-16 2022-02-15 西南科技大学 一种脆性材料多道次纳米划痕深度预测的方法
CN112836445B (zh) * 2021-01-18 2023-06-16 中车大连机车研究所有限公司 一种用于求解滚动轴承的动态接触混合润滑参数的方法
CN113503197B (zh) * 2021-07-12 2022-11-18 哈尔滨工程大学 一种考虑结构振动的船用凸轮-挺柱副弹流润滑分析方法
CN113779826B (zh) * 2021-08-17 2024-02-02 北京化工大学 一种动载滑动轴承瞬态与时滞变形交替耦合作用的混合润滑模型建模方法
CN114611433B (zh) * 2022-03-22 2024-02-23 郑州大学 一种耦合流态与粗糙度的动静压浮环轴承模型计算方法
CN114781158B (zh) * 2022-04-21 2023-05-05 西南交通大学 一种简化轮轨非赫兹切向接触计算的等效方法
CN116933510B (zh) * 2023-07-10 2024-06-25 哈尔滨工业大学 一种轴承打滑蹭伤失效行为预测分析方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1414194A (zh) * 2002-11-21 2003-04-30 西安理工大学 含多柔性结合部梁的拉、压、弯扭耦合边界元法
WO2007145712A2 (en) * 2006-06-07 2007-12-21 Fei Company Slider bearing for use with an apparatus comprising a vacuum chamber
CN106441196A (zh) * 2016-08-30 2017-02-22 北京理工大学 一种基于摩擦力的轴孔配合间隙测量装置及方法
CN106649982A (zh) * 2016-11-08 2017-05-10 大连工业大学 风力发电机大锥角圆锥主轴承摩擦力矩计算方法
CN107506562A (zh) * 2017-09-29 2017-12-22 西安科技大学 一种水润滑橡胶轴承双向热流固耦合计算方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1414194A (zh) * 2002-11-21 2003-04-30 西安理工大学 含多柔性结合部梁的拉、压、弯扭耦合边界元法
WO2007145712A2 (en) * 2006-06-07 2007-12-21 Fei Company Slider bearing for use with an apparatus comprising a vacuum chamber
CN106441196A (zh) * 2016-08-30 2017-02-22 北京理工大学 一种基于摩擦力的轴孔配合间隙测量装置及方法
CN106649982A (zh) * 2016-11-08 2017-05-10 大连工业大学 风力发电机大锥角圆锥主轴承摩擦力矩计算方法
CN107506562A (zh) * 2017-09-29 2017-12-22 西安科技大学 一种水润滑橡胶轴承双向热流固耦合计算方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Auxiliary Bearing Life Prediction Using Hertzian Contact Bearing Model;Sun, G;《ASME. J. Vib. Acoust》;20060430;第128卷(第02期);203–209 *
赫兹接触理论在采用边界元法分析轧机轴承载荷中的应用;肖宏 等;《中国机械工程》;20110413;第21卷(第21期);2532-2535,2540 *

Also Published As

Publication number Publication date
CN108984933A (zh) 2018-12-11

Similar Documents

Publication Publication Date Title
CN108984933B (zh) 弹流润滑条件下计算滚动轴承载荷和压力的边界元法
Jin et al. Heat generation modeling of ball bearing based on internal load distribution
CN106649982B (zh) 风力发电机大锥角圆锥主轴承摩擦力矩计算方法
Noel et al. Complete analytical expression of the stiffness matrix of angular contact ball bearings
Arghir et al. A simplified structural model of bump-type foil bearings based on contact mechanics including gaps and friction
Kania et al. A catalogue capacity of slewing bearings
CN109753723B (zh) 一种向心滚动轴承疲劳寿命计算方法
Liu et al. A novel method to model effects of natural defect on roller bearing
Lacroix et al. Four-point contact ball bearing model with deformable rings
Śpiewak Methodology for calculating the complete static carrying capacity of twin slewing bearing
CN111177898A (zh) 基于bp神经网络的滚动轴承-转子系统耦合性能求解方法
CN114580292A (zh) 一种滚动轴承的塑性变形的可靠性和灵敏度分析方法
Xing et al. A novel wear prediction model for planetary roller screw based on universal sliding distance model
Kang et al. Stiffness determination of angular-contact ball bearings by using neural network
Liu et al. An analysis of the load distribution characteristics of a cylindrical roller bearing including the component deformation and waviness
Zhang et al. An efficient multi-scale modeling method that reveals coupled effects between surface roughness and roll-stack deformation in cold sheet rolling
CN115470584B (zh) 应用于行星齿轮与滚动轴承耦合系统的动力学建模方法
CN111898216B (zh) 一种螺旋曲面的磨损预测方法
Zhilnikov et al. A method of calculating the friction moment in cageless bearings
Li et al. Dynamic characteristics of joint surface considering friction and vibration factors based on fractal theory
CN110929350B (zh) 一种考虑齿面磨损的摩擦热流量预测方法
CN107391893B (zh) 基于空心圆柱滚子接触变形量和载荷分布的轴承参数优化方法
Scurria et al. An advanced modeling technique for rolling element bearings in elastohydrodynamic field
Yang et al. The prediction of load distribution in four-row cylindrical roller bearings under the elastohydrodynamic lubrication based on boundary element method
Karpenko et al. A numerical method for analysis of extended rough wavy surfaces in contact

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
EE01 Entry into force of recordation of patent licensing contract

Application publication date: 20181211

Assignee: Shanxi advanced forming intelligent equipment Co.,Ltd.

Assignor: TAIYUAN University OF SCIENCE AND TECHNOLOGY

Contract record no.: X2023980040616

Denomination of invention: Boundary Element Method for Calculating Load and Pressure of Rolling Bearings under Elastohydrodynamic Lubrication Conditions

Granted publication date: 20220520

License type: Common License

Record date: 20230825

EE01 Entry into force of recordation of patent licensing contract