考虑基体变形的结合面连续光滑接触刚度模型

2021-12-01王晶晶史小辉杨亚兰阮晓光

工程力学 2021年12期

李玲，王晶晶，史小辉，杨亚兰，阮晓光

(1.西安建筑科技大学机电工程学院，西安 710055；2.西安陕鼓动力股份有限公司，西安 710075)

结合面在机床运行过程中起着传递运动、载荷和能量的重要作用，同时也是引起机床非线性和不确定性的根源[1−3]。研究表明，结合面的刚度约占机床中总刚度的60%～80%，其刚度值的大小和变化规律，直接影响着机床的静、动态力学性能和精度保持性水平[4−5]。因此，开展粗糙表面接触刚度特性的研究，对整机性能的评估与动态优化具有重要的理论指导意义。

Greenwood和Williamson[6]是较早开展接触理论研究的先驱者。他们将粗糙表面上微凸体高度分布看成随机函数，并基于Hertz理论提出统计学微观弹性接触模型(GW模型)，奠定了接触力学基础。为进一步完善粗糙表面的接触理论，Chang等[7]考虑微凸体变形过程中体积守恒原则，建立了粗糙表面的弹-塑性接触模型(CEB模型)。但该模型却忽略了微凸体由弹性到塑性变形的中间过程，使得在临界屈服点处接触载荷出现不连续现象。Zhao等[8]为弥补CEB模型的不足，将微凸体的变形分为弹性变形、弹塑性变形和完全塑性变形阶段，并利用多项式数学插值法建立了粗糙表面的弹塑性接触模型(ZMC模型)。该模型增强了接触载荷和平均接触压力在弹性屈服点与塑性临界接触点的连续性，但却忽略了光滑性，导致建立的接触刚度模型在2个临界点处均是不连续的。Kogut和Estion[9−10]利用有限元法，将微凸体的弹塑性接触变形区细分为2个子区域，建立出更高精度的弹塑性接触模型(KE模型)。但在这两个子区域的转换点处和弹塑性变形的开始与结束点处，接触载荷也是不连续的，从而引起其刚度模型的不连续性。Brake[11]采用三次Hermite多项式插值函数，提出了一种在不同接触状态之间，接触载荷光滑且连续变化的粗糙表面弹塑性接触模型，但由于其多项式阶次较高，获得的接触载荷与变形量之间的关系曲线出现振荡，从而导致使用其所建立的接触刚度模型也会出现振荡，甚至出现负刚度的情况。

上述模型除了各自所存在的缺陷外，还忽略了微凸体所附着的基体的变形情况。对于一些含有硬涂层的粗糙表面，其接触特性主要由软基体上的硬涂层材料决定，此时较软基体的变形对粗糙表面的接触特性有着显著的影响[1,12 −13]。在之前的研究中[1]，仅对纯弹性接触下的基体变形进行了分析，但未包含微凸体的弹塑性和完全塑性变形阶段。为此，本文基于微凸体在弹塑性接触变形过程中，其接触载荷和接触刚度具有连续且光滑的特性，首先对Hermite多项式直接插值函数进行了改进，提出了一种在微凸体的整个变形过程中，其接触刚度满足连续且光滑变化的模型。其次，引入了基体的弹性变形，提高了传统接触刚度建模的精度。然后，基于统计学方法，将单个微凸体-基体系统模型拓展到整个粗糙表面上，建立了考虑基体变形的连续光滑结合面刚度模型。最后，与GW、ZMC、KE、Brake模型进行对比分析，验证了本文接触刚度模型的正确性。

1 单个微凸体-基体系统的接触刚度

基于以下4点假设开展结合面刚度建模研究。

1)将图1(a)两个粗糙表面的接触，等效为图1(b)所示的光滑刚性平面与等效粗糙表面的接触。

图1 粗糙表面接触模型Fig.1 Rough surface contact model

2)微凸体发生接触变形前后的体积保持不变，其变形刚度为ka；基体的变形为弹性变形，其变形刚度为kb，如图1(c)所示。

3)微凸体的分布是各向同性的，且各个微凸体间的变形相互独立。

4)微凸体高度服从高斯分布，所有微凸体的峰顶都等效成半球体，且具有相同的峰顶曲率半径。

图1(c)为单个微凸体-基体系统的受力分析。当刚性平面被施加一法向载荷f时，微凸体产生δa的变形，同时基体也发生δb宏观变形，z为微凸体变形前的高度，d为基体平均平面与刚性平面间的距离。考虑基体变形影响下系统总变形量δ可表示为：

基于Chang等[7]的研究，微凸体由弹性变形转变为弹塑性变形时的临界接触变形量可表示为：

微凸体开始发生完全塑性变形时的临界变形量δp与初始屈服点临界变形量δc的关系为[15]：

