李文婧,王泓晖,刘贵杰,田晓洁,李亚楠

(中国海洋大学 工程学院,山东 青岛 266000)

0 引 言

随着科学技术的快速发展,现代动车组装备逐渐趋于多元化、复杂化。高速动车组服役环境变化多样,常处于高速重载、剧烈冲击振动等恶劣工况下,其一旦出现故障,势必影响行车平稳性和安全性[1],严重时将造成不可估量的经济损失及人员伤亡。

轴箱轴承作为高速动车组走行部的关键部件之一,其往往处于高速甚至超高速的运行状态,此时轴承的磨损是一个非常突出的问题,有必要对轴箱轴承开展磨损损伤状态分析与研究[2]。但目前,获取磨损状态信息不充分、难度大、成本高等因素是转向架轴承磨损退化分析时存在的主要不足[3]。

因此,建立转向架轴承的磨损模型,模拟在轨道激励载荷因素作用下滚动轴承的磨损状态,同时开展轴承的磨损损伤分析具有重要的研究价值和实际意义。

现阶段,围绕着模拟轴承的磨损退化研究,国内外学者已经展开了许多研究工作。HWANG S Y等人[4]利用网格删除方式来复现轴承实际的磨损形貌,根据轴承几何形状和跳动值的变化,预测了轴承的磨损程度。SCHMIDT A A等人[5]通过实施节点位移偏移,得到了轴承接触表面复杂的动态磨损行为,最后检查了所有网格的畸变尺寸,得到的磨损量符合实际验证结果。张金萍等人[6]使用ADAMS建立了轴承转动系统虚拟样机,通过改变内部间隙来获得不同磨损状态间的差异,最终得到了能够有效反映轴承磨损状态的特征变量。ASHRAF M A等人[7]利用有限元分析方法,探究了聚合物表面接触的滑动磨损过程,从而确定了磨损体积随时间的变化,并预测了磨损的几何结构。MUKRAS S等人[8]运用非线性有限元分析方法,模拟了轴承磨损形貌变化,同时提出了数值积分方案,以及并行计算两种不同的最小化磨损计算成本方法。

由于磨损的模拟方法通常需要多次迭代循环,计算耗时长、成本高,尤其对于轴箱轴承这类尺寸大且工况复杂的零部件,上述问题则尤为突出。为此,非常有必要探索一种在保证精度可靠的前提下能够有效降低仿真时间成本的数值模拟方法。

依据以上论述,笔者以转向架轴箱轴承为研究对象,结合Archard磨损模型[9],提出一种切片式半解耦磨损损伤计算方法,对三维轴承模型进行多组二维切片磨损状态的仿真,利用高斯过程拟合切片的仿真结果,探究轴箱轴承内圈磨损分布规律。

在保证结果有效性的条件下,笔者最小化有限元磨损分析方法的计算时间,以期为转向架轴承的磨损损伤分析方法提供依据和参考。

1 轴箱轴承磨损有限元分析方法

1.1 研究对象

为了研究轴箱轴承的磨损损伤分析方法,笔者选取CRH380A高速动车组转向架轴承作为具体的研究对象。由于高铁轴箱轴承基本使用双列圆锥滚子轴承,其中两列轴承的材料属性以及几何参数完全一致,所以笔者仅对其中的一列轴承进行分析。

单列圆锥滚子轴承如图1所示。

图1 单列圆锥滚子轴承

对于高速动车组而言,其运行条件往往较为严苛。其内部轴箱轴承主要受到来自列车车体的垂向载荷、速度变化带来的冲击力和轮轨激励等,轴承内部部件间的接触状态十分复杂[10-12]。

由于轴承内滚道外径小,单位面积的接触应力大,使得轴承内圈成为磨损的重灾区。因此,笔者的研究侧重于轴箱轴承内圈的磨损损伤分析,探究的重点也聚焦于其内滚道。

轴承各部件相关的材料参数如表1所示。

表1 轴承各部件材料参数

1.2 Archard磨损模型

