基于多尺度下塑性指数模型的结合面接触刚度计算方法

2022-02-22赵永胜牛娜娜杨聪彬刘志峰孟令军

振动与冲击 2022年3期

赵永胜，牛娜娜，杨聪彬，刘志峰，姜凯，孟令军

(1.北京工业大学先进制造与智能技术研究所，北京 100124；2.北京工业大学先进制造技术北京市重点实验室，北京 100124；3.北京工业大学机械工业重型机床数字化设计与测试技术重点实验室，北京 100124)

实际工程设备中的栓接接触表面在微观上是由无数个微凸体组成的，真实接触面积比名义接触面积小得多。因此准确地描述结合面的接触特性对研究密封、摩擦、磨损和导电导热有着十分重要的作用[1]。

早期的解析法是Greenwood等[2-3]基于统计学建立的G-W(Greenwood-Williamson)接触模型，G-W模型中将接触表面简化为粗糙表面和一个刚性平面的接触，并假设表面微凸体在粗糙表面的水平方向上均匀分布，法向高度服从正态分布且具有相同的曲率半径。Whitehouse等[4]研究了微凸体峰高与曲率半径的相关性，引入指数形式的自相关函数系数，得到了表征粗糙表面的二维W-A(Whitehouse-Archard)接触模型。Kadin等[5]基于Etsion粗糙表面单个微凸体的有限元模型，建立了加载-卸载统计模型。微凸体在加载和卸载过程中，发生的塑性变形会导致残余应变和微凸体高度分布函数的变化。但是统计学参数的获取容易受到样本长度以及仪器分辨率的影响[6]，且有试验发现当结合面受到较大载荷时，因为微凸体附着在一个较大半径的微凸体上，微凸体的边缘可能与相邻微凸体融合以及一个微凸体的受压可能会导致其邻近的微凸体的整体向下位移。但单一尺度下的统计学研究并没有考虑到相邻微凸体的相互作用，因此统计学无法较好地反映出大载荷情况下的真实接触特性。

由于单一尺度下的建模存在不足，考虑粗糙表面的多尺度特性对接触分析仍是有必要的。为此，Ciavarella等[7]采用Weierstrass级数离散表征二维分形曲面，与原G-W模型相比，获得了较好的曲面表示效果。Persson[8]提出了另一种解决粗糙表面多尺度问题的方法，通过引入了“放大”的概念，使观察到的表面轮廓逐渐演变为精细的尺度。放大由参数量化ζ表示为:在最低倍率放大(即ζ=1)时,粗糙表面看起来是平的，没有观察到微凸体；随着放大倍数ζ的增加，较小尺度的微凸体逐渐被检测到。借助于多尺度的思想，可以从一个新的角度研究许多微凸体的接触问题，如表面粗糙度对密封泄漏率的影响。但多尺度思想的一个假设是表面之间完全接触。因此多尺度适合用于有较大载荷、精度要求较高的接触建模，但对于较大的载荷或较小的接触面积而言，精度反而会变差。

有学者研究微观接触面发现，在不同的测量尺度下，粗糙表面具有自相似性和自仿射性，由此Mandelbrot[9]提出了用分形几何的概念来模拟自然界中不规则和破碎的表面。因为分形参数不受尺度影响，所以分形理论被引入并广泛应用于结合面的接触特性研究。Berry等[10]基于Weierstrass-Mandelbrot函数表征了加工表面微凸体的二维轮廓。Majumdar等[11]首次结合Hertz接触理论，提出了一种基于分形几何的二维粗糙表面M-B(A.Majumder-B.Bhushan)接触模型。与统计学不同，M-B模型中W-M(Weierstrass-Mandelbrot)函数描述了不同尺度下微凸体的二维轮廓曲线，考虑了微凸体曲率半径对接触变形的影响。而统计学中认为微凸体具有相同的曲率半径，随着载荷的增大，微凸体的变形遵循弹性-弹塑性-塑性变形的规律。实际上微凸体发生何种变形与微凸体的曲率半径存在很大关系，与法向载荷无关。因此，分形理论的临界参数考虑到了不同尺度下微凸体的曲率半径，认为接触面积较小的微凸体发生塑性变形，接触面积较大的微凸体发生弹性变形。Zhao等[12]建立了一个修正的三维分形模型来分析栓接结构粗糙接触面的刚度和阻尼。Yuan等[13]基于分形理论，建立了三维粗糙表面之间的加载-卸载接触模型。得到了无量纲的总接触载荷与总真实接触面积在加卸载过程中的函数。Jalali等[14]提出了一种新的通用接触模型来模拟结合面。研究了不同预紧力和表面粗糙度与结合面接触特性的线性关系。最后得出预紧力和表面粗糙度分别影响接触模型的法向刚度和剪切刚度的结论。Liu等[15]提出了基于分形理论的法向动态接触刚度模型。推导了弹性变形和弹塑性变形与刚度的线性关系，建立了结合面的法向刚度模型。同时，得到了接触面动刚度与静刚度之间的关系。分析结果表明，当激励频率远小于固有频率时，接触动刚度近似等于静刚度，当激励频率远大于固有频率时，系统的动刚度远大于静刚度，同时存在使动刚度最小的共振频率。Wang等[16]提出了考虑微凸体相互作用影响的分形接触刚度模型。在弹性理论的基础上，对给定粗糙度的接触表面进行了变形推导。通过与分形模型的比较，研究了粗糙相互作用对接触刚度的影响。Liang等[17]考虑了微凸体的相对滑动对材料硬度的影响，建立了表面接触的预测模型。分析了总的接触载荷与总的真实接触面积之间的关系。但是没有考虑微凸体的弹塑性变形阶段。Liu等[18-19]考虑了重型机械受不均匀载荷的情况，根据应力特性将宏观接触角转换为微观接触角，最后，建立了一个在接触表面之间具有倾斜度的试验系统，分析了不同的分形参数，倾角值及其他重要参数对接触表面特性的影响。以上通过分形理论建立的接触模型未考虑微凸体弹塑性变形阶段，建模精度有待提高。

