粗糙面滑动摩擦热-应力-磨损耦合仿真分析

2022-10-21安先龙王佳赵诚高超

装备制造技术 2022年7期

安先龙，王佳，赵诚，高超

（中机科（北京）车辆检测工程研究院有限公司，北京 102100）

0 引言

工程表面微观上是凹凸不平的，其摩擦过程实质为摩擦副粗糙表面上微凸体间的相互摩擦过程[1]。摩擦过程中，摩擦功产生大量的热量，改变摩擦副摩擦磨损特性，同时由于压力的作用，摩擦副间磨损的产生也会改变接触表面的压力、温度场的分布状态。因此，摩擦过程涉及复杂的应力、温度和磨损耦合现象[2，3]。

摩擦副粗糙面具有分形特征，Mandelbort[4]创立分形几何学，并在Weierstrass函数的基础上，提出了分形W-M函数，因其具有自相似性、自仿射性，因此，可被用来描述粗糙表面形态。后来有不少多学者对分形理论进行了研究，黄健萌等[5]通过建立粗糙面滑动接触模型，仿真分析粗糙实体的热、应力分布规律。丁雪兴[6]通过分形理论建立了磨粒磨损的预测模型，并推导出磨损率和分形模型参数之间的关系。Wang等[7]通过分形理论表征分形粗糙表面，分析了粗糙表面微凸体处的摩擦温升。Gong和Komvopoulos[8]通过建立球体弹塑性有限元模型，研究了球体几何形状对接触压力分布和表面应力-应变场的影响规律。

依据W-M分形函数，在Matlab、Catia软件中联合建立接触式干滑动模型，然后基于应力、温度和磨损耦合方法，通过编写Umeshmotion子程序代码，在ABAQUS软件中建立微凸体表面干滑动磨损仿真模型，探究滑动过程中粗糙面微凸体的接触压力、磨损和温度场的分布影响规律。

1 计算模型建立

1.1 粗糙面模型建立

W-M分形函数模型[9]，其三维函数表达式：

式中，L为取样长度；G为特征尺度系数；cosφm，n表示（0，2）间均匀分布的随机相位；D为分形维数；M为重构表面时叠加峰值的数量，一般选为10；γ为随机轮廓的空间频率；n为频率系数，nmax为频率的上限。分形维数和特征尺度系数是决定粗糙面形貌的主要参数[10]，为减少有限元仿真分析占用过多计算资源，仅对摩擦副上粗糙面采用分形理论建模，下粗糙面则采用理想平面实体。

分形接触模型如图1所示。摩擦块尺寸：内半径（半径）270 μm，外径（半径）420 μm，厚度80 μm，摩擦盘尺寸：内径（半径）270 μm，外径（半径）430 μm，厚度为100 μm。

图1 销-盘模型

1.2 计算模型假设条件

此计算模型的假设条件如下：

（1）摩擦副材料（摩擦块和摩擦盘）各向同性。

（2）摩擦过程符合库伦定律，摩擦盘在转动过程中，由于时间较短，默认摩擦系数为定值。

（3）忽略摩擦过程中惯性力和振动造成的影响。

（4）忽略热辐射。

1.3 计算模型热边界条件

热边界条件是温度场周围的表面热交换条件，共有3种[11]。如下所示：

（1）第一类边界条件表明边界上的温度值是已知的，故不存在第一类边界条件。

（2）第二类边界条件规定了边界上热流密度，其方程表示为：

式中，q为边界热流密度，μ摩擦系数，p为接触压力，ϑγ为相对滑动速度。

默认接触面产生热量在摩擦块和摩擦盘间自由分配，根据能量守恒和热传导的傅立叶定律，摩擦副热传导方程表示为：

式中，λ为热导率，ρ为密度，c为比热容。

（3）第三类边界条件定义了流体温度以及物体与流体之间的传热系数，流体与固体表面间的对流换热[12]方程表示为：

式中，q为散热的热流密度；h为对流换热系数；Ts为固体温度；Tf为流体的温度。

1.4 计算模型Archard磨损理论

使用Archard理论模拟表面磨损损伤[13]，其广义增量磨损高度定义如下：

式中，△h为磨损高度增量，k为无量纲常数，是表征材料耐磨性的磨损系数，可由试验获得；p为接触压力；△s为相对滑动位移增量。

假设摩擦块有限元模型接触节点为I，积分步骤为j，则磨损增量表示为：

累计磨损表示为：

2 热力耦合和磨损仿真分析

2.1 有限元模型建立

摩擦块材料为机械化合金，摩擦盘材料为铸钢，摩擦系数设为0.3，材料性能参数见表1。摩擦块和摩擦盘六面体单元属性均为C3D8T，摩擦块单元数为54000，摩擦盘单元数为19136。为提高仿真结果准确性，摩擦块接触区域网格做加密处理，划分后模型网格如图2所示。

