一种基于有限元的热-机械-磨损耦合分析数值模拟方法

摘要

本发明涉及一种基于有限元的热-机械-磨损耦合分析数值模拟方法，其包括以下步骤：建立有限元模型；将热-机械-磨损耦合过程离散化；利用商用有限元软件进行热-机械耦合分析；输出接触对的温度场和接触压力场；计算磨损量，确定磨损方向；输出节点磨损量和磨损方向；修正节点位移，更新有限元模型；判断所有增量步是否均已完成，如果未完成则返回循环操作，直到所有增量步均已完成为止，如果完成则输出所有增量步的接触节点所对应的接触压力场分布、温度场分布和磨损量分布。由于本发明将商用有限元软件作为平台，因此简单而可靠，磨损量的计算在有限元软件内部完成，快速而高效。本发明可以广泛应用于涉及干滑动摩擦的结构分析中。

说明书

技术领域

本发明涉及一种数值模拟方法，特别是关于一种适用于零件干滑动摩擦过程中模 拟接触对之间接触压力、摩擦生热和磨损过程的基于有限元的热-机械-磨损耦合分析 数值模拟方法。

背景技术

机械零部件的损坏往往与零件之间摩擦磨损状况相关。与其他类型的摩擦磨损现 象相比，干滑动摩擦磨损会对零件造成更大的损害。现有技术中，机械零部件的生产 通常首先根据实际需要设计出产品，对产品加上载荷进行实验，如果实验过程中产品 达不到设计要求重新设计产品，然后再进行实验，直到产品合格为止，此过程不仅耗 时较长，而且返工率比较高，导致生产成本较大和效率较低。利用数值模拟方法准确 地模拟干滑动过程中出现的现象可以预测产品的使用寿命，改进零部件的设计方案， 节省研发成本。

接触对之间的干滑动摩擦是一种强烈的热-机械-摩擦磨损耦合现象。摩擦热主要 产生于两个接触面之间的摩擦接触表层，接触表层局部产生的热流密度与切向摩擦应 力和相对滑移速度有关。热流密度通常很高，会导致接触面的温度迅速升高，并通过 热传导的方式迅速传播开去。在较高的温度和温度梯度下，接触面特别是摩擦表面由 于热变形而发生翘曲，引起接触压力的重新分布。在接触压力和相对滑动下，接触表 面逐渐磨损。由于接触对之间接触压力分布不均匀以及不同位置处相对滑移速度存在 差异，因此接触面不同位置处的磨损量也不相等。这种接触面上磨损的不均匀性会导 致接触压力的重新分布，进而影响切向摩擦应力的分布，并且导致热流密度发生改变。 另一方面，温度的改变也会显著改变接触对材料的摩擦特性，例如改变摩擦系数和磨 损系数，还会改变材料的本构特性，导致材料屈服极限、硬度下降等。因此，接触对 之间的温度场、接触压力场和磨损量是强烈相互耦合的。因此，只有将它们耦合在一 起研究，才能准确地预测接触对之间的摩擦磨损现象。

接触问题本身就是强非线性问题，又由于存在上述多场耦合现象，因此几乎不可 能求得该问题的解析解。在目前的工业应用中，往往忽略多场耦合的影响以及非线性 因素，分别单独求解接触状况、温度分布和磨损量，这使得计算结果十分不准确。工 业应用中的另外一种手段是通过实验获得相关数据，然而实验成本较高，而且实验过 程中只能获得某些点的信息，无法得到数据在整个区域的分布。随着有限元技术的进 步，尤其是大型商用有限元软件突飞猛进的发展，强非线性问题和多场耦合问题可以 利用这些软件得到很满意的数值解答。然而遗憾的是，目前的有限元软件只能解决热- 机械耦合问题，却无法模拟接触对之间的磨损过程。

发明内容

针对上述问题，本发明的目的是提供一种基于有限元的热-机械-磨损耦合分析数 值模拟方法，能够更加完整地模拟干滑动磨损过程，精确得到干滑动摩擦过程中接触 对之间接触压力场分布、温度场分布和磨损量分布。

