用有限元的方法模拟滑动摩擦磨损

用有限元的方法模拟滑动摩擦磨损
摘要
磨损往往是影响产品寿命的一个主要因素。

因此磨损预测就成为工程的一个重要部分。

这篇论文介绍了用有限元软件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），r0为接触半径（m），F N为载荷，H为接触对中较软材料的硬度（P a），v为相对滑动速度（m/s）。

表一列出了Lim和Ashby使用的方程和参数。

图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]。

可以同时使
用不同的单元类型，复杂的载荷和边界条件。

在结构分析中，把自由的程度定义为节点位移。

把各个单元按原来的结构重新连接起来，形成整体的有限元方程[D]{μ}={F}
D为结构刚度阵或总体刚度阵，{μ}为结构的节点位移列阵或变形矢量，{F}为结构节点的载荷列阵。

节点的应力通过变形得出。

商用有限元软件ANSYS可以处理几种材料和非线性问题，例如塑性，粘弹性，摩擦等[6]。

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

在此情况下，使用了点对面单元。

有限元软件配备了多种工具，提高了非线性分析程序，但选择参数时要格外细心。

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

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

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

每个几何和载荷的情况都要很好的离散化。

模型中采用较多的单元会得到更精确的结果，但是会增加计算时间和占用硬盘空间。

当得到了FEA迭代应力后，就可以确定接触区域。

然后可以确定每个接触单元的状态。

用接触单元的结点坐标确定接触区域。

根据接触区域节点应力来确定压力的分配。

用欧拉方法求磨损相对于时间的积分。

系统参数在磨损模拟的每一个载荷步中被假定为常数，根据下面所描述的磨损模型，在每个单元节点上对磨损深度产生影响。

（4）
Δh j,n是在节点n,j上的磨损增量。

根据已知的应力分布，可以估计出节点的磨损增量Δh j,n（m）。

如果在一个磨损步中某个节点处的磨损增量过大，模拟的结果就可能变得不确定。

接触区域就可能出现间隙。

因此就根据以往的经验事先假定和引入了一个允许的最大磨损增量。

进行时间较短的模拟实验来调整这个值使它尽可能大。

最初的节点磨损增量经过一个连续的时间增量Δt（s）计算出来。

每个求解步j中的时间因数M j可以通过下式估测
Δh in,j,max为节点磨损增量Δh in,j，n的最大值。

那个求解步j中实际的时间间隙被定为M jΔt。

然后模型的几何参数改变，根据方程（4）把节点移到新的位置。

这种方法，没有用连续的时间步长，提高了FEA计算速度。

在每个求解步后保存输出数据是很重要的。

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

3.3 FEA 结果确认
也许确认FEA结果最可靠的方法是将它与已知的分析结果相对比。

ANSYS软件也配备了能量误差估计方法，建立在FEA结构分析结果是一个从一个单元到另一个单元连续的位移域[6]。

为了得到更满意的应力，单元节点的应力被平均。

节点应力误差矢量相应的被估测出来，作为单元和整个模型能量误差估计的基础。

当每一个单元能量误差都相等时，这个离散化的模型才是最有效的。

3.4 Sphere-on-plane FE模型
Pin-on-disk结构如图3所示，用上面所示的FEA方法作了分析。

在这种情况下塑性变形和摩擦对压力的分配被看作是可以忽略的。

一个顶端为球形半径R=5mm的pin通过一个轴向均质的Sphere-on-plane FE模型来表达。

这个结构采用了两维实体结构单元PlANE42。

接触面采用二维点对面接触单元CONTACT48。

Pin和disk都被认为是弹性模量E=210GPa，泊松比µ=0.3的钢做成的。

两个负载为F N=21N和F N=50N。

接触区单元的X方向的尺寸根据载荷不同分别为25μm和32.5μm。

ANSYS接触刚度参数定为KN=5X107N/m。

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

（图4所示）
E*=E/2（1-µ2）为常态弹性模量。

这个模型中忽略了弹性变形和摩擦。

用FE和用Hertz 分别计算的结果相差不超过5%。

4.实验程序
一个半径R=5mm 的pin在disk上滑动，压力分别为F N=21N和F N=50N。

disk和pin分别硬化到刚度分别为HV=4.6GPa和HV=3GPa。

Hertz计算的最大的接触应力被假定为在弹性范围内。

实验设备可以直线测量磨损深度和摩擦力矩（图5）。

实验中的滑动速度为v=25mm/s。

实验结果图6所示，做了两个载荷下的实验。

图6（C）所示的磨损率由图6所示两个实验的的平均磨损深度决定。

