CN104239627B - 一种干滑动摩擦热‑应力‑磨损分步耦合的模拟方法 - Google Patents

Abstract

本发明涉及一种干滑动摩擦热‑应力‑磨损分步耦合的模拟方法，步骤为：将摩擦热‑应力‑磨损耦合过程分步；对热‑应力耦合分析进行初始化设置；采用瞬态热传导分析程序依次完成N个增量步的热传导分析，获得输出温度场；对应力‑磨损分析进行初始化设置；设置有限元模型应力分析的单元类型，进行应力分析，获得输出接触压力场；根据温度场、接触压力场及两接触面接触节点的相对滑移速率，计算磨损量并确定其空间方向；计算接触压力偏差；根据磨损量和方向，修正接触节点的位移，更新有限元模型；判断接触压力的偏差是否超过容差，更新接触压力；判断模拟过程是否结束。本发明计算精度和计算效率较高且适用于干滑动摩擦热‑应力‑磨损强耦合问题的模拟。

Description

一种干滑动摩擦热-应力-磨损分步耦合的模拟方法

技术领域

本发明涉及一种干滑动摩擦耦合模拟方法，特别是关于一种干滑动摩擦热-应力-磨损分步耦合的模拟方法。

背景技术

干滑动摩擦界面存在摩擦热-应力-摩擦磨损耦合现象。摩擦热主要产生于两个接触面之间的摩擦接触表层，摩擦热的热流密度与切向摩擦应力和相对滑移速度有关。摩擦热会直接导致接触面的温度迅速升高，迅速升高的温度致使接触面发生热变形，从而最终使得接触压力在接触面上重新不均匀的分布，在接触压力和相对滑移速度的综合作用下，接触面逐渐磨损。由于接触面上接触压力分布不均以及接触面不同位置处相对滑移速度存在差异，因此，接触面不同位置处的磨损量不同。这种接触面上磨损量的不均匀性导致接触压力的再次重新分布，进而改变切向摩擦应力的分布和热流密度的大小。另一方面，接触面温度的改变还会显著改变两接触材料的摩擦特性，比如摩擦系数和磨损系数，除此之外，还会改变材料的本构特性，比如屈服极限和硬度等。因此，接触面的摩擦热(温度场)、应力(接触压力场)和磨损(磨损量)之间是相互耦合的。

这种耦合行为难以求得解析解。在采用数值方法进行求解数值解时，如果采用完全耦合分析方法，每个时间增量步都考虑摩擦热-应力-磨损三者之间的双向相互耦合效应，则计算效率极低，无法求解较长时间的摩擦热-应力-磨损耦合问题，但如果采用顺序耦合方法，即先分析温度场，再进行应力-磨损耦合分析，忽略应力和磨损对热传导单方向的耦合项，则牺牲了计算精度，仅适用于摩擦热-应力-磨损的弱耦合问题。

发明内容

针对上述问题，本发明的目的是提供一种计算精度和计算效率较高，且适用于强耦合问题的干滑动摩擦热-应力-磨损分步耦合的模拟方法。

