线接触滚动轴承弹塑性变形计算方法与流程

本发明属于滚动轴承技术领域，涉及线接触滚动轴承弹塑性变形计算方法。

背景技术：

线接触滚动轴承(圆柱滚子轴承和滚针轴承)在工程机械、风力发电、轨道交通及航空航天等领域有着广泛的应用，它依靠元件间的滚动接触来支撑转动零件，是现代机器中重要的动力传递部件。轴承钢材料在外载荷作用下存在一个应变极限，如果超过此极限，在卸去载荷后轴承便不可能完全恢复其原有尺寸，滚道上可能会留下压痕，而滚子上可能会出现长条状斑痕。通常，这种永久变形很小，不至于引起正常工作载荷下的轴承发生断裂失效。经验表明，如果任一点上的永久变形量限制在最大为0.0001d时(d为滚子直径)，它对轴承稳定运转的影响很小。但是，如果永久变形量很大，则会在滚道上形成凹坑，引起轴承振动，摩擦加剧。当压痕与边界润滑状态同时出现时，有可能导致轴承表面首先疲劳失效。

传统的轴承设计方法均针对光滑轴承，无法用于分析加工精度对轴承零件弹塑性变形(也即永久变形量)的影响，无法准确计算轴承在多大径向载荷下会发生0.0001d的可能导致轴承运转失稳的弹塑性变形，也即无法准确指导线接触滚动轴承径向载荷的施加。所以，传统的线接触滚动轴承弹塑性变形计算方法更多的是具有统计与经验意义，其考察因素不够全面、计算精度上不够准确。

技术实现要素：

有鉴于此，本发明的目的在于提供一种线接触滚动轴承弹塑性变形计算方法，该方法能够考虑加工精度对轴承零件弹塑性变形的影响，可准确、快速地获得计算区域内的轴承弹塑性变形，计算出的弹塑性变形结果可用于指导控制线接触滚动轴承的径向载荷在不至于使轴承发生永久变形量过大、摩擦振动加剧的合理范围之内。

为达到上述目的，本发明提供如下技术方案：

一种线接触滚动轴承弹塑性变形计算方法，包括以下步骤：

s1：用控制误差的试解法求解线接触滚动轴承各滚子承受的法向载荷qψj；

s2：对各方位角处受载的滚子，将其承受的法向载荷qψj分为nl个准静态载荷步进行加载，将第n个准静态载荷施加给所分析的滚子，n＝1,…,nl；

s3：计算含表面粗糙度hr和滚子弹塑性变形u3的接触间隙h；

s4：采用cgm(conjugategradientmethod，共轭梯度法)计算表面法向接触应力p，其中法向弹性变形的计算需用到fft；

s5：采用icm(influencecoefficientmethod，影响系数法)及fft计算次表面应力σ；

s6：根据vonmises屈服准则判断塑性区域；

s7：采用newton-raphson法计算塑性应变增量δε^p；

s8：计算塑性应变ε^p所引起的残余应力σ⁽²⁾；

s9：判断塑性应变增量δε^p是否收敛，若不收敛，则更新塑性应变增量δε^p后返回步骤s7；若收敛，则转入步骤s10；

s10：计算表面残余变形

s11：判断表面残余变形是否收敛，若不收敛，则返回步骤s3更新接触间隙h；若收敛，则转入步骤s12；

s12：判断加载是否结束，即准静态载荷步数n是否达到nl，若未结束，则增大准静态载荷，即令n＝n+1后返回步骤s2；若加载已结束，则终止循环。

进一步，所述步骤s1具体为：

s101：给定/修正轴承载荷分布范围参数ε；

s102：查询线接触载荷分布积分表并通过线性插值得到载荷分布积分jr(ε)；

s103：由式δr＝pd/[2(1-2ε)]计算得到内、外圈在径向载荷方向上的相对位移δr，其中pd表示径向游隙；

s104：由式计算得到计算外载荷其中z表示滚子数目，kn表示滚子与内、外圈之间总的负荷-变形常数；

s105：判断计算外载荷与给定的径向载荷fr之间的关系是否满足收敛条件，若则按照式修正ε并重复步骤s101-步骤s104，其中，外载荷收敛精度ξ取0.001，外载荷修正系数β在0.001与0.01之间取任意值；若则依据式计算各滚子的载荷其中，ψj表示各滚子的方位角，j＝1,…,z。

