CN104200034B - 一种干滑动摩擦热‑应力‑磨损顺序耦合的模拟方法 - Google Patents

Abstract

本发明涉及一种干滑动摩擦热‑应力‑磨损顺序耦合的模拟方法，步骤为：在摩擦热‑应力‑磨损耦合过程的时间轴上将该过程离散为N个增量步，瞬态热传导分析、应力分析和磨损分析采用相同时间轴和增量步；建立干滑动摩擦的有限元模型，进行瞬态热传导分析，获得输出温度场；设置有限元模型应力分析的单元类型，进行应力分析，获得输出接触压力场；根据温度场、接触压力场及两接触面接触节点的相对滑移速率，计算两接触面接触节点的磨损量并确定其空间方向；根据磨损量和方向，修正接触节点的位移，更新有限元模型；判断应力‑磨损耦合过程的模拟是否完成。本发明计算精度损失小且计算效率较高，能广泛应用于各类干滑动摩擦的摩擦热‑应力‑磨损弱耦合问题的模拟。

Description

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

技术领域

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

背景技术

干滑动摩擦界面存在摩擦热-应力-磨损耦合现象，其中，摩擦热主要产生于两个接触面之间的摩擦接触表层，摩擦热的热流密度与切向摩擦应力和相对滑移速度有关。摩擦热会直接导致接触面的温度迅速升高，迅速升高的温度致使接触面发生热变形，从而最终使得接触压力在接触面上重新不均匀的分布，在接触压力和相对滑移速度的综合作用下，接触面逐渐磨损。由于接触面上接触压力分布不均以及接触面不同位置处相对滑移速度存在差异，因此，接触面不同位置处的磨损量不同。这种接触面上磨损量的不均匀性导致接触压力的再次重新分布，进而改变切向摩擦应力的分布和热流密度的大小。另一方面，接触面温度的改变还会显著改变两接触材料的摩擦特性，比如摩擦系数和磨损系数，除此之外，还会改变材料的本构特性，比如屈服极限和硬度等。因此，接触面的摩擦热(温度场)、应力(接触压力场)和磨损(磨损量)之间是相互耦合的，这种耦合行为几乎不可能求得解析解，而采用数值方法进行求解数值解时，如果考虑摩擦热-应力-磨损三者之间的双向耦合效应，那么计算的效率极低，无法求解较长时间的摩擦热-应力-磨损耦合问题。

发明内容

针对上述问题，本发明的目的是提供一种计算精度损失小且计算效率高的干滑动摩擦热-应力-磨损顺序耦合的模拟方法。

为实现上述目的，本发明采取以下技术方案：一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，包括以下步骤：1)时域离散：假设摩擦热-应力-磨损耦合过程的时间轴为(t₀,t_N)，在该时间轴上将摩擦热-应力-磨损耦合过程离散为N个增量步，瞬态热传导分析、应力分析和磨损分析采用相同的时间轴和增量步；2)瞬态热传导分析：建立干滑动摩擦的有限元模型，在步骤1)的时间轴(t₀,t_N)上采用瞬态热传导分析程序进行瞬态热传导分析，其步骤为：①初始化有限元模型的温度，选定瞬态热传导分析程序，设置有限元模型中两构件的单元类型、生热热流边界和散热边界；②对每个增量步进行瞬态热传导分析；③获得并输出温度场；3)应力分析：设置有限元模型应力分析的单元类型，定义边界条件、载荷和接触条件，将步骤2)输出的温度场作为热载荷施加在有限元模型中，并对有限元模型进行应力分析，获得并输出接触压力场；4)确定磨损量：根据热传导分析得到的温度场、应力分析得到的接触压力场以及两接触面接触节点的相对滑移速率，计算两接触面接触节点的磨损量，并且确定其空间方向；5)更新网格：根据两接触面接触节点的磨损量和方向，修正该接触节点的位移，更新有限元模型；6)判断应力-磨损耦合过程的模拟是否完成，若增量步达到N，则模拟过程结束；否则，进入下一个增量步，回到步骤2)。

所述步骤1)中的瞬态热传导分析中的数值积分方法采用向后欧拉方法，应力分析中的数值积分方法采用Newmark-β方法，磨损分析中的数值积分方法采用显式欧拉算法。

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