为实现上述目的，本发明采取以下技术方案：一种基于有限元的热-机械-磨损耦 合分析数值模拟方法，其包括以下步骤：1）对构成接触对的两个零件分别建立有限元 模型；2）将热-机械-磨损耦合过程所持续的时间历程离散化为N个增量步；3）对每 个增量步分别进入4）～8）进行操作；4）在每一增量步对接触对的有限元模型进行 热-机械耦合分析，得到接触对的所有接触节点所对应的接触压力和温度，所述接触节 点包括边界节点和非边界节点；5）输出每一增量步所对应接触对的所有接触节点的接 触压力和温度，用以绘制相应云图；6）根据接触压力场分布及接触节点的相对滑移速 度，计算相互接触零件各自接触面上的各个接触节点的磨损量，并确定各个接触节点 的磨损方向；7）输出步骤6中得到的所有接触节点的磨损量和磨损方向，用以绘制磨 损量云图；8）对接触压力与相对切向滑移均不为零的接触节点的位移进行修正，更新 有限元模型；9）判断所有增量步是否均已完成，如果未完成则返回步骤4），循环操 作，直到所有增量步均已完成为止；如果完成则输出所有增量步的接触节点所对应的 接触压力场分布、温度场分布和磨损量分布。

所述步骤2）采用显式欧拉向前积分法则对热-机械-磨损耦合过程所持续的时间 历程进行离散化，离散化过程中增量步数目和增量步大小满足以下公式：

$\underset{l=1}{\overset{N}{\Sigma }}\Delta t_l=T$

式中，Δt_l为第l个增量步时间大小，N为整个耦合过程增量步数目，T为耦合过程的 总时间长度。

所述步骤6）中磨损量的计算采用Archard磨损法则，其公式为：

Δh＝k_d·p_c·Δs

式中，Δh为节点磨损量，k_d为材料磨损系数，p_c为节点接触压力，Δs为节点相对滑 移位移。

所述步骤6）中磨损方向的确定方式：对于非边界节点的，磨损方向为其所在表 面的法方向；对于边界节点，磨损方向为其与厚度方向对应点的连线方向。

所述步骤8）对接触压力与相对切向滑移均不为零的接触节点的位移进行修正采 用ALE方法，ALE的作用域是所有包含需要修改位移的节点的单元，ALE的约束范围为 所有需要修改位移的节点，ALE的约束类型为位移约束或者速度约束。

本发明由于采取以上技术方案，其具有以下优点：1、本发明将零件的干滑动摩擦 过程中接触对的压力、温度和磨损耦合在一起进行计算，能够更加准确地预测接触对 之间的摩擦磨损现象。2、本发明将热-机械-磨损耦合过程所持续的时间历程离散化为 N个增量步，并分别对每一个增量步进行接触压力场、温度场和磨损量的计算，可以 根据实际需要得到很多点的信息，进而得到数据在整个区域的分布，因此根据预先得 到零件的接触压力场分布、温度场分布和磨损量分布，预测产品的使用寿命，改进零 部件的设计方案，节省研发成本，减低返工率，有效提高生产效率。3、由于本发明使 用了ALE方法，根据接触面节点的磨损量和磨损方向修正节点位移，实现了热-机械- 磨损的完全耦合，并且可以很好的控制显式欧拉方法的增量步长，因此能够提高仿真 的计算效率和计算精度。4、本发明由于利用商用有限元软件计算热-机械耦合过程， 设置简单，容易上手，避免了繁杂的计算，具有较高的精度和准确性，而且磨损量的 计算也可以利用有限元软件的子程序计算，不需要第三方软件，免去了数据传输和交 互，加快了访问结果数据的速度，提高了效率。本发明使用范围广，可以广泛应用于 各类涉及接触的磨损仿真中，特别是适用于指导各种机械零件生产制造或汽车制动磨 损分析中。

附图说明

图1是本发明立方体垫片与环形旋转盘接触磨损过程示意图；

图2是本发明模拟热-机械-磨损耦合过程的流程示意图；

图3是本发明热-机械-磨损耦合仿真效果示意图，图3（a）是立方体垫片Pad.N25 接触压力随时间变化示意图，横坐标为时间Time，单位为s，纵坐标为压力CPress， 单位为MPa；图3（b）是立方体垫片Pad.N25温度随时间变化示意图，横坐标为时间 Time，单位为s，纵坐标为温度Temperature，单位为℃；图3（c）是环形旋转盘Dish.N11 温度随时间变化示意图，横坐标为时间Time，单位为s，纵坐标为温度Temperature， 单位为℃；图3（d）是立方体垫片Pad.N25累积磨损量随时间变化示意图，横坐标为 时间Time，单位为s，纵坐标为磨损量Wear，单位为um。

具体实施方式

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