为实现上述目的，本发明采取以下技术方案：一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，包括以下步骤：1)时域分步：将摩擦热-应力-磨损耦合过程分步，初始化步数计数器i＝1和仿真总时间计时器t＝0，计算初始接触压力p₀＝p_i；2)进入第i步热-应力-磨损顺序耦合，其包括以下步骤：(1)热-应力耦合分析初始化：预设第i步耦合过程的总时间步ΔT，增量步的步数为N，ΔT为耦合过程的总时间T′与增量步的步数N之比，初始化第i步耦合过程的当前时间增量Δt_i＝0和增量步计数器j＝1；(2)热传导分析：在已有的瞬态热传导分析程序ABAQUS中依次完成N个增量步的热传导分析，每个增量步对应输出一个温度场，用于应力-磨损分析；(3)应力-磨损分析初始化：开始进入应力-磨损耦合分析之前，再次初始化增量步计数器j＝1，将步骤2)结束时的状态作为本步应力-磨损分析的初始状态，并记录初始接触压力p₀；(4)应力分析：进入第j个增量步的应力分析，设置有限元模型应力分析的单元类型，定义边界条件、载荷和接触条件，将热传导分析的温度场作为热载荷施加在有限元模型中，应用瞬态热传导分析程序ABAQUS对有限元模型进行应力分析，获得并输出接触压力场；(5)确定磨损量：根据热传导分析得到的温度场、应力分析得到的接触压力场以及两接触面接触节点的相对滑移速率，计算两接触面接触节点的磨损量，并且确定其空间方向；(6)接触压力偏差计算：计算当前增量步的接触压力p_j与初始接触压力p₀之间的相对偏差e(p_j,p₀)；(7)更新网格：根据两接触面接触节点的磨损量和方向，修正该节点的位移，更新有限元模型；(8)判断相对偏差e(p_j,p₀)是否超过容差：考察接触压力的相对偏差e(p_j,p₀)是否达到容许值TOL，如果偏差没有超过容许值TOL，转到步骤(4)，进入下一个增量步的应力分析；否则，抛出异常，终止本步耦合过程，并更新本步耦合过程的当前时间增量Δt_i＝jΔt，进入步骤3)；如果所有增量步都没有抛出异常，则在最后一个增量步结束后更新总时间增量Δt_i＝NΔt，进入步骤3)；3)更新总时间t：根据第i步耦合过程的当前时间增量Δt_i更新总时间t＝t+Δt_i，判断当前总时间t和耦合过程总时间T′的大小，如果t＜T′，则更新计数器i＝i+1和初始接触压力p₀＝p_j，进入步骤2)，开始下一步摩擦热-应力-磨损耦合分析过程，否则，模拟过程结束。

所述步骤2)中，瞬态热传导分析的生热热流边界由摩擦功率引起，其摩擦热流q表达式为：式中，q为接触节点的摩擦热流；μ为接触节点的摩擦系数；p₀为接触面的初始接触压力；为两接触表面在接触节点的相对滑动速率；摩擦系数μ为接触面温度T和初始接触压力p₀等因素的函数：μ＝μ(T,p₀,…)；根据两构件的具体结构和相对运动求取；摩擦热流q通过热流用户子程序DFLUX施加在其中一构件的底面上，两接触面接触节点温度T、接触节点到构件边沿的距离均从计算结果文件中读取，用户子程序DFLUX包括以下步骤：①从瞬态热传导分析程序ABAQUS中跳转至用户子程序DFLUX；②初始化常量和系数；③根据两接触面接触节点坐标计算回转半径根据两接触面接触节点温度T和平均接触压力计算摩擦系数μ；④根据回转半径和摩擦系数μ计算两接触面接触节点处的摩擦热流⑤返回瞬态热传导分析程序ABAQUS。

所述步骤2)中磨损量的计算采用广义非线性Archard磨损公式：Δh＝κ·p·Δs；式中，Δh为两接触面接触节点的磨损增量；κ为广义Archard磨损公式的磨损系数；p为两接触面接触节点的接触压力；Δs为两接触表面在接触节点的相对滑动位移增量；磨损系数κ是温度T、接触压力p和接触节点的相对滑动速率等因素的函数：

所述步骤2)中磨损量的空间方向的确定方法如下：两接触面接触节点分为两类，一类是位于接触面边缘的节点，简称边界节点；另一类是位于接触面内部的节点，简称非边界节点，对于非边界节点的磨损量方向由节点法方向确定，节点法方向直接从有限元分析结果文件中读取。

所述步骤2)中磨损量的确定在瞬态热传导分析程序ABAQUS中通过网格移动用户子程序UMESHMOTION实现，两接触面接触节点的温度T、接触压力p、节点到构件边缘的距离和增量步时间长度Δt均从计算结果文件中读取，用户子程序UMESHMOTION包括以下步骤：①初始化常量和系数；②访问计算结果文件，读取两接触面接触节点的接触压力p和温度T；③计算两接触面接触节点的磨损量，确定两接触面接触节点的磨损方向；④出两接触面接触节点磨损信息，用于绘制磨损量云图。

所述步骤2)中，接触节点位移的修正和有限元模型的更新采用任意拉格朗日欧拉法实现；需要修正位移的接触节点为接触压力与相对切向滑移均不为零的接触节点。