表1 材料属性

图2 有限元模型网格划分

为明显区分磨损所造成的影响，因此选取工况为：第一步为压力分析步，即在40 μs内将20 MPa法向压力通过摩擦盘施加到摩擦块上；第二步为滑动分析步，即在法向压力的作用下摩擦盘以10 m/s的线速度单向匀速旋转40 μs。在摩擦盘中心创建参考点RP-1，由RP-1来控制摩擦盘的转动和载荷的施加。仿真过程中，摩擦盘仅保留Z方向的自由度，其他方向位移施加约束限制。

有限元仿真简化流程如图3所示，在ABAQUS软件中完成材料赋予、接触属性的创建等。仿真过程中，设置摩擦块网格为ALE自适应网格区域，通过提取ABAQUS接触压力和滑动位移，并根据式（5）编写的磨损子程序代码，计算节点的磨损增量，并按照计算的偏移方向，不断调整表面网格，仿真模拟磨损深度处于不断变化的磨损过程，计算摩擦块的累计磨损量。基于不断更新变化的表面形貌，在温度-位移耦合分析步中，得到不断更新的表面压力、温度场的分布状态。

图3 有限元仿真简化流程

2.2 仿真分析结果与讨论

图4为粗糙实体摩擦块接触面不同时刻的接触压力分布，从图中可以看出，滑动分析步初期，摩擦块端面接触压力分布不均匀，微凸体尖端处产生了巨大的应力集中，最大接触压力均达到1853 MPa，是法向压力的92.65倍。

图4 滑动过程接触压力分布

滑动初期，粗糙实体摩擦块表面微凸体与摩擦盘接触后相互摩擦，此时实际摩擦面积为微凸体接触面积，随着摩擦盘的转动，摩擦时间的增长，在法向压力持续作用下，摩擦表面微凸体逐渐磨损变形，摩擦表面微凸体不断被磨平，摩擦面积不断增大，从而使得接触面整体的接触压力逐渐下降，并逐渐在微凸体的边缘呈现出较高压力值。滑动后最高接触压力下降到419.7 MPa，滑动后最高接触压力下降为288.9 MPa。

图5为滑动摩擦块的节点磨损分布云图。从图中可以看出，接触压力分布云图与磨损分布云图具有较高的一致性，此现象符合Archard磨损定律，即磨损与接触压力成正相关。滑动后，摩擦块形成了大量块状磨损区，微凸体中心处具有较大的磨损高度。最高微凸体节点，在滑动后磨损高度为6.51。选取磨损最为严重的A、B、C、D节分析点节点磨损、温度随时间的变化规律，所在位置如图5所示。

图5 滑动后摩擦块磨损深度分布

图6为A、B、C、D四个节点磨损高度随滑动时间变化的曲线。四个节点的磨损高度变化规律相似性较高，可分为两个阶段。第一个阶段为压力分析步阶段（0～40μs），微凸体节点经过快速偏移到缓慢偏移的过程。这是因为接触初期，接触面积为微凸体接触面积，压力的增加致使微凸体节点快速沿z方向偏移，但随着接触面积的增大，节点应力降低，节点偏移速率减缓。除此之外，接触初期，A节点为最高微凸体，接触后迅速向z方向偏移，但B、C、D节点在未和摩擦盘接触时，仍发生轻微节点偏移，考虑为A节点受到压力的作用，致使摩擦块材料发生弹塑性变形，从而使得B、C、D节点沿压力方向发生轻微偏移所致。第二个阶段滑动分析步（40 μs～80 μs），此时表面微凸体节点的偏移速率由公式（5）archard磨损定律计算所得，滑动初期，磨损迅速增加，随着接触面积的进一步增大，接触压力进一步下降，磨损速率变得缓慢。

图6 节点磨损高度随时间变化曲线

滑动后模型的温度场分布如图7所示。温度场分布以块状热区为主，滑动后最高温度为76.34℃，A、B、C、D 4个节点的温度随时间变化曲线如图8所示，滑动初期，微凸体节点摩擦产生的热量来不及传递，温度聚集在微凸体内部[14]，节点温度迅速增加。以温升较大的节点B为例，除接触压力较大影响外，其位置靠近摩擦块外径，相对滑动速度较大也是造成温升较大的原因。节点B在3.18内迅速温升至171.847℃，随着滑动的进行，节点B微凸体被磨平，其他微凸体陆续进入摩擦过程，接触面积不断增大，接触压力不断下降，磨损消耗了部分热量导致节点产生的热量减少，同时加上节点热传导的散热作用，温度大幅下降后趋于稳定。