进一步，步骤s3中，所述滚子弹塑性变形u3由两部分组成：(1)由表面法向接触应力p引起的法向弹性变形(2)由塑性应变ε^p引起的残余变形

进一步，步骤s4中，所述表面法向接触应力p的求解包括以下步骤：

s401：给定初始接触应力分布p0；

s402：利用fft计算接触体的法向弹性变形

s403：确定第k步搜索方向也即共轭梯度方向；

s404：确定第k步沿方向修正应力p所用的步长βp；

s405：按照式修正应力p；

s406：判断法向接触应力p是否收敛，即判断应力在计算域上的积分是否与滚子的准静态载荷qm,n相等，若不收敛，则按照式修正应力p后返回步骤v2；若收敛，则结束循环，其中，m为受载滚子序号，假定有nk个滚子受载，则m＝1,…,nk，n为准静态载荷步数，n＝1,…,nl。

进一步，步骤s5中，所述次表面应力σ由两部分组成：(1)由表面法向接触应力p引起的次表面弹性应力σ⁽¹⁾；(2)由塑性应变ε^p引起的次表面残余应力σ⁽²⁾。

进一步，在所述步骤s6中，在根据vonmises屈服准则判断塑性区域时，塑性应变的变化服从塑性流动准则，材料的硬化规律服从双线性硬化准则。在求得塑性应变增量后，残余应力根据塑性应变增量所引起的应力扰动进行修正；在求得塑性应变后，求解表面残余变形；表面残余变形改变表面形貌，并进一步改变表面法向接触应力；对更新形貌后的表面接触问题进行重新计算；至此，建立了关于两连续迭代步间接触应力增量、塑性应变增量及表面形貌变化量的闭环迭代方案；当相邻两迭代步残余变形的误差小于既定误差时，迭代停止，否则进入下一迭代步；当前载荷步计算收敛后，则增大准静态载荷，进入下一载荷步。

本发明的有益效果在于：本方法可用于同时计算线接触滚动轴承各滚子与内圈或外圈之间的弹塑性变形量及各滚子与内外圈总的弹塑性变形量，可在考虑加工精度影响的前提下准确计算轴承在多大载荷下会发生0.0001d的弹塑性变形，计算结果可用于指导控制线接触滚动轴承的径向载荷在不至于使轴承发生永久变形量过大、摩擦振动加剧的合理范围之内，本发明提出的方法可为工程实践中线接触滚动轴承的使用提供理论指导，对延长轴承使用寿命、提高机械系统可靠性是有益的。本发明的技术内容是线接触弹塑性变形，主要针对光滑或粗糙圆柱滚子轴承、滚针轴承的弹塑性变形计算；本发明除了公开弹塑性变形计算方法外，还将求得的弹塑性变形应用在指导线接触滚动轴承径向载荷加载上，具有一定的实践指导意义。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为圆柱滚子轴承二维图；

图2为光滑当量圆柱体-半空间平面接触示意图；

图3为圆柱滚子轴承弹塑性变形数值求解流程图；

图4为粗糙当量圆柱体-半空间平面接触示意图；

图5为表面计算域的离散化示意图；

图6为计算域网格划分示意图；

图7为法向接触应力分布计算流程图。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

典型的线接触滚动轴承(以圆柱滚子轴承为例)的二维图如图1所示，其内圈沟顶直径为di，外圈沟底直径为do，径向游隙为pd，滚子直径为d，滚子数为z。圆柱滚子轴承的接触问题可简单划分为滚子与内圈接触、滚子与外圈接触两种情况，它们均属于两圆柱体间的接触问题。通常，两个圆柱体的接触可以等效为一个半径为re的当量圆柱体与一个半空间平面的接触，如图2所示。

对于滚子与内圈的接触，两圆柱的中心位于接触线的两侧，当量圆柱体的半径为：

对于滚子与外圈的接触，两圆柱的中心位于接触线的同侧，当量圆柱体的半径为：

根据轴承刚性套圈假设，轴承套圈是不发生弹性变形的刚体，滚子与套圈的接触只会使滚子产生局部的弹塑性接触变形而不会改变套圈的整体形状和尺寸。据此，在等效模型中，也假定弹塑性变形只发生在当量圆柱体上，半空间平面不发生变形。

