磨损率

1/11

有限元模拟滑动摩擦磨损

摘要

磨损往往是影响产品寿命的一个主要因素。因此磨损预测就成为工程的一个重要部分。

这篇论文介绍了用有限元软件ANSYS来模拟磨损的方法。用线性磨损定律和欧拉解析积分

提出了一个模型化的模拟程序。然而，还要考虑保证模型的正确性和数学方法的收敛性。

分别用实验和有限元的方法分析了球形pin-on–disk系统在没有润滑条件下的接触问题，使

用了Lim和Ashby磨损图来区分磨损机理。在给定几何尺寸和载荷的条件下，可以用有限

元的方法模拟磨损，得到磨损率对滑动距离的对应关系。有限元软件ANSYS非常适合解决

接触问题和磨损模拟。实际磨损率的分布范围在40-60%的界限内会导致磨损模拟结果相

当大的偏离。因此这些结果必须在一个相对的值上进行估测，从而比较不同的设计。

关键词：磨损模拟；FEA；磨损试验；接触温度

1．绪论

摩擦副之间最可靠的摩擦学行为的知识可以通过做磨损实验来获得。然而，当特别是设

计改变时需要在日常的内部程序基础上进行迅速的估测。已经进行了大量的研究工作来帮助

设计者实现这一步。

已经证实一个给定系统滑动磨损的主要参数是接触载荷和相对滑动速度。速度由机构运

动来决定。系统载荷怎么影响接触应力是很复杂的一个问题。第一个分析两个弹性实体接触

应力的人是赫兹。他认为接触体是弹性的，接触部分为椭圆形，而且没有摩擦的。这些假设

被用在接触应力的计算中。磨损发生在机械构件相互接触时。一个重要的实际问题是在给定

的时间里有多少的材料损失。由于功能和加工误差等表面的形状是不同的。而且会因为磨损

和弹性变形而改变。因此压力的分配就依赖于这些条件。有限元的方法是一个通用的工具来

解决应力应变的问题。这篇论文使用有限元软件ANSYS5.0A分析了接触压力和磨损模拟。

2.磨损模型

磨损过程可以认为是动态的，由许多参数决定，这个过程的预测可以看作是一个初始值

的问题。从而磨损率就可以由一个总的方程来描述。

dh/ds=f(载荷，速度，温度，材料参数，润滑，….)

h为磨损深度（m），s为滑动距离（m）。文献中可以查到很多磨损模型。这些数学表达多种

多样，从简单的经验式到复杂的依赖于物理概念和定义的方程[1]。常常包括特定的参数和

变量，只有在特定的情况下是正确的，在手册中是查不到的。因此这些模型很少被用于实际

中磨损的估测。

Lim和Ashby给出了钢板一个更广泛的在大范围载我的爷爷 荷和速度下的磨损分类方法。他们的

工作建立在简化磨损方程并在很多pin-on-disk实验数据的基础上进行调整。通过这个工作

得到了一个磨损图，图1所示。给出了磨损机理的轮廓和无量纲化磨损率，作为无量纲

化压力p和无量纲化速度v的作用结果，定义如下

（1）

V为磨损体积（m3），A为接触面积（m），r

0

为接触半径（m），F

N

为载荷，H为接触对中

较软材料的硬度（P

a

），v为首尾相应 相对滑动速度（m/s）。表一列出了Lim和Ashby使用的方程和

参数。

2/11

图1所示磨损图中的温度的分析是假定一个简单的温度沿一维方向流动的基础上进行的。更

进一步的说就是迅速传播的温度对磨损起了一个很重要作用，热量分配系数12=0.5，如果

接触的温度达到7000C，就会发生严重的氧化磨损。在这个温度一下，磨损与载荷成线性关

系，不受速度影响。

最常用的模型是线性磨损方程Q=kp，磨损体积率与载荷成比例关系。这个模型被认为

是Archard的磨损定律，尽管它的基本形式首先由Holm发表。这个模型是建立在实验观察