由于轴箱轴承内部承受较大的接触载荷,滚子的自转运动导致滚动体与滚道间极易发生磨损失效。而磨损会改变滚道的几何外形,从而改变滚子与滚道间的接触应力,加速轴承的磨损失效进程。

Archard磨损模型目前被广泛应用于材料磨损的计算场合,依据磨损体积的改变映射出材料磨损程度。

在该模型中,材料的磨损量与法向接触力及滑动距离成正比,与材料的硬度成反比,即滚动轴承内圈的磨损量可表示为:

(1)

式中:Vm为磨损量,mm3;W为接触面法向力,N;Lm为相对滑移距离,mm;H为材料硬度,N/mm2;K为磨损系数。

对于滚子轴承而言,内滚道与滚动体的接触、外滚道与滚动体的接触以及滚子与保持架的接触为线或面接触,接触区域属于小面积接触。

其公式可改写为:

(2)

式中:hs为磨损深度,mm;k为接触面无量纲磨损系数与接触面材料表面硬度的比值,即有量纲磨损系数,k=K/H;p为接触面正压力,MPa;s为接触滑移距离,mm。

笔者利用有限元方法分析轴承磨损过程,将整个磨损过程离散成多个增量步的磨损增量堆积。

由于轴承内部的接触情况是瞬时发生的,即滚子与内外滚道的接触位置、接触压力以及滑移率都随时间的变化而发生快速变化,并且也会随磨损而变化;所以,此时笔者假设在每个非常小的增量步内,滚子与轴承内圈接触点处的压力和磨损系数是一个固定值。

其接触点的磨损深度表示为:

(3)

式中:hn为在第n个增量步下的磨损深度;hn-1为第n-1个增量步下的磨损总深度;pn为第n个增量步下的接触应力。

1.3 切片式半解耦损伤分析方法

磨损分析每个周期均需要考虑非线性接触问题,若每个增量步都进行迭代计算,需通过大量迭代循环来模拟磨损演化过程;再加之接触点发生磨损后,每个增量步对应的网格需要重绘,相应的计算成本大幅增加,并且模型计算效率也十分低下。

为应对此问题,笔者提出一种切片式半解耦损伤分析法。该损伤分析方法的关键步骤如图2所示。

图2 切片式半解耦损伤分析法关键步骤

在图2中,切片式半解耦损伤分析方法的具体实施过程如下:

1)选取合适间隔,将轴承三维有限元模型沿径向等距切分为多组,根据每组切面几何特征,构建二维有限元模型,并将其作为后续分析的基本单位(后文称之为二维切片模型);

2)将若干次磨损循环看作一个磨损周期Δq,该周期内的磨损体积变化很小,即可以认为在每个磨损周期内的应力场不因磨损而改变,即为定值,该过程称为“半解耦”;同时设定该圈数Δq=10为磨损分析步长,并通过计算得到该步长下每个二维切片模型的磨损增量(深度),即:

(4)

3)根据该磨损增量完成二维切片模型的几何重构;

4)根据所有单元磨损增量更新轴承磨损状态,重新求解新应力场,并进行下一段循环,直至达到预设的磨损仿真增量步或者磨损深度阈值;

5)期间,考虑到切片的原因,使三维轴承离散化,仅通过有限数量的二维切片模型难以直接给出轴承的磨损状态,需借助高斯过程非线性问题的强大拟合能力,对切片模型之间的磨损状态进行拟合。

假设多组二维切片计算得到的磨损深度数据集N={(xi,yi)|i=1,2,…,n},使用高斯函数表示:

(5)

式中:i为第i个接触对;xi为在第i个接触对下接触点距轴承端面的距离;yi为对应接触节点的磨损深度;ymax,xmax为待估计参数,S分别为函数曲线的最大值、最大值位置以及高斯曲线半宽度信息。

对上式等号两侧同时取对数,可以得到:

(6)

将右侧平方公式展开可得:

(7)

(8)

将其简化为:

Z=XA

(9)

再依据最小二乘法,矩阵A的解为:

A=(XTX)-1XTZ

(10)

联立式(5)可以求解出待估计参数ymax,xmax以及S。

依据此高斯拟合式,笔者拟合多组切片的磨损分析结果,复现三维磨损状态。

1.4 轴承磨损有限元分析

1.4.1 磨损仿真分析流程

由于磨损是一个增量过程,一个完整的磨损周期仿真包括以下几个步骤:

1)结合有限元理论,利用ABAQUS求解轴承滚子与内外圈的非线性接触问题,获得节点接触力、相对位移和节点坐标等信息;

2)利用Archard磨损公式求解节点的偏移量,偏移方向为接触点的法向,二者均通过UMESHMOTION用户子程序[13,14]来实现,将磨损深度计算结果反映于节点的移动。

轴承内圈发生一定程度的磨损退化后,会改变内圈接触表面的几何形状,为此,笔者采用ALE网格偏移技术,对轴承内圈磨损退化后的几何形状进行更新,如图3所示。

图3 ALE网格自动偏移技术

在每个增量步开始时,笔者重新求解非线性接触问题,同时计算该增量步下的磨损增量,不断循环往复,直到完成全部增量步计算。

USESHMOTION磨损计算子程序使用Fortran语言编写,在ABAQUS二次开发中实现。

据此得到磨损仿真分析具体流程,如图4所示。

图4 磨损流程图

1.4.2 磨损有限元模型

接下来,笔者采用ABAQUS隐式方法对二维轴承切片磨损状态进行分析。

为了便于计算收敛,笔者对整个加载过程进行了细化:1)内圈全部节点耦合在圆心处,在初始分析步中,仅放开沿y轴的垂向自由度;2)保持架全部节点耦合在圆心处,在初始分析步中,仅放开绕中心轴线的旋转自由度;3)外圈始终添加全约束。

全程设置2个分析步step1至step2:step1在内圈耦合点添加垂向载荷25 000 N;step2放开内圈绕轴线转动的自由度,同时添加内圈转速1 000 r/min。滚子与内外圈的接触形式采用表面与表面接触,以滚子表面作为主面,内滚道外表面作为从面的接触方式;在磨损分析之前,笔者对内圈磨损区域预设ALE网格自适应区域,考虑到网格单元的适用性,划分时单元类型选择CPS4(四节点双线性平面应力四边形单元),同时对滚子与内外滚道的边缘区域进行网格细化处理。

构建的二维切片有限元模型如图5所示。

各部件均设置为弹性材料。接触面设置为表面与表面(surface to surface)接触,接触属性分为切向和法向两部分,切向行为采用静摩擦-动摩擦指数衰减:静摩擦系数为0.1,动摩擦指数为0.05,衰减指数为0.01;法向行为采用硬接触。

由于模型尺寸较大,为了保证计算的可靠性,需先探讨有限元模拟计算中轴承内部接触状态的网格有效性。笔者共设计5组网格数量,针对滚子及内外滚道采取不同程度的网格细化处理,以此来完成对网格有效性的验证。

5组网格的具体划分尺寸及数量如表2所示。

表2 网格具体参数

模拟进行0.1 s后,笔者观察其模拟接触力结果,如图6所示。

图6 内滚道平均接触应力

由表2及图6可知:网格数量对有限元模拟仿真结果有较大影响,当网格划分尺寸为2 mm时,此时网格数量为49 795,模拟结果开始趋于稳定。

依据上述分析,笔者选用尺寸为2 mm的网格对轴承接触区域进行细化处理。

2 磨损结果分析与讨论

2.1 方法有效性验证

2.1.1 模型高效性与收敛性对比

不同切片数量会影响轴承内圈摩擦磨损分析的时间差异。

笔者以5为初始切片数量,对三维模型进行切片,并行研究数量为5片、10片、20片以及30片,共4组切片式来验证半解耦分析法的计算效率,将计算结果与传统三维磨损仿真分析的结果进行比较,并且将该对比实验在同一台计算机上执行(用于比较时间)。