为解决统计学理论中存在的尺度单一、关键参数受仪器分辨率影响的问题，本文将分形参数引入统计学的塑性指数表达式，提出用塑性指数作为临界变形参数的判断依据，建立考虑微凸体弹塑性变形阶段的预测模型。最后将仿真结果与试验结果进行对比，误差在8.7%以内，验证了理论模型的正确性。

1 基于塑性指数求取临界变形参数

Mandelbrot根据布朗函数的形式构造出了一个分形函数，它满足连续性、处处不可微和自仿射的数学特性，该函数称为Weierstrass-Mandelbrot函数[20]

(1)

式中：D为分形维数；γ为频率相关常数，取γ=1.5；x为x方向坐标值；m为微凸体频率序数；mmin为最低截止频率序数；φm为相位角。

取Weierstrass-Mandelbrot函数的实部，并设相位为0，可得到现用于表面轮廓模拟的W-M函数，其表达式为

(2)

式中：G为表面粗糙度系数；z为微凸体轮廓高度。

根据W-M函数得到与尺度相关的微凸体曲率半径及微凸体变形量的表达式为

(3)

(4)

式中：R为微凸体曲率半径；δ为微凸体的变形量。

瓦伦丁[21]根据G-W模型得出实际接触面积正比于法向力，而比例系数几乎与法向力无关的结论。并且给出了接触表面平均接触压力的表达式式(5)，且该式已通过精确的数值计算验证，适用于各种粗糙表面，包括分形特征表面。

(5)

根据式(5)可将实际接触面上的压力表示为等效模量E*与表面轮廓平均梯度∇z乘积的一半

(6)

如果实际接触面上的压力大于材料屈服应力σ0，即

(7)

则微凸体处于完全塑性变形阶段，ψ被Greenwood等称为塑性指数，且当ψ<2/3时，接触表面表现为纯弹性。因此，由式(7)可知微凸体表现为弹性还是塑性与法向力无关。且表面轮廓梯度∇z与表面的检测分辨率有关，即尺度相关。如果在不同尺度，梯度不同，那么表面只在符合∇z>2σ0/E*条件的尺度发生塑性变形。

与传统的统计学参数不同，W-M函数的分形参数D、G均与表面微凸体的频率序数无关。Liu等给出了统计学表面均方差与分形参数之间的函数关系式为

(8)

式中：ωl为由采样长度决定的最低频率；ωh为由仪器分辨率和滤波决定的最高频率；S(ω)为W-M函数的功率谱。

因为统计学中假设表面微凸体服从正态分布，所以表面微凸体的均方差就等于均方根，即J=l。将式(8)代入式(7)即可得到与尺度相关的塑性指数表达式

(9)

由塑性指数的定义可知，微凸体发生何种变形与表面参数、材料参数以及微凸体的频率序数有关，而与法向载荷无关。那么根据与尺度相关的塑性指数表达式可得到弹性-弹塑性-完全塑性变形的临界频率序数m1、m2。当m

2 理论建模

本文将两个接触表面简化为一个粗糙表面和一个刚性平面，不考虑微凸体之间相互作用。在M-B模型中，假设粗糙表面与平面接触时，接触点面积分布规律与海洋面岛屿面积的分布规律相似，且遵循以下幂函数变化规律[22]