本发明由于采取以上技术方案，其具有以下优点：1、本发明由于是在顺序耦合方法的基础上，当接触压力场的变化达到一定量级时，修正接触压力，重新计算摩擦热流，然后继续开始顺序耦合分析过程，依次反复进行直至完成整个耦合过程，因此，这种摩擦热-应力-磨损分步耦合的模拟方法兼顾了计算精度和计算效率，也适合于摩擦热-应力-磨损的强耦合过程。2、本发明由于是基于现有的商用有限元软件平台，因此，设置简单、容易上手且便于推广。3、本发明由于磨损计算的程序在商用有限元软件的应力分析算法基础上二次开发完成，完全嵌入商用有限元软件之中，免去了数据传输和交互，因此，大大提高了计算效率。综上所述，本方法具有非常广泛的适用范围，可以用于各类干滑动摩擦的摩擦热-应力-磨损耦合问题的模拟。

附图说明

图1是本发明的总体流程示意图；

图2是本发明的第i步摩擦热-应力-磨损耦合的计算框图；

图3是本发明的用户子程序UMASFL的流程示意图；

图4是本发明的用户子程序DFLUX的流程示意图；

图5是本发明的用户子程序UMESHMOTION的流程示意图；

图6是本发明的销-盘有限元模型示意图；

图7是本发明的销接触节点示意图；

图8是本发明中销-盘磨损试验的温度分布示意图；

图9是本发明与现有完全耦合方法和顺序耦合方法得到的销-盘磨损试验的接触面温度对比示意图，其中，实线表示节点N101501采用完全耦合方法得到的销-盘磨损试验的接触面温度，虚线表示节点N101501采用顺序耦合方法得到的销-盘磨损试验的接触面温度，点划线表示节点N101501采用分步耦合方法得到的销-盘磨损试验的接触面温度；

图10是本发明与现有完全耦合方法和顺序耦合方法得到的销的接触压力对比示意图，其中，实线表示节点N101501采用完全耦合方法得到的销的接触压力，虚线表示节点N101501采用顺序耦合方法得到的销的接触压力，点划线表示节点N101501采用分步耦合方法得到的销的接触压力；

图11是本发明与现有完全耦合方法和顺序耦合方法得到的销的磨损量对比示意图，其中，实线表示节点N101501采用完全耦合方法得到的销的磨损量，虚线表示节点N101501采用顺序耦合方法得到的销的磨损量，点划线表示节点N101501采用分步耦合方法得到的销的磨损量。

具体实施方式

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

本发明提供一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，实施本发明方法之前，首先建立干滑动摩擦的有限元模型，有限元模型包括两个具有相对干滑动摩擦的构件，约束不能使两构件产生相对干滑动摩擦的自由度，初始化两构件之间的作用载荷，初始化两构件的材料参数，材料参数包括不同温度对应的比热容、导热系数、线膨胀系数、杨氏模量、密度和泊松比。本发明方法包括以下步骤(如图1、图2所示)：

1)时域分步：将摩擦热-应力-磨损耦合过程分步，初始化步数计数器i＝1和仿真总时间计时器t＝0，计算初始接触压力p₀＝p_i。

2)进入第i步热-应力-磨损顺序耦合，其包括以下步骤(如图2所示)：

(1)热-应力耦合分析初始化：预设第i步耦合过程的总时间ΔT，增量步的步数为N，ΔT为耦合过程的总时间T′与增量步的步数N之比，初始化第i步耦合过程的当前时间增量Δt_i＝0和增量步计数器j＝1。

(2)热传导分析：在已有的瞬态热传导分析程序ABAQUS中依次完成N个增量步的热传导分析，每个增量步对应输出一个温度场，用于应力-磨损分析；热传导分析步骤为：

①初始化有限元模型的温度，设置有限元模型中两构件的单元类型、生热热流边界和散热边界；

②对每个增量步进行瞬态热传导分析；

③获得并输出温度场。

如果其中一单元类型允许其对应的构件的材料脱离网格运动，那么材料的运动利用用户子程序UMASFL实现(如图3所示)，以销-盘磨损实验为例，用户子程序UMASFL包括以下步骤：

①从瞬态热传导分析程序ABAQUS中跳转至用户子程序UMASFL；

②初始化构件的常量和系数，如角速度材料密度ρ等；

③根据销-盘表面接触节点坐标计算回转半径如果不是旋转件，直接读取该节点速度即可；对于销-盘实验，知道销-盘表面接触节点坐标和旋转角速度即可计算出节点速度；

④根据回转半径计算销-盘表面接触节点处的质量流速

⑤返回瞬态热传导分析程序ABAQUS。

瞬态热传导分析的生热热流边界由摩擦功率引起，其摩擦热流q表达式为：