为计算圆柱滚子轴承各滚子的弹塑性变形，首先需要确定其载荷分布。对于给定了径向游隙pd和径向载荷fr的圆柱滚子轴承来说，圆柱滚子轴承的载荷分布qψj可用图3左半部分所示的控制误差的试解法进行求解。具体求解流程如下：

首先给定/修正轴承载荷分布范围参数ε，查询载荷分布积分表并通过线性插值得到载荷分布积分jr(ε)。然后由式δr＝pd/[2(1-2ε)]计算得到内、外圈在径向载荷方向上的相对位移δr，由式计算得到计算外载荷最后，判断计算外载荷与给定的径向载荷fr之间的关系是否满足若不满足则按照式修正ε并重复上述过程，若满足则依据式计算各滚子载荷

需要注意的是，对于圆柱滚子轴承，各滚子与内、外圈的接触角同为0°，以上计算出的滚子载荷既表示滚子与内圈间的接触载荷也表示滚子与外圈间的接触载荷。

各种滚动轴承产品必然存在一定的加工精度(表面粗糙度)，传统的轴承设计方法通常不考虑其对轴承性能的影响而简单地当作光滑表面来处理，这虽然可以简化分析过程，但无疑不能用于准确分析滚子的变形，特别是塑性变形，也就无法准确指导线接触滚动轴承径向载荷的施加。当考虑轴承表面的加工精度时，图2中的当量圆柱体-半空间平面接触示意图应用更为详细的图4表示。

圆柱滚子轴承在工作过程中，对单一滚子而言，它承受着由小增大又由大变小的循环载荷，在当量圆柱体-半空间平面接触模型中，为求解当量圆柱体的弹塑性变形，可只取上述载荷变化中的“由小增大”过程，也即当量圆柱体被逐渐增大的法向载荷压入半空间平面。这是由于塑性问题是非线性的，并且与加载路径有关，因而需采用增量加载方法来描述塑性应变的累积过程。上述载荷“由小增大”的过程是通过将各受载滚子的法向载荷qψ分为nl个准静态载荷逐步加载而实现的，载荷在每一载荷步中保持不变。受载滚子的准静态载荷用qm,n表示，m为受载滚子序号，假定有nk个滚子受载，则m＝1,…,nk；n为准静态载荷步数，n＝1,…,nl。

以各方位角处受载的滚子的准静态载荷qm,n作为输入量，一种线接触滚动轴承弹塑性变形计算方法可通过以下步骤实现：

s1：对各方位角处受载的滚子施加准静态载荷qm,n；

s2：计算含表面粗糙度hr和滚子弹塑性变形u3的接触间隙h；

s3：采用cgm计算表面法向接触应力p，其中法向弹性变形的计算需用到fft；

s4：采用icm及fft计算次表面应力σ；

s5：根据vonmises屈服准则判断塑性区域；

s6：采用newton-raphson法计算塑性应变增量δε^p；

s7：计算塑性应变ε^p所引起的残余应力σ⁽²⁾；

s8：判断塑性应变增量δε^p是否收敛，若不收敛，则更新塑性应变增量δε^p后返回步骤s6；若收敛，则转入步骤s9；

s9：计算表面残余变形

s10：判断表面残余变形是否收敛，若不收敛，则返回步骤s2更新接触间隙h；若收敛，则转入步骤s11；

s11：判断加载是否结束，即准静态载荷步数n是否达到nl，若未结束，则增大准静态载荷(令n＝n+1)后返回步骤s1；若加载已结束，则终止循环。

上述过程可用图3右半部分所示的弹塑性变形数值求解流程图表示。

在步骤s2中，表面弹塑性变形u3的求解过程如下：

计算出法向接触应力分布p之后，即可求解当量圆柱体的表面法向变形u3(也即弹塑性变形)。表面法向变形u3由两部分组成：(1)由表面法向接触应力p引起的法向弹性变形(2)由次表面塑性应变ε^p引起的残余变形

半空间表面点(x,y)处的弹性变形可根据分布的应力并利用以下积分得到：

式中，g3^p是表面法向弹性变形的格林函数：