基础上的，用公式表示

（2）

引入磨损率K使理论和实验的结果相吻合。Holm把它作为一个常量，来表述耐磨原子的数

量。在Archard’s工作中，提出了微凸体相互作用导致磨损粒子产生的可能性[4]。然而，那

并不是唯一可能的解释。Lim和Ashby通过计算认为分层和塑性磨损机理是主要的。对于刚，

他们建议用以下的值：

然而，对于一个接触实际的K值需要通过实验的方法得到。对于工程应用来说，相对于磨

损体积，对磨损深度更感兴趣。这里Archard想通过接触面积A[4]把方程（2）都分开，给

出

H为磨损深度（m），k为空间磨损率，p为接触压力。磨损过程可以看作是个动态的过程，

它的预测是个初始值的问题。从而磨损模型可以通过一个不同的方程来表述，在线性的情况

下，方程（2）可以用下面公式表示

（3）

3.有限元磨损模拟程序

3.1有限元理论

用有限元方法计算磨损的主要工作是计算接触应力。工程结构被离散成许多单元，单

元与单元之间用节点连接起来。在有限元中，可以把单元中一些物理量(位移，

温度等)通过多项式分段拟合来近似描述，通过结点位古老肉的做法 移来表示[5]。可以同时使

3/11

用不同的单元类型，复杂的载荷和边界条件。在结构分析中，把自由的程度定义为节点位移。

把各个单元按原来的结构重新连接起来，形成整体的有限元方程[D]{}={F}

D为结构刚度阵或总体刚度阵，{}为结构的节点位移列阵或变形矢量，{F}为结构节点的

载荷列阵。节点的应力通过变形得出。商用有限元软件ANSYS可以处理几种材料和非线性

问题，例如塑性，粘弹性，摩擦等[6]。

有限元磨损计算包括解决全面的接触问题，实际接触面积是未知的，所以分析为非线性。

在此情况下，使用了点对面单元。有限元软件配备了多种工具，提高了非线性分析程序，

但选择参数时要格外细心。

3.2磨损模拟程序

磨损模拟程序的流程图包括一系列的结构的解决步骤和额外的计算，如图2所示。

初始的参数用来定义模型的几何尺寸，载荷，约束和磨损模型参数以及单元和材料数据。

并开发了一些特殊的小程序生成FE模型和自动的定义载荷和约束。每个几何和载荷的情况

都要很好的离散化。模型中采用较多的单元会得到更精确的结果，但是会增加计算时间和占

用硬盘空间。

4/11

当得到了FEA迭代应力后，就心理素质情况 可以确定接触区域。然后可以确定每个接触单元的状态。

用接触单元的结点坐标确定接触区域。根据接触区域节点应力来确定压力的分配。

用欧拉方法求磨损相对于时间的积分。系统参数在磨损模拟的每一个载荷步中被假定为

常数，根据下面所描述的磨损模型，在每个单元节点上对磨损深度产生影响。

（4）

h

j,n

是在节点n,j上的磨损增量。根据已知的应力分布，可以估计出节点的磨损增量

h

j,n

（m）。如果在一个磨损步中某个节点处的磨损增量过大，模拟的结果就可能变得不确定。主题班会方案

接触区域就可能出现间隙。因此就根据以往的经验事先假定和引入了一个允许的最大磨损增

量。进行时间较短的模拟实验来调整这个值使它尽可能大。最初的节点磨损增量经过一个连

续的时间增量t（s）计算出来。每个求解步j中的时间因数Mj可以通过下式估测

h

in,j,max

为节点磨损增量h

in,j，n

的最大值。那个求解步j中实际的时间间隙被定为M

j

t。

然后模型的几何参数改变，根据方程（4）把节点移到新的位置。这种方法，没有用连续的

时间步长，提高了FEA计算速度。在每个求解步后保存输出数据是很重要的。这就要求如果