如图1所示，本发明实施例中构成接触对的两个零件包括一个立方体垫片1（Pad） 和一个环形旋转盘2（Dish），立方体垫片1在设定的1Mpa压力作用下压紧环形旋转 盘2，环形旋转盘2以设定的恒定转速15r/s进行旋转，两零件接触面之间无润滑， 摩擦系数为0.37。本实例中在环形旋转盘2旋转一周的周期内，立方体垫片1始终处 于接触磨损状态，而环形旋转盘2处于间歇接触磨损状态，因此在整个过程中立方体 垫片1的磨损量远大于环形旋转盘2，所以本实例中只考虑立方体垫片1的磨损，而 将环形旋转盘2当作刚体处理。如果相互接触的两个零件均有较大磨损量，则需要同 时考虑两者的磨损，每个零件磨损的设置方法和立方体垫片1的磨损设置方法相同。 本发明将以模拟此工况下立方体垫片1和环形旋转盘2这对接触对在1s时间内的热- 机械-磨损耦合行为为实施例对本发明进行详细说明。

如图2所示，本发明基于有限元的热-机械-磨损耦合行为的数值模拟方法，包括 以下步骤：

1、对构成接触对的两个零件分别建立有限元模型。

如图1所示，采用现有的商用有限元软件例如ABAQUS（以此为例，但不限于此， 可以根据实际需要选择使用其它的有限元软件）分别对立方体垫片1和环形旋转盘2 建立有限元模型，具体建模方法为现有技术在此不再赘述。本发明实施例中对立方体 垫片1和环形旋转盘2建立有限元模型所选择的材料参数如表1所示，两个有限元分 析模型的接触特性如表2所示，表1和表2中的单位制采用国际单位制。

表1材料参数

表2接触特性

2、将热-机械-磨损耦合过程所持续的时间历程离散化为N个增量步。

本发明采用显式欧拉向前积分法则对热-机械-磨损耦合过程所持续的时间历程进 行离散化，离散化后认为每个增量步中热、机械和磨损三者的场变量保持不变。离散 化过程中增量步数目和增量步大小满足以下公式：

$\underset{l=1}{\overset{N}{\Sigma }}\Delta t_l=T---\left(1\right)$

式中，Δt_l为第l个增量步时间大小，N为整个耦合过程增量步数目，T为耦合过 程的总时间长度。Δt_l的设置可以根据所需要分析问题进行确定，对于耦合过程中场变 量变化剧烈的时间段可以适当增加增量步数目并且减少增量步大小，以便提高仿真结 果的精度。本实施例的整个时间历程为0～2s，其中0～1s为预加载过程，1～2s为热 -机械-磨损耦合过程。本实例的热-机械-磨损耦合总时间为T＝1s，假设总增量步数量 N＝400，且每相邻两个增量步时间间隔相等，则Δt_l＝Δt＝0.0025s,l＝1,2,…N。

3、对每个增量步分别进入4～8进行操作。

4、在每一增量步对接触对的有限元模型进行热-机械耦合分析，得到接触对的所 有接触节点所对应的接触压力和温度；其中，接触节点可以分为两类，一类是位于接 触面边缘的节点（简称边界节点），另一类是位于接触面内部的节点（简称非边界节点）。

5、输出每一增量步所对应接触对的所有接触节点的接触压力和温度，用以绘制相 应云图。

6、根据接触压力场分布及接触节点的相对滑移速度，计算相互接触零件各自接触 面上的各个接触节点的磨损量，并确定各个接触节点的磨损方向。本发明实施例中立 方体垫片1的磨损量远大于环形旋转盘2，因此本实施例中只计算立方体垫片1各个 接触节点的磨损量，但是不限于此，可以根据实际需要选择某个零件或者所有零件进 行磨损量的计算。

磨损量的计算采用Archard磨损法则，其公式为：

Δh＝k_d·p_c·Δs   （2）

式中，Δh为节点磨损量，k_d为材料磨损系数，p_c为节点接触压力，Δs为节点相 对滑移位移。

公式（2）两边同时处以增量步时间，则为：

$\stackrel{·}{h}=k_d·p_c·\stackrel{·}{\gamma }---\left(3\right)$

式中，为节点的磨损速率，为节点相对滑移速率。

公式（2）和公式（3）中，k_d为已知参数可以由实验得到，本实施例中k_d取1×10^-11Pa^-1； p_c和Δs均可以从有限元分析结果文件中读取；可以计算得到，计算公式为：Δt为增量步时间大小。根据公式（2）或（3）可以求得一个增量步中所有接触节点处 的磨损量或者磨损速率。

磨损方向的确定方式：对于非边界节点的，磨损方向为其所在表面的法方向，节 点法方向可以直接从有限元分析结果文件中读取；对于边界节点，由于不存在法方向， 磨损方向为其与厚度方向对应点的连线方向。

从图1可以获得本发明实施例中立方体垫片1的节点编号的规律（以此为例，不 限于此，可以根据实际需要进行编号），编号为｛21～25，46，50，71，75，96，100， 121～125｝(部分节点编号图中未显示)的16个节点属于边界节点，其磨损方向如表3 所示。编号为｛47—49，72—74，97—99｝(部分节点编号图中未显示)的9个节点属 于非边界节点，其磨损方向直接取节点所在表面的法方向即可。