基于式(1)～式(3)，在考虑基体变形的同时，对微凸体的弹性、弹塑性以及完全塑性接触变形进行逐一分析。为便于描述，下文依次按弹性变形、塑性变形和弹塑性变形3个阶段建模分析。

1.1 弹性接触变形阶段

当微凸体变形量0＜δa＜δc时，根据Hertz接触理论，微凸体的弹性接触载荷fae可表示为：

式中，δae为微凸体弹性阶段的接触变形。由刚度的定义，可得单个微凸体发生弹性变形时的接触刚度：

当微凸体受到载荷时，微凸体和基体的接触区是一个半径为rb的圆形区域(如图1(c)所示)，且在该区域内各点受到的压力各不相同，最大接触压力P0发生在接触中心，可表示为[1]：

距离接触区域中心的距离r不同，基体表面发生的法向位移不同，当基体受到的最大压力P0不断变化时，基体表面的法向位移可表示为：

联立方程式(6)和式(8)，可获得基体的弹性刚度：

根据力的相互作用原理和叠加原理，可得：

将式(4)、式(9)和式(11)代入式(10)，可得微凸体变形δae和系统总变形δ之间的关系：

由式(11)可知，微凸体变形δae是系统总变形δ的函数，令δae=g(δ)，则式(13)可重新表示为：

为获得δae与δ的直接关系，以便下文进行基体变形影响下粗糙表面的接触建模与分析，在之前的工作中已采用不动点迭代法对式(14)进行了数值求解，获得弹性接触时，考虑基体变形影响下微凸体接触变形δae的计算表达式[16]：

1.2 完全塑性接触变形阶段

当微凸体变形量δa＞δp时，微凸体将发生完全塑性接触变形，此时微凸体的接触载荷fap和接触刚度kap可分别表示为[17]：

式中，δap表示完全塑性阶段的微凸体接触变形。同理，根据力的相互作用原理，可得：

联立式(9)、式(11)、式(17)和式(18)，可得微凸体发生完全塑性变形时的变形量δap与系统总变形δ之间的关系为：

1.3 混合弹塑性接触变形阶段

当微凸体变形量δc≤δa＜δp时，微凸体变形处于混合弹塑性阶段。通过研究发现，改进Brake模型[11]的Hermite多项式直接插值函数，来弥补插值函数阶次较高而引起微凸体接触参数振荡行为，可实现微凸体的接触载荷连续光滑变化。利用对数变换可将混合弹塑性区域插值间区间(1，110)减小为(0，4.7≈ln110)，进而很大程度降低了接触行为的振荡性，最终可得弹塑性接触变形阶段的法向接触载荷和法向接触刚度分别为：

式中：δaep为混合弹塑性阶段的微凸体接触变形；C1、C2和C3为Hermite多项式插值函数的系数项，分别表示为：其中，q=δp/δc=110，表示塑性临界接触变形与弹性临界变形之比。

同理，在微凸体的弹塑性接触变形阶段中引入基体变形的影响。根据力的相互作用原理，有：

联立式(9)、式(11)、式(20)和式(22)，可得：

为简化式(23)的数值求解过程，忽略参数值较小的C3项，式(23)可重新表示为：

进一步化简为：

式中，λ=R/4rb为确定微凸体几何形状的基本参数。

类似地，令δaep=g(δ)，则式(25)可表示为：

采用不动点迭代法可得式(26)的近似解，迭代公式为：

式中，n=0,1,···,∞，迭代的次数越多，得到的结果越精确。经过数值计算发现，只需迭代2次得到的近似解就非常接近于数值解，且两者之间的误差不超过5%。因此微凸体的变形量δa可以近似地表示为一个关于总变形δ的函数，即：

图2显示了微凸体变形量δa数值解与近似解的对比。由图可知，对于不同的微凸体几何参数λ，数值解与近似解之间的误差很小，从而证明用近似解代替数值解的有效性。从图中可以发现，对于相同的系统变形量δ，微凸体变形量δa随着λ的增大而减小，而基体的变形逐渐增大。此时，基体的接触变形在微凸体-基体系统中起主导作用。

图2 微凸体变形量数值解与近似解的对比Fig.2 Comparison of numerical solution and approximate solution of asperity deformation

2 考虑基体变形的接触刚度模型

为研究结合面的微观接触特性，需将单个微凸体-基体系统模型扩展到整个粗糙表面上，从而建立一种具有连续光滑特性并考虑基体变形的粗糙表面接触新模型。根据统计学方法和微凸体高度服从高斯分布，假设粗糙表面名义接触面积上微凸体的总数为N，则不同接触间隙d下发生接触的微凸体数量可表示为[14]：