分析由于某种原因中断时，快速数据回顾和保存前一个步中的数据。

3.3FEA结果确认

也许确认FEA结果最可靠的方法是将它与已知的分析结果相对比。ANSYS软件也配备了

能量误差估计方法，建立在FEA结构分析结果是一个从一个单元到另一个单元连续的位移域

[6]。为了得到更满意的应力，单元节点的应力被平均。节点应力误差矢量相应的被估测出

来，作为单元和整个模型能量误差估计的基础。当每一个单元能量误差都相等时，这个离散

化的模型才是最有效的。

3.4Sphere-on-planeFE模型

Pin-on-disk结构如图3所示，用上面所示的FEA方法作了分析。在这种情况下塑性变

形和摩擦对压力的分配被看作是可以忽略的。一个顶端为球形半径R=5mm的pin通过一个轴

向均质的Sphere-on-planeFE模型来表达。这个结构采用了两维实体结构单元PlANE42。

接触面采用二维点对面接触单元CONTACT48。Pin和disk都被认为是弹性模量E=210GPa，

泊松比=0.3的钢做成的。两个负载为FN=21N和FN=50N。接触区单元的X方向的尺寸根据

载荷不同分别为25m和32.5m。ANSYS接触刚度参数定为KN=5X107N/m。

为了检验模型的正确性，接触应力的分配分别通过FEA和赫兹公式进行了计算。（图4

所示）

5/11

E*=E/2（1-2）为常态弹性模量。这个模型中忽略了弹性变形和摩擦。用FE和用Hertz

分别计算的结果相差不超过5%。

4.实验程序

一个半径R=5mm的pin在disk上滑动，压力分别为FN=21N和FN=50N。disk和pin分

别硬化到刚度分别为HV=4.6GPa和HV=3GPa。Hertz计算的最大的接触应力被假定为在弹性

范围内。实验设备可以直线测量磨损深度和摩擦力矩（图5）。实验中的滑动速度为v=25mm/s。

实验结果图6所示，做了两个载荷下的实验。图6（C）所示的磨损率由图6所示两个

实验的的平均磨损深度决定。使用了如下方程

6/11

体积磨损增量通过公式[7]计算

i≥1是样品点的数量，s=0.15m为滑动距离增量。Disk的硬度比pin的硬度大。平均磨

损率通过试验数据估测，滑动距离为s=3m和s=4.5m[图6（c）]。当滑动距离分别为s=3m

和s=5m，这些值和偏差分别为在F=21N时K=（1.250.44）*10-13Pa-1和K=（2.261.44）

*10-13Pa-1，在F=50N时K=（1.330.54）*10-13Pa-1和K=（2.011.21）*10-13Pa-1。所测得

7/11

的摩擦系数的平均值为fFr=0.70.2。

通过Archard建议的方法分析了接触产生的温度[8]。假定摩擦热流量qn（W/m2）没

有流入接触体内部。qn=fFrpv。P为接触区域的平均接触压力，计算时没有考虑系统变形

通过一个圆形的热源加热，计算出了两个接触体的平均和最大温度。没有考虑辐射和转化。

使用了下面的公式[9]

接触温度由[8]计算

无量纲参数Pe=0.5v=vr0/2a0为数。计算所得温度没有超过6K[图6（d）].在载

荷F=21N和F=50N时，实验速度和压力分别为v=0.2,…,1.6,p=0.007,…,0.27和

v=0.3,…,2.0,p=0.009,…0.37。在平均赫兹压力的基础上，计算相对于初始情况更高

8/11

的无量纲压力。认为磨损机理为分层或黏着磨损，图1，得到了线形磨损规律，公式（3）。

5.磨损模拟结果

假定符合线性磨损定律，FEA磨损模拟结果可以用磨损率相对于磨损距离的关系来表

示。给定几何尺寸和载荷时，如果ks不变，磨损的深度不会改变。

5.1Sphere-on-plane滑动接触

