CALCULATION METHOD FOR VARIABLE FRICTION CONTACT WEAR BASED ON ADAPTIVE PENALTY FACTOR
-
摘要: 磨损广泛出现在各种物体的相互接触中, 准确且高效地模拟物体的接触磨损行为对研究磨损危害及其预防至关重要. 文章基于罚函数接触算法和Archard磨损模型, 结合考虑常数、线性和指数压力相关摩擦系数的库伦摩擦模型, 构建了一个包含几何非线性、材料非线性和边界非线性等复杂非线性因素的接触磨损数值模拟方法. 为了提高变摩擦系数模型的效率和精度, 提出了变摩擦自适应罚因子(VFAPF)算法, 并将其应用于线弹性材料和Mooney-Rivlin超弹性材料的接触磨损行为仿真分析中. 研究结果表明: 相比于自适应罚因子(APF)算法, VFAPF算法在变摩擦大滑移问题中具备更好的精度稳定性, 当弹性滑移量占主导时, 虽效率有小幅下降(约12%), 但精度的提升效果更明显(约62%). 针对线弹性材料的小滑移接触磨损问题, 尽管3种摩擦系数模型在磨损量和应力分布上存在差异, 但对磨损位置和接触压力的影响较小. 相比之下, 在Mooney-Rivlin超弹性材料的大滑移接触磨损问题中, 3种摩擦系数模型在磨损量、磨损位置、接触压力和应力分布方面均表现出显著差异, 其中线性摩擦系数模型的差异性最突出.
-
关键词:
- 接触磨损 /
- 压力相关变摩擦系数 /
- VFAPF算法 /
- Mooney-Rivlin超弹性材料 /
- 大滑移问题
Abstract: Wear is commonly observed in the interaction between various objects, and accurately and efficiently simulating the contact wear behavior of these objects is crucial for studying the hazards of wear and its prevention. Based on a penalty function contact algorithm and the Archard wear model, combined with the Coulomb friction model that considers constant, linear and exponential pressure-dependent friction coefficients, incorporating complex nonlinear factors such as geometric nonlinearity, material nonlinearity, and boundary nonlinearity into the numerical simulation method for contact wear. To improve the efficiency and accuracy of the variable friction coefficient model, this paper proposes a variable friction adaptive penalty factor (VFAPF) algorithm and applies it to the simulation and analysis of contact wear behavior in linear elastic materials and Mooney-Rivlin hyperelastic materials. The results show that compared with the adaptive penalty factor (APF) algorithm, the VFAPF algorithm possesses better accuracy stability in the variable friction large slip problem, and when the elastic slip is dominant, although there is a small decrease in the efficiency (about 12%), the improvement of the accuracy is more obvious (about 62%). For the small-sliding contact wear problem of linear elastic material, despite differences in wear gap and stress distribution among the three friction coefficient models, their effects on the wear location and contact pressure are minimal. In contrast, in the large-sliding contact wear problem for Mooney-Rivlin hyperelastic material, the three friction coefficient models show significant differences in wear gap, wear location, contact pressure, and stress distribution, with the linear friction coefficient model exhibiting the most pronounced variation. -
引 言
在现实生活中, 物体磨损是一种常见现象, 普遍出现在轮胎、刹车盘和齿轮等多种零部件中. 磨损不仅影响零部件的功能和性能, 如振动、噪音和可控性, 还关系到使用寿命和安全[1]. 因此, 利用磨损机制对物体磨损进行模拟预测, 可以显著优化零部件的设计, 提升性能并延长使用寿命, 从而提高整体运行的效率和安全性[2].
一般来说, 磨损可以分为磨料磨损、黏着磨损、腐蚀磨损和表面疲劳磨损[3-4]. 磨损的实验建模方法有很多学者进行了研究, Vieira等[5]利用实验测试对多种橡胶材料的摩擦和磨损进行评估, 揭示了橡胶的黏弹性特性对摩擦和磨损率的影响. 邹龙庆等[6]研究了金属-橡胶的粗糙面磨损接触, 验证了磨料磨损是刚柔接触界面的主要磨损形式. Kahms等[7]对轮胎花纹块样本进行磨损测试, 结合摩擦系数、接触压力和滑动速度, 成功模拟了不同条件下的飞机轮胎磨损行为. Sakhnevych等[8]和Huang等[9]考虑轮胎材料的黏弹性特性、路面粗糙度和环境, 从动力学和热力学的角度分析了轮胎的磨损特性. 付豪等[10]研究了SiC/AZ91D镁基复合材料在活塞制造中的磨损特性, 主要磨损机制为磨料和剥层磨损, 同时磨损深度随载荷增加而增加.
从数值模拟的角度来看, 常见的磨损模型有Archard磨损模型[11]和能量磨损模型[12]. 磨损问题包含着接触非线性, 而接触是计算力学中最难的问题之一. 接触计算的难点主要体现在两个方面, 一方面是接触区域的不确定性, 且接触可能涉及到一对一或者一对多的情况. 另一方面是接触边界上具有高度非线性, 当接触未发生时接触力为0, 当发生接触的一瞬间接触力开始突变, 并且接触过程还需要考虑摩擦力的影响[13]. 接触区域的不确定性和接触界面的高度非线性给求解接触问题带来了许多麻烦. 随着技术的不断发展, 目前已经发展了许多用于求解接触问题的数值方法, 如罚函数法[14]、拉格朗日乘子法[15]、增广拉格朗日乘子法[16]、扰动拉格朗日乘子法[17]和双势方法[18-19]等.
接触磨损过程还需要考虑摩擦力的影响, 常见的摩擦模型有Coulomb摩擦[20]、广义黏滞摩擦[21]和Amontons摩擦[22]等. 在Coulomb摩擦模型中, 切向接触力不能超过极限摩擦力, 极限摩擦力可以用摩擦系数与法向接触力的乘积表示[20]. 对于同一种材料, 摩擦系数通常近似为一个恒定的常数, 但实际上摩擦系数会与物体的运动速度、加速度、温度和接触压力等内变量相关联[23]. Luo[24]在研究TiAlN/VN涂层的氧化铝无润滑滑动磨损实验中发现, 从室温到700 °C内涂层的摩擦系数和磨损特性显著受温度的影响. Mäntylä等[25]考虑了切向牵引力相关的摩擦系数, 将Archard磨损模型应用于柱面磨损实验. 将实验结果与模拟结果比较后发现, Archard 磨损模型具有一定局限性, 没有考虑到接触中磨损碎片的夹带从而高估了磨损量. Yue等[26]研究了摩擦系数随磨损循环次数变化的情况, 发现在部分滑动和黏着阶段时, 采用变摩擦系数的有限元模型能够更准确地预测实验结果. Ning等[27]基于双势接触算法, 采用压力相关摩擦系数模型对线弹性材料的接触磨损问题进行了数值求解.
现有研究缺乏对复杂接触磨损问题的完整数值模拟方法, 并且未考虑压力相关摩擦系数对超弹性材料大滑移接触磨损特性的影响. 因此, 本文基于罚函数接触算法, 采用不同的压力相关摩擦系数模型计算接触磨损问题. 同时, 考虑变摩擦系数的自适应罚因子(VFAPF)算法. 基于此算法, 分别讨论变摩擦系数模型对小滑移和大滑移接触磨损过程的影响, 探究接触磨损问题中压力相关摩擦系数对物体的接触压力、磨损量、磨损位置和应力分布的影响, 以期为零部件的复杂接触磨损模拟提供参考方法.
1. 基于完全拉格朗日格式的非线性有限元
非线性有限元是线性有限元的拓展, 在分析过程中包含几何非线性、材料非线性和边界非线性等非线性因素. 非线性有限元包含两种格式, 第一种是以当前构型为参考构型的更新拉格朗日格式; 第二种是以初始构型为参考构型的完全拉格朗日格式, 这两种格式都是采用增量分析方法.
基于初始构型, 在完全拉格朗日格式下的虚功方程可以表示为[28]
$$ \int_{{{\varOmega }^0}} {{{({\delta }{{\boldsymbol{E}}})}^{\text{T}}}{{\boldsymbol{S}}}{\text{d}}{\varOmega }} - \int_{{{\varOmega }^0}} {{{({\delta }{{\boldsymbol{u}}})}^{\text{T}}}{{\boldsymbol{f}}}{\text{d}}{\varOmega }} - \int_{{\varGamma }_{^{\sigma }}^0} {{{({\delta }{{\boldsymbol{u}}})}^{\text{T}}}{{\boldsymbol{T}}}{\text{d}}{\varGamma }} = 0 $$ (1) 式中, E表示拉格朗日应变张量, S表示第二类Piola-Kirchhoff应力张量(PK2应力张量), f表示体力, T表示面力, ${{\varOmega }^0}$和${\varGamma }_\sigma ^0$分别表示初始构形对应的体积和力边界.
对方程(1)进行增量线性化得到
$$ \begin{split} & \int_{{{\varOmega }^0}} {{{\left( {{\delta }\Delta {{\boldsymbol{\varepsilon}} }} \right)}^{\text{T}}}{{\boldsymbol{D}}}\Delta {{\boldsymbol{\varepsilon}} }{\text{d}}{\varOmega }} + \int_{{{\varOmega }^0}} {{{\left( {{\delta }\Delta {{\boldsymbol{\varphi}} }} \right)}^{\text{T}}}{{\boldsymbol{S}}}' {\text{d}}{\varOmega }} = \\ &\qquad \int_{{{\varOmega }^0}} {{{\left( {{\delta }\Delta {{\boldsymbol{u}}}} \right)}^{\text{T}}}{{\boldsymbol{f}}}{\text{d}}{\varOmega }} + \int_{{\varGamma }_{\sigma }^0} {{{\left( {{\delta }\Delta {{\boldsymbol{u}}}} \right)}^{\text{T}}}{{\boldsymbol{T}}}{\text{d}}{\varGamma }} - \int_{{{\varOmega }^0}} {{{\left( {{\delta }\Delta {{\boldsymbol{\varepsilon}} }} \right)}^{\text{T}}}{{\boldsymbol{S}}}' {\text{d}}{\varOmega }} \end{split}$$ (2) 式中, $ {{\boldsymbol{D}}} $表示材料刚度矩阵, $ {{\boldsymbol{S}}}' $表示上一时刻的PK2应力张量, $\Delta {{\boldsymbol{\varepsilon}} }$和$\Delta {{\boldsymbol{\varphi}} }$分别表示增量拉格朗日应变$ \Delta {E} $的线性项和非线性项, 具体形式为
$$ \left. \begin{aligned} & \Delta {{\boldsymbol{\varepsilon}} } = \frac{1}{2}\left( {\frac{{\partial \Delta {{\boldsymbol{u}}}}}{{\partial {{\boldsymbol{X}}}}} + \frac{{\partial \Delta {{{\boldsymbol{u}}}^{\text{T}}}}}{{\partial {{\boldsymbol{X}}}}} + \frac{{\partial \Delta {{{{\boldsymbol{u}}'}}^{\text{T}}}}}{{\partial {{\boldsymbol{X}}}}}\frac{{\partial {{\boldsymbol{u}}}}}{{\partial {{\boldsymbol{X}}}}} + \frac{{\partial {{{\boldsymbol{u}}}^{\text{T}}}}}{{\partial {{\boldsymbol{X}}}}}\frac{{\partial \Delta {{\boldsymbol{u}}}' }}{{\partial {{\boldsymbol{X}}}}}} \right) \\ & \Delta {{\boldsymbol{\varphi}} } = \frac{1}{2}\frac{{\partial \Delta {{{\boldsymbol{u}}}^{\text{T}}}}}{{\partial {{\boldsymbol{X}}}}}\frac{{\partial \Delta {{\boldsymbol{u}}}}}{{\partial {{\boldsymbol{X}}}}} \end{aligned}\right\} $$ (3) 式中, $ {{\boldsymbol{u}}} $和$ \Delta {{\boldsymbol{u}}} $分别表示当前时刻的位移和位移增量, $ {{\boldsymbol{u}}}' $和$ \Delta {{\boldsymbol{u}}}' $分别表示上一时刻的位移和位移增量, $ {{\boldsymbol{X}}} $表示初始构型的质点坐标.
采用八节点六面体单元对方程(2)进行离散化, 则方程(2)的左端项可以离散为[29]
$$\begin{split} & \int_{{{\varOmega }^0}} {{{\left( {{\delta }\Delta {{\boldsymbol{\varepsilon}} }} \right)}^{\text{T}}}{{\boldsymbol{D}}}\Delta {{\boldsymbol{\varepsilon}} }{\text{d}}{\varOmega }} + \int_{{{\varOmega }^0}} {{{\left( {{\delta }\Delta {{\boldsymbol{\varphi}} }} \right)}^{\text{T}}}{{\boldsymbol{S}}}' {\text{d}}{\varOmega }} = \\ &\quad \sum\limits_e {{{\left( {{\delta }\Delta {{{\boldsymbol{a}}}^e}} \right)}^{\text{T}}}\left( {\int_{{\varOmega }_e^0} {\boldsymbol{B}_L^{\text{T}}{{\boldsymbol{D}}}{\boldsymbol{B}_L}{\text{d}}{\varOmega }} {\text{ + }}\int_{{\varOmega }_e^0} {{{\boldsymbol{B}}}_{NL}^{\text{T}}{\hat {\boldsymbol{S}}}{{{\boldsymbol{B}}}_{NL}}{\text{d}}{\varOmega }} } \right)\Delta {{{\boldsymbol{a}}}^e}} = \\ &\quad \sum\limits_e {{{\left( {{\delta }\Delta {{{\boldsymbol{a}}}^e}} \right)}^{\text{T}}}\left( {{{\boldsymbol{K}}}_L^e + {{\boldsymbol{K}}}_{NL}^e} \right)\Delta {{{\boldsymbol{a}}}^e}}\\[-1pt]\end{split} $$ (4) 式中, ${{{\boldsymbol{B}}}_L}$和${{{\boldsymbol{B}}}_{NL}}$分别表示线性应变矩阵和非线性应变矩阵, ${\hat {\boldsymbol{S}}}$表示PK2应力矩阵, 具体形式见Kim的著作[29]. 方程(2)的右端项可以离散为
$$ \begin{split} & \int_{{{\varOmega }^0}} {{{\left( {{\delta }\Delta {{\boldsymbol{u}}}} \right)}^{\text{T}}}{{\boldsymbol{f}}}{\text{d}}{\varOmega }} + \int_{{\varGamma }_{\sigma }^0} {{{\left( {{\delta }\Delta {{\boldsymbol{u}}}} \right)}^{\text{T}}}{{\boldsymbol{T}}}{\text{d}}{\varGamma }} - \int_{{{\varOmega }^0}} {{{\left( {{\delta }\Delta {{\boldsymbol{\varepsilon}} }} \right)}^{\text{T}}}{{\boldsymbol{S}}}' {\text{d}}{\varOmega }} = \\ &\;\; \sum\limits_e {{{\left( {{\delta }\Delta {{{\boldsymbol{a}}}^e}} \right)}^{\text{T}}}\left( {\int_{{\varOmega }_e^0} {{{{\boldsymbol{N}}}^{\text{T}}}{{\boldsymbol{f}}}{\text{d}}V} + \int_{{\varGamma }_{{\sigma ,e}}^0} {{{{\boldsymbol{N}}}^{\text{T}}}{{\boldsymbol{T}}}{\text{d}}S} - \int_{{\varOmega }_e^0} {{{\boldsymbol{B}}}_L^{\text{T}}\bar {{\boldsymbol{S}}} {\text{d}}V} } \right)} = \\ &\;\;\sum\limits_e {{{\left( {{\delta }\Delta {{{\boldsymbol{a}}}^e}} \right)}^{\text{T}}}\left( {{{{\boldsymbol{P}}}^e} - {{\boldsymbol{F}}}_{\text{int} }^e} \right)}\\[-1pt] \end{split} $$ (5) 式中, ${{{\boldsymbol{P}}}^e}$和${{\boldsymbol{F}}}_{\text{int} }^e$分别表示单元外力和内力向量, $\bar {{\boldsymbol{S}}} $表示PK2应力向量. 将方程(4)和(5)代入方程(2)得
$$ \sum\limits_e {{{\left( {{\delta }\Delta {{{\boldsymbol{a}}}^e}} \right)}^{\text{T}}}\left( {{{\boldsymbol{K}}}_L^e + {{\boldsymbol{K}}}_{NL}^e} \right)\Delta {{{\boldsymbol{a}}}^e}} = \sum\limits_e {{{\left( {{\delta }\Delta {{{\boldsymbol{a}}}^e}} \right)}^{\text{T}}}\left( {{{{\boldsymbol{P}}}^e} - {{\boldsymbol{F}}}_{\text{int} }^e} \right)} $$ (6) 对每个单元的刚度矩阵和力向量组装得到系统增量平衡方程, 即
$$ \left( {{{\boldsymbol{K}}}_L^{} + {{\boldsymbol{K}}}_{NL}^{}} \right)\Delta {{{\boldsymbol{a}}}^{}} = {{{\boldsymbol{P}}}^{}} - {{\boldsymbol{F}}}_{\text{int} }^{} $$ (7) 式中, ${{\boldsymbol{K}}}_L^{},{{\boldsymbol{K}}}_{NL}^{},\Delta {{{\boldsymbol{a}}}^{}},{{{\boldsymbol{P}}}^{}}$和${{\boldsymbol{F}}}_{\text{int} }^{}$分别表示系统线性刚度矩阵、非线性刚度矩阵, 节点位移增量、内力向量和外载荷向量.
2. 罚函数接触算法
2.1 接触界面本构关系
如图1所示, 物体${{\varOmega }_1}$和${{\varOmega }_2}$相互接触, ${\varGamma }_{C}^1$和${\varGamma }_{C}^2$分别表示两个物体的接触边界, 物体${{\varOmega }_1}$和边界${\varGamma }_{C}^1$分别称为从接触体和从接触面, 物体${{\varOmega }_2}$和边界${\varGamma }_{C}^2$分别称为主接触体和主接触面. 在边界${\varGamma }_{C}^2$上建立局部坐标系${\xi ^\alpha }\left( {\alpha = 1,2} \right)$, $ {\bar {\boldsymbol{n}}} $和${{\bar {\boldsymbol{a}}}_\alpha }$表示局部坐标系的法向和切向基矢量. 边界${\varGamma }_{C}^1$和${\varGamma }_{C}^2$上分别存在点${{{\boldsymbol{x}}}_s}$和${\bar {\boldsymbol{x}}}$, 点${\bar {\boldsymbol{x}}}$是边界${\varGamma }_{C}^2$上距离点${{{\boldsymbol{x}}}_s}$最近的点, 通过正交投影法则进行确定[13,28,30]
$$ \frac{{{{{\boldsymbol{x}}}_s} - {\bar {\boldsymbol{x}}}\left( {{{\bar \xi }^1},{{\bar \xi }^2}} \right)}}{{\left\| {{{{\boldsymbol{x}}}_s} - {\bar {\boldsymbol{x}}}\left( {{{\bar \xi }^1},{{\bar \xi }^2}} \right)} \right\|}} \cdot {{\bar {\boldsymbol{a}}}_\alpha }\left( {{{\bar \xi }^1},{{\bar \xi }^2}} \right) = 0 $$ (8) 式中, $ {\bar \xi ^\alpha } $表示主接触点${\bar {\boldsymbol{x}}}$的局部坐标.
法向基矢量$ {\bar {\boldsymbol{n}}} $可以通过切向基矢量的外积进行确定
$$ {\bar {\boldsymbol{n}}} = \frac{{{{{\bar {\boldsymbol{a}}}}_1} \times {{{\bar {\boldsymbol{a}}}}_2}}}{{\left\| {{{{\bar {\boldsymbol{a}}}}_1} \times {{{\bar {\boldsymbol{a}}}}_2}} \right\|}} $$ (9) 法向间隙${g_N}$由点${\boldsymbol{x}_s}$到${\bar {\boldsymbol{x}}}$距离的法向投影计算得到
$$ {g_N} = {\left( {{\boldsymbol{x}_s} - {\bar {\boldsymbol{x}}}} \right)^{\text{T}}}{\bar {\boldsymbol{n}}} $$ (10) 在考虑摩擦接触后, 需要引入点${\boldsymbol{x}_s}$在主接触面${\varGamma }_{C}^2$上的相对切向滑移
$$ {g_T} = \int_{{t_0}}^t {\left\| {{{{\dot {\boldsymbol{g}}}}_T}} \right\|} {\text{ d}}t = \int_{{t_0}}^t {\left\| {{{\dot {\bar {\xi}} }^\alpha }{{{\bar {\boldsymbol{a}}}}_\alpha }} \right\|} {\text{ d}}t = \int_{{t_0}}^t {\sqrt {{{\dot {\bar {\xi}} }^\alpha }{{\dot {\bar {\xi}} }^\beta }{{\bar a}_{\alpha \beta }}} } {\text{ d}}t $$ (11) 式中, $ {\bar a_{\alpha \beta }} = {{\bar {\boldsymbol{a}}}_\alpha } \cdot {{\bar {\boldsymbol{a}}}_\beta } $表示主接触面${\varGamma }_{C}^2$在点${\bar {\boldsymbol{x}}}$处的度量张量. 方程(11)中局部坐标$ {\bar \xi ^\beta } $的时间导数$ {\dot {\bar {\xi}} ^\beta } $可以由方程(8)得到
$$ {\dot {\bar {\xi}} ^\beta } = {\bar a^{\alpha \beta }}\left( {{{{\dot {\boldsymbol{x}}}}_s} - {\dot {\bar {\boldsymbol{x}}}}} \right) \cdot {{\bar {\boldsymbol{a}}}_\alpha } = {\dot g_{{T_\alpha }}} $$ (12) 上式只取了线性项, 且$ {\bar a^{\alpha \beta }} = {\left( {{{\bar a}_{\alpha \beta }}} \right)^{ - 1}} $. 在主接触面${\varGamma }_{C}^2$上点${\bar {\boldsymbol{x}}}$处的接触力可以分解为法向部分和切向部分
$$ {{\boldsymbol{t}}} = {{{\boldsymbol{t}}}_N} + {{{\boldsymbol{t}}}_T} = {{{t}}_N}{\bar {\boldsymbol{n}}} + {{{t}}_{{T_\alpha }}}{{\bar {\boldsymbol{a}}}^\alpha } $$ (13) 式中, $ {t_N} $和$ {t_{{T_\alpha }}} $表示局部法向和切向接触力, $ {{\bar {\boldsymbol{a}}}^\alpha } = {\bar a_{\alpha \beta }}{{\bar {\boldsymbol{a}}}_\beta } $表示逆变基向量. 为防止从节点${{x}_s}$的穿透, 法向接触部分需满足Kuhn-Tucker条件
$$ {g_N} \geqslant 0,\quad{\text{ }}{t_N} \leqslant 0, \quad{\text{ }}{t_N}{g_N} = 0 $$ (14) Kuhn-Tucker条件的第一式表示从节点允许分离但不允许穿透主接触面; 第二式表示法向接触力始终为压力; 第三式表示接触发生时才会产生接触力.
在切向方向上, 需要区分黏着(stick)和滑移(slip)两种情况. 根据库伦摩擦模型, 可以定义接触状态函数
$$ \varPhi = \left\| {{{{\boldsymbol{t}}}_T}} \right\| - \mu \left| {{{\boldsymbol{t}}_N}} \right| $$ (15) 式中, μ表示摩擦系数. 当$ \varPhi \leqslant 0 $时, 处于黏着状态; 当$ \varPhi > 0 $时, 处于滑移状态.
采用罚函数法计算接触力, 则法向接触力可以表示为
$$ {t_N} = {\omega _N}\left| {{g_N}} \right| $$ (16) 式中, $ {\omega _N} $表示法向罚因子. 切向接触力需要区分黏着和滑移两种状态
$$ {t_{{T_\alpha }}} = \left\{\begin{aligned} & t_{{T_\alpha }}^{tr}{\text{ , }}\quad \varPhi \leqslant 0 \\ & \mu \left| {{t_N}} \right|{p_{{T_\alpha }}}{\text{ }}{ }{\text{,}}\quad { }\varPhi > 0 \end{aligned} \right. $$ (17) 式中, $ {p_{{T_\alpha }}} = t_{{T_\alpha }}^{tr}/\left\| {{{\boldsymbol{t}}}_T^{tr}} \right\| $表示切向接触力的方向. 假设第$n$步的切向接触力$ {t_{T_\alpha ^n}} $已知, 则第$n + 1$步的预估切向接触力$ t_{T_\alpha ^{n + 1}}^{tr} $可以表示为
$$ t_{T_\alpha ^{n + 1}}^{tr} = t_{T_\alpha ^n}^{tr} + {\omega _T}{\bar a_{\alpha \beta }}\left( {\bar \xi _{n + 1}^\beta - \bar \xi _n^\beta } \right) $$ (18) 式中, $ {\omega _T} $表示切向罚因子.
2.2 压力相关摩擦系数
实际情况中, 摩擦系数μ并不是一个恒定的常数, 可能与接触压力、滑动位移和滑动速度等变量相关. 本文采用接触压力相关摩擦系数模型, 并使用常数、线性和指数等不同模型进行对比分析.
在常数模型中, 摩擦系数μ等于一个常数${\mu _c}$. 在线性模型中, 接触压力与摩擦系数保持着线性关系, 即[27]
$$ \mu = {k_\mu }{P_N} + b $$ (19) 式中, ${k_\mu }$表示比例常数, $b$表示线性常数. 指数模型是基于Vollertsen等[31]的实验结果得到的, 可以表示为
$$ \mu = {C_1} + {C_2}{{\mathrm{e}}^{ - {P_N}{C_4}}} + {C_3}{{\mathrm{e}}^{ - {P_N}{C_5}}} $$ (20) 式中, $ {C}_{1} \sim {C}_{5} $均为材料常数, 通过实验进行确定.
3种摩擦系数模型的参数如表1所示[27], 摩擦系数随接触压力的变化如图2所示. 可以看到, 线性模型的摩擦系数随接触压力的增加而增加, 指数模型的摩擦系数随接触压力的增加而减小, 最后摩擦系数收敛于0.25.
Parameter Value Parameter Value ${\mu _c}$ 0.3 ${C_2}$ 0.12 $ {k_\mu } $ 0.000 2 ${C_3}$ 0.16 $b$ 0.2 ${C_4}$ 0.006 ${C_1}$ 0.25 ${C_5}$ 0.007 ![]()
图 2 3种模型的摩擦系数随接触压力变化的情况Figure 2. The variation of the friction coefficient with contact pressure for three models2.3 基于点到面的接触离散化
接触问题相比于非接触问题只是增加了接触面上的约束条件, 将物体${{\varOmega }_1}$和物体${{\varOmega }_2}$作为求解区域, 接触界面${\varGamma }_{C}^1$和${\varGamma }_{C}^2$可以视为给定面力边界. 在几何非线性情况下, 与平衡条件等效的虚功方程可以表示为[13, 28]
$$ \begin{split} &\int_{{{\varOmega }^0}} {{{({\delta }{{\boldsymbol{E}}})}^{\text{T}}}{{\boldsymbol{S}}}{\text{d}}{\varOmega }} - \int_{{{\varOmega }^0}} {{{({\delta }{{\boldsymbol{u}}})}^{\text{T}}}{{\boldsymbol{f}}}{\text{d}}{\varOmega }} -\\ &\qquad \int_{{\varGamma }_{^{\sigma }}^0} {{{({\delta }{{\boldsymbol{u}}})}^{\text{T}}}{{\boldsymbol{T}}}{\text{d}}{\varGamma }} - {W_C} = 0 \end{split} $$ (21) 式中, ${W_C}$表示接触力虚功, $ \varGamma _C^0 = \varGamma _C^1 + \varGamma _C^2 $表示初始构型对应的接触力边界. ${W_C}$可以进一步表示为
$$ {W_C} = \int_{{\varGamma }_C^0} {\left( {{t_N}\delta {g_N} + {{{\boldsymbol{t}}}_T} \cdot \delta {{{\boldsymbol{g}}}_T}} \right){\text{d}}{\varGamma }} $$ (22) 为将接触力虚功方程(22)与非线性迭代格式(7)进行统一, 需要对方程(22)进行增量线性化
$$\begin{split} & \Delta {W_C} = \int_{{\varGamma }_C^0} {\Delta \left( {{t_N}\delta {g_N}} \right){\text{d}}{\varGamma }} + \int_{{\varGamma }_C^0} {\Delta \left( {{{{\boldsymbol{t}}}_T} \cdot \delta {{{\boldsymbol{g}}}_T}} \right){\text{d}}{\varGamma }} =\\ &\qquad \Delta W_{CN} + \Delta W_{CT}^{}\end{split} $$ (23) 增量方程(23)的具体形式略有复杂, 不是本文的研究重点, 具体细节见Wriggers[13]的著作.
在第1节中, 已经得到了非线性有限元的求解格式(7). 在罚函数接触算法中, 可以直接将其拓展为接触非线性有限元的求解格式
$$ {\text{(}}{{\boldsymbol{K}}}_L^{} + {{\boldsymbol{K}}}_{NL}^{}{\text{ + }}{{{\boldsymbol{K}}}_C}{\text{)}}\Delta {{{\boldsymbol{a}}}^{}} = {{{\boldsymbol{P}}}^{}} - {{\boldsymbol{F}}}_{\text{int} }^{} + {{{\boldsymbol{F}}}_C} $$ (24) 式中, ${{{\boldsymbol{K}}}_C}$表示系统罚函数接触刚度, ${\boldsymbol{F}_C}$表示系统罚函数接触力向量. 采用如图1所示的点到面离散化形式, 定义投影向量[32]
$$ {{\boldsymbol{N}}} = \left[ {\begin{array}{*{20}{c}} {{\bar {\boldsymbol{n}}}} \\ { - {N_1}{\bar {\boldsymbol{n}}}} \\ { - {N_2}{\bar {\boldsymbol{n}}}} \\ { - {N_3}{\bar {\boldsymbol{n}}}} \\ { - {N_4}{\bar {\boldsymbol{n}}}} \end{array}} \right],{\text{ }}{{{\boldsymbol{T}}}_\alpha } = \left[ {\begin{array}{*{20}{c}} {{{{\bar {\boldsymbol{a}}}}_\alpha }} \\ { - {N_1}{{{\bar {\boldsymbol{a}}}}_\alpha }} \\ { - {N_2}{{{\bar {\boldsymbol{a}}}}_\alpha }} \\ { - {N_3}{{{\bar {\boldsymbol{a}}}}_\alpha }} \\ { - {N_4}{{{\bar {\boldsymbol{a}}}}_\alpha }} \end{array}} \right],{\text{ }}{{{\boldsymbol{D}}}^\alpha } = {\bar a^{\alpha \beta }}{{{\boldsymbol{T}}}^\alpha } $$ (25) 式中, ${N_k}$表示接触面节点$k$的形函数, 则接触力${\boldsymbol{F}_C}$可以表示为
$$ {\boldsymbol{F}_C} = {t_N}{{\boldsymbol{N}} + }{t_{{T_\alpha }}}{\boldsymbol{D}^\alpha } $$ (26) 刚度矩阵${\boldsymbol{K}_C}$可以分解为法向${\boldsymbol{K}_{CN}}$和切向${\boldsymbol{K}_{CT}}$两个部分
$$ {\boldsymbol{K}_C} = {\boldsymbol{K}_{CN}} + {\boldsymbol{K}_{CT}} $$ (27) 法向接触刚度可以表示为
$$ {\boldsymbol{K}_{CN}} = {\omega _N}{{\boldsymbol{N}}}{{{\boldsymbol{N}}}^{\text{T}}} $$ (28) 对于切向接触刚度${\boldsymbol{K}_{CT}}$, 需要区分黏着和滑移两种状态, 对于黏着状态
$$ {{\boldsymbol{K}}}_{CT}^{} = {\omega _T}{\bar a_{\alpha \beta }}{\boldsymbol{D}^\alpha }{\boldsymbol{D}^{\beta {\text{T}}}} $$ (29) 对于滑移状态
$$ {{\boldsymbol{K}}}_{CT}^{} = \mu {\omega _N}{p_{{T_\alpha }}}{\boldsymbol{D}^\alpha }{\boldsymbol{N}^{\text{T}}} + \frac{{\mu {\omega _T}{t_N}}}{{\left\| {{{\boldsymbol{t}}}_T^{tr}} \right\|}}{\bar a_{\beta \gamma }}\left( {{\delta _{\alpha \beta }} - {p_{{T_\alpha }}}p_T^\beta } \right){\boldsymbol{D}^\alpha }{\boldsymbol{D}^{\gamma {\text{T}}}} $$ (30) 式中, $ {\delta _{\alpha \beta }} $表示Kronecker函数. 由于库伦摩擦模型是非关联本构方程, 故滑移状态的接触刚度矩阵式(30)是非对称的.
3. 考虑变摩擦系数的自适应罚因子(VFAPF)算法
在罚函数接触算法中, 由法向罚因子${\omega _N}$和切向罚因子${\omega _T}$控制算法的效率和精度, 考虑自适应罚因子可以更好地维持效率和精度之间的平衡. Lee[33]考虑了一种基于法向穿透量的自适应罚因子(APF)算法, 设罚因子${{\boldsymbol{\omega}} } = {\left[ {{\omega _N},{\omega _T}} \right]^{\text{T}}}$, 则第$n + 1$增量步的罚因子${{{\boldsymbol{\omega}} }^{n + 1}}$可以通过第$ n $步的罚因子${\boldsymbol{\omega }^n}$计算得到
$$ \left. \begin{aligned} & {\boldsymbol{\omega }^{n + 1}} = f(R){\boldsymbol{\omega }^n} \\ & f(R) = {\alpha ^{{\text{INT(}}{{\log }_\alpha }R{\text{)}}}} \end{aligned} \right\} $$ (31) 式中, $f(R)$表示以$\alpha $为底的指数函数, INT表示取整函数. 在迭代过程中, 罚因子比值$\lambda = {\omega _N}/{\omega _T}$始终保持着同一个常数, 保证了两个连续解之间的摩擦兼容性.
在Lee[33]的自适应方法中, 自适应参数$R = {R_N}$为法向穿透量的函数
$$ {R_N} = \left\{\begin{aligned} & {\frac{{{g_N}}}{{g_N^{\max }}}{\mkern 1mu} {\text{ }}{\mkern 1mu} ,{\text{ }}{g_N} > g_N^{\max }} \\ & \frac{{{g_N}}}{{g_N^{\min }}}{\mkern 1mu} {\mkern 1mu} {\text{ }},{\text{ }}{g_N} < g_N^{\min } \\ &1{\text{ }},{\text{ otherwise}} \end{aligned} \right. $$ (32) 式中, $ g_N^{\min } $和$ g_N^{\max } $分别表示最小和最大法向穿透量.
下面将进一步考虑切向弹性(黏着)滑移和变摩擦系数对自适应算法的影响. 参考式(32), 引入切向自适应参数${R_T}$
$$ {R_T} = \left\{\begin{aligned} & {\frac{{g_T^e}}{{g_T^{\max } + g_T^\mu }}{\mkern 1mu} {\text{ }}{\mkern 1mu} ,{\text{ }}g_T^e > g_T^{\max } + g_T^\mu } \\ & \frac{{g_T^e}}{{g_T^{\min } + g_T^\mu }}{\mkern 1mu} {\mkern 1mu} {\text{ }},{\text{ }}g_T^e < g_T^{\min } + g_T^\mu \\ &1{\text{ }},{\text{ otherwise}} \end{aligned} \right. $$ (33) 式中, $ g_T^{\min } $和$ g_T^{\max } $分别表示最小和最大弹性滑移量, $ g_T^e $表示当前弹性滑移量, $ g_T^\mu $表示摩擦系数变化引起的弹性滑移量.
当前弹性滑移量$ g_T^e $可以由当前的切向接触力和切向罚因子$ {\omega _T} $进行确定
$$ g_T^e = \frac{{\left\| {{{{\boldsymbol{t}}}_T}} \right\|}}{{{\omega _T}}} $$ (34) 在变摩擦系数模型中, 同一个法向接触力的极限摩擦力会随着摩擦系数的变化而变化, 因此临界弹性滑移量$ g_T^{\min } $和$ g_T^{\max } $需要适应这种摩擦系数的变化关系.
如图3所示, ${\mu _0}$和${\mu _1}$分别表示参考摩擦系数和当前摩擦系数, ${t_{{N_0}}}$和${t_{{N_1}}}$分别表示最小和最大弹性滑移量对应的法向接触力. 随着摩擦系数从${\mu _0}$变化到${\mu _1}$, 临界弹性滑移量也会从$ g_{{N_0}}^{\min } $和$ g_{{N_0}}^{\max } $平移至$ g_{{N_1}}^{\min } $和$ g_{{N_1}}^{\max } $, 平移量为$ g_T^\mu $, 由几何关系可知
![]()
图 3 临界弹性滑移量随摩擦系数变化的几何关系Figure 3. The geometric relationship between the critical elastic slip and the variation of the friction coefficient$$ g_T^\mu = \Delta \mu \lambda \left| {{g_N}} \right| $$ (35) 式中, $ \Delta \mu = {\mu _1} - {\mu _0} $表示摩擦系数变化量.
综合考虑法向穿透量和弹性滑移量的影响, 自适应参数$R$可以表示为
$$ R = \max \{ {R_N},{R_T}\} $$ (36) 上式考虑了精度优先, 而考虑效率优先或者综合影响可以采用$ R = \min \{ {R_N},{R_T}\} $或 $ R = ({R_N} + {R_T})/2 $的形式.
同时, 还需要设置一个调整率$\theta $, 避免罚因子过大或过小
$$ \frac{{{{{\boldsymbol{\omega}} }^0}}}{\theta } \leqslant {{{\boldsymbol{\omega}} }^n} \leqslant \theta {{{\boldsymbol{\omega}} }^0} $$ (37) 式中, $ {{{\boldsymbol{\omega}} }^0} $表示初始罚因子. 在本文中调整率$\theta $ = 100, 参考摩擦系数${\mu _0}$ = 0.3.
4. Archard磨损模型
在接触磨损过程中, 材料的不可逆损伤源于微观尺度的多重机制协同作用, 如硬质磨粒的犁削和黏着结点的剪切失效、循环载荷下的疲劳裂纹扩展, 以及环境介质的化学腐蚀等. 这些机制共同导致表面材料剥离、承载能力退化, 并形成自加速的磨损循环. Archard[11]基于对宏观磨损现象的统计描述, 发现磨损速率${\dot g_w}$与相对滑动速度${{\dot {\boldsymbol{g}}}_T}$和接触压力${P_N}$的乘积成正比
$$ {\dot g_w} = {k_{\text{A}}}{P_N}\left\| {{{{\dot {\boldsymbol{g}}}}_T}} \right\| $$ (38) 式中, ${k_{\text{A}}}$表示Archard磨损系数, 该系数隐含材料硬度、韧性及环境响应的综合影响.
为了将Archard磨损模型应用于非线性有限元的求解中, 需要将式(38)进行增量线性化, 磨损间隙增量$ \Delta g_w^{} $可以近似表示为
$$ \Delta g_w^{} = {k_{\text{A}}}P_N^{}\sqrt {{{\left( {\Delta {g_{{T_1}}}} \right)}^2} + {{\left( {\Delta {g_{{T_2}}}} \right)}^2}} $$ (39) 式中, $ \Delta {{g}_{{{T}_{\alpha }}}} $表示滑移增量. 在每次增量步迭代结束后对上一时刻的磨损间隙$ g'_w $进行更新
$$ g_w^{} = g'_w + \Delta g_w^{} $$ (40) 同时, 两个接触体之间法向间隙${g_N}$应考虑磨损间隙$ g_w^{} $的影响, 即
$$ {g_N} \leftarrow {g_N} + {g_w} $$ (41) 5. 数值算例
5.1 半圆平板模型的小滑移接触磨损
如图4(a)所示, 一个半圆平板模型, 半圆半径$r = 6{\text{ mm}}$, 平板高$h = 4{\text{ mm}}$, 长$l = 12{\text{ mm}}$. 平板下端固定, 在半圆上端竖直方向施加$s = $$ - 0.008{\text{ mm}}$的预位移, 水平方向施加$U = $$ \pm 0.005{\text{ mm}}$的交变位移. 采用线弹性材料, 杨氏模量$E = 2.1 \times {10^5}{\text{ MPa}}$, 泊松比$\nu = 0.3$, Archard磨损系数${k_{\text{A}}} = 8.5 \times {10^{ - 9}}{\text{ MP}}{{\text{a}}^{ - 1}}$[34], 摩擦系数采用第2.2节的模型. 考虑平面应变问题, 网格划分参考了Ning等[35]的模型, 如图4(b)所示, 全局网格尺寸为$0.5{\text{ mm}}$, 局部网格尺寸为$0.01{\text{ mm}}$.
本算例是一个Hertz接触问题, 将程序计算结果(FEM)与Hertz解[4,36]进行对比, 如图5所示. 可以看到, FEM解与Hertz解吻合良好, 验证了结果的准确性.
![]()
图 5 半圆平板模型的接触压力对比Figure 5. Comparison of contact pressure in the semicircular plate model将半圆平板模型100次循环载荷结束时3种摩擦系数模型的磨损间隙与Ning等[27]使用双势接触算法的计算结果进行对比, 如图6所示. 可以看到, 虽然中间部分的磨损间隙有着略微的误差, 但是整体上的计算结果与Ning的计算结果基本保持一致. 此外, 可以发现3种模型的磨损间隙分布有着明显的差异性, 但磨损位置始终保持在同一处. 常数和指数模型的磨损间隙分布呈现着“拱形”, 而线性模型的接触面中间区域几乎没有磨损. 图7是半圆平板模型100次循环载荷结束时3种摩擦系数模型的接触压力. 可以看到, 不同摩擦模型对接触压力的影响比较小, 几乎可以忽略.
![]()
图 6 半圆平板模型100次循环载荷结束时3种压力相关摩擦系数模型的磨损量对比Figure 6. Comparison of wear for three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the semicircular plate model![]()
图 7 半圆平板模型100次循环载荷结束时3种压力相关摩擦系数模型的接触压力对比Figure 7. Comparison of contact pressure of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the semicircular plate model图8是半圆平板模型100次循环载荷结束时3种摩擦系数模型的Mises应力分布图及其局部放大图. 可以发现, 常数模型和指数模型的应力分布基本保持一致, 最大应力分别为8.408 × 102和8.411 × 102 MPa, 指数模型最大应力只比常数模型大了0.036%, 差别几乎可以忽略. 线性模型的应力分布比较集中且最大应力区的分布与其他两种模型有着较大差别, 最大应力为8.823 × 102 MPa, 比常数模型大了4.936%.
![]()
图 8 半圆平板模型100次循环载荷结束时3种压力相关摩擦系数模型的Mises应力对比Figure 8. Comparison of Mises Stress of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the semicircular plate5.2 平面橡胶圈模型的大滑移接触磨损
超弹性材料是一类具有特殊力学性质的材料, 其特点是在应力作用下可以发生大变形, 但在卸载后能够完全恢复到初始形状[37]. Mooney-Rivlin材料是典型的超弹性材料之一, 常用于模拟橡胶、轮胎等. 在近似不可压缩性条件下, Mooney-Rivlin的应变能密度函数可以定义为[38]
$$ W({J_1},{J_2},{J_3}) = {C_{10}}({J_1} - 3) + {C_{01}}({J_2} - 3) + \frac{K}{2}{({J_3} - 1)^2} $$ (42) 式中, $ {C_{10}} $和$ {C_{01}} $表示材料常数, $ K $表示体积模量.
如图9所示的橡胶圈, 外圈半径${r_1} = 10{\text{ mm}}$, 内圈半径${r_2} = 8{\text{ mm}}$, 上端竖直方向受到$s = 4{\text{ mm}}$的预位移, 水平方向受到$U = 0\sim 20{\text{ mm}}$的交变位移, 橡胶圈与刚性地面接触. 考虑平面应变问题, 采用Mooney-Rivlin超弹性材料, 材料参数${C_{10}} = 2000{\text{ MPa}}$, ${C_{01}} = 500{\text{ MPa}}$, $K = 1.0 \times {10^7}{\text{ MPa}}$. Archard磨损系数${k_{\text{A}}} = 8.5 \times {10^{ - 9}}{\text{ MP}}{{\text{a}}^{ - 1}}$, 摩擦系数采用第2.2节的模型.
为了比较VFAPF算法与APF算法的效率和精度, 考虑如表2所示的模型参数. 径向单元尺寸统一为0.4 mm, 考虑不同的周向网格尺寸lm、一次循环载荷的加载步数Nc、初始法向罚因子${\omega _N}$和罚因子比${\omega _T}/{\omega _N}$. 对照参数为周向网格尺寸lm = 0.2 mm, 加载步数Nc = 400, ${\omega _N} = 1.0 \times {10^5}$以及${\omega _T}/{\omega _N} = 0.1$.
表 2 橡胶圈模型的参数设置Table 2. Parameter settings of the rubber ring modelParameter Value lm/mm {M1, M2, M3} = {0.1, 0.2, 0.4} Nc {S1, S2, S3} = {200, 400, 800} ${\omega _N}$ {P1, P2, P3} = {104, 105, 106} ${\omega _T}/{\omega _N}$ {R1, R2, R3} = {0.1, 0.3, 1.0} 首先考虑无磨损常摩擦系数的橡胶圈模型, 将水平位移U = 0, 1, 2, 20 mm时的接触压力(FEM)与ABAQUS的计算结果(Ref)进行对比, 如图10所示. 可以看到, FEM的接触压力与ABAQUS的接触压力吻合良好, 验证了计算结果的准确性.
在方程(31)中, 取底数$\alpha = 2$的指数函数, 临界穿透量$ g_N^{\max } = g_T^{\max } = 5.0 \times {10^{ - 4}} {l_c} $ 和 $ g_N^{\min } = g_T^{\min } $$ = 5.0 \times {10^{ - 5}} {l_c} $, ${l_c}$表示接触面的特征长度, ${l_c}$取为从接触面的平均长度, 其他模型参数设置如表2所示. 图11是橡胶圈模型不同参数下VFAPF算法和APF算法的迭代次数、法向穿透量和弹性滑移量的对比. 从图11(a)可以看到, VFAPF算法在线性摩擦模型中的迭代次数较APF算法有一定程度的增加(约12%), VFAPF算法在指数摩擦模型中与APF算法的迭代次数基本一致. 此外, 将拉格朗日乘子接触算法的迭代次数与这两种算法进行了对比, 发现拉格朗日乘子法的迭代次数都高于这两种算法, 但拉格朗日乘子接触算法可以使接触约束精确满足.
![]()
图 11 橡胶圈模型不同参数下VFAPF算法和APF算法的迭代次数、法向穿透量和弹性滑移量对比Figure 11. Comparison of the VFAPF algorithm and APF algorithm in the rubber ring model under different parameters regarding number of iterations, normal penetration, and elastic slip在图11(b)和图11(c)中, 当罚因子比${\omega _T}/{\omega _N}$较小(R1)时, 弹性滑移量占主导地位(弹性滑移量大于法向穿透量), VFAPF算法的法向穿透量和弹性滑移量较APF算法有大约59% 的提升. 随着罚因子比${\omega _T}/{\omega _N}$的提高(R1 ~ R3), 法向穿透量占主导地位, VFAPF算法与APF算法精度基本保持一致.
图12是橡胶圈模型100次循环载荷结束时3种摩擦系数模型的磨损间隙对比. 可以发现, 3种模型的磨损间隙分布均不相同, 线性模型两端略小, 中间略大; 常数和指数模型比较接近, 但接触面两侧的“折点”位置也不相同, 该位置对应接触面左端接触压力的峰值点(见图13).
![]()
图 12 橡胶圈模型100次循环载荷结束时3种压力相关摩擦系数模型的磨损量对比Figure 12. Comparison of wear gap of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber ring model![]()
图 13 橡胶圈模型100次循环载荷结束时3种压力相关摩擦系数模型的接触压力对比Figure 13. Comparison of contact pressure of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber ring model将橡胶圈模型100次循环载荷结束时3种摩擦系数模型的接触压力进行对比, 如图13所示. 在接触面右端3种模型的最大接触压力相差不超过1%, 但是在接触面左端3种模型均呈现出了不一样的压力分布, 峰值点和最大接触压力均相差较大. 线性模型最大接触压力较小且峰值点更靠近接触面中段, 指数模型最大接触压力较大且峰值点更远离接触面中段. 从图14可以发现线性和指数模型接触界面的摩擦系数呈现出高度非线性, 最大摩擦系数分别达到0.507和0.53, 最小摩擦系数分别为0.201和0.25, 差值在0.3左右.
![]()
图 14 橡胶圈模型100次循环载荷结束时3种压力相关摩擦系数模型的摩擦系数对比Figure 14. Comparison of friction coefficient of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber ring model图15是橡胶圈模型100次循环载荷结束时3种压力相关摩擦系数模型的Mises应力分布. 可以发现, 常数模型和指数模型的应力分布基本保持一致, 最大应力分别为1.835 × 104和1.820 × 104 MPa, 指数模型最大应力只比常数模型小了0.817%, 差别几乎可以忽略. 从图15(b)中发现, 线性模型的应力分布相较于常数和线性模型有应力偏大区(region of excessive stress)和应力偏小区(region of low stress), 最大应力为1.777 × 104 MPa, 比常数模型大了3.161%.
![]()
图 15 橡胶圈模型100次循环载荷结束时3种压力相关摩擦系数模型的接触压力对比Figure 15. Comparison of Mises stress of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber ring model5.3 三维橡胶块模型的大滑移接触磨损
考虑如图16所示的橡胶块模型, 上端的橡胶块边长$l = 10{\text{ mm}}$, 下端为刚性支撑. 橡胶块上端面z方向受到${U_z} = - 0.25{\text{ mm}}$的预位移, 在x-y方向施加${U} = \left( { \pm 3, \pm 2,0} \right){\text{ mm}}$的交变位移, 材料参数和磨损参数与5.2节的橡胶圈模型一致. 考虑如表3所示的参数设置, 对照参数为网格尺寸lm = 2.0 mm, 加载步数Nc = 400, 初始法向罚因子${\omega _N} = 1.0 \times {10^5}$以及罚因子比${\omega _T}/{\omega _N} = 0.1$.
表 3 橡胶块模型的参数设置Table 3. Parameter settings of the rubber block modelParameter Value lm/mm {M1, M2, M3} = {1.0, 2.0, 3.3} Nc {S1, S2, S3} = {200, 400, 800} ${\omega _N}$ {P1, P2, P3} = {104, 105, 106} ${\omega _T}/{\omega _N}$ {R1, R2, R3} = {0.1, 0.3, 1.0} 将常系数摩擦的橡胶块模型1次循环载荷结束时的接触压力和磨损分布与ABAQUS的计算结果进行对比, 如图17和图18所示. FEM最大接触压力和磨损间隙分别为1.174 × 103 MPa和6.670 × 10−5 mm, ABAQUS最大接触压力和磨损间隙分别为1.163 × 103 MPa和7.059 × 10−5 mm, 最大接触压力相对误差小于1%, 最大磨损间隙相对误差约为5%, 且分布形状几乎一致, 验证了结果的准确性.
为比较VFAPF算法和APF算法在三维大滑移问题的效率与精度, 分析了橡胶块模型在不同参数(见表3)下两种算法的迭代次数、法向穿透量和弹性滑移量, 结果如图19所示. 可以看到, 不同参数下VFAPF算法较APF算法迭代次数均有一定程度的提高, 大约为12%. VFAPF算法较APF算法的法向穿透量和弹性滑移量下降了64%左右. 随罚因子比${\omega _T}/{\omega _N}$的增加(R1 ~ R3), VFAPF算法与APF算法的结果基本一致, 说明了当弹性滑移量占主导地位时, VFAPF算法有着更好的精度稳定性.
![]()
图 19 橡胶块模型不同参数下VFAPF算法和APF算法的迭代次数、法向穿透量和弹性滑移量对比Figure 19. Comparison of the VFAPF algorithm and APF algorithm in the rubber block model under different parameters regarding number of iterations, normal penetration, and elastic slip图20和图21分别是橡胶块100次循环载荷结束时, 3种系数模型的接触压力和磨损分布情况. 常数、线性和指数模型的最大接触压力分别为1.134 × 103, 1.378 × 103和1.003 × 103 MPa, 线性模型比常数模型高了21.517%, 而指数模型比常数模型低了11.552%. 常数、线性和指数模型的最大磨损量分别为6.570 × 10−3, 7.941 × 10−3和5.992 × 10−3 mm, 线性模型比常数模型高了20.868%, 而指数模型比常数模型低了8.798%. 在图21(b)中, 还可以发现线性模型的磨损分布较常数、指数模型有着较大的差异性, 线性模型的磨损分布在接触面中间区域呈现出偏小的趋势, 常数模型和指数模型磨损分布几乎一致. 这些差别说明了不同的摩擦系数模型对物体的接触压力、磨损量和磨损位置有着显著的影响.
![]()
图 20 橡胶块模型100次循环载荷结束时3种压力相关摩擦系数模型的接触压力对比Figure 20. Comparison of contact pressure of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber block model![]()
图 21 橡胶块模型100次循环载荷结束时3种压力相关摩擦系数模型的磨损分布对比Figure 21. Comparison of wear gap of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber block model6. 结论
本文基于罚函数接触算法和Archard磨损模型, 考虑非线性有限元格式, 构造了常数、线性和指数的压力相关摩擦系数模型来计算接触磨损问题. 将变摩擦系数引入到自适应罚因子算法中, 构建了复杂接触磨损问题的数值模拟方法, 并将线弹性材料和Mooney-Rivlin超弹性材料分别用于小滑移和大滑移的接触磨损问题求解中. 根据求解结果, 探讨了VFAPF算法对压力相关摩擦系数模型求解过程中效率和精度的影响效果, 还分析了压力相关摩擦系数模型对物体接触过程中磨损量、磨损位置、接触压力和应力分布的影响. 本文得到主要结论如下.
(1) 相比于APF算法, VFAPF算法在变摩擦大滑移问题中有着更好的精度稳定性. 弹性滑移量占主导地位(弹性滑移量大于法向穿透量)时, 12%的效率下降带来了约62%的精度提升; 随着罚因子比${\omega _T}/{\omega _N}$的增加, VFAPF算法与APF算法基本一致.
(2) 对于线弹性材料的小滑移接触磨损问题, 尽管3种压力相关摩擦系数模型在磨损量和应力分布上存在差异, 但对磨损位置和接触压力的影响较小.
(3) 对于超弹性材料的大滑移接触磨损问题, 3种压力相关摩擦系数模型的磨损量、磨损位置、接触压力和应力分布均表现出了显著差异性. 在平面问题中, 最突出的是对接触压力的影响, 磨损量和磨损位置有略微的影响, 且线性模型的应力分布差异最大; 在三维问题中, 线性模型的接触压力和磨损量较常数模型有20%左右的差异, 指数模型较常数模型有10%左右的差异, 且线性模型的磨损分布较常数和指数模型有着显著区别.
-
![]()
图 2 3种模型的摩擦系数随接触压力变化的情况
Figure 2. The variation of the friction coefficient with contact pressure for three models
![]()
图 3 临界弹性滑移量随摩擦系数变化的几何关系
Figure 3. The geometric relationship between the critical elastic slip and the variation of the friction coefficient
![]()
图 5 半圆平板模型的接触压力对比
Figure 5. Comparison of contact pressure in the semicircular plate model
![]()
图 6 半圆平板模型100次循环载荷结束时3种压力相关摩擦系数模型的磨损量对比
Figure 6. Comparison of wear for three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the semicircular plate model
![]()
图 7 半圆平板模型100次循环载荷结束时3种压力相关摩擦系数模型的接触压力对比
Figure 7. Comparison of contact pressure of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the semicircular plate model
![]()
图 8 半圆平板模型100次循环载荷结束时3种压力相关摩擦系数模型的Mises应力对比
Figure 8. Comparison of Mises Stress of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the semicircular plate
![]()
图 11 橡胶圈模型不同参数下VFAPF算法和APF算法的迭代次数、法向穿透量和弹性滑移量对比
Figure 11. Comparison of the VFAPF algorithm and APF algorithm in the rubber ring model under different parameters regarding number of iterations, normal penetration, and elastic slip
![]()
图 12 橡胶圈模型100次循环载荷结束时3种压力相关摩擦系数模型的磨损量对比
Figure 12. Comparison of wear gap of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber ring model
![]()
图 13 橡胶圈模型100次循环载荷结束时3种压力相关摩擦系数模型的接触压力对比
Figure 13. Comparison of contact pressure of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber ring model
![]()
图 14 橡胶圈模型100次循环载荷结束时3种压力相关摩擦系数模型的摩擦系数对比
Figure 14. Comparison of friction coefficient of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber ring model
![]()
图 15 橡胶圈模型100次循环载荷结束时3种压力相关摩擦系数模型的接触压力对比
Figure 15. Comparison of Mises stress of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber ring model
![]()
图 19 橡胶块模型不同参数下VFAPF算法和APF算法的迭代次数、法向穿透量和弹性滑移量对比
Figure 19. Comparison of the VFAPF algorithm and APF algorithm in the rubber block model under different parameters regarding number of iterations, normal penetration, and elastic slip
![]()
图 20 橡胶块模型100次循环载荷结束时3种压力相关摩擦系数模型的接触压力对比
Figure 20. Comparison of contact pressure of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber block model
![]()
图 21 橡胶块模型100次循环载荷结束时3种压力相关摩擦系数模型的磨损分布对比
Figure 21. Comparison of wear gap of three pressure-dependent friction coefficient models at the end of 100 cycles of loading in the rubber block model
Parameter Value Parameter Value ${\mu _c}$ 0.3 ${C_2}$ 0.12 $ {k_\mu } $ 0.000 2 ${C_3}$ 0.16 $b$ 0.2 ${C_4}$ 0.006 ${C_1}$ 0.25 ${C_5}$ 0.007 
下载: 导出CSV
表 2 橡胶圈模型的参数设置
Table 2 Parameter settings of the rubber ring model
Parameter Value lm/mm {M1, M2, M3} = {0.1, 0.2, 0.4} Nc {S1, S2, S3} = {200, 400, 800} ${\omega _N}$ {P1, P2, P3} = {104, 105, 106} ${\omega _T}/{\omega _N}$ {R1, R2, R3} = {0.1, 0.3, 1.0}
下载: 导出CSV
表 3 橡胶块模型的参数设置
Table 3 Parameter settings of the rubber block model
Parameter Value lm/mm {M1, M2, M3} = {1.0, 2.0, 3.3} Nc {S1, S2, S3} = {200, 400, 800} ${\omega _N}$ {P1, P2, P3} = {104, 105, 106} ${\omega _T}/{\omega _N}$ {R1, R2, R3} = {0.1, 0.3, 1.0}
下载: 导出CSV
-
[1] Li Y, Zuo SG, Lei L, et al. Analysis of impact factors of tire wear. Journal of Vibration and Control, 2011, 18(6): 833-840
[2] Kragelsky IV, Dobychin MN, Kombalov VS. Friction and Wear: Calculation Methods. Pergamon, 1982
[3] Nguyen VH, Zheng D, Schmerwitz F, et al. An advanced abrasion model for tire wear. Wear, 2018, 396: 75-85
[4] Popov VL. Contact Mechanics and Friction. Springer, 2010
[5] Vieira T, Ferreira RP, Kuchiishi AK, et al. Evaluation of friction mechanisms and wear rates on rubber tire materials by low-cost laboratory tests. Wear, 2015, 328: 556-562
[6] 邹龙庆, 黄聪聪, 付海龙等. 表面粗糙峰坐标点云重构的金属-橡胶接触分析. 表面技术, 2021, 50(10): 255-262 (Zou Longqing, Huang Congcong, Fu Hailong, et al. Metal-rubber rigid soft contact analysis based on gaussian rough surface. Surface Technology, 2021, 50(10): 255-262 (in Chinese) Zou Longqing, Huang Congcong, Fu Hailong, et al. Metal-rubber rigid soft contact analysis based on gaussian rough surface. Surface Technology, 2021, 50(10): 255-262 (in Chinese)
[7] Kahms S, Wangenheim M. Experimental investigation and simulation of aircraft tire wear. Tire Science and Technology, 2021, 49(1): 55-74 doi: 10.2346/tire.20.180201
[8] Sakhnevych A, Genovese A. Tyre wear model: A fusion of rubber viscoelasticity, road roughness, and thermodynamic state. Wear, 2024, 542: 205291
[9] Huang M, Guibert M, Thévenet J, et al. A new test method to simulate low-severity wear conditions experienced by rubber tire materials. Wear, 2018, 410: 72-82
[10] 付豪, 尧军平, 梁超群等. 基于Archard磨损模型研究SiC/AZ91D复合材料干摩擦磨损特性. 复合材料学报, 2024, 41(11): 6206-6214 (Fu Hao, Yao Junping, Liang Chaoqun, et al. Research on dry friction and wear characteristics of SiC/AZ91D composites based on Archard wear mode. Acta Materiae Compositae Sinica, 2024, 41(11): 6206-6214 (in Chinese) Fu Hao, Yao Junping, Liang Chaoqun, et al. Research on dry friction and wear characteristics of SiC/AZ91D composites based on Archard wear mode. Acta Materiae Compositae Sinica, 2024, 41(11): 6206-6214 (in Chinese)
[11] Archard JF. Contact and rubbing of flat surfaces. Journal of Applied Physics, 1953, 24(8): 981-988 doi: 10.1063/1.1721448
[12] Fouvry S, Kapsa P, Vincent L. Quantification of fretting damage. Wear, 1996, 200(1-2): 186-205 doi: 10.1016/S0043-1648(96)07306-1
[13] Wriggers P. Computational Contact Mechanics. Springer, 2006
[14] Weyler R, Oliver J, Sain T, et al. On the contact domain method: A comparison of penalty and Lagrange multiplier implementations. Computer Methods in Applied Mechanics and Engineering, 2012, 205(15): 68-82
[15] Papadopoulos P, Solberg JM. A Lagrange multiplier method for the finite element solution of frictionless contact problems. Mathematical and Computer Modelling, 1998, 28(4-8): 373-384 doi: 10.1016/S0895-7177(98)00128-9
[16] Pietrzak G, Curnier A. Large deformation frictional contact mechanics: continuum formulation and augmented Lagrangian treatment. Computer Methods in Applied Mechanics and Engineering, 1999, 177(3-4): 351-381 doi: 10.1016/S0045-7825(98)00388-0
[17] Simo JC, Wriggers P, Taylor RL. A perturbed Lagrangian formulation for the finite element solution of contact problems. Computer Methods in Applied Mechanics and Engineering, 1985, 50(2): 163-180 doi: 10.1016/0045-7825(85)90088-X
[18] De Saxcé G, Feng ZQ. New inequality and functional for contact with friction: the implicit standard material approach. Mechanics of Structures and Machines, 1991, 19(3): 301-325 doi: 10.1080/08905459108905146
[19] De Saxcé G, Feng ZQ. The bipotential method: a constructive approach to design the complete contact law with friction and improved numerical algorithms. Mathematical and Computer Modelling, 1998, 28(4-8): 225-245 doi: 10.1016/S0895-7177(98)00119-8
[20] Bernard JE. The simulation of coulomb friction in mechanical systems. Simulation, 1980, 34(1): 11-16 doi: 10.1177/003754978003400104
[21] Quiñones-Cisneros SE, Deiters UK. Generalization of the friction theory for viscosity modeling. The Journal of Physical Chemistry B, 2006, 110(25): 12820-12834 doi: 10.1021/jp0618577
[22] Gao JP, Luedtke WD, Gourdon D, et al. Frictional forces and Amontons' law: from the molecular to the macroscopic scale. The Journal of Physical Chemistry B, 2004, 108(11): 3410-3425 doi: 10.1021/jp036362l
[23] Lischinsky P, Canudas-De-Wit C, Morel G. Friction compensation for an industrial hydraulic robot. IEEE Control Systems Magazine, 1999, 19(1): 25-32 doi: 10.1109/37.745763
[24] Luo Q. Temperature dependent friction and wear of magnetron sputtered coating TiAlN/VN. Wear, 2011, 271(9-10): 2058-2066 doi: 10.1016/j.wear.2011.01.054
[25] Mäntylä A, Hintikka J, Frondelius T, et al. Prediction of contact condition and surface damage by simulating variable friction coefficient and wear. Tribology International, 2020, 143: 106054 doi: 10.1016/j.triboint.2019.106054
[26] Yue TY, Wahab MA. Finite element analysis of fretting wear under variable coefficient of friction and different contact regimes. Tribology International, 2017, 107: 274-282 doi: 10.1016/j.triboint.2016.11.044
[27] Ning P, Li Y, Feng Z. A Newton-like algorithm to solve contact and wear problems with pressure-dependent friction coefficients. Communications in Nonlinear Science and Numerical Simulation, 2020, 85: 105216 doi: 10.1016/j.cnsns.2020.105216
[28] 王勖成. 有限单元法. 北京: 清华大学出版社, 2003 (Wang Xucheng. Finite Element Method. Beijing: Tsinghua University Press, 2003 (in Chinese) Wang Xucheng. Finite Element Method. Beijing: Tsinghua University Press, 2003 (in Chinese)
[29] Kim NH. Introduction to Nonlinear Finite Element Analysis. Springer Science & Business Media, 2014
[30] Vulovic S, Zivkovic M, Grujovic N, et al. A comparative study of contact problems solution based on the penalty and Lagrange multiplier approaches. Journal of the Serbian Society for Computational Mechanics, 2007, 1(1): 174-183
[31] Vollertsen F, Hu ZY. Determination of size-dependent friction functions in sheet metal forming with respect to the distribution of the contact pressure. Production Engineering, 2008, 2(4): 345-350 doi: 10.1007/s11740-008-0130-4
[32] Laursen TA, Simo JC. A continuum-based finite element formulation for the implicit solution of multibody, large deformation-frictional contact problems. International Journal for Numerical Methods in Engineering, 1993, 36(20): 3451-3485 doi: 10.1002/nme.1620362005
[33] Lee SH. Rudimentary considerations for adaptive gap/friction element based on the penalty method. Computers & Structures, 1993, 47(6): 1043-1056
[34] Stromberg N. An augmented Lagrangian method for fretting problems. European Journal of Mechanics, A/Solids, 1997, 16(4): 573-593
[35] Ning P, Feng ZQ, Quintero JAR, et al. Uzawa algorithm to solve elastic and elastic-plastic fretting wear problems within the bi-potential framework. Computational Mechanics, 2018, 62(6): 1327-1341 doi: 10.1007/s00466-018-1567-8
[36] McColl IR, Ding J, Leen SB. Finite element simulation and experimental validation of fretting wear. Wear, 2004, 256(11-12): 1114-1127 doi: 10.1016/j.wear.2003.07.001
[37] 彭向峰, 李录贤. 超弹性材料本构关系的最新研究进展. 力学学报, 2020, 52(5): 1221-1234 (Peng Xiangfeng, Li Luxian. State of the art of constitutive relations of hyperelastic materials. Chinese Journal of Theoretical and Applied Mechanics, 2020, 52(5): 1221-1232 (in Chinese) Peng Xiangfeng, Li Luxian. State of the art of constitutive relations of hyperelastic materials. Chinese Journal of Theoretical and Applied Mechanics, 2020, 52(5): 1221-1232 (in Chinese)
[38] Mooney M. A theory of large elastic deformation. Journal of Applied Physics, 1940, 11(9): 582-592 doi: 10.1063/1.1712836
计量
- 文章访问数: 67
- HTML全文浏览量: 10
- PDF下载量: 18