式中，q为接触节点的摩擦热流；μ为接触节点的摩擦系数；p₀为接触面的初始接触压力；为两接触表面在接触节点的相对滑动速率。

摩擦系数μ为接触面温度T和初始接触压力p₀等因素的函数：

μ＝μ(T,p₀,…) (2)

根据两构件的具体结构和相对运动求取。

摩擦热流q通过热流用户子程序DFLUX施加在其中一构件的底面上，节点温度T、回转半径均可以从计算结果文件中读取，以销-盘实验为例，用户子程序DFLUX包括以下步骤(如图4所示)：

①从瞬态热传导分析程序ABAQUS中跳转至用户子程序DFLUX；

②初始化常量和系数，如常量π和平均接触压力等；

③根据销-盘表面接触节点坐标计算回转半径根据销-盘表面接触节点温度T和平均接触压力计算摩擦系数μ；

④根据回转半径和摩擦系数μ计算销-盘表面接触节点处的摩擦热流

⑤返回瞬态热传导分析程序ABAQUS。

(3)应力-磨损分析初始化：开始进入应力-磨损耦合分析之前，再次初始化增量步计数器j＝1，将步骤2)结束时的状态作为本步应力-磨损分析的初始状态，并记录初始接触压力p₀。

(4)应力分析：进入第j个增量步的应力分析，设置有限元模型应力分析的单元类型，定义边界条件、载荷和接触条件，将热传导分析的温度场作为热载荷施加在有限元模型中，应用瞬态热传导分析程序ABAQUS对有限元模型进行应力分析，获得并输出接触压力场。

(5)确定磨损量：根据热传导分析得到的温度场、应力分析得到的接触压力场以及两接触面接触节点的相对滑移速率，计算两接触面接触节点的磨损量，并且确定其空间方向。其中，磨损量的计算采用广义非线性Archard磨损公式：

Δh＝κ·p·Δs (3)

式中，Δh为两接触面接触节点的磨损增量；κ为广义Archard磨损公式的磨损系数；p为两接触面接触节点的接触压力；Δs为两接触表面在接触节点的相对滑动位移增量。

磨损系数κ是温度T、接触压力p和接触节点的相对滑动速率等因素的函数：

两接触面接触节点磨损量的空间方向的确定方法如下：两接触面接触节点可以分为两类，一类是位于接触面边缘的节点，简称边界节点；另一类是位于接触面内部的节点，简称非边界节点。对于非边界节点的磨损量方向可以由节点法方向确定，节点法方向可以直接从有限元分析结果文件中读取。而对于边界节点，由于不存在法方向，因此，采用边界节点与其厚度方向对应点的连线方向作为磨损量方向。

磨损量的确定在瞬态热传导分析程序ABAQUS中通过网格移动用户子程序UMESHMOTION实现，两接触面接触节点的温度T、接触压力p、节点到构件边缘的距离和增量步时间长度Δt均可以从计算结果文件中读取，以销-盘实验为例，用户子程序UMESHMOTION包括以下步骤(如图5所示)：

①初始化常量和系数，如常量π和角速度等；

②访问计算结果文件，读取销-盘表面接触节点的接触压力p和温度T；

③计算销-盘表面接触节点的磨损量，确定销-盘表面接触节点的磨损方向；

④输出销-盘表面接触节点磨损信息，用于绘制磨损量云图；

(6)接触压力偏差计算：计算当前增量步的接触压力p_j与初始接触压力p₀之间的相对偏差e(p_j,p₀)。

(7)更新网格：根据两接触面接触节点的磨损量和方向，修正接触压力与相对切向滑移均不为零的两接触面接触节点的位移，更新有限元模型，更新网格是通过更新两接触面接触节点位移，然后重画磨损区域网格实现，本发明采用了任意拉格朗日欧拉法(ALE)。ALE需要设置ALE作用域和ALE约束条件，其中，ALE的作用域是所有包含需要修正位移的接触节点的单元，ALE的约束范围为所有需要修正位移的接触节点，ALE的约束类型为位移约束或者速度约束。

以速度约束为例，若接触节点i的磨损速率为其在全局坐标系下的方向向量为接触节点i的局部坐标系三个轴的归一化方向矢量分别为则该接触节点的ALE速度约束矢量在局部坐标系中的分量{v_xi,v_yi,v_zi}^T，由下面公式计算：