用以上的模型进行FE磨损模拟，假定符合公式（3）的线性磨损定律。模拟磨损时估测

空间磨损率K=（1.330.54）*10-13Pa-1。最大的磨损增量定为hlim=0.1m。求解步磨损增

量通过公式（3）h=kps计算出来。Disk在实验中被假定为硬度较高，因此只有pin发

生磨损。

磨损率做为常数来处理。对比实验结果和FEA磨损模拟结果，如图9所示。粗线表示在

平均k值下的磨损曲线，细线表示磨损误差的影响。在摩擦过程中表面状况在不断地发生变

化，从而会影响实际的磨损过程和磨损率。

5.2锥体对锥体有变形时旋转接触

分析了锥体和座旋转接触并且锥体和座都有磨损[11]。在磨损过程中，接触面大小不发

生变化，磨损特性接近于直线。引用了空间磨损率K=（1.330.54）*10-13Pa-1。锥角=60，

轴向载荷F=40N，接触体的材料为纲。

5.3锥体对锥体并且没有变形旋转接触

在这种情况下，随着磨损的进行，接触面积在不断增加。在不同的角=10和=20

时进行了分析。载荷和磨损率和上面相同，锥角=60（图11）。

9/11

5.4锥体对圆环面旋转接触

和上面相同的锥体和圆环形的座旋转接触，接触面半径R1=37.5mm，同样用FEM进行了

分析，假定相同的材料与载荷，图12所示。

磨损曲线由在磨损过程中给定几何参数和载荷情况下的实际接触面积的变化来决定。

一个相关的问题是磨损模拟结果的正确性需要验证。每次试验所得的磨损率的值是行政部职责 不同

的[11]。目前的pin-on-disk磨损试验得到的磨损率的值有41%的偏差。因此值得怀疑的

是这种模拟的结果能否预测特定接触系统的寿命。结果证明，是可以的。可以对比不同的设

计解决方法和可选择的方案。锥体和环形座的角度查=20对磨损的影响，可以通过对比

有变形的环形接触和锥体对圆环接触的摩擦行为，考虑磨损率的偏差，如图13所示。

10/11

6.讨论与结论

FEA的正确性依赖于模型的离散化。细的结点划分能得到更好的结果，但是会增加计算

时间和占用更大的硬盘空间。经常用一些额外的方法使结果的正确性更高。接触问题在有限

元中是一个非线性问题。有限元模型的离散化和给定载荷和约束下的刚度直接与迭代过程的

收敛有关。建立较好的模型结构需要建立在经验的基础上。

时间步长积分是模拟结果可靠的一个重要参数。太长的时间步长会导致不确定的结果

和FEA程序不能收敛。太少的时间间隔占用太多的计算时间。因此开发了一个简单的模拟时

间步的最优化程序，在确定最大磨损增量的基础上估计每个求解步的持续时间。

在模拟过程中，必须考虑磨损的机理和预知它的改变。对于刚可以使用Lim和Ashby

磨损图。

假定线性磨损定律是正确的，一个给定几何尺寸和载荷的FEA模拟结果可以通过磨损

率-滑动距离来表示。

因为模型简化和输入数据的偏差，FEA磨损模拟的结果需要在一个相对的值上来估算，

对比不同的设计选择方案，而不是用来预测绝对的磨损寿命。

线性磨损定律：磨损率与所加压力成线性关系

接触问题：接触面积变化而产生的非线性以及由接触压力分布变化而产生的非线性，也

有摩擦作用产生的非线性。

载荷步：作用在给定时间间隔内的一系列载荷

子步：载荷步中的时间点，在这些点，可以求得中间解。

时间步长：两个连续的子步间的时间差称为时间步长或时间增量。

塑性：在某种载荷下，材料产生永久变形的材料特性，对大多数工程材料来说，当其盈

利低于比例极限时，应力-应变关系是线性的。

11/11

（范文素材和资料部分来自网络，供参考。可复制、编制，期待你的好评与关注）