式中：N=ηAn，η为微凸体密度，An表示名义接触面积；φ(z)为微凸体高度分布的概率密度函数。对于某一给定的表面接触间隙d，整个变形区内粗糙表面的接触载荷Fnew可表示为：

式中，Fe、Fp、Fep分别表示弹性接触载荷、塑性接触载荷及弹塑性接触载荷，分别表示如下：

将式(4)和式(15)代入式(32)，式(16)和式(19)代入式(33)，式(20)和式(28)代入式(34)，即可分别获得基体变形影响下粗糙表面的弹性接触载荷、塑性接触载荷及弹塑性接触载荷，从而将式(31)进一步表示为：

在微凸体的三个变形阶段中，考虑基体变形δb影响的无量纲微凸体接触变形可根据式(15)、式(19)和式(28)分别表示为：

无量纲概率密度函数φ(z)，可表示为：

式中，z*=z/σ，表示无量纲的微凸体高度参数，σs为微凸体高度的标准差，与粗糙表面高度的标准差σ的关系为[14]：

式中，γ=ηRσ为粗糙表面的综合粗糙度参数。

同理，整个粗糙表面的法向接触刚度Knew可表示为：

类似地，将式(5)、式(15)、式(17)、式(19)、式(21)和式(28)分别代入式(40)即可获得考虑基体变形影响下的粗糙表面接触刚度，则有：

3 分析与讨论

3.1 接触刚度模型对比分析

表1 四种不同粗糙表面形貌参数[18]Table1 Four different rough surface topography parameters

利用式(34)和式(38)，可得不同粗糙表面下接触刚度与接触载荷之间的关系，如图3所示。由图3可知，本文模型与GW、ZMC、KE、Brake模型计算的接触刚度曲线均一致，且均随着接触压力的增大而增大，从而表明本文刚度模型的合理性，可适用于不同粗糙形貌表面的接触刚度特性分析。

通过对比图3(a)～图3(d)可知：随着表面粗糙度的增大，不同接触模型在相同载荷下的接触刚度增大，会引起粗糙表面间真实接触面积减少，从而使得粗糙表面的接触刚度增大。因此，可以通过降低表面粗糙度来改善粗糙表面的接触刚度特性。

此外，从图3还发现：①对于光滑接触表面，Brake模型预测的接触刚度最大，且不同模型之间的差距较小，尤其是当接触压力较小时，它们之间的差距很小；②对于中等粗糙接触表面，依然是Brake模型的结果最大，但此时本文未考虑基体变形的模型最接近于GW模型结果；③对于粗糙和非常粗糙的表面，GW模型预测的接触刚度最大，且其与其他模型的计算结果存在较大的差异。

图3 不同粗糙表面模型下接触刚度与法向接触压力的关系Fig.3 Relationship between contact stiffness and normal contact pressure for different rough surface models

3.2 基体变形对接触刚度的影响

图4显示了不同粗糙表面形貌下基体变形对表面接触刚度的影响。由图可知，对于相同的无量纲表面接触间隙d*，与未考虑基体变形下的接触刚度值相比，考虑基体变形的接触刚度较小。这是因为，在未考虑基体变形的接触模型中，认为微凸体的变形等于系统的总变形，此时获得的微凸体变形量不是其真实的变形量，包含基体的变形量。换言之，未考虑基体的宏观变形时，计算得到的微凸体变形量大于其实际的变形，从而导致计算的接触刚度偏大，而考虑基体变形时预测的接触刚度更准确。

图4 不同粗糙表面形貌下基体变形对接触刚度的影响Fig.4 The influence of matrix deformation with different rough surface morphologies on contact stiffness

对比图4(a)～图4(d)可知，随着表面粗糙度的增大，两者之间的差距逐渐减小，即基体变形的影响减弱了，直到表面非常粗糙时，两者之间的差距很小，这与实际情况是一致的。从图中还可以发现，随着表面无量纲接触间隙的减小，两者之间的差距递增，这是因为表面接触间隙的减小意味着接触载荷的增大，此时基体变形越大。

此外，从图4中还可以获得，表面接触刚度随表面间隙的增大而减小，且其减小的速率呈先快后慢的变化趋势，这是因为当表面粗糙度与外载荷一定时，随着表面接触间隙的增大，表面间发生接触微凸体的数目减小，引起粗糙表面真实接触面积的减小，从而导致接触刚度的减小。此外，当无量纲表面接触间隙d*=d/σ增大到3时，此时粗糙表面的接触刚度非常小，趋近于0，这是因为对于微凸体高度服从高斯分布的粗糙表面，微凸体高度分布在[−3σ，3σ]的概率为99.73%，即几乎所有的微凸体高度z均小于3σ，所以当d*=3时，粗糙表面间接触的微凸体数目极少，从而导致表面接触刚度很小，几乎趋近于0。