(8)判断相对偏差e(p_j,p₀)是否超过容差：考察接触压力的相对偏差e(p_j,p₀)是否达到容许值TOL，如果偏差没有超过容许值TOL，转到步骤(4)，进入下一个增量步的应力分析；否则，抛出异常，终止本步耦合过程，并更新本步耦合过程的当前时间增量Δt_i＝jΔt，进入步骤3)；如果所有增量步都没有抛出异常，则在最后一个增量步结束后更新当前时间增量Δt_i＝NΔt，进入步骤3)。

3)更新总时间t：根据第i步耦合过程的当前时间增量Δt_i更新总时间t＝t+Δt_i，判断当前总时间t和耦合过程总时间T′的大小，如果t＜T′，则更新计数器i＝i+1和初始接触压力p₀＝p_j，进入步骤2)，开始下一步摩擦热-应力-磨损耦合分析过程，否则，模拟过程结束。

实施例：以销-盘磨损试验为例，采用本发明对销-盘之间发生的干滑动摩擦的摩擦热-应力-磨损分布耦合进行模拟。

如图6所示，实施本发明之前，首先建立销-盘的有限元模型，包括一销夹具1，一夹持在销夹具1的正下方且其下端面与销夹具1下端面平齐的销2，一转动连接在一轴承3上的盘夹具4，盘夹具4的正上方可拆卸的紧固连接一盘5，销2的下端面和盘5的上端面摩擦接触，轴承3坐落在一轴承座6上。销2和盘5之间沿y轴方向施加的载荷为500N，盘绕y轴负方向的转速为600rpm，耦合全过程的时间设定为10s，销2和盘5的材料参数如表1和表2所示，模拟过程包括以下步骤(如图1、图2所示)：

1)时域分步：将摩擦热-应力-磨损耦合过程分步，初始化步数计数器i＝1和仿真总时间计时器t＝0，利用已有的瞬态热传导分析程序ABAQUS计算初始接触压力p₀＝p_i。

2)进入第i步热-应力-磨损顺序耦合，其包括以下步骤(如图2所示)：

(1)热-应力耦合分析初始化：预设第i步耦合过程的总时间ΔT＝1s，然后设置增量步的步数N＝1000，初始化第i步耦合过程的当前时间增量Δt_i＝0，增量步计数器j＝1。

(2)热传导分析：在已有的瞬态热传导分析程序ABAQUS中依次完成N个增量步的热传导分析，每个增量步对应输出一个温度场，用于应力-磨损分析。

销-盘有限元模型的初始温度为25℃，销的单元类型采用3维热传导单元DC3D8，热传导分析不允许网格空间运动，为了模拟盘的转动，盘的单元类型采用欧拉单元DCC3D8D，该单元可以允许材料脱离网格运动，材料的运动可以利用用户子程序UMASFL实现(如图3所示)，用户子程序UMASFL的步骤如下：

①瞬态热传导分析程序ABAQUS中跳转至用户子程序UMASFL；

②初始化构件的常量和系数，如角速度材料密度ρ等；

③根据销-盘表面接触节点坐标计算节点的回转半径

④根据回转半径计算销-盘表面接触节点处的质量流速

⑤返回瞬态热传导分析程序ABAQUS。

瞬态热传导分析的热流边界由摩擦功率引起，其摩擦热流q表达式为：

摩擦系数μ为温度T和初始接触压力p₀的函数：

相对滑动速率可以近似表示为：

式中，n为转速10rps，r为接触点到回转轴线的距离。

摩擦热流q通过热流用户子程序DFLUX施加在销的底面上，节点温度T、回转半径均可以在子程序计算结果文件中读取，用户子程序DFLUX包括以下步骤(如图4所示)：

①从瞬态热传导分析程序ABAQUS中跳转至用户子程序UMASFL；

②初始化常量和系数，如常量π和平均接触压力等；

③根据销-盘表面接触节点坐标计算回转半径根据销-盘表面接触节点温度T和平均接触压力计算摩擦系数μ；

④根据回转半径和摩擦系数μ计算销-盘表面接触节点处的摩擦热流

⑤返回瞬态热传导分析程序ABAQUS。

(3)应力-磨损分析初始化：开始进入应力-磨损耦合分析之前，再次初始化增量步计数器j＝1，将前一步步骤2)结束时的状态作为本步应力-磨损分析的初始状态，并记录初始接触压力p₀。

