机械结合面切向接触刚度的三维分形理论建模
潘五九 李小彭 李木岩 王雪 高建卓 闻邦椿
摘要： 将微凸体的弹塑性变形区进一步划分，并考虑三维结合面形貌的WM函数，推导了三维机械结合面切向分形接触刚度的理论模型。数值模拟了结合面的三维切向接触刚度随着分形维数D、分形尺度系数G、材料特征参数间的变化趋势，以及二维分形和三维分形间的对比分析。仿真结果显示：机械结合面的三维切向接触刚度与法向载荷和材料特征参数成单调递增关系，与分形尺度系数成单调递减关系；而其与分形维数之间以D=2.5为界，依次成递增与递减关系；三维分形下的结合面切向接触刚度大于二维分形下的结合面切向接触刚度。切向接触刚度模型的构建可为后续粗糙表面接触非线性动力学及整机动力学模型的建立提供基础。关键词： 结合面； 三维分形； 弹塑性区再划分； 切向刚度； 分形理论
中图分类号： TH113.1； TB123文献标志码： A文章编号： 10044523（2017）04057710
DOI：10.16385/j.cnki.issn.10044523.2017.04.008
引言
机械结合面大量存在于机械系统中，如高精密机床导轨间的结合面、刹车制动盘间的结合面、精密减速器中齿轮间的结合面、微纳器件间的结合面等。对机械结合面研究的核心本质问题，就是对两个粗糙表面间的相互接触进行分析计算。经过机械加工的零部件表面看似光滑，实际上从微观来看表面存在着类似于山峰高低起伏般的粗糙度。当机械进行作业任务，传递运动载荷或其受到外部激励时，结合面间微凸体形貌的变化将影响系统的振动、摩擦、磨耗、润滑、热阻、电阻及接触刚度和接触阻尼等。而这些表征量又直接影响一个系统的静、动态特性，进而宏观表现为影响系统的加工精度，运行精度以及振动噪声等现象。
因此，有必要基于Hertz理论[1]和分形理论[2]对微凸体接触进行分析，进而推导出机械结合面的接触刚度模型，这是研究系统静、动态特性的核心问题之一，也是学术界和企业界一直研究的热点。为此很多国内外学者对粗糙结合面接触特性进行了深入的探讨。Bemporad等[3]给出了一种可用于自仿射分形表面法向接触下的优化算法，该法可有效应用于工程粗糙表面间接触分析。Pohrt等[4]利用边界有限元法来分析法向和切向接触问题。Jourani[5]在进行加载接触分析时考虑微凸体间相互作用影响，并分析了三维分形维数对接触面积的影响。温淑花和张学良等[67]基于MB模型[8]建立了考虑域扩展因子的结合面切向刚度模型，并分析了切向刚度和很多结合面间参数的影响关系；田红亮等[910]在改进分形理论和严格应用Hertz接触力学基础上建立切向接触刚度模型，并进行了法、切向刚度的相关验证。李小彭等[11]建立了考虑摩擦因素的结合面法向刚度分形模型；Jiang等[12]通过实验给出了不同加工表面下的刚度值，并将它们与理论值进行对比；You [13]给出了结合面上法向和切向的统计模型。Zhao等[14]建立了结合面上微凸体从弹性到完全塑性变形的微接触模型。
经过对以上文献中理论模型分析后发现，以上理论建模均是在二维接触曲线下，即在分形维数为1 1机械结合面三维分形接触理论基础 根据赫兹理论知，两微凸体相互接触可等效为一刚性光滑平面和一等效微凸体相接触。简化之后的模型如下图1所示。 图1光滑刚性平面与等效微凸体接触 Fig.1Contact between smooth rigid plane and equivalent asperity 图1中直角三角形obc由勾股定理可得R2=r′2+（R-ω）2（1）式（1）变形可得R=ω2+r′22ω（2）由于变形量ω远小于R，故假设Rω2，可得近似式为R≈r′22ω（3）第4期潘五九，等： 机械结合面切向接触刚度的三维分形理论建模振 动 工 程 学 报第30卷图1中没有变形的等效微凸体和刚性光滑平面接触的横截面积是A=πr′2=π[R2-（R-ω）2]≈2πRω（4）微凸体赫兹接触时，实际接触面积的半径[13]r为：r=3peR4E13（5） E=1-ν21E1+1-ν22E2-1（6）式中E表示等效弹性模量；E1和E2表示两相互接触的微凸体的弹性模量；ν1和ν2则表示它们的泊松比。 式（5）中pe为单个微凸体处于弹性变形阶段时所承受的法向载荷[15]pe=43ER0.5ω1.5（7）将式（7）代入式（5）可得r=R0.5ω0.5（8）將式（3）代入式（8）可得r=22r′（9）根据式（4）和（9），图1中等效微凸体和刚性光滑平面实际接触面积为a=πr2=12πr′2=πRω（10）Yan等[16]改进WM函数，得到了能更加准确模拟三维表面形貌的曲面函数式。修改后的WM函数关系式为 z（x，y）=LGL（D-2）lnγM1/2∑Mm=1∑nmaxn=0γ（D-3）n× cosφm，n-cos2πγn（x2+y2）1/2L× cos（tan-1yx-πmM+φm，n（11） 式中L为表面形貌取样长度；D为三维表面形貌分形维数（2
由式（13）～（15）能推导出微凸体的临界弹性变形面积
ac=23D-112-D33kμ〖〗4022-Dπ4-D2-D（lnγ）1D-2G2（16）
由文献[1718]可知，对于微凸体的弹塑性区仍可再划分为两个区域，即弹塑性第一区域和弹塑性第二区域，且给出这两个分区的各自临界面积为aep1=11012-Dac（17）
aep2=612-Dac（18）2.2建立微凸体在各变形阶段的法向接触载荷
当单个微凸体受载处于弹性变形阶段时，其变形量与受载的关系见式（7）。根据式（7），（13）与（14）可得表示成载荷是接触面积的函数p5（a）=13EπD-42215-3D2（lnγ）12GD-2a4-D2（19）单个微凸体塑性变形阶段法向载荷[11]pp=λσya（20）式中λ=H/σy为定义的系数，H为较软材料的硬度。
根据文献[17]，弹塑性接触第一、二区的法向接触载荷和表面微凸体变形量之间存在关系有：
当1≤ω/ωc≤6时，法向接触载荷p与微凸体变形量间关系为ppc=1.03（ω〖〗ωc）1.425（21）当6≤ω/ωc≤110时，法向接触载荷p与微凸体变形量间关系为ppc=1.40（ωωc）1.263（22）式（21）和（22）中，pc为ω=ωc时的法向接触载荷，pc=43ER0.5ω1.5c。
联合式（13）～（15）和式（21）～（22），可推导出微凸体弹塑性第一，第二区域的法向接触载荷与接触面积的关系：
当1≤ω/ωc≤6时，有
pep1=1.033E（33kμ40）0.15（lnγ）0.425×
π0.425D-1.726.75-1.275DG0.85D-1.7a1.85-0.425D（23）
当6≤ω/ωc≤110时，有
pep2=1.403E（33kμ40）0.474（lnγ）0.263×
π0.263D-1.05224.893-0.789DG0.526D-1.052a1.526-0.263D（24）
2.3微凸體的切向载荷与切向接触刚度
据文献[19]的研究，单个微凸体所承受的切向载荷为Q=8aσy6-3νπ+82ν-1π6-3νp（25）式中σy为材料屈服强度；ν为泊松比。
微凸体在结合面上的面积分布函数[16]n（a）和最大微凸体接触面积al间的关系为na=D-12aD-12la-D+12，
0
（1）D≠2.5时
Te=∫alac8aσy6-3νπ+82ν-1π6-3νpen（a）da=
8σyD-1π6-3ν3-Da0.5D-0.5l（a1.5-0.5Dl-
a1.5-0.5Dc）+2ν-1E3π6-3ν5-2D·
π0.5D-2lnγ0.5GD-2×211.5-1.5Da0.5D-0.5l·
a2.5-Dl-a2.5-Dc（27）
（2）D=2.5时
Te=∫al ac 8aσy6-3νπ+82ν-1π6-3νpe nada=
24σy π6-3νa0.75l（a0.25l-a0.25c）+
2ν-1Eπ6-3νπ-0.75（lnγ）0.5G0.524.75a0.75llnalac（28）
对于所有处于弹塑性阶段的微凸体切向载荷为：
当1≤ω/ωc≤6时，根据式（23），（25）和（26）有
Tep1=∫aep2aep18aσy6-3νπ+82ν-1π6-3νpep1nada=
8σyD-1π6-3ν3-Da0.5D-0.5l（a1.5-0.5Dep2-
a1.5-0.5Dep1）+2.06E2ν-136-3ν2.35-0.925D·
33kμ400.15×（lnγ）0.425π0.425D-2.728.75-1.275D·
G0.85D-1.7a0.5D-0.5l×（a2.35-0.925Dep2-
a2.35-0.925Dep1）（29）
当6≤ω/ωc≤110时，根据式（24） ～（26）有
Tep2=∫acaep28aσy6-3νπ+82ν-1π6-3νpep2nada=
8σyD-1〖〗π6-3ν3-Da0.5D-0.5l（a1.5-0.5Dc-
a1.5-0.5Dep2）+1.4E2ν-136-3ν2.026-0.763D·
33kμ400.474×（lnγ）0.263π0.263D-2.05227.893-0.789D·
G0.526D-1.052×a0.5D-0.5l（a2.026-0.763Dc-
a2.026-0.763Dep2） （30）
单个微凸体的切向接触刚度kt表示为[6]kt=8a0.52-νπ1-1μTP13（31）式中为结合面材料的等效剪切弹性模量。
3建立机械结合面三维分形下的接触状态〖*2〗3.1机械结合面上的法向载荷对面积分布函数在各个变形区间进行连续积分，可计算出整个结合面上的实际接触面积Ar，它应包含微凸体经历从弹性变形，弹塑性第一、二区域变形和完全塑性变形各个阶段的接触面积总和。
Ar=∫aep10n（a）ada+∫acaep2n（a）ada+
∫alacn（a）ada=D-13-Dal（32）
（1）D≠2.5时，总法向载荷由式（17） ～（20），（23），（24）和（26）得
P=∫aep10ppn（a）da+∫aep2aep1pep2n（a）da+
∫acaep2pep1n（a）da+∫alacpen（a）da=
D-13-Dλσy1103-D4-2Da0.5D-0.5la1.5-0.5Dc+
1.403E33kμ400.474D-12.026-0.763D·
23.893-0.789Dπ0.263D-1.052（lnγ）0.263·
G0.526D-1.052a0.5D-0.5la2.026-0.763Dc[62.026-0.763D2-D-
1102.026-0.763D2-D]+1.03〖〗3E33kμ400.15·
D-1〖〗2.35-0.925D25.75-1.275Dπ0.425D-1.7·
（lnγ）0.425G0.85D-1.7a0.5D-0.5la2.35-0.925Dc[1-
62.35-0.925D2-D]+（D-1）3（2.5-D）·
（26.5-1.5DEπ0.5D-2lnγ）0.5GD-2a0.5D-0.5l·
（a2.5-Dl-a2.5-Dc）（33）
对式（33）进行无量纲化
P* = λh1 1103-D4-2Da（0.5D-0.5）*l a（1.5-0.5D）*c +
1.403（33kμ 40）0.474h2 23.393-0.289Dπ0.263D-1.052·
（lnγ）0.263G（0.526D-1.052）*×a（0.5D-0.5）*l ·
a（2.026-0.763D）*c [62.026-0.763D2-D-
1102.026-0.763D2-D] +1.033（33kμ 40）0.15h3 ·
25.25-0.775Dπ0.425D-1.7（lnγ）0.425G（0.85D-1.7）*·
a（0.5D-0.5）*l a（2.35-0.925D）*c [1-62.35-0.925D2-D] +
13π0.5D-226.5-1.5Dh4 （lnγ）0.5G（D-2）*·
a（0.5D-0.5）*l [a（2.5-D）*l -a（2.5-D）*c ] （34）
式中P=P/（EAa），al=al/Aa，ac=ac/Aa，G=G/Aa，h1=D-13-D，h2=D-12.026-0.763D，h3=D-1〖〗2.35-0.925D，h4=D-12.5-D。
（2）D=2.5時，总法向载荷由式（17） ～（20），（23），（24）和（26）得
P=∫aep10ppn（a）da+∫aep2aep1pep2n（a）da+
∫acaep2pep1n（a）da+∫alacpen（a）da=
3λσy110-0.5a0.75la0.25c+0.7E33kμ400.474·
21.9205π-0.3945（lnγ）0.263×G0.263a0.75la0.1185c（6-0.237-
110-0.237）+0.515E33kμ400.1522.5625π-0.6375·
（lnγ）0.425G0.425a0.75la0.0375c（1-6-0.075）+0.125E·
π-0.7524.75（lnγ）0.5G0.5a0.75l（lnal-lnac） （35）
对式（35）进行无量纲化
P* = 3λ110-0.5a0.75*l a0.25*c + 0.733kμ 400.474·
21.9205π-0.3945（lnγ）0.263G0.263*a0.75*l a0.1185*c ·
（6-0.237-110-0.237） + 0.51533kμ 400.15·
22.5625π-0.6375（lnγ）0.425G0.425*a0.75*l a0.0375*c ·
（1-6-0.075） + 0.125π-0.7524.75（lnγ）0.5·
G0.5*a0.75*l lna*la*c（36）
3.2结合面的三维分形切向接触刚度
综合式（17）～（19），（23），（24）和（26）～（30）可以得到三维分形切向刚度模型
Kt=∫aep2aep1ktep2n（a）da+∫acaep2ktep1n（a）da+
∫alackten（a）da=8G2-νπ1-1〖〗μTep2Pep213·
D-1〖〗2-D（60.5-1100.5）a0.5D-0.5la1-0.5Dc+1-1μ·
Tep1Pep113D-12-D（1-60.5）a0.5D-0.5la1-0.5Dc+1-1μTePe13D-12-Da0.5D-0.5l（a1-0.5Dl-a1-0.5Dc）（37）
对式（37）进行无量纲化
Kt=∫aep2aep1ktep2n（a）da+∫acaep2ktep1n（a）da+
∫alackten（a）da=82-νπ1-1μTep2Pep213·
D-12-D（60.5-1100.5）a（0.5D-0.5）*la（1-0.5D）*c+
1-1μTep1Pep113D-12-D（1-60.5）·
a（0.5D-0.5）*la（1-0.5D）*c+1-1μTePe13·
D-12-Da（0.5D-0.5）*l（a（1-0.5D）*l-a（1-0.5D）*c）（38）
式中Kt=Kt/Aa，al=al/Aa，ac=ac/Aa。
4数值仿真4.1切向刚度和分形维数D之间的关系根据文中式（32），（34），（36）和（38），给定无量纲参数G=10-10 ，=2.5 ，ν=0.3，μ=0.19。研究结合面切向接触刚度在不同的D下的变化关系，如图2所示。当2.1≤D≤2.5时，结合面切向刚度随分形维数D成递增关系；当2.5≤D≤2.9时，结合面切向刚度随分形维数D成递减关系。这是因为当D大于等于2.5时，结合面上微凸体变小，且单一微凸体的完全塑性变形比例提高，导致单一微凸体刚度减小，宏观表现为三维结合面形貌的切向接触刚度减小。且知，D对切向刚度的影响很大，表现为图2 （a），（b）和（c）的纵坐标数值出现数量级上的差异。
4.2切向刚度与分形尺度系数G间的关系
根据文中式（32），（34），（36）和（38），分别给定无量綱参数G=10-9，10-10，10-11在D为2.2，2.5和2.8时的变化关系，如图3所示。易见，结合面切向刚度随无量纲分形尺度系数G增大而相应的减小。原因为G越大，结合面的宏观表面形貌越不光滑，整个表面上的微凸体的弹性变形能下降，宏观表现为结合面的切向接触刚度减小。
图2无量纲切向刚度在不同分形维数下的变化关系
Fig.2The change relation of dimensionless tangential stiffness under different fractal dimension4.3切向刚度和材料特性参数间的关系
根据式（32），（34），（36）和（38），参数选取同4.1节，来研究结合面切向接触刚度与材料特性参数=σy/E间的变化关系，如图4所示。随材料特性参数增大，结合面的切向刚度相应增大。且在确定的下，随着分形维数的增大曲线由非线性逐渐趋于线性化。
图3无量纲切向刚度在不同分形尺度系数下的变化关系
Fig.3The change relation of dimensionless tangential stiffness under different topothesy图4无量纲切向刚度在不同材料特性参数下的变化关系
Fig.4The change relation of dimensionless tangential stiffness under different material parameters
图5二维分形与三维分形切向刚度的关系
Fig.5The relationship between two dimensional fractal and three dimensional fractal tangential stiffness4.4对比结合面二维分形切向刚度和三维分形切向刚度三维形貌分形与二维曲线分形间的关系为D=Ds+1[18]，其中Ds表示二维分形维数。如图5所示。由三维分形模型计算得出的切向刚度要比二维分形模型[9]计算得出的切向刚度要大，两者间的差值随结合面上法向载荷的增大而增大；随着分形维数的增大，两种模型计算得出的切向刚度差值在缩小。因此，在这两种不同模型计算下的结合面切向刚度是有所差异的。
5试验验证
本节将对两块由16个M6螺栓连接而成的45号钢板进行试验模态分析，并将其结果与有限元结果进行对比，来验证文中推导的结合面间切向刚度模型的合理性。两块钢板尺寸均为400 mm×50 mm×6 mm。据文献[20]，可对机械结合面进行“固隙固”等效层处理。如图6（b）所示。图中h为等效层厚度，此处取1 mm[21]。由下式[22]，可分别求得等效层的弹性模量E，剪切模量G和泊松比ν。E=Knh/（Aa）（39）
G=Kth/（Aa）（40）
ν=E/2G-1 （41）图6试验用钢板及等效处理示意图
Fig.6Steel plate for test and schematic diagram of equivalent treatment
结合面中两个重要分形参数D和G的获取是通过T1000型轮廓仪，通过测算分析得到。经计算得D=2.427，G=1.342×10-5 m。假定每个螺栓的拧紧力矩是3 N·m，则每个的预紧力是2500 N，结合面上总法向载荷是40009 N。根据表1中钢板材料常数，再通过文中式（37）求得切向刚度Kt=5.5002×1010 N/m，由文献[23]求得法向刚度Kn=1.328×1011 N/m。Aa为名义接触面积，等于2×10-2 m2。将这些数值代入公式（39），（40）和（41），可求得等效层的E=6.64×109 Pa， ν=0.2072。等效层的参数将用于有限元模态分析。为了验证有限元模态分析的结果，下面进行模态试验分析。试验设备及过程分别如表2和图7所示。
实际测得的幅频图和两种方法下得到的结果对比分别如图8和表3所示。可知，将结合面作等效处理得到的结果和试验结果较为接近，误差均在20%以内。误差原因主要为利用等效层法做有限元模态分析时，所建的物理模型是对实际的螺栓连接试验钢板做了简化处理的，且不考虑其转动惯量和剪切变形的影响。进行试验自由模态分析时，图7中用悬索悬挂也是自由状态的一种近似。 因此， 上
表1钢板材料常数和表面分形参数
Tab.1Material constants and fractal parameters of steel plate
E/Paνρ/（kg·m-3）σy/PaDG/m2.1×10110.378503.55×1082.4271.342×10-5
表2试验设备
Tab.2The experiment instrument
序号设备名称1 DH5956数据采集分析系统24508B加速度传感器3DHL050型脉冲响应力锤
图7现场测试图
Fig.7Experiment process图8幅频曲线
Fig.8Curves of amplitude frequency
表3有限元结果与试验结果对比
Tab.3Comparison of finite element results and test results
模态
阶数模态试验
结果/Hz三维分形
有限元
结果/Hz二维分形
有限元
结果/Hz三维分形
误差二维分形
误差1365395.5412.6+8.36%+13.04%28971070.71071.6+19.36%+19.46%31690.61513.91872.8-10.45%+10.78%424302045.62891.0-15.82%+18.97%530303090.63910.1+2.00%+29.05%
述误差在可接受范围之内，同时可验证文中所建切向刚度模型的合理性。为了进一步分析基于二维分形曲线得到的二维切向接触刚度[24]相较于基于三维分形形貌得到的三维切向接触刚度，与模态试验结果的差值大小，仍然可以将二维的切向接触刚度按照本节的方法进行计算。从表3中，可看出二维分形下的相对误差要更大些。
6结论
（1）三维分形模型下的结合面切向接触刚度与法向载荷成单调递增关系。分形维数D对三维结合面切向接触刚度的影响较复杂，当2.1≤D≤2.5时，结合面切向刚度随D成递增关系；当2.5≤D≤2.9时，结合面切向刚度随D成递减关系。
（2）三维结合面切向接触刚度与分形尺度系数G成单调递减关系；与材料特性参数成单调递增关系。
（3）研究对比了三维分形和二维分形下的切向接触刚度；对比发现三维分形下的结合面切向接触刚度大于二维分形下的结合面切向接触刚度。
参考文献：
[1]Johnson K L. Contact Mechanics[M]. Cambridge： Cambridge University Press， 1985.
[2]Mandelbort B B. The Fractal Geometry of Nature[M]. New York： W H Freeman，1982.
[3]Bemporad A， Paggi M. Optimization algorithms for the solution of the frictionless normal contact between rough surfaces[J]. International Journal of Solids & Structures， 2015，S6970（3）：94—105.
[4]Pohrt R， Li Q. Complete boundary element formulation for normal and tangential contact problems[J]. Physical Mesomechanics， 2014，17（4）：334—340.
[5]Jourani A. Effect of 3D fractal dimension on contact area and asperity interactions in elastoplastic contact[J]. Aip Advances， 2016，6（5）：5799—76.
[6]温淑花，张学良，文晓光，等.结合面切向接触刚度分形模型建立与仿真[J].农业机械学报，2009，40（12）：223—227.
Wen Shuhua， Zhang Xueliang， Wen Xiaoguang， et al. Fractal model of tangential contact stiffness of joint interfaces and its simulation[J]. Transactions of the Chinese Society for Agricultural Machinery， 2009，40（12）：223—227.
[7]温淑花，张学良，武美先，等.结合面法向接触刚度分形模型建立与仿真[J].農业机械学报，2009，40（11）：197—202.
Wen Shuhua， Zhang Xueliang， Wu Meixian， et al. Fractal model and simulation of normal contact stiffness of joint interfaces and its simulation[J]. Transactions of the Chinese Society for Agricultural Machinery， 2009，40（11）：197—202.
[8]Majumdar A， Bhushan B. Fractal model of elasticplastic contact between rough surfaces[J]. ASME Journal of Tribology， 1991，113（1）：1—11.
[9]田红亮， 陈从平，方子帆，等.应用改进分形几何理论的结合部切向刚度模型[J].西安交通大学学报，2014，48（7）：46—52.
Tian Hongliang， Chen Congping， Fang Zifan， et al. Tangential stiffness model for joint interface adopting the revised fractal geometric theory[J]. Journal of Xi′an Jiaotong University，2014，48（7）：46—52.
[10]田红亮，赵春华，朱大林，等.金属材料结合部法切向刚度修正与实验验证[J].农业机械学报，2012，43（6）：207—214.
Tian Hongliang， Zhao Chunhua， Zhu Dalin， et al. Modification of normal and tangential stiffness for joint interface with metallic material and experimental validation[J]. Transactions of the Chinese Society for Agricultural Machinery， 2012，43（6）：207—214.
[11]李小彭，郭浩，刘井年，等.考虑摩擦的结合面法向刚度分形模型及仿真[J].振动、测试与诊断，2013，33（2）：210—213.
Li Xiaopeng， Guo Hao， Liu Jingnian， et al. Normal stiffness fractal model of joint surface considering friction factors and its simulation[J]. Journal of Vibration，Measurement&Diagnosis;， 2013，33（2）：210—213.
[12]Jiang S， Zheng Y， Zhu H. A contact stiffness model of machined plane joint based on fractal theory[J]. Journal of Tribology， 2010，132（1）：01140110114017.
[13]You J M， Chen T N. Statistical model for normal and tangential contact parameters of rough surfaces[J]. ARCHIVE Proceedings of the Institution of Mechanical Engineers Part C： Journal of Mechanical Engineering Science， 2010，1（1）：1—15.
[14]Zhao Y， Maietta D M， Chang L. An asperity microcontact model incorporating the transition from elastic deformation to fully plastic flow[J]. Journal of Tribology， 2000，122（1）：86—93.
[15]Chang W R， Etsion I， Bogy D B. Static friction coefficient model for metallic rough surfaces[J]. ASME Journal of Tribology， 1988，110（1）：57—63.
[16]Yan W， Komvopoulos K. Contact analysis of elasticplastic fractal surfaces[J]. Journal of Applied Physics， 1998，84（7）：3617—3624.
[17]Luen L J， Fin L J. A new microcontact model developed for variable fractal dimension， topothesy， density of asperity， and probability density function of asperity heights[J]. Journal of Applied Mechanics， 2007，74（4）：603—613.
[18]Liou J L， Lin J F. A modified fractal microcontact model developed for asperity heights with variable morphology parameters[J]. Wear， 2010， 268（12）：133—144.
[19]盛選禹，雒建斌， 温诗铸.基于分形接触的静摩擦系数预测[J].中国机械工程，1998，9（7）：16—18.
Sheng Xuanyu， Luo Jianbin， Wen Shizhu. Prediction of static friction coefficient based on fractal contact[J]. China Mechanical Engineering， 1998，9（7）：16—18.
[20]李小彭，梁亚敏，郭浩，等.结合面广义间隙的等效模型研究[J].振动工程学报，2014，27（1）：25—32.
Li Xiaopeng， Liang Yaming， GUO Hao， et al. Study on equivalent model of generalized clearance of joint surface[J]. Journal of Vibration Engineering， 2014，27（1）：25—32.
[21]田红亮，刘芙蓉，方子帆，等.引入各向同性虚拟材料的固定结合部模型[J].振动工程学报，2013，26（4）：561—573.
Tian Hongliang， Liu Furong， Fang Zifang， et al. Immovable joint surface′s model using isotropic virtual material[J]. Journal of Vibration Engineering， 2013，26（4）：561—573.
[22]Li X， Liang Y， Zhao G， et al. Dynamic characteristics of joint surface considering friction and vibration factors based on fractal theory[J]. Journal of Vibroengineering， 2013， 15（2）：872—883.
[23]李小彭，王雪，运海萌，等.三维分形固定结合面法向接触刚度的研究[J].华南理工大学学报（自然科学版），2016，44（1）：114—122.
Li Xiaopeng， Wang Xue， Yun Haimeng， et al. Investigation into normal contact stiffness of fixed joint surface with threedimensional fractal[J]. Journal of South China University of Technology （Natural Science Edition）， 2016，44（1）：114—122.
[24]赵光辉.固定结合面接触刚度建模与广义间隙等效方法研究[D].沈阳：东北大学，2014.
Zhao Guanghui. Study on contact stiffness modeling and equivalent method of generalized clearance of fixed joint surface[D]. Shenyang： Northeastern University， 2014.
Threedimensional fractal theory modeling of tangential contact
stiffness of mechanized joint surfaces
PAN Wujiu， LI Xiaopeng， LI Muyan， WANG Xue， GAO Jianzhuo， WEN Bangchun
（School of Mechanical Engineering & Automation， Northeastern University， Shenyang 110819， China）
Abstract： Establishing the stiffness model of joint surfaces is one of the key problems to study the static and the dynamic characteristics of the mechanical systems. The elasticplastic deformation zone of asperity is further divided and the modified WM function of threedimensional surface topography is considered， then the threedimensional fractal tangential contact stiffness model of mechanical joint surfaces is established. Through numerical simulation the relationship between the tangential contact stiffness and the fractal dimension， the topothesy and the material characteristic parameters are studied. Comparative analysis between two dimensional fractal and three dimensional fractal is carried out. Results show： the tangential contact stiffness increases with the normal load and the characteristic parameters of material increase， and decreases with the topothesy increases； the relationship between fractal dimension and tangential contact stiffness is more complicated； the threedimensional fractal tangential stiffness is greater than the twodimensional case. The construction of the tangential contact stiffness model can provide the foundation for the establishment of the nonlinear dynamics of rough contact surface and the dynamic model of whole machine in future.Key words： joint surfaces； threedimensional fractal； elasticplastic zone subdivide； tangential contact stiffness； fractal theory作者简介： 潘五九（1986—），男，博士研究生。电话： 15940263160； Email： panspace@sina.cn
通讯作者： 李小彭（1976—），男，教授，博士研究生导师。电话： 13940029225； Email： xpli@me.neu.edu.cn
聞邦椿（1930—），男，中国科学院院士，教授，博士研究生导师。 Email： bcwen1930@vip.sina.com