所述步骤4)中的磨损量大小的计算采用广义非线性Archard磨损公式：Δh＝κ·p·Δs，式中，Δh为两接触面接触节点的磨损增量；κ为广义Archard磨损公式的磨损系数；p为两接触面接触节点的接触压力；Δs为两接触表面在接触点的相对滑动位移增量；上式的左右两边同时除以增量步时间，则得到如下公式：式中，为两接触面接触节点的磨损速率；对于参数κ，为温度和接触压力等因素的函数，其表达式为：由实验测量得到；p和Δs从有限元分析结果的文件中读取；从下式计算得到：

式中，Δt为增量步时间大小。

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

所述步骤4)中的磨损量大小和方向的确定通过网格移动用户子程序UMESHMOTION实现，包括以下步骤：①初始化常量和系数；②读取两接触面接触节点的接触压力p和温度T；③计算两接触面接触节点的磨损量，确定两接触面接触节点的磨损方向；④输出两接触面接触节点磨损信息，用于绘制磨损量云图。

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

本发明由于采取以上技术方案，其具有以下优点：1、本发明由于在实际工程应用，如制动过程中，摩擦热-应力属于弱耦合问题，可以忽略应力和磨损对摩擦热传导单方向的耦合项，本发明采用的顺序耦合方法模拟摩擦热-应力-磨损耦合行为，即将摩擦热-应力-磨损弱耦合问题分解为瞬态热传导分析和应力-磨损耦合分析，因此，在计算精度损失较小的情况下极大地提高计算效率。2、本发明由于基于现有商用的有限元软件平台，因此设置简单、容易上手且便于推广。3、本发明由于采用的瞬态热传导分析可以直接使用商用软件的算法，应力-磨损耦合分析在商用有限元软件的应力分析算法的基础上二次开发完成，磨损的计算完全嵌入有限元软件之中，免去了数据传输和交互，因此，使得计算效率大大提高。本发明适用范围广，可以广泛应用于各类干滑动摩擦的摩擦热-应力-磨损弱耦合问题的模拟。

附图说明

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

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

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

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

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

图6是本发明的销底部的局部放大示意图；

图7是图2中Ⅰ处局部放大示意图；

图8是采用顺序耦合方法和完全耦合方法得到的销-盘磨损试验的温度分布对比示意图，其中图(a)为采用完全耦合方法得到的销-盘磨损试验的温度分布，图(b)为采用顺序耦合方法得到的销-盘磨损试验的温度分布；

图9是本发明采用顺序耦合方法和完全耦合方法得到的销-盘磨损试验的接触节点温度分布对比示意图，其中，表示节点N101101采用完全耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101105采用完全耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101109采用完全耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101113采用完全耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101501采用完全耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101101采用顺序耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101105采用顺序耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101109采用顺序耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101113采用顺序耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101501采用顺序耦合方法得到的销-盘磨损试验的接触节点温度；

图10是本发明采用顺序耦合方法和完全耦合方法得到的销-盘磨损试验的接触节点接触压力对比示意图，其中，表示节点N101101采用完全耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点N101105采用完全耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点N101109采用完全耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点N101113采用完全耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点N101501采用完全耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点N101101采用顺序耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点N101105采用顺序耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点N101109采用顺序耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点N101113采用顺序耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点N101501采用顺序耦合方法得到的销-盘磨损试验的接触节点接触压力；

图11是本发明采用顺序耦合方法和完全耦合方法得到的销-盘磨损试验的接触节点磨损量对比示意图，其中，表示节点N101101采用完全耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点N101105采用完全耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点N101109采用完全耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点N101113采用完全耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点N101501采用完全耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点N101101采用顺序耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点N101105采用顺序耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点N101109采用顺序耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点N101113采用顺序耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点N101501采用顺序耦合方法得到的销-盘磨损试验的接触节点磨损量。

具体实施方式

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

如图1所示，本发明提供一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，包括以下步骤：

1)时域离散：假设摩擦热-应力-磨损耦合过程的时间轴为(t₀,t_N)，在该时间轴上将摩擦热-应力-磨损耦合过程离散为N个增量步，对瞬态热传导分析、应力分析和磨损分析采用相同的时间轴和增量步，增量步的设定需要满足：