(10)

这里将两个粗糙表面的接触简化为一个刚性平面与一个粗糙平面的接触，忽略微凸体之间的相互作用。根据微凸体变形后的几何关系，可得到不同尺度下的微凸体横截面积为

a′=2πRδ=γ-2m

(11)

则根据式(10)、式(11)可以得到不同尺度下接触微凸体的分布函数为

(12)

式中：n(m)为不同频率下微凸体的分布函数;mL为接触微凸体的最小频率序数。

根据赫兹接触理论，可以得到弹性变形区间、弹塑性变形区间和塑性变形区间的接触载荷与频率序数的关系分别为[23]

(13)

(14)

fp=Ha′=Hγ-2m

(15)

式中：H为较软材料的硬度，可用材料的屈服强度表示为H=2.8Y；k为与泊松比有关的系数，k=0.454+0.41v；δ1为弹性变形到弹塑性变形的临界变形量，δ1=GD-1γ-m1(2-D)。

将单个微凸体的接触载荷在其变形的频率区间内积分可得整个结合面弹性接触载荷、弹塑性接触载荷和塑性接触载荷分别为

(16)

(17)

(18)

则结合面总的接触载荷为

F=Fe+Fep+FP

(19)

弹性、弹塑性和塑性变形微凸体总的真实接触面积分别为

(20)

(21)

(22)

则结合面总的真实接触面积为

A=Ae+Aep+Ap

(23)

单个弹性和弹塑性变形微凸体的法向接触刚度分别为

(24)

(25)

因为发生完全塑性变形的微凸体接触刚度为0，结合面上总的接触刚度为弹性变形微凸体接触刚度和弹塑性变形微凸体接触刚度的总和。

发生弹性变形的频率区间上微凸体总的接触刚度为

(26)

发生弹塑性变形频率区间上的微凸体总接触刚度为

[γ(0.38D-0.76)m2-γ(0.38D-0.76)m1]

(27)

总的法向接触刚度为

Kn=Kne+Knep

(28)

根据Wang等的研究可知：法向载荷作用下，接触表面处于塑性变形的微凸体由于受到局部接触载荷作用而发生塑性流动，将不能继续承受切向载荷。所以切向接触刚度只包括弹性变形微凸体和弹塑性变形微凸体。考虑滑动的切向接触问题，对于法向力和切向力同时作用时的组合接触问题，根据库伦摩擦定律得到的黏着条件

τ≤μp

(29)

由图1可知，黏着面边界的剪切应力接近无穷大，而赫兹接触的边界处法向应力为零，这意味着，接触边界不可能满足黏着的条件，即使在很小的切向载荷下接触边界总会有滑动。

血清IL-2、IFN-γ、TNF-α与IGF-1水平：2组患者于治疗前和治疗4周期后空腹采静脉血5 mL，离心取血清于-80 ℃储存，采用酶联免疫吸附法(ELISA)测量患者血清IL-2、IFN-γ、TNF-α与IGF-1水平。

图1 法向应力和切向应力

因此将总接触面分为一个内无滑动区(黏着)和一个外滑动区，如图2所示。用条件τ=μp可以得到与法向载荷和切向载荷的比值有关的黏着区域半径

图2 圆形切向接触中的黏着区和滑动区

(30)

式中：FT为切向载荷；c为黏着区域半径；r为微凸体截面半径。

又已知黏着区域的位移和分布载荷

(31)

(32)

根据式(31)、式(32)可得到黏着区域的接触刚度为

(33)

(34)

(35)

则可得单个弹性变形和弹塑性变形微凸体切向接触刚度分别为

(36)

(37)

(38)

结合面弹塑性变形微凸体的切向接触刚度

(39)

总的切向接触刚度

Kt=Kte+Ktep

(40)

3 有限元仿真及试验验证

本章设计了一个框形试件，用来验证多螺栓连接的装配体表面接触刚度的理论计算结果。因为用直接测量栓接结构接触刚度的试验方法误差较大，且难以表征接触刚度在接触表面分布的不均匀性，所以本次考虑到栓接表面接触应力分布的不均匀性，通过静力学分析提取接触表面节点的接触应力，再根据式(19)、式(28)、式(40)得到法向接触刚度及切向接触刚度与接触载荷的非线性关系，从而计算得到各节点的接触刚度，通过ANSYS有限元软件中的Matrix27单元建立节点接触对，并定义节点刚度，通过模态分析得到装配体的固有频率，并与振动测试试验的结果进行对比。具体实施步骤如图3所示。

图3 实施流程图