(4)应力分析：进入第j个增量步的应力分析，设置销-盘有限元模型应力分析的单元类型为C3D8R，按照实例工况定义边界条件、载荷和接触条件，将热传导分析的温度场作为热载荷施加在销-盘有限元模型中，应用瞬态热传导分析程序ABAQUS对销-盘有限元模型进行应力分析，获得并输出接触压力场。

(5)确定磨损量：根据热传导分析得到的温度场、应力分析得到的接触压力场以及销-盘表面接触节点的相对滑移速率，计算销-盘表面接触节点的磨损量，并且确定其空间方向。其中，销-盘表面接触节点的磨损量Δh的计算采用广义非线性Archard磨损公式：

Δh＝κ·p·Δs，

材料的磨损系数κ是温度和接触压力函数：

κ＝1908+115.4T-1662p

-0.3834T²-3.762T·p+360.7p²；

相对滑动位移增量可以近似表示为：

Δs＝2πnrΔt；

磨损量的确定在瞬态热传导分析程序ABAQUS中通过网格移动用户子程序UMESHMOTION实现，销-盘表面接触节点的温度T、接触压力p、回转半径和增量步时间长度Δt均可以从计算结果文件中读取，用户子程序UMESHMOTION包括以下步骤(如图5所示)：

①初始化常量和系数，如常量π和角速度等；

②问计算结果文件，读取销-盘表面接触节点的接触压力p和温度T；

③算销-盘表面接触节点的磨损量，确定销-盘表面接触节点的磨损方向；

④输出销-盘表面接触节点磨损信息，用于绘制磨损量云图；

销-盘表面接触节点磨损量的空间方向的确定方法如下：销-盘表面接触节点可以分为两类，一类是位于接触面边缘的节点，简称边界节点；另一类是位于接触面内部的节点，简称非边界节点。对于非边界节点的磨损量方向可以由节点法方向确定，节点法方向可以直接从有限元分析结果文件中读取。而对于边界节点，由于不存在法方向，因此，采用边界节点与其厚度方向对应点的连线方向作为磨损量方向。

如图7所示，图中销截面的节点分为边界节点7和非边界节点8。位于截面圆周上的节点为边界节点7，而位于圆周以内的节点为非边界节点8。对于非边界节点8，其磨损量的方向为内法线方向。而对于边界节点7，以节点N101114为例，由于不存在内法线，可以定义其磨损方向为节点N101114与节点N102114的连线方向。如果忽略变形影响，所有节点的磨损方向都可以定义为y轴正方向。

(6)接触压力偏差计算：计算当前增量步的接触压力p_j与初始接触压力p₀之间的相对偏差e(p_j,p₀)：

式中，ci为销截面的接触节点编号；cN为销截面接触节点总数。

(7)更新网格：根据销-盘表面接触节点的磨损量和方向，修正该节点的位置，并利用瞬态热传导分析程序ABAQUS的任意拉格朗日欧拉法ALE技术实现有限元模型的更新。在本实施例中，由于盘5表面接触节点的磨损是间歇的，并且盘5的磨损系数远小于销2，因此仿真时将其忽略。设置任意拉格朗日欧拉法ALE的作用域为销2底层的3层单元，ALE约束选择为位移约束，作用在销2底面所有节点上，用以模拟这些接触节点的磨损。

(8)判断相对偏差e(p_j,p₀)是否超过容差：设置容许值TOL＝10％。考察接触压力的相对偏差e(p_j,p₀)是否达到容许值TOL，如果偏差没有超过容许值TOL，转到步骤(4)，进入下一个增量步的应力分析；否则，抛出异常，终止本步耦合过程，并更新本步耦合过程的总时间增量Δt_i＝jΔt，进入步骤3)；如果所有增量步都没有抛出异常，则在最后一个增量步结束后更新总时间增量Δt_i＝NΔt，同样进入步骤3)。

3)更新总时间t：根据第i步耦合过程的当前时间增量Δt_i更新总时间t＝t+Δt_i，判断当前总时间t和耦合过程总时间T′的大小，如果t＜T′，则更新计数器i＝i+1和初始接触压力p₀＝p_j，进入步骤2)，开始下一步摩擦热-应力-磨损耦合分析过程，否则，模拟过程结束。

该方法模拟得到的10s时温度场分布如图8所示，呈现近似轴对称特点，最高温度约为63℃，出现在销盘接触的环形区域。