式中，Δt_l为第l个增量步在时间轴对应的时间增量，l∈(0,N),Δt_l＝t_l-t_l-1；N为增量步总数；T′为耦合过程总时间。

需要说明的是：各增量步的大小可以相同，也可以不同，但是要同时兼顾瞬态热传导分析、应力分析和磨损分析的收敛性和精度，即各增量步大小的设置需要具体问题具体分析。对于耦合过程中场变量剧烈变化、强非线性、变量之间强烈耦合的时间段应适当增加增量步数目并且减少增量步大小，以便提高仿真结果的精度。

在上述实施例中，瞬态热传导分析中的数值积分方法采用向后欧拉方法，应力分析中的数值积分方法采用Newmark-β方法，磨损分析中的数值积分方法采用显式欧拉算法。

2)瞬态热传导分析：建立干滑动摩擦的有限元模型，有限元模型包括两具有相对干滑动运动的构件，约束不能使两构件产生相对干滑动运动的自由度，初始化两构件之间的作用载荷，初始化两构件的材料参数，材料参数包括不同温度对应的比热容、导热系数、线膨胀系数、杨氏模量、密度和泊松比。

在步骤1)的时间轴(t₀,t_N)上采用瞬态热传导分析程序进行瞬态热传导分析，其步骤为：

①初始化有限元模型的温度，选定瞬态热传导分析程序，设置有限元模型中两构件的单元类型、生热热流边界和散热边界；

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

③获得并输出温度场。

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

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

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

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

④计算销-盘表面接触节点处的质量流速

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

瞬态热传导分析程序步骤①中的生热热流边界由摩擦功率引起，其热流表达式为：

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

摩擦系数μ为两接触面接触节点的温度T、平均接触压力等因素的函数，其表达式为：

平均接触压力和两接触表面在接触点的相对滑动速率根据两构件的具体结构和相对运动求解。

摩擦热流通过热流用户子程序DFLUX施加在其中一构件的底面上，其中DFLUX为程序名，两接触面接触节点温度T、接触节点到构件边沿的距离均可以从计算结果文件中读取。以销-盘实验为例，用户子程序DFLUX包括以下步骤(如图3所示)：

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

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

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

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

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

3)应力分析：设置有限元模型应力分析的单元类型，定义边界条件、载荷和接触条件，将步骤2)输出的温度场作为热载荷施加在有限元模型中，并对有限元模型进行应力分析，获得并输出接触压力场。

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

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

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

公式(4)的左右两边同时除以增量步时间，则得到如下公式：

式中，为两接触面接触节点的磨损速率。

对于参数κ，为温度和接触压力等因素的函数，其表达式为：可由实验测量得到；p和Δs均可从有限元分析结果的文件中读取；可以由下式计算得到：

式中，Δt为增量步时间大小。

利用公式(4)～公式(6)可以求得一个增量步中所有两接触面接触节点处的磨损量或者磨损速率。

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

磨损量的确定通过网格移动用户子程序UMESHMOTION(用户子程序名称)实现(如图4所示)，包括以下步骤：

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

②读取两接触面接触节点的接触压力p和温度T；

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

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

5)更新网格：根据两接触面接触节点的磨损量和方向，修正接触压力与相对切向滑移均不为零的两接触面接触节点的位移，更新有限元模型；更新网格是通过更新两接触面接触节点位移，然后重画磨损区域网格实现的。磨损仿真要求更新两接触面接触节点位移而不引入附加节点力，否则两接触面接触节点位移的更新相当于只改变了接触材料的变形，而非磨损。为了实现磨损的仿真而不引入附加节点力，本发明采用任意拉格朗日欧拉法(ALE)。通过ALE方法，可以实现在不改变有限元分析结果中其他变量的前提下只改变两接触面接触节点的位移。ALE方法需要设置ALE作用域和ALE约束条件。ALE的作用域是所有包含需要修改位移的两接触面接触节点的单元，ALE的约束范围为所有需要修改位移的两接触面接触节点，ALE的约束类型为位移约束或者速度约束。以速度约束为例，若两接触面接触节点i的磨损速率为其在全局坐标系下的方向向量为节点i的局部坐标系三个轴的归一化方向矢量分别为则该节点的ALE速度约束矢量在局部坐标系中的分量{v_xi,v_yi,v_zi}^T，由下公式计算：

6)判断应力-磨损耦合过程的模拟是否完成，若增量步达到N，则模拟过程结束；否则，进入下一个增量步，回到步骤2)。

实施例：以销-盘模型为例，采用本发明方法对干滑动摩擦热-应力-磨损顺序耦合进行模拟，包括以下步骤：

1)时域离散：假设摩擦热-应力-磨损耦合过程的时间轴为(0,90s)，则耦合过程总时间T′＝90s。考虑到本实施例数值积分的收敛性和准确性，选取增量步数量为N＝40000，每个增量步时间相等，则Δt_l＝2.25×10^-3s,l＝1,2,…N。

2)瞬态热传导分析：如图5所示，建立销-盘的有限元模型，销上端固定约束，盘除了沿y方向平移和绕y轴转动自由度外的其它自由度均被约束。沿y轴正方向，销对盘施加50N载荷，盘绕y轴负方向的转速为600rpm，销和盘的材料参数如表1和表2所示。

表1

表2

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

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

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

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

④计算销-盘表面接触节点处的质量流速

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

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

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

摩擦系数μ为两接触面接触节点的温度T、平均接触压力等因素的函数，其表达式为：

平均接触压力通过下式求取：

式中，F为载荷，F＝50N，R为销半径，R＝3.9mm。

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

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

摩擦热流通过热流用户子程序DFLUX(如图3所示)施加在销的底面上，销-盘表面接触节点温度T、回转半径均可以从子程序中计算结果文件中读取。

用户子程序DFLUX包括以下步骤：

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

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

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

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

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

销-盘接触界面的热阻设置为0，盘外表面设置为对流散热边界，其对流换热系数为100W/m².K，环境温度为25℃。

3)应力分析，所有增量步完成瞬态热传导分析后，开始应力-磨损耦合分析，由于每个增量步内温升和磨损均很小，接触压力几乎不变，为了进一步增加仿真速度，每40个增量步内使用相同的接触压力，因此，每隔40个增量步重新计算接触状况。设置有限元模型应力分析的单元类型为C3D8R，按照实例工况定义边界条件、载荷、接触条件，在ABAQUS中将热传导分析的温度场作为热载荷施加在模型中，并对模型进行应力分析，获得并输出接触压力场。

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

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

本实例中，磨损系数κ是销-盘表面接触节点温度T和销-盘表面接触节点压力p的函数，其表达式为：

κ＝5.300+3.206×10^-1T-4.616p (13)

-1.065×10^-3T²-1.045×10^-2T·p+1.002p²

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

Δs＝2πnrΔt (14)

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

如图6和图7所示，销截面的节点分为边界节点和非边界节点，位于截面圆周上的节点为边界节点，而位于圆周以内的节点为非边界节点。对于非边界节点，其磨损量的方向为内法线方向，如果忽略变形则为y轴正方向。而对于边界节点，以节点N101114为例，其磨损方向为节点N101114与节点N102114的连线方向。输出销-盘表面接触节点的磨损量和方向，用于绘制磨损云图。

磨损量的确定通过网格移动用户子程序UMESHMOTION(用户子程序名称)实现(如图4所示)，包括以下步骤：

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

②读取销-盘表面接触节点的接触压力p和温度T；

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

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

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

6)判断应力-磨损耦合过程的模拟是否完成，若增量步达到N，则模拟过程结束；否则，进入下一个增量步，回到步骤2)。

为了说明本发明方法对于热-应力-磨损弱耦合问题的适用性，将上述实例的仿真结果与利用全耦合方法计算得到的结果对比如下：

全耦合方法考虑了热传导、应力以及磨损两两之间的双向耦合效应，本发明忽略了应力和磨损对热传导单方向的耦合作用，保留了其他所有耦合关系。

图8是采用本发明得到的销-盘磨损试验的温度分布与全耦合方法的对比，两者分布形式十分接近。

图9～图11分别是采用本发明得到的销-盘磨损试验的接触节点的接触面温度、接触压力、磨损量与全耦合方法的对比，图中，这两种方法得到的接触面5节点的曲线均十分接近，最大相对误差在5％左右。这说明了本发明对于热-应力-磨损弱耦合问题的适用性。