使用了如下方程
体积磨损增量通过公式[7]计算
i≥1是样品点的数量，Δs=0.15m为滑动距离增量。

Disk的硬度比pin的硬度大。

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

当滑动距离分别为s=3m 和s=5m，这些值和偏差分别为在F=21N时K=（1.25±0.44）*10-13Pa-1和K=（2.26±1.44）*10-13Pa-1，在F=50N时K=（1.33±0.54）*10-13Pa-1和K=（2.01±1.21）*10-13Pa-1。

所测得
的摩擦系数的平均值为f Fr=0.7±0.2。

通过Archard建议的方法分析了接触产生的温度[8]。

假定摩擦热流量q n（W/m2）没有流入接触体内部。

q n=f Fr pv。

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。

在平均赫兹压力的基础上，计算相对于初始情况更高
的无量纲压力。

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

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

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

5.1 Sphere-on-plane滑动接触
用以上的模型进行FE磨损模拟，假定符合公式（3）的线性磨损定律。

模拟磨损时估测空间磨损率K=（1.33±0.54）*10-13Pa-1。

最大的磨损增量定为Δh lim=0.1µm。

求解步磨损增量通过公式（3）Δh=kpΔs计算出来。

Disk在实验中被假定为硬度较高，因此只有pin发生磨损。

磨损率做为常数来处理。

对比实验结果和FEA磨损模拟结果，如图9所示。

粗线表示在平均k值下的磨损曲线，细线表示磨损误差的影响。

在摩擦过程中表面状况在不断地发生变化，从而会影响实际的磨损过程和磨损率。

5.2锥体对锥体有变形时旋转接触
分析了锥体和座旋转接触并且锥体和座都有磨损[11]。

在磨损过程中，接触面大小不发生变化，磨损特性接近于直线。

引用了空间磨损率K=（1.33±0.54）*10-13Pa-1。

锥角α=60°，轴向载荷F=40N，接触体的材料为纲。

5.3 锥体对锥体并且没有变形旋转接触
在这种情况下，随着磨损的进行，接触面积在不断增加。

在不同的角Δα=10´和Δα=20´时进行了分析。

载荷和磨损率和上面相同，锥角α=60°（图11）。

5.4锥体对圆环面旋转接触
和上面相同的锥体和圆环形的座旋转接触，接触面半径R1=37.5mm，同样用FEM进行了分析，假定相同的材料与载荷，图12所示。

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

一个相关的问题是磨损模拟结果的正确性需要验证。

每次试验所得的磨损率的值是不同的[11]。

目前的pin-on-disk磨损试验得到的磨损率的值有±41%的偏差。

因此值得怀疑的是这种模拟的结果能否预测特定接触系统的寿命。

结果证明，是可以的。

可以对比不同的设计解决方法和可选择的方案。

锥体和环形座的角度查Δα=20´对磨损的影响，可以通过对比有变形的环形接触和锥体对圆环接触的摩擦行为，考虑磨损率的偏差，如图13所示。

6.讨论与结论
FEA的正确性依赖于模型的离散化。

细的结点划分能得到更好的结果，但是会增加计算时间和占用更大的硬盘空间。

经常用一些额外的方法使结果的正确性更高。

接触问题在有限元中是一个非线性问题。

有限元模型的离散化和给定载荷和约束下的刚度直接与迭代过程的收敛有关。

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

时间步长积分是模拟结果可靠的一个重要参数。

太长的时间步长会导致不确定的结果和FEA程序不能收敛。

太少的时间间隔占用太多的计算时间。

因此开发了一个简单的模拟时间步的最优化程序，在确定最大磨损增量的基础上估计每个求解步的持续时间。

在模拟过程中，必须考虑磨损的机理和预知它的改变。

对于刚可以使用Lim和Ashby 磨损图。

假定线性磨损定律是正确的，一个给定几何尺寸和载荷的FEA模拟结果可以通过磨损率-滑动距离来表示。

因为模型简化和输入数据的偏差，FEA磨损模拟的结果需要在一个相对的值上来估算，对比不同的设计选择方案，而不是用来预测绝对的磨损寿命。

线性磨损定律：磨损率与所加压力成线性关系
接触问题：接触面积变化而产生的非线性以及由接触压力分布变化而产生的非线性，也有摩擦作用产生的非线性。

载荷步：作用在给定时间间隔内的一系列载荷
子步：载荷步中的时间点，在这些点，可以求得中间解。

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

塑性：在某种载荷下，材料产生永久变形的材料特性，对大多数工程材料来说，当其盈利低于比例极限时，应力-应变关系是线性的。