假定滚子材料的弹性模量和泊松比分别为e1和v1，轴承内外圈使用相同的材料，其弹性模量和泊松比分别为e2和v2，当量圆柱体的弹性模量用当量弹性模量e^*来表示，

因此，由表面法向接触应力p产生的法向弹性变形可由以下公式求得：

如图5，假定将待计算表面离散成n1×n2个矩形单元，每个单元的大小为2δ1×2δ2，且计算点位于各单元的中心。u3[α,β]和p[α,β]分别是以(2αδ1,2βδ2)为中心的单元的法向变形和应力。该单元的应力分布遵循p[α,β]·y(x-2αδ1,y-2βδ2)，其中y(x,y)为形函数。式(5)中的为离散影响系数，定义为单位应力作用于以原点为中心的单元时(即p[0,0]＝1)，点(2αδ1,2βδ2)处所产生的响应：

根据求解塑性应变问题的半解析方法，本征应变在半空间表面点(x,y)处引起的残余变形可通过以下体积分求解：

式中，为作用在表面点(x,y)处的集中法向应力p所引起的(x′,y′,z′)处的弹性应力，为半空间表面原点的单位法向力所产生的弹性应力，其封闭解如下：

其中，r²＝x²+y²，γx＝ln(x+ρ)，γy＝ln(y+ρ)，γz＝ln(z+ρ)，

如图6，采用大小为2δ1×2δ2×2δ3的均匀微小立方单元对求解区域进行网格划分。n1，n2，n3分别为沿x，y，z方向的单元数。将每个表面下单元中的应力和应变视作常数，例如以点(2αδ1,2βδ2,2γδ3)为中心的单元的塑性应变可采用该点处的塑性应变表示。相应地，为以(2αδ1,2βδ2)为中心的表面离散区域的法向残余变形。

因此，式(7)的离散形式可表达为：

式中，为残余法向变形的影响系数，的三重不定积分可由式(8)-(13)中得到。

弹塑性变形为弹性变形与残余变形之和：

在步骤s3中，表面法向接触应力p的求解过程如下：

如图4，以初始接触点为原点建立笛卡尔坐标系，且z轴指向当量圆柱体内部。计入边界约束的通用干接触模型如下：

h(x,y)＝h0(x,y)+u3(x,y)+hr(x,y)-ω(17)

其中，式(16)为载荷平衡条件，p为法向接触应力。如图4，h为接触间隙，h0为初始接触间隙(对当量半径为re的圆柱体，h0＝x²/2re)，u3为法向弹塑性变形，hr为接触体表面粗糙度，ω为表面法向刚体趋近量。

接触产生了表面接触区域ac，在接触区内，需满足：

p(x,y)＞0,h(x,y)＝0(x,y)∈ac(18)

在非接触区，需满足：

可以看出，上述方程所描述的接触问题构成了一个线性补余问题，可用矩阵形式表示为：

h(x,y)≥0,p(x,y)≥0,h(x,y)p(x,y)＝0(21)

式中，为常数矩阵，k为弹性变形影响系数矩阵。

根据变分原理，式(20)和(21)描述的补余问题等价于二次函数的条件极值问题。也即寻找目标函数：

在约束条件(21)下的极小值。

二次函数的条件极值问题通常可采用cgm求解，表面法向接触应力p的具体求解过程如下：

v1：给定初始接触应力估计p0；

v2：利用fft计算接触体的法向弹性变形

v3：确定第k步搜索方向也即共轭梯度方向：

式中，wp为目标函数w(p)关于接触应力p的偏导数，k表示第k步迭代。

在接触区ac内先求得接触间隙的平均值然后修正接触间隙则：

v4：确定第k步沿方向修正应力所用的步长βp：

为得到式中的先利用fft计算中间变量

然后在接触区ac内求得的平均值同时修正

v5：修正接触应力p：

v6：计算应力迭代误差ep和载荷平衡误差el：

判断计算所得外载荷是否与准静态载荷qm,n平衡，是否需要修正应力：

式(29)和(30)中，若ep＞εp或el＞εl，则按照式(32)修正应力p后返回步骤v2；若ep≤εp且el≤εl，则停止迭代，得到法向接触应力分布p，其中，εp、εl表示迭代精度。

图7给出了上述表面法向接触应力分布计算流程。