在仿真时,笔者将2种轴承磨损损伤的分析方法均设置为105个增量步。切片计算结果的高斯拟合过程在MATLAB中完成。由于该拟合过程运行时间较短,笔者主要对比有限元分析完整过程的时间损耗。

相关计算结果,即磨损方法计算时间对比如表3所示。

表3 磨损方法计算时间对比

需要注意的是,在切片式半解耦损伤计算方法中,每一个二维切片模型的计算时长较为固定,通常在6.5 h～7.2 h之间,随着切片数量的增加,相关模型的计算时间损耗呈现线性增长趋势;当切片数量控制在30片以下时,半解耦损伤分析方法可以有效缩短仿真的时间损耗;若切片数量大于30,该方法的高效性难以体现。

由于磨损分析采用隐式分析步进行,在非线性条件下极易出现模型不收敛,导致仿真终止的现象(比如:欠约束、复杂接触以及模型变形过大等,都可能使计算结果不收敛)。

笔者在确保材料属性与边界条件准确施加的前提下,讨论切片式半解耦损伤分析方法与传统三维分析方法的迭代收敛情况。

在计算中,迭代是在每个增量步中进行的。在一个增量步中,迭代的最高次数为16次,若在16次迭代中ABAQUS求得平衡解,就会进入下一个增量步,若未获得平衡解,软件将会减小当前增量步时间,继续迭代;若增量步时间折减次数超过5次,则停止计算,同时当增量步总数超过设定的最大增量步数目时,停止计算。

此处,笔者设定了105个增量步,用于模拟轴承内圈的磨损损伤演化过程,历史收敛图(收敛情况对比分析)如图7所示。

图7 收敛情况对比分析

图7为三维模型与多组切片式半耦合模型的收敛情况,图中X轴表示对应的不同增量步,Y轴表示对应增量步下的最大连续迭代不收敛次数。

笔者分析图7中曲线得到以下结论:

1)无论切片数量为多少,其对应的二维磨损模型的收敛情况均相似;在30 000个增量步之前,模型的最大连续不收敛次数呈下降趋势;在连续计算40 000个增量步之后,切片模型的磨损分析会偶尔出现1次迭代不收敛;

2)与切片式半解耦损伤分析法相比,三维磨损模型的最大连续迭代不收敛次数非常不稳定,主要集中在3次～4次之间,该次数濒临模型计算的崩溃阶段。

2.1.2 模型准确性验证

验证上述方法有效性的一般手段为将切片式半解耦损伤分析方法的结果与传统的有限元分析结果进行比较。笔者将该方法与轴承三维有限元磨损仿真分析方法进行对比,对切片式半解耦磨损分析法的模拟结果进行了验证。

传统的三维建模包括轴承内外圈、滚动体及保持架四部分,建模时不考虑密封与倒油孔的影响;在三维磨损分析中,笔者同样设置2个分析步(step1和step2),施加法向载荷25 000 N,添加转速为1 000 r/min;接触形式、增量步及边界条件设置与上述一致。

为了方便探讨三维磨损模型中的最大磨损深度,笔者在轴承内滚道处列出了检测磨损增量的测点位置标识,如图8所示。

图8 轴承内滚道磨损测点位置

随着磨损仿真增量步的增加,磨损增量测点的深度也发生相应的变化。

在三维模型完成全部增量步计算后,笔者统计每个测点的磨损深度,得到的结果如表4所示。

表4 对应测点磨损深度

通过分析表4中数据可知:在105个增量步仿真时长结束后,三维模型轴承内圈仿真计算的最大磨损深度为2.460 mm,最小磨损深度为0.072 mm。

笔者依旧采用切片数量为5片、10片、20片以及30片的半解耦分析方法的磨损计算结果与之进行比较。

相关磨损对比结果如表5所示。

表5 磨损结果对比

分析表5可知:

5切片的二维模型的最大磨损深度拟合结果为1.907 mm,与三维模型的误差达22.5%;当切片增加至10片时,误差下降到10%以内,为9.1%,可信度提高;随着切片数量的逐渐增加,仿真计算时间成本提升的同时,误差也逐渐降低;当切片数量达30片时,拟合最大磨损深度为2.426 mm,误差仅有1.4%,但此时该组耗时与三维模型接近。

2.2 磨损结果分析

根据上述分析结果,笔者提取20片切片的半解耦损伤仿真计算的接触力及磨损深度等数据,以探讨轴承内圈磨损情况。

轴承内圈接触力情况,即轴箱轴承滚动体与内圈接触应力曲线,如图9所示。

图9 轴承内圈接触压力

图9中,在滚子与内圈的接触起点,据内圈中心的转动角度为0度的接触对中,处于承载区的滚动体由于径向载荷作用受到来自轴承内外圈的挤压作用,使得轴承内部径向间隙减小至0,滚动体与内圈接触连续且稳定,可以得出内圈接触区域的平均接触应力为703.2 MPa,然后其沿着滚动体中心线与内圈中心线夹角的增大而逐渐降低;当滚动体运动到非承载区时,滚动体所受内外圈作用力、振动及撞击力较小且运动状态较不稳定,内部径向间隙较大,接触并不完全,使得内圈非承载区接触区域的接触应力接近0 MPa。

利用二维拟合思路,笔者将轴承内圈不同位置节点对应的磨损量进行拟合,得到了轴承内圈磨损分布结果,如图10所示。

图10 轴承内圈磨损分布

由图10可知:由于边缘压膜效应[15],在接触发生的边缘区域产生应力集中现象,使得滚动体与内滚道接触面的最大磨损深度出现在滚动体边缘位置;而在内圈中心区域,由于滚子与内圈的接触状态稳定,磨损也较为平均;在105个增量步仿真结束后,内圈接触边缘区域的最大磨损深度达2.38 mm。

在接触区域的中间位置,滚动体自转运动状态较稳定,接触应力分布比较均匀,其平均磨损深度为1.44 mm。

轴承内滚道最大磨损深度变化,即随着增量步的增加,轴承内圈最大磨损深度的变化情况,如图11所示。

图11 轴承内圈最大磨损量变化

图11中,在磨损系数不变的情况下,磨损深度主要受接触应力及接触位移的影响[16]。

从图11中曲线可以看出:最大磨损深度并非随着轴承转动增量步的增加而产生线性变化,其斜率(磨损率)随着磨损循环次数的增加而逐渐减小;导致轴承内圈最大磨损深度非线性增加的重要因素是随着磨损深度的增加使得内部径向间隙增大,其直接造成了内圈接触区域的接触应力非线性变化[17,18]。

3 结束语

笔者以转向架轴箱轴承为研究对象,基于数值模拟的方式,提出了一种切片式半解耦磨损分析方法,构建了多组轴承内圈磨损损伤的二维有限元仿真模型,并在给定工况下,对该轴承内圈的磨损状态进行了模拟及结果拟合,最终得到如下结论:

1)笔者使用切片式半解耦磨损损伤分析法、高斯拟合以及ABAQUS中强大的ALE网格偏移及UMESHMOTION子程序技术,其可替代三维模型磨损分析方法,在保证误差在10%以下的同时,该模型的计算效率比三维模型提高了约3倍。针对转向架轴箱轴承及其他结构尺寸较大、仿真分析时间较长的模拟模型,该方法不失为一种高效的仿真研究方法;

2)研究发现,轴承内圈表面磨损最为严重的地方出现在滚动体与内滚道接触的边缘区域,而接触的中心区域磨损程度较轻(由于其接触状态稳定)。同时,内圈最大磨损深度并不是随着增量步的变化而线性增加,这是由于轴承的磨损致使部件间隙发生变化,导致滚动体与内圈接触应力发生非线性变化。

目前,笔者仅对轴承内圈的磨损损伤进行了数值模拟分析。在后续研究中,笔者将开展其他零部件磨损损伤研究,以进一步了解轴承的磨损特性,同时也拟通过实验,对磨损的数值模拟结果进行验证。