为了说明该方法的特点，将上述实施例的模拟结果与利用全耦合方法、顺序耦合方法计算得到的结果对比如下：

如图9～图11所示，采用本发明方法与现有完全耦合方法和顺序耦合方法得到的接触面中心节点的节点温度、接触压力和磨损量均十分接近，本发明得到的曲线位于完全耦合与顺序耦合曲线之间。完全耦合方法考虑了热传导、应力以及磨损两两之间的双向耦合效应，顺序耦合方法忽略了应力和磨损对热传导单方向的耦合作用，本发明通过反馈机制改善了顺序耦合方法的缺点，计算精度介于完全耦合方法和顺序耦合方法之间。

在相同的商用机HP Compaq Elite 8300 Convertible Microtower上，本发明计算用时147min，完全耦合方法计算用时1318min，顺序耦合方法计算用时106min。本发明计算用时介于完全耦合方法和顺序耦合方法之间，更接近于顺序耦合方法，因此，使用本发明保留了顺序耦合方法计算高效的优点。

上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (10)

1.一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，包括以下步骤：

1)时域分步：将摩擦热-应力-磨损耦合过程分步，初始化步数计数器i＝1和仿真总时间计时器t＝0，计算初始接触压力p₀＝p_i；

2)进入第i步热-应力-磨损顺序耦合，其包括以下步骤：

(1)热-应力耦合分析初始化：预设第i步耦合过程的总时间步ΔT，增量步的步数为N，ΔT为耦合过程的总时间T′与增量步的步数N之比，初始化第i步耦合过程的当前时间增量Δt_i＝0和增量步计数器j＝1；

(2)热传导分析：在已有的瞬态热传导分析程序ABAQUS中依次完成N个增量步的热传导分析，每个增量步对应输出一个温度场，用于应力-磨损分析；

(3)应力-磨损分析初始化：开始进入应力-磨损耦合分析之前，再次初始化增量步计数器j＝1，将步骤2)结束时的状态作为本步应力-磨损分析的初始状态，并记录初始接触压力p₀；

(4)应力分析：进入第j个增量步的应力分析，设置有限元模型应力分析的单元类型，定义边界条件、载荷和接触条件，将热传导分析的温度场作为热载荷施加在有限元模型中，应用瞬态热传导分析程序ABAQUS对有限元模型进行应力分析，获得并输出接触压力场；

(5)确定磨损量：根据热传导分析得到的温度场、应力分析得到的接触压力场以及两接触面接触节点的相对滑移速率，计算两接触面接触节点的磨损量，并且确定其空间方向；

(6)接触压力偏差计算：计算当前增量步的接触压力p_j与初始接触压力p₀之间的相对偏差e(p_j,p₀)；

(7)更新网格：根据两接触面接触节点的磨损量和方向，修正该节点的位移，更新有限元模型；

(8)判断相对偏差e(p_j,p₀)是否超过容差：考察接触压力的相对偏差e(p_j,p₀)是否达到容许值TOL，如果偏差没有超过容许值TOL，转到步骤(4)，进入下一个增量步的应力分析；否则，抛出异常，终止本步耦合过程，并更新本步耦合过程的当前时间增量Δt_i＝jΔt，进入步骤3)；如果所有增量步都没有抛出异常，则在最后一个增量步结束后更新总时间增量Δt_i＝NΔt，进入步骤3)；

3)更新总时间t：根据第i步耦合过程的当前时间增量Δt_i更新总时间t＝t+Δt_i，判断当前总时间t和耦合过程总时间T′的大小，如果t＜T′，则更新计数器i＝i+1和初始接触压力p₀＝p_j，进入步骤2)，开始下一步摩擦热-应力-磨损耦合分析过程，否则，模拟过程结束。

2.如权利要求1所述的一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，其特征在于：所述步骤2)中，瞬态热传导分析的生热热流边界由摩擦功率引起，其摩擦热流q表达式为：

$q={\mathrm{\μp}}_0\stackrel{\·}{\mathrm{\γ}},$

式中，q为接触节点的摩擦热流；μ为接触节点的摩擦系数；p₀为接触面的初始接触压力；为两接触表面在接触节点的相对滑动速率；

摩擦系数μ为接触面温度T和初始接触压力p₀因素的函数：

μ＝μ(T,p₀)；

根据两构件的具体结构和相对运动求取；