在步骤s4中，次表面应力σ的求解过程如下：

当量圆柱体次表面应力的求解是塑性变形演化问题中的关键步骤。次表面应力σ由两部分组成：(1)由表面法向接触应力p引起的次表面弹性应力σ⁽¹⁾；(2)由次表面塑性应变ε^p所引起的次表面残余应力σ⁽²⁾。

当量圆柱体点(x,y,z)处由表面法向接触应力p引起的次表面弹性应力σ⁽¹⁾可采用如下的二维卷积公式计算得到：

采用与图6相同的立方单元系统，在每个离散表面的小区域内的应力为常数，进而，以坐标(2αδ1,2βδ2)为中心的表面单元区域的应力采用p[α,β]表示。是在以点(2αδ1,2βδ2,2γδ3)为中心的单元上的均匀弹性应力。任意离散单元处的次表面弹性应力为：

式中，是由表面法向接触应力对应力的影响系数。的二重不定积分可从式(8)-(13)中得到。

以坐标(2αδ1,2βδ2,2γδ3)为中心的单元的残余应力可表示为：

式中，为影响系数。

将弹性应力和残余应力叠加，即可得到总的次表面应力：

在步骤s5中，根据vonmises屈服准则判断塑性区域的方法如下：

塑性是材料对所施加载荷表现出的一种不可逆行为。如式(37)，采用适用于大多数金属材料的vonmises屈服准则对材料屈服进行判断。当vonmises等效应力超过材料局部屈服强度时(也即f＞0)，材料发生塑性变形：

式中，σvm是vonmises等效应力，“:”表示矩阵内积运算，sij＝σij-σkkδij/3为偏应力，g(λ)是屈服强度函数(g(0)等于初始屈服强度σy)，为有效累计塑性应变。加载过程中，有效塑性应变的增量dλ和屈服函数f满足以下条件：

f≤0,dλ≥0,f·dλ＝0(38)

当采用vonmises屈服准则时，塑性应变的变化服从以下塑性流动准则：

屈服强度函数g(λ)服从双线性硬化准则，在屈服阶段，g(λ)可表示为：

式中，et是材料的切向模量，e是材料的弹性模量。

在步骤s6中，塑性应变增量δε^p的求解过程如下：

塑性应变的增量可通过上述应力场和塑性模型计算求得。当f(λ)＞0时，材料发生塑性屈服。在塑性区域，实际的有效塑性应变增量δλ应满足f(λ+δλ)＝0。f(λ)是关于δλ的非线性表达式。因而，可采用newton-raphson法求解塑性应变增量。屈服函数可近似扩展为：

f⁽ⁿ⁺¹⁾＝f⁽ⁿ⁾+δλ⁽ⁿ⁾f,λ⁽ⁿ⁾＝0(41)

两连续迭代步之间的有效塑性应变修正量δλ⁽ⁿ⁾可表示为：

式中，

式中，μ是材料的剪切模量，kp、γ是运动硬化常数。

所有相关变量根据下列式子进行更新：

式中，λ⁽¹⁾、σvm⁽¹⁾和σij⁽¹⁾分别表示初始有效塑性应变、背应力、等效vonmises应力和柯西应力分量。迭代在满足以下收敛条件和收敛精度εg时停止：

若未收敛，则重复执行式(42)-(45)，直至收敛。塑性应变增量可从如下塑性流动准则中得出：

在求得塑性应变增量后，残余应力需根据塑性应变增量所引起的应力扰动进行修正，见式(35)。在求得塑性应变后，表面残余变形可通过由式(14)求解。表面残余变形将改变表面形貌，而进一步改变表面法向接触应力。因此，需对更新形貌后的表面接触问题进行重新计算。至此，建立了关于两连续迭代步间接触应力增量、塑性应变增量及表面形貌变化量的闭环迭代方案。当相邻两迭代步残余变形的误差小于既定误差时，迭代停止，否则进入下一迭代步。当前载荷步计算收敛后，则增大准静态载荷，进入下一载荷步。完整的圆柱滚子轴承弹塑性变形数值求解流程图如图3所示。计算结果可为工程实践中线接触滚动轴承的使用提供理论指导，控制轴承的径向载荷在不至于使轴承发生永久变形量过大、摩擦振动加剧的合理范围之内。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。