表3边界节点磨损方向

7、输出步骤6中得到的所有接触节点的磨损量和磨损方向，用以绘制磨损量云图。

8、对接触压力与相对切向滑移均不为零的接触节点的位移进行修正，更新有限元 模型。

磨损仿真过程中要求节点位移的更新不得引入附加节点力，否则节点位移的更新 只相当于附加了接触材料的变形，而不是接触对的磨损。为了实现接触材料的磨损而 不引入附加节点力，本发明采用了任意拉格朗日欧拉法（ALE）。通过ALE方法，可以 实现在不改变有限元分析结果中其他变量的前提下只改变接触节点的位移。ALE的作 用域是所有包含需要修改位移的节点的单元，ALE的约束范围为所有需要修改位移的 节点，ALE的约束类型为位移约束或者速度约束。本发明以速度约束为例进行说明。

本发明的实施例中，由于环形旋转盘2表面的磨损是间歇的，因此只需计算立方 体垫片1的磨损，而不需要考虑环形旋转盘2的磨损。设置ALE的作用域为立方体垫 片1表层与环形旋转盘2接触的16个单元，ALE约束作用的节点为立方体垫片1与环 形旋转盘2接触的25个节点｛21～25，46～50，71～75，96～100，121～125｝，并称 该集合为NSet_ALE，ALE约束类型设置为速度约束。

节点速度约束的大小取决于节点磨损速率和磨损方向。对于节点i∈NSet_ALE， 其磨损速率为归一化磨损方向向量为局部坐标系三个轴的归一化方向矢量分 别为则该节点的ALE速度约束{v_xi,v_yi,v_zi}^T由公式（4）确定：

$v_{\mathrm{xi}}={\stackrel{·}{h}}_i{\overrightarrow{l}}_i·{\overrightarrow{x}}_i$

$v_{\mathrm{yi}}={\stackrel{·}{h}}_i{\overrightarrow{l}}_i·{\overrightarrow{y}}_i---\left(4\right)$

$v_{\mathrm{zi}}={\stackrel{·}{h}}_i{\overrightarrow{l}}_i·{\overrightarrow{z}}_i$

式中，和可以从有限元计算结果文件中读取。

9、判断所有增量步是否均已完成，如果未完成则返回步骤4，循环操作，直到所 有增量步均已完成为止；如果完成则输出所有增量步的接触节点所对应的接触压力场 分布、温度场分布和磨损量分布。根据所有增量步的接触节点所对应的接触压力场、 温度场分布和磨损量分布可以用于对机械零件的生产进行指导或者对汽车制动磨损分 析中。

如图3（a）～（d）所示，Pad.N25表示立方体垫片1的25节点，此节点在整个 仿真过程中接触应力始终最大，温度始终最高，磨损量也最大。观察Pad.N25的接触 压力和温度随时间变化曲线可以发现，刚出现滑动摩擦时，接触压力分布很不均匀， 因此Pad.N25接触压力大于背面的1MPa压力。由于摩擦生热，Pad.N25温度急速上升， 热膨胀加剧了接触压力的不均匀性，Pad.N25的接触压力升高。随着时间的推移，热 传导和边界散热作用逐渐明显，Pad.N25的温度逐渐趋于稳定，这时磨损量对接触压 力分布的作用越来越明显，整个接触面上接触压力的分布随着磨损的增加而逐渐趋于 均匀，因此Pad.N25处的接触压力逐渐减小。Pad.N25处累积磨损量主要取决于接触 压力的大小，随着接触压力的减小，磨损量增加速度减缓。

Dish.N11为环形旋转盘2上的一个节点，其初始时刻和立方体垫片1接触，当N11 和立方体垫片1接触滑动时，其温度升高，当它脱离接触时，温度下降。在仿真过程 中，Dish.N11周期性与立方体垫片1接触和分离，因此其温度呈周期性上升。由于Dish 转速为15r/s，因此如图3（c）所示，在1s的仿真过程中Dish.N11温升出现15个阶 梯。

上述实施例中，接触对的温度场和接触压力场的云图可以直接利用有限元软件绘 制，而接触对的磨损量的分布需要根据输出结果利用Matlab软件绘制。累积磨损量云 图的绘制需要对增量步磨损量累加，其公式为：其中，h(i,j)为接 触节点i在第j个增量步时的累积磨损量，Δh(i,l)为接触节点i在第l个增量步的磨损量。

上述各实施例仅用于说明本发明，其中实施方法的各步骤等都是可以有所变化的， 凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护 范围之外。