摩擦热流q通过热流用户子程序DFLUX施加在其中一构件的底面上，两接触面接触节点温度T、接触节点到构件边沿的距离均从计算结果文件中读取，用户子程序DFLUX包括以下步骤：

①从瞬态热传导分析程序ABAQUS中跳转至用户子程序DFLUX；

②初始化常量和系数；

③根据两接触面接触节点坐标计算回转半径根据两接触面接触节点温度T和平均接触压力计算摩擦系数μ；

④根据回转半径和摩擦系数μ计算两接触面接触节点处的摩擦热流

⑤返回瞬态热传导分析程序ABAQUS。

3.如权利要求1所述的一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，其特征在于：所述步骤2)中磨损量的计算采用广义非线性Archard磨损公式：

Δh＝κ·p·Δs；

式中，Δh为两接触面接触节点的磨损增量；κ为广义Archard磨损公式的磨损系数；p为两接触面接触节点的接触压力；Δs为两接触表面在接触节点的相对滑动位移增量；

磨损系数κ是温度T、接触压力p和接触节点的相对滑动速率因素的函数：

$\mathrm{\κ}=\mathrm{\κ}(T,p,\stackrel{\·}{\mathrm{\γ}}).$

4.如权利要求2所述的一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，其特征在于：所述步骤2)中磨损量的计算采用广义非线性Archard磨损公式：

Δh＝κ·p·Δs；

式中，Δh为两接触面接触节点的磨损增量；κ为广义Archard磨损公式的磨损系数；p为两接触面接触节点的接触压力；Δs为两接触表面在接触节点的相对滑动位移增量；

磨损系数κ是温度T、接触压力p和接触节点的相对滑动速率因素的函数：

$\mathrm{\κ}=\mathrm{\κ}(T,p,\stackrel{\·}{\mathrm{\γ}}).$

5.如权利要求1或2或3或4所述的一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，其特征在于：所述步骤2)中磨损量的空间方向的确定方法如下：两接触面接触节点分为两类，一类是位于接触面边缘的节点，简称边界节点；另一类是位于接触面内部的节点，简称非边界节点，对于非边界节点的磨损量方向由节点法方向确定，节点法方向直接从有限元分析结果文件中读取。

6.如权利要求1或2或3或4所述的一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，其特征在于：所述步骤2)中磨损量的确定在瞬态热传导分析程序ABAQUS中通过网格移动用户子程序UMESHMOTION实现，两接触面接触节点的温度T、接触压力p、节点到构件边缘的距离和增量步时间长度Δt均从计算结果文件中读取，用户子程序UMESHMOTION包括以下步骤：

①初始化常量和系数；

②访问计算结果文件，读取两接触面接触节点的接触压力p和温度T；

③计算两接触面接触节点的磨损量，确定两接触面接触节点的磨损方向；

④出两接触面接触节点磨损信息，用于绘制磨损量云图。

7.如权利要求5所述的一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，其特征在于：所述步骤2)中磨损量的确定在瞬态热传导分析程序ABAQUS中通过网格移动用户子程序UMESHMOTION实现，两接触面接触节点的温度T、接触压力p、节点到构件边缘的距离和增量步时间长度Δt均从计算结果文件中读取，用户子程序UMESHMOTION包括以下步骤：

①初始化常量和系数；

②访问计算结果文件，读取两接触面接触节点的接触压力p和温度T；

③计算两接触面接触节点的磨损量，确定两接触面接触节点的磨损方向；

④出两接触面接触节点磨损信息，用于绘制磨损量云图。

8.如权利要求1或2或3或4所述的一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，其特征在于：所述步骤2)中，接触节点位移的修正和有限元模型的更新采用任意拉格朗日欧拉法实现；需要修正位移的接触节点为接触压力与相对切向滑移均不为零的接触节点。

9.如权利要求5所述的一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，其特征在于：所述步骤2)中，接触节点位移的修正和有限元模型的更新采用任意拉格朗日欧拉法实现；需要修正位移的接触节点为接触压力与相对切向滑移均不为零的接触节点。

10.如权利要求6所述的一种干滑动摩擦热-应力-磨损分步耦合的模拟方法，其特征在于：所述步骤2)中，接触节点位移的修正和有限元模型的更新采用任意拉格朗日欧拉法实现；需要修正位移的接触节点为接触压力与相对切向滑移均不为零的接触节点。