在相同的商用机HP Compaq Elite8300Convertible Microtower上，本发明计算用时121min，全耦合方法计算用时867min，两者相差7倍以上。因此，使用本发明可以大大提高计算效率。

上述各实施例仅用于说明本发明，各部件的连接和结构都是可以有所变化的，在本发明技术方案的基础上，凡根据本发明原理对个别部件的连接和结构进行的改进和等同变换，均不应排除在本发明的保护范围之外。

Claims (10)

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

1)时域离散：假设摩擦热-应力-磨损耦合过程的时间轴为(t₀,t_N)，在该时间轴上将摩擦热-应力-磨损耦合过程离散为N个增量步，瞬态热传导分析、应力分析和磨损分析采用相同的时间轴和增量步；

2)瞬态热传导分析：建立干滑动摩擦的有限元模型，在步骤1)的时间轴(t₀,t_N)上采用瞬态热传导分析程序进行瞬态热传导分析，其步骤为：

①初始化有限元模型的温度，选定瞬态热传导分析程序，设置有限元模型中两构件的单元类型、生热热流边界和散热边界；生热热流边界由摩擦功率引起，则两接触面接触节点的摩擦热流通过热流自定义子程序DFLUX施加在构件的底面上，两接触面接触节点温度T、接触节点到构件边沿的距离均从计算结果文件中读取：

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

ⅱ)初始化常量和系数；

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

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

ⅴ)返回瞬态热传导分析程序；

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

③获得并输出温度场；

如果其中一单元类型允许其对应的构件的材料脱离网格运动，那么材料的运动利用用户子程序UMASFL实现：

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

②初始化构件的常量和系数；

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

④计算销-盘表面接触节点处的质量流速 为角速度，ρ为材料密度；

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

3)应力分析：设置有限元模型应力分析的单元类型，定义边界条件、载荷和接触条件，将步骤2)输出的温度场作为热载荷施加在有限元模型中，并对有限元模型进行应力分析，获得并输出接触压力场；

4)确定磨损量：根据热传导分析得到的温度场、应力分析得到的接触压力场以及两接触面接触节点的相对滑移速率，计算两接触面接触节点的磨损量，并且确定其空间方向；磨损量大小和方向的确定通过网格移动用户子程序UMESHMOTION实现，包括以下步骤：

①初始化常量和系数；

②读取两接触面接触节点的接触压力p和温度T；

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

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

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

6)判断应力-磨损耦合过程的模拟是否完成，若增量步达到N，则模拟过程结束；否则，进入下一个增量步，回到步骤2)。

2.如权利要求1所述的一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，其特征在于：所述步骤1)中的瞬态热传导分析中的数值积分方法采用向后欧拉方法，应力分析中的数值积分方法采用Newmark-β方法，磨损分析中的数值积分方法采用显式欧拉算法。

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

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

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

摩擦系数μ为两接触面接触节点的温度T、平均接触压力因素的函数，其表达式为：

$\mathrm{\μ}=\mathrm{\μ}(T,\stackrel{\‾}{p});$

平均接触压力和两接触表面在接触点的相对滑动速率根据有限元模型的具体结构和相对运动求解。

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

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

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

摩擦系数μ为两接触面接触节点的温度T、平均接触压力因素的函数，其表达式为：

$\mathrm{\μ}=\mathrm{\μ}(T,\stackrel{\‾}{p});$

平均接触压力和两接触表面在接触点的相对滑动速率根据有限元模型的具体结构和相对运动求解。

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

Δh＝κ·p·Δs，

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

上式的左右两边同时除以增量步时间，则得到如下公式：

$\stackrel{\·}{h}=\mathrm{\κ}\·p\·\stackrel{\·}{\mathrm{\γ}},$

式中，为两接触面接触节点的磨损速率；

对于参数κ，为温度和接触压力因素的函数，其表达式为：由实验测量得到；p和Δs从有限元分析结果的文件中读取；从下式计算得到：

$\stackrel{\·}{\mathrm{\γ}}=\frac{\mathrm{\Δ}s}{\mathrm{\Δ}t};$

式中，Δt为增量步时间大小。

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

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

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

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

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