一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法
【专利摘要】本发明涉及一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，步骤为：在摩擦热-应力-磨损耦合过程的时间轴上将该过程离散为N个增量步，瞬态热传导分析、应力分析和磨损分析采用相同时间轴和增量步；建立干滑动摩擦的有限元模型，进行瞬态热传导分析，获得输出温度场；设置有限元模型应力分析的单元类型，进行应力分析，获得输出接触压力场；根据温度场、接触压力场及两接触面接触节点的相对滑移速率，计算两接触面接触节点的磨损量并确定其空间方向；根据磨损量和方向，修正接触节点的位移，更新有限元模型；判断应力-磨损耦合过程的模拟是否完成。本发明计算精度损失小且计算效率较高，能广泛应用于各类干滑动摩擦的摩擦热-应力-磨损弱耦合问题的模拟。
【专利说明】一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法
【技术领域】
[0001]本发明涉及一种干滑动摩擦耦合模拟方法，特别是关于一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法。
【背景技术】
[0002]干滑动摩擦界面存在摩擦热-应力-磨损耦合现象，其中，摩擦热主要产生于两个接触面之间的摩擦接触表层，摩擦热的热流密度与切向摩擦应力和相对滑移速度有关。摩擦热会直接导致接触面的温度迅速升高，迅速升高的温度致使接触面发生热变形，从而最终使得接触压力在接触面上重新不均匀的分布，在接触压力和相对滑移速度的综合作用下，接触面逐渐磨损。由于接触面上接触压力分布不均以及接触面不同位置处相对滑移速度存在差异，因此，接触面不同位置处的磨损量不同。这种接触面上磨损量的不均匀性导致接触压力的再次重新分布，进而改变切向摩擦应力的分布和热流密度的大小。另一方面，接触面温度的改变还会显著改变两接触材料的摩擦特性，比如摩擦系数和磨损系数,除此之夕卜，还会改变材料的本构特性，比如屈服极限和硬度等。因此，接触面的摩擦热(温度场)、应力(接触压力场)和磨损(磨损量)之间是相互耦合的，这种耦合行为几乎不可能求得解析解，而采用数值方法进行求解数值解时，如果考虑摩擦热-应力-磨损三者之间的双向耦合效应，那么计算的效率极低，无法求解较长时间的摩擦热-应力-磨损耦合问题。
【发明内容】
[0003]针对上述问题，本发明的目的是提供一种计算精度损失小且计算效率高的干滑动摩擦热-应力-磨损顺序耦合的模拟方法。
[0004]为实现上述目的，本发明采取以下技术方案:一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，包括以下步骤:1)时域离散:假设摩擦热-应力-磨损耦合过程的时间轴为(V tN)，在该时间轴上将摩擦热-应力-磨损耦合过程离散为N个增量步，瞬态热传导分析、应力分析和磨损分析采用相同的时间轴和增量步；2)瞬态热传导分析:建立干滑动摩擦的有限元模型，在步骤I)的时间轴Uci^n)上采用瞬态热传导分析程序进行瞬态热传导分析，其步骤为:①初始化有限元模型的温度，选定瞬态热传导分析程序，设置有限元模型中两构件的单元类型、生热热流边界和散热边界对每个增量步进行瞬态热传导分析；
③获得并输出温度场；3)应力分析:设置有限元模型应力分析的单元类型，定义边界条件、载荷和接触条件，将步骤2)输出的温度场作为热载荷施加在有限元模型中，并对有限元模型进行应力分析，获得并输出接触压力场；4)确定磨损量:根据热传导分析得到的温度场、应力分析得到的接触压力场以及两接触面接触节点的相对滑移速率，计算两接触面接触节点的磨损量，并且确定其空间方向；5)更新网格:根据两接触面接触节点的磨损量和方向，修正该接触节点的位移，更新有限元模型；6)判断应力-磨损耦合过程的模拟是否完成，若增量步达到N，则模拟过程结束；否则，进入下一个增量步，回到步骤2)。
[0005]所述步骤I)中的瞬态热传导分析中的数值积分方法采用向后欧拉方法，应力分析中的数值积分方法采用Newmark- β方法，磨损分析中的数值积分方法采用显式欧拉算法。
[0006]所述步骤2)中，瞬态热传导分析程序中的生热热流边界由摩擦功率引起，其热流密度表达式为.Λ = μΤ,式中，f为两接触面接触节点的摩擦热流；μ为两接触面接触节点的摩擦系数为接触面的平均接触压力；f为两接触表面在接触点的相对滑动速率；摩擦系数μ为两接触面接触节点的温度Τ、平均接触压力P等因素的函数，其表达式为://二/I,...):平均接触压力$和两接触表面在接触点的相对滑动速率^根据有限元模型的具体结构和相对运动求解；摩擦热流f通过热流自定义子程序DFLUX施加在构件的底面上，两接触面接触节点温度T、接触节点到构件边沿的距离均从计算结果文件中读取；用户子程序DFLUX包括以下步骤:①从瞬态热传导分析程序中跳转至用户子程序DFLUX ;②初始化常量和系数；③根据两接触面接触节点坐标计算回转半径？，根据两接触面接触节点温度T和平均接触压力]H十算摩擦系数μ ;④根据回转半径？和摩擦系数μ计算两接触面接触节点处的摩擦热流^ ；⑤返回瞬态热传导分析程序。
[0007]所述步骤4)中的磨损量大小的计算采用广义非线性Archard磨损公式:Ah =K.ρ.As,式中，Ah为两接触面接触节点的磨损增量；K为广义Archard磨损公式的磨损系数；Ρ为两接触面接触节点的接触压力；△ s为两接触表面在接触点的相对滑动位移增量；上式的左右两边同时除以增量步时间，则得到如下公式:A = ，式中，A为两接触面接触节点的磨损速率；对于参数K，为温度和接触压力等因素的函数，其表达式为:
=，由实验测量得到；1)和As从有限元分析结果的文件中读取彳从下式计算得到:
'—tsS
[0_] r=-_
[0009]式中，At为增量步时间大小。
[0010]所述步骤4)中的磨损量的方向的确定方法如下:两接触面接触节点分为两类，一类是位于接触面边缘的节点，简称边界节点；另一类是位于接触面内部的节点，简称非边界节点；对于非边界节点的磨损量方向是由节点法方向确定，节点法方向是直接从有限元分析结果文件中读取。
[0011]所述步骤4)中的磨损量大小和方向的确定通过网格移动用户子程序UMESHMOT1N实现，包括以下步骤:①初始化常量和系数；②读取两接触面接触节点的接触压力P和温度T ;③计算两接触面接触节点的磨损量，确定两接触面接触节点的磨损方向；
④输出两接触面接触节点磨损信息，用于绘制磨损量云图。
[0012]所述步骤5)中，两接触面接触节点位移的修正和有限元模型的更新采用任意拉格朗日欧拉法实现；需要修正位移的两接触面接触节点为接触压力与相对切向滑移均不为零的接触节点。
[0013]本发明由于采取以上技术方案，其具有以下优点:1、本发明由于在实际工程应用，如制动过程中，摩擦热-应力属于弱耦合问题，可以忽略应力和磨损对摩擦热传导单方向的耦合项，本发明采用的顺序耦合方法模拟摩擦热-应力-磨损耦合行为，即将摩擦热-应力-磨损弱耦合问题分解为瞬态热传导分析和应力-磨损耦合分析，因此，在计算精度损失较小的情况下极大地提高计算效率。2、本发明由于基于现有商用的有限元软件平台，因此设置简单、容易上手且便于推广。3、本发明由于采用的瞬态热传导分析可以直接使用商用软件的算法，应力-磨损耦合分析在商用有限元软件的应力分析算法的基础上二次开发完成，磨损的计算完全嵌入有限元软件之中，免去了数据传输和交互，因此，使得计算效率大大提高。本发明适用范围广，可以广泛应用于各类干滑动摩擦的摩擦热-应力-磨损弱耦合问题的模拟。
【专利附图】
【附图说明】
[0014]图1是本发明的流程示意图；
[0015]图2是本发明的用户子程序UMASFL流程示意图；
[0016]图3是本发明的用户子程序DFLUX流程示意图；
[0017]图4是本发明的用户子程序UMESHMOT1N流程示意图；
[0018]图5是本发明的销-盘的有限元模型示意图；
[0019]图6是本发明的销底部的局部放大示意图；
[0020]图7是图2中I处局部放大示意图；
[0021]图8是采用顺序耦合方法和完全耦合方法得到的销-盘磨损试验的温度分布对比示意图，其中图(a)为采用完全耦合方法得到的销-盘磨损试验的温度分布，图(b)为采用顺序耦合方法得到的销-盘磨损试验的温度分布；
[0022]图9是本发明采用顺序耦合方法和完全耦合方法得到的销-盘磨损试验的接触节点温度分布对比示意图，其中，表示节点N101101采用完全耦合方法得到的销-盘磨损试验的接触节点温度，O表示节点N101105采用完全耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101109采用完全耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101113采用完全耦合方法得到的销-盘磨损试验的接触节点温度，
表示节点N101501采用完全耦合方法得到的销-盘磨损试验的接触节点温度，* 表示节点N101101采用顺序稱合方法得到的销-盘磨损试验的接触节点温度， 表示节点NlOl 105采用顺序耦合方法得到的销-盘磨损试验的接触节点温度，表示节点NlOl 109采用顺序耦合方法得到的销-盘磨损试验的接触节点温度， 表示节点N101113采用顺序耦合方法得到的销-盘磨损试验的接触节点温度，表示节点N101501采用顺序耦合方法得到的销-盘磨损试验的接触节点温度；
[0023]图10是本发明采用顺序耦合方法和完全耦合方法得到的销-盘磨损试验的接触节点接触压力对比示意图，其中， 表示节点N101101采用完全耦合方法得到的销-盘磨损试验的接触节点接触压力，-表示节点NlOl 105采用完全耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点NlOl 109采用完全耦合方法得到的销-盘磨损试验的接触节点接触压力， ' 表示节点N101113采用完全耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点N101501采用完全耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点Ν101101采用顺序耦合方法得到的销-盘磨损试验的接触节点接触压力，O*表示节点Ν101105采用顺序耦合方法得到的销-盘磨损试验的接触节点接触压力，表示节点Ν101109采用顺序稱合方法得到的销-盘磨损试验的接触节点接触压力，表示节点Ν101113采用顺序耦合方法得到的销-盘磨损试验的接触节点接触压力，' ^表示节点Ν101501采用顺序稱合方法得到的销-盘磨损试验的接触节点接触压力；
[0024]图11是本发明采用顺序耦合方法和完全耦合方法得到的销-盘磨损试验的接触节点磨损量对比示意图，其中，?~表示节点Ν101101采用完全稱合方法得到的销-盘磨损试验的接触节点磨损量， 表示节点Ν101105采用完全耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点Ν101109采用完全耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点Ν101113采用完全耦合方法得到的销-盘磨损试验的接触节点磨损量， 表示节点Ν101501采用完全耦合方法得到的销-盘磨损试验的接触节点磨损量，.0表示节点Ν101101采用顺序耦合方法得到的销-盘磨损试验的接触节点磨损量，
表示节点Ν101105采用顺序耦合方法得到的销-盘磨损试验的接触节点磨损量，
表示节点Ν101109采用顺序耦合方法得到的销-盘磨损试验的接触节点磨损量，表示节点Ν101113采用顺序稱合方法得到的销-盘磨损试验的接触节点磨损量，表示节点Ν101501采用顺序耦合方法得到的销-盘磨损试验的接触节点磨损量。
【具体实施方式】
[0025]下面结合附图和实施例对本发明进行详细的描述。
[0026]如图1所示，本发明提供一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，包括以下步骤:
[0027]I)时域离散:假设摩擦热-应力-磨损耦合过程的时间轴为(V tN)，在该时间轴上将摩擦热-应力-磨损耦合过程离散为N个增量步，对瞬态热传导分析、应力分析和磨损分析采用相同的时间轴和增量步，增量步的设定需要满足:
【权利要求】
1.一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，包括以下步骤: 1)时域离散:假设摩擦热-应力-磨损耦合过程的时间轴为tN)，在该时间轴上将摩擦热-应力-磨损耦合过程离散为N个增量步，瞬态热传导分析、应力分析和磨损分析采用相同的时间轴和增量步； 2)瞬态热传导分析:建立干滑动摩擦的有限元模型，在步骤I)的时间轴(^t,)上采用瞬态热传导分析程序进行瞬态热传导分析，其步骤为: ①初始化有限元模型的温度，选定瞬态热传导分析程序，设置有限元模型中两构件的单元类型、生热热流边界和散热边界； ②对每个增量步进行瞬态热传导分析； ③获得并输出温度场； 3)应力分析:设置有限元模型应力分析的单元类型，定义边界条件、载荷和接触条件，将步骤2)输出的温度场作为热载荷施加在有限元模型中，并对有限元模型进行应力分析，获得并输出接触压力场； 4)确定磨损量:根据热传导分析得到的温度场、应力分析得到的接触压力场以及两接触面接触节点的相对滑移速率，计算两接触面接触节点的磨损量，并且确定其空间方向； 5)更新网格:根据两接触面接触节点的磨损量和方向，修正该接触节点的位移，更新有限元模型； 6)判断应力-磨损耦合过程的模拟是否完成，若增量步达到N，则模拟过程结束；否则，进入下一个增量步，回到步骤2)。
2.如权利要求1所述的一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，其特征在于:所述步骤I)中的瞬态热传导分析中的数值积分方法采用向后欧拉方法，应力分析中的数值积分方法采用Newmark-β方法，磨损分析中的数值积分方法采用显式欧拉算法。
3.如权利要求1所述的一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，其特征在于:所述步骤2)中，瞬态热传导分析程序中的生热热流边界由摩擦功率引起，其热流密度表达式为:
Q = MPY， 式中，q为两接触面接触节点的摩擦热流；μ为两接触面接触节点的摩擦系数；尹为接触面的平均接触压力^为两接触表面在接触点的相对滑动速率； 摩擦系数μ为两接触面接触节点的温度Τ、平均接触压力P等因素的函数，其表达式为:
μ = μ(Τ,ρ,...)., 平均接触压力P和两接触表面在接触点的相对滑动速率^根据有限元模型的具体结构和相对运动求解； 摩擦热流通过热流自定义子程序DFLUX施加在构件的底面上，两接触面接触节点温度Τ、接触节点到构件边沿的距离均从计算结果文件中读取；用户子程序DFLUX包括以下步骤: ①从瞬态热传导分析程序中跳转至用户子程序DFLUX ； ②初始化常量和系数； ③根据两接触面接触节点坐标计算回转半径？，根据两接触面接触节点温度T和平均接触压力声计算摩擦系数μ ； ④根据回转半径？和摩擦系数μ计算两接触面接触节点处的摩擦热^ ⑤返回瞬态热传导分析程序。
4.如权利要求2所述的一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，其特征在于:所述步骤2)中，瞬态热传导分析程序中的生热热流边界由摩擦功率引起，其热流密度表达式为:
q = μΤπ， 式中，H为两接触面接触节点的摩擦热流；μ为两接触面接触节点的摩擦系数；戶为接触面的平均接触压力；^为两接触表面在接触点的相对滑动速率； 摩擦系数μ为两接触面接触节点的温度Τ、平均接触压力#等因素的函数，其表达式为:
μ = μ{Τ,ρ,...)., 平均接触压力P和两接触表面在接触点的相对滑动速率#根据有限元模型的具体结构和相对运动求解； 摩擦热流^通过热流自定义子程序DFLUX施加在构件的底面上，两接触面接触节点温度Τ、接触节点到构件边沿的距离均从计算结果文件中读取；用户子程序DFLUX包括以下步骤: ①从瞬态热传导分析程序中跳转至用户子程序DFLUX； ②初始化常量和系数； ③根据两接触面接触节点坐标计算回转半径？，根据两接触面接触节点温度T和平均接触压力多计算摩擦系数μ ； ④根据回转半径？和摩擦系数μ计算两接触面接触节点处的摩擦热流歹； ⑤返回瞬态热传导分析程序。
5.如权利要求1或2或3或4所述的一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，其特征在于:所述步骤4)中的磨损量大小的计算采用广义非线性Archard磨损公式: Ah= K.P.As， 式中，Ah为两接触面接触节点的磨损增量；K为广义Archard磨损公式的磨损系数；P为两接触面接触节点的接触压力；△ s为两接触表面在接触点的相对滑动位移增量； 上式的左右两边同时除以增量步时间，则得到如下公式:
h-K' p-f j 式中，《为两接触面接触节点的磨损速率； 对于参数K，为温度和接触压力等因素的函数，其表达式为= 六…)，由实验测量得到；p和△ S从有限元分析结果的文件中读取 >从下式计算得到:.Av Y=—
At.1 式中，At为增量步时间大小。
6.如权利要求1或2或3或4所述的一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，其特征在于:所述步骤4)中的磨损量的方向的确定方法如下:两接触面接触节点分为两类，一类是位于接触面边缘的节点，简称边界节点；另一类是位于接触面内部的节点，简称非边界节点；对于非边界节点的磨损量方向是由节点法方向确定，节点法方向是直接从有限元分析结果文件中读取。
7.如权利要求5所述的一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，其特征在于:所述步骤4)中的磨损量的方向的确定方法如下:两接触面接触节点分为两类，一类是位于接触面边缘的节点，简称边界节点；另一类是位于接触面内部的节点，简称非边界节点；对于非边界节点的磨损量方向是由节点法方向确定，节点法方向是直接从有限元分析结果文件中读取。
8.如权利要求1到7任一项所述的一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，其特征在于:所述步骤4)中的磨损量大小和方向的确定通过网格移动用户子程序UMESHMOT1N实现，包括以下步骤: ①初始化常量和系数； ②读取两接触面接触节点的接触压力P和温度T； ③计算两接触面接触节点的磨损量，确定两接触面接触节点的磨损方向； ④输出两接触面接触节点磨损信息，用于绘制磨损量云图。
9.如权利要求1到8任一项所述的一种干滑动摩擦热-应力-磨损顺序耦合的模拟方法，其特征在于:所述步骤5)中，两接触面接触节点位移的修正和有限元模型的更新采用任意拉格朗日欧拉法实现；需要修正位移的两接触面接触节点为接触压力与相对切向滑移均不为零的接触节点。
【文档编号】G06F17/50GK104200034SQ201410457959
【公开日】2014年12月10日 申请日期:2014年9月10日 优先权日:2014年9月10日
【发明者】桂良进, 张方宇, 范子杰 申请人:清华大学