框型试件选用的材料是球墨铸铁QT600-3，表面是粗糙度Ra=1.6 μm的铣削加工平面，尺寸如图4所示。材料具体参数如表1所示。

图4 试验件尺寸示意图(mm)

表1 材料参数表

通过三维形貌仪获得表面轮廓参数，如图5、图6所示。并运用功率谱密度法得到表面分形参数D=1.42，G=1.21×10-8m，由此可得到弹性到弹塑性变形的临界频率序数m1=18，弹塑性到塑性变形的临界频率序数m2=22。

图5 试件粗糙表面微观形貌图

图6 功率谱密度

考虑到栓接结合面接触应力分布不均匀的情况，如图7所示。用通过静力学分析得到的接触单元的应力分布来表征接触刚度分布的不均匀性，将提取的接触表面应力导入MATLAB程序来计算各节点刚度。通过有限元软件中的Matrix27单元建立各节点之间的接触，最终对装配体进行模态分析求得装配体的固有频率。

图7 试件接触应力分布示意图

试件用M16的螺栓装配起来，螺栓材料为45钢。试验施加的螺栓预紧力为18 000 N。为了准确取得栓接试件结合部的动态特性，用绳子将框型试件吊起，并在栓接装配体上共布置了10个加速度传感器，如图8所示。利用力锤连续敲击10次以提供激振力，利用LMS装置采集并处理振动数据，得到试验件的各阶固有频率及振型,如表2所示。

图8 试验装置示意图

表2 仿真与试验所得固有频率结果对比表

通过理论结果与试验结果的对比，误差在8.7%以内，且与经典的M-B模型结果进行了对比，为结合面的刚度建模提供了更加精确的方法。

4 讨论

根据推导出的塑性指数表达式可以得到不同尺度下的塑性指数，塑性随频率序数的变化规律如图9、图10所示。

由图9可知，随着频率序数的增大，塑性指数呈指数增大。这是因为高频的微凸体更易发生塑性变形，而低频的微凸体有着较大的曲率半径，其真实接触面积较大，更容易发生弹性变形。所以随着频率序数的增加微凸体依次发生弹性、弹塑性、塑性变形。根据接触微凸体分布函数可知，在小载荷情况下，高频微凸体最先进入接触，然后随着载荷的增大，就会有低频的微凸体发生弹塑性变形以及弹性变形。且同一频率序数下塑性指数随着材料屈服强度的增大而减小，这说明屈服强度较低的材料表面上的微凸体更容易发生塑性变形。由图10可知，当材料参数确定时，粗糙度越大的加工表面，塑性指数越大，更多的微凸体发生塑性变形，整个结合面上的接触刚度越小。

图9 塑性指数与频率序数函数关系图

图10 塑性指数与表面粗糙度的关系图

接触载荷对真实接触面积的影响,如图11所示。M-B模型与本文真实接触面积结果对比,如图12所示。

图11 载荷与真实接触面积的关系图

由图11可知,结合面的真实接触面积随接触载荷的增大而增大，真实接触面积只占名义接触面积的千分之一左右。这是因为随着载荷的增大，更多的微凸体发生接触和形变，真实接触面积随之增大。由W-M函数对表面轮廓的描述可知，随着D的增大和G的减小，表面会趋于平滑，导致真实接触面积变大。由图12可知，塑性变形的微凸体在相同载荷下，真实接触面积会更大，这是因为弹塑性变形的微凸体相比于完全塑性变形的微凸体曲率半径越更大，曲率半径较大的微凸体接触面积更大。由于真实接触面积占名义接触面积的比值很小，所以栓接结构因螺栓松动而产生的二次以及多次装配过程中，微凸体的接触是一个随机的过程，并不会对表面的接触刚度有显著影响。

图12 M-B模型与本文真实接触面积结果对比

接触载荷对接触刚度的影响,如图13、图14所示。

图13 法向接触刚度与载荷的关系图

图14 切向接触刚度与载荷的关系图

由图13、图14可知，法向和切向接触刚度均随着载荷F的增大而增大。这是因为随着载荷F的增大，更多的微凸体参与接触，真实接触面积变大，所以接触刚度随之增大。由图13、图14还可知:随着D的增大，G减小，接触刚度随之增大，这是因为趋于平滑的表面会使真实接触面积变大，接触刚度随之变大；并且随着D增大，G减小，微凸体的曲率半径变大，更多的微凸体会发生弹性以及弹塑性变形，因此接触刚度也会越大。现有研究给出了实际工程当中几种常见的加工方式和表面分形参数，如表3所示[25]。合理选择装配表面的加工方式，可以有效提高栓接装配体的接触特性。