机床结构中接合面的摩擦特性分析

Frictional Behavior Analysis of Contact Interfaces in Machine Tool Structures

摘要: 本文针对机床制造有重大作用的结合面，提出了一种基于有限元模拟的摩擦分析模型，对结合面的摩擦特性进行了研究。采取一定的假设条件对表面轮廓完成了数值简化模拟，并将其导入ABAQUS软件，进行建模，基于此模型分析其结合面性质。结果表明，纯弹性结合面的切向力与切向位移成正比，并且其比值随着摩擦系数的增加而增加。塑性材料的切向力的大小同样随切向位移增加而增加，但两者之间不再互成正比，切向力–切向位移曲线的斜率随着切向位移的增加而逐渐减小，最终趋于零。切向力随着屈服极限以及摩擦系数增加而增加。当模型塑性变形较大时，其切向力的大小将与屈服极限更加相关。

Abstract: This study proposes a finite element-based friction analysis model to investigate the tribological behavior of contact interfaces, which play a critical role in machine tool manufacturing. By applying a series of well-defined assumptions, the surface topography was numerically simplified and subsequently imported into ABAQUS for model construction and analysis. The results demonstrate that, for purely elastic contact interfaces, the tangential force is proportional to the tangential displacement, with the proportionality constant increasing as the coefficient of friction rises. In contrast, for plastic materials, the tangential force also increases with tangential displacement; however, the relationship is no longer linear. The slope of the tangential force-displacement curve gradually decreases with increasing displacement and ultimately approaches zero. Moreover, the tangential force increases with both the yield strength and the coefficient of friction. In scenarios involving significant plastic deformation, the magnitude of the tangential force becomes more strongly correlated with the yield strength of the material.

1. 引言

在工业生产过程中，“制造”是其中的关键所在[1]，除此之外，摩擦与机器的正常运行也同样息息相关。而在对高速高精度机床的开发中，摩擦分析无疑是其中重要一环，摩擦在低速时表现出强烈的非线性特性，机械结构分别在其启动、停止以及速度的转向时受到摩擦的影响最严重。为了进一步提高机械系统的效率，我们需要想办法消除或者减少摩擦带来的影响，针对这个问题，人们在各个方面都做出了一定的尝试。例如提高零件的加工精度、改进润滑条件以及采用高性能轴承(如气浮轴承、磁悬浮轴承等)等等都能够有效地减少摩擦的影响。

然而上述的解决方法中，机械的价格一般都极为昂贵，导致性价比极低，甚至有些精度很难通过现有的技术实现。因此对于摩擦力在现有的机械结构之中进行补偿和修正来消除其不利的影响是十分有必要的。摩擦非线性几乎无处不在[2]。很多日常的案例之中，摩擦扮演着对人们有利的角色，行走、自行车的运动、各种加持、切割和洗刷等都受到了摩擦力的帮助；除此之外，在我们的机械结构中，同样有不少依托摩擦力的原理工作，如摩擦压力机、摩擦离合器、摩擦传动机构、摩擦制动装置和螺栓连接等。然而在大多数的制造情况之下，摩擦是有害的。其在高精度伺服机构位置控制和低速控制中存在感很强，是提高性能最主要的障碍之一。摩擦不仅仅会使系统在位置伺服方式时产生死区或极限环，从而引起稳态误差，而且会使系统进入无序运动的混沌状态，甚至会使系统在速度伺服时产生爬行、振荡现象。这一点在精密加工、超精密加工、纳米加工等技术体现的尤为明显，这些加工方式都要求伺服系统具有极高的定位精度和低速性能。因此合理地解决机械系统中存在的摩擦环节尤其是非线性摩擦环节的制约问题成为当前研究的热点[3]。

2. 建模及加载

2.1. 简化平面模型

在本文的研究中，鉴了对于机床结合面的摩擦特性研究的仿真模型的结合面简化方式，其根据实验的结果对结合面进行了一定量的简化假设：(1) 假设式样沿垂直方向带有一定的周期性，将其简化为周期分布，从而选取其中一个单元进行简化分析；(2) 假设任取式样沿垂直进给方向的一组截面，其微观形貌相同，即将结合面简化为平面应变的问题来进行建模。继而进一步对结合面的形貌进行分析。根据实验结果可以得到，接合面的形貌峰谷呈现类正态分布的情况，更据这个结论，可将结合面简化为一个半圆弧状的周期性单元，圆弧的半径R即为平均粗糙度$R\_a$ 。并在一定量的实验分析后进一步近似为正弦曲线分布。而本文所研究的简化结合面即参考了此种周期性的正弦曲线分布形[4]。

在简化模型参数设定中，结合面粗糙度采用正弦分布，如图1所示，其模型来源于前期实验中通过激光扫描法测量得到的表面轮廓数据。周期长度、波峰高度等参数则根据该平均粗糙度进行数学拟合，确保简化模型能反映实际接触形貌的统计特征。摩擦系数的范围选择基于已有文献对钢材或常见机床结合材料的干摩擦测试值，同时考虑了不同加工润滑状态下的可能取值区间。

Figure 1. Interface surface topography scan

图1. 结合面扫描图

具体建模步骤：首先在外部制作python设定所需的简化平面的点集，导入abaqus中，此处突起的长度D_C = 10 mm，高度H_C = 1 mm，外部Part为一正方体，边长为4D_C = 40 mm。其次，在结合面上方建立一个刚性的平面，设立刚性平面与接合面的面面接触关系，在接触性质处设定为“硬接触”，这种接触方式虽然会在接触面分离时产生剧烈的压力变化，但能够有效限制计算中可能发生的穿透现象，而本文所研究的模型中并不存在接触分离的情况，因此此种接触形式非常适合此模型。在纯弹性变形下，分别设置了μ_int = 0.1、0.2、0.5、0.8、15共5种情况；在塑性材料变形情况下，分别设置了μ_int = 0.1、15共2种情况。此外，在接触界面同样设置了不同的摩擦系数来模拟计算各种不同情况下的结合面摩擦特性[5]。

Figure 2. Contact settings for interface simulation

图2. 结合面接触设置

如图2所示为abaqus软件中的结合面接触设置，其中$IntProp-1$ 即为前文所提到的摩擦系数设置以及硬接触设置。

2.2. 加载示意

(a)

(b)

Figure 3. (a) Loading conditions; (b) loading step configuration

图3. (a) 载荷设置；(b) 加载步骤设置

如图3所示加载时，本文参考文献的加载方式，选择并计算了 $U_{11}=\frac{1}{4}*H_C=0.25\text{ mm}、U_{12}=\frac{1}{2}*H_C=0.5\text{ mm}$ 两种情况。在加载法向$\left(Y\right)$ 位移后，切向$\left(X\right)$ 位移选择

U₂ = 0.1 mm。切向加载时，参考了滑动起始准则约束方面，可变形的平面在底部受到法向方向的约束，而两侧面受到切向方向的约束。在abaqus中，刚性平面需要通过一个RP点来对其进行处理。仿真中刚体平面的设定是为了近似等效于真实结构中具有显著刚度的机床滑轨或支承件，从而能够真实反映机床中的相对运动。

本文在此处选择了刚体平面的中点作为参考点(RP)，通过此参考点可以对刚体平面进行位移加载。本文材料本构关系采用理想弹塑性模型，其中弹性阶段满足胡克定律，屈服后进入理想塑性区。该模型适用于结构分析中对金属接触面发生局部塑性屈服的情境，能够较好反映局部压入引发的材料非线性响应。

2.3. 网格划分及分析步设置

本文的模型网格如图4所示，模型中的元素总数为29,992，节点数为30,753。运用了从结合面突起部分至主体部分由密到疏的划分方式。这种方式可以在保证计算的准确性的前提之下，尽量减少模型的复杂程度，从而一定程度上加快模型的运算速度。

Figure 4. Mesh generation

图4. 网格划分

在分析步的设定方面，由于Abaqus本身自带且必须存在的intial分析步并不能添加类似接触，位移载荷等具体计算步骤，因而本文的模型均从第一分析步开始添加对模型的具体设定：在第一步添加约束以及接触、在第二步添加小量的刚体平面法向位移载荷、在第三步扩展这个位移载荷至所需大小、第四步添加切向位移载荷。整个过程为线性过程，前后并不相互干扰，同时在后处理阶段可以分别查看每一个分析步的情况并根据需要进行调整。

这样的一种分析步的处理方式也同样可以在计算完成后的提取数据进行分析的阶段有所帮助，abaqus的软件输出的曲线是关于每一个step的时间和所提取的变量拟合而成的一个曲线，有序的分析步设置可以让提取的过程更加有序以及直观。

3. 结果与讨论

3.1. 计算结果

Table 1. The purely elastic deformation data were obtained at an indentation depth of 0.25 mm

表1. 纯弹性变形数据，压入深度0.25 mm

位移(mm)	切向力(N)	位移(mm)	切向力(N)	位移(mm)	切向力(N)	位移(mm)	切向力(N)
0.1	0.0990932	0.1	0.374317	0.1	0.190231	0.1	1.49195
0.3	1.19087	0.3	0.588846	0.3	5.07307	0.3	7.19695
0.4	1.63825	0.4	−0.632631	0.4	3.13409	0.4	5.84563
0.5	1.23368	0.5	−0.146674	0.5	6.62809	0.5	9.62272
0.6	1.12943	0.6	−0.643686	0.6	7.30214	0.6	9.92175
0.7	1.26703	0.7	−0.371441	0.7	8.74079	0.7	10.4491

续表

0.8	0.999673	0.8	−0.818716	0.8	8.42823	0.8	10.5779
0.9	0.993278	0.9	−0.191602	0.9	10.0938	0.9	14.6068
1	1.27424	1	−1.55052	1	11.5299	1	18.3475
1.1	597.147	1.1	613.59	1.1	685.968	1.1	705.451
1.3	1726.7	1.3	1807.55	1.3	2015.24	1.3	2066
1.4	1834.17	1.4	2383.52	1.4	2667.95	1.4	2743.54
1.5	1834.22	1.5	2938.01	1.5	3306.83	1.5	3408.33
1.6	1834.17	1.6	3456.32	1.6	3938.5	1.6	4073.35
1.7	1834.17	1.7	3659.3	1.7	4557.31	1.7	4737.88
1.8	1834.17	1.8	3659.63	1.8	5178.6	1.8	5390.84
1.9	1834.17	1.9	3659.29	1.9	5783.31	1.9	6044.02
2	1834.17	2	3659.32	2	6383.44	2	6682.07

Table 2. The purely elastic deformation data were obtained at an indentation depth of 0.5 mm

表2. 纯弹性变形数据，压入深度0.5 mm

位移(mm)	切向力(N)	位移(mm)	切向力(N)	位移(mm)	切向力(N)	位移(mm)	切向力(N)
0.1	2.62064	0.1	0.19265	0.1	1.27617	0.1	−1.0293
0.3	0.554199	0.3	−0.459147	0.3	9.0835	0.3	10.0312
0.4	−0.281142	0.4	−0.685041	0.4	9.39094	0.4	14.7936
0.5	−1.99673	0.5	−1.45569	0.5	10.2413	0.5	21.172
0.6	−2.28773	0.6	−0.844133	0.6	13.0753	0.6	25.3416
0.7	−2.14659	0.7	0.118503	0.7	13.6118	0.7	26.377
0.8	0.301895	0.8	0.923786	0.8	13.5656	0.8	22.4385
0.9	−3.34389	0.9	−1.61808	0.9	13.7785	0.9	22.3383
1	−4.52395	1	−0.840369	1	14.3568	1	24.5442
1.1	651.309	1.1	683.571	1.1	770.866	1.1	792.46
1.3	1949.43	1.3	2037.4	1.3	2263.47	1.3	2317.94
1.4	2584.27	1.4	2702.39	1.4	2998.11	1.4	3079.77
1.5	3196.36	1.5	3363.17	1.5	3729.61	1.5	3834.16
1.6	3769.22	1.6	4013.24	1.6	4453.27	1.6	4584.13
1.7	3999.92	1.7	4653.29	1.7	5177.8	1.7	5332.81
1.8	3999.82	1.8	5284.58	1.8	5894.07	1.8	6081.14
1.9	3999.64	1.9	5902.35	1.9	6606.99	1.9	6822.12
2	3999.64	2	6506.31	2	7308.03	2	7558.09

如表1，表2所示为本文计算的数据，提取数据需要在step步骤对计算的历史输出进行定义，虽然将结果全部输出在操作方面更加便捷，但会大大加长运算时间，所以在定义输出时，本文仅选择了必要的变量导出。此处同样使用到了上文所提到的刚体平面的中点，即其参考点(RP)，输出其计算过程中的法向力、切向力、法向位移、切向位移。此处刚体平面的性质保证了RP点上的力与位移即结合面的合力与位移。

3.2. 结果分析

如图5所示，所示为摩擦系数μ_int = 0.1、0.2、0.5、0.8、15，结合面压入深度为0.5 mm时切向力Q与切向位移$U_x$ 的图表，在此图中，可以看出当压入深度与材料弹性模量一致时，在到达峰值前，纯弹性简化结合面的切向力与切向位移成正比，并且其比值随着摩擦系数的增加而增加。这个增加的幅度与摩擦系数的大小成反比，即摩擦系数越大，其增长对于这个比值的影响会越小。

Figure 5. Purely elastic deformation—indentation depth: 0.5 mm

图5. 纯弹性变形–压入深度0.5 mm

Figure 6. Purely elastic deformation—indentation depth: 0.25 mm

图6. 纯弹性变形–压入深度0.25 mm

同时也可以看出在摩擦系数为0.1的曲线上，其摩擦力到达峰值前，这个比值将会产生一段突变减小，根据分析比较发现，这是突变是缘于软件的分析时对于这个分析步所设置的分析步过少导致的，并不影响其峰值大小。

如图6所示，所示为摩擦系数μ_int = 0.1、0.2、0.5、0.8、1.5，结合面压入深度为0.25 mm时切向力Q与切向位移$U_x$ 的图表。其与图7大致一致，都可以得出切向力与切向位移成正比的结论。而在此图中，作为比较修正了上文所提到的峰值前曲线的斜率有突变的情况，并且纵向的对比两张图表的具体数据可以发现，当压入深度更小时，相同摩擦系数以及弹性模量的材料，在法向载荷下承受切向载荷时的切向力峰值相较而言更小。经过分析，本文认为这是由于压入深度较小时，模型的结合面接触面积更小，从而导致的摩擦力峰值更小。

Figure 7. Plastic deformation—indentation depth: 0.5 mm

图7. 塑性变形–压入深度0.5 mm

如图7所示，所示为摩擦系数$\mu \_int$ = 0.1、15，屈服极限$\sigma \_\max$ 分别为200 MPa与600 MPa，压入深度为0.5 mm的各个模型计算结果制成的图表曲线。可以看出，在纯弹性的模型加入了屈服极限之后，切向力的大小同样随切向位移增加而增加，然而与上述的纯弹性变形的模型结果有所不同的是，切向力与切向位移之间不再互成正比，可以看出，图表上的曲线的斜率随着切向位移的增加而逐渐减小，并且摩擦系数相对较小的模型其切向力最终停止增长且随着切向位移的增大趋于平稳，即标志这结合面之间的摩擦力达到峰值，结合面与结合面直接开始滑移，在这种条件下，不难发现，摩擦系数较大的材料，即使在切向力停止增长之前，相同切向位移的条件下，其切向力也更加大，这与图7纯弹性材料的结论保持一致。除此之外，还可以看出屈服极限$\sigma \_\max$ 较大的模型产生切向载荷同样较大，并且其峰值也高出许多，这可能与材料的塑性变形所产生的残余应力有所关联。

如图8所示，所示为摩擦系数$\mu \_int$ = 0.1、15，屈服极限${\sigma }_{\max }$ 分别为200 MPa与600 MPa，压入深度为0.25 mm的各个模型计算结果制成的图表曲线。本组分析结果在曲线形态上基本与上文提到的图7所式结果一致。即切向力的大小同样随切向位移增加而增加，而摩擦系数与屈服极限增大时，切向载荷也会随之增大。相互作为对照时，可以看出，在相同简化模型的以及相同材料的前提下，压入深度更大的结合面在相同切向位移的条件下产生的切向力更大，其切向力的峰值也更大，这个结论与纯弹性型得到的结论一致。同时从绿色以及橙色两条达到峰值的曲线具体数值中还看出，相同材料参数的模型在压入深度更小时其切向载荷达到峰值的力更小。这是由于在压入深度更小时，其模型的法相压力相对而言也更小，即可得切向力大小的峰值与法向压力存在关联。但其增长的速率，即曲线的斜率没有明显得变化，可能仅与切向载荷的大小有关。

Figure 8. Plastic deformation—indentation depth: 0.25 mm

图8. 塑性变形–压入深度0.25 mm

Figure 9. Tangential/normal force ratio in the plastic regime (0.25 mm depth)

图9. 切向力/法向力，塑性变形–压入深度(0.25 mm)

Figure 10. Tangential/normal force ratio in the plastic regime (0.5 mm depth)

图10. 切向力/法向力，塑性变形–压入深度(0.5 mm)

Figure 11. Contact area-plastic deformation at 0.25 mm indentation depth

图11. 接触面积，塑性变形–压入深度0.25 mm

如图9、图10所示为摩擦系数μ_int = 0.1、1.5，屈服极限${\sigma }_{\max }$ 分别为200 MPa与600 MPa，压入深度分别为0.25 mm与0.5 mm时的模型计算数据曲线。在Q/P，即切向力/法向力比切向位移的图形中，其单条曲线在材料相同的情况下变化趋势与上文切向力比切向位移的图形十分相似，这是由于模型的分析步被设置为先下压，再平移，继而此处的法向压力P基本上为一定值，即对应材料结合面下压分别为0.5 mm、0.25 mm时结合面的法向压力。从此图中还可以得出，当摩擦系数相同时，屈服极限的大小虽然与切向力的大小成正比，但当切向力达到峰值时，屈服强度并不影响最终Q/P的数值，在目前本文的分析范围内，简化平面Q/P的峰值仅与摩擦系数有关。且最终值与摩擦系数的值几乎相同，在${\mu }_{int}=0.1$ 时能得到Q/P的峰值，即为0.1左右。

如图11所示为摩擦系数${\mu }_{int}=0.1$ ，屈服极限${\sigma }_{\max }$ 分别为200 MPa与600 MPa，压入深度为0.5 mm的模型计算结果。对于本图像所示的接触面积以及法向位移的关系曲线来说，其余模型提取的数据曲线中，摩擦系数对此图表几乎不产生影响；同时，法向载荷$U_y$ 的改变仅仅演唱了曲线的长度，其总体趋势和性质并没有变化，因此本图表仅列出所示的两组数据。如图10所示，随着法向位移的增大，结合面的接触面积呈现震荡上升的趋势，且纯弹性材料的结果上同样发现了类似的情况。可以推测这并不是由于屈服极限导致的震荡现象，推测这可能与隐式计算的收敛方式与网格划分有所关联。除此之外，可以看出屈服极限${\sigma }_{\max }$ 为200 MPa的模型在同一法向位移的情况下接触面积更大，这是由于其接触面在受压的情况下更加容易产生屈服，即其表面显得更“软”。

3.3. 云图结果

除开上述的计算结果曲线之外，本文得到了较多mises应力云图。如图12所示，这份mises应力云图是${\sigma }_{\max }$ = 200 MPa，${\mu }_{int}=0.1$ ，$U_y=0.25\text{ mm}$ 的模型的计算结果，在图中可以看出，其模型中心应力较大的区域主要是呈现一个新月形的分布情况，围绕突起区域呈现一个非轴对称的分布。而在突起区域两侧则转变为较为对称形状，可以看出距离中心区域越远的单元体所受的应力越小，最后趋向于零，且越远越趋向于轴对称分布。除此之外，由于切向载荷是向右方加载的，因此即使两侧应力大致相同，中点右边的单元体承受的应力仍然是大于左侧单元体的。

Figure 12. Von mises stress contour plot-elastic deformation

图12. 弹性变形mises应力云图

在研究的过程之中，我们发现模型的屈服极限对结果曲线的形态与数值影响很大，因此在图13、图14中列出了模型的塑性应变云图，可以看出塑性变形的区域同样集中在模型中间的突起区域，且其变形较大的区域形成了一个倒三角形的包络线，变形最大的区域在结合面的左侧收到切向载荷拉伸的区域。而在结合面的中间区域塑性变形为零，且向包络线慢慢过渡变大。在倒三角形包络线的外部塑性变形的区域类似一个长方体，再向外扩展直至减少为零。值得一提的是，与实现设想的不同，在结合面的塑性应变方面，从文本的计算结果来看，材料的屈服极限的改变所引起的塑性应变的变化远远不如结合面的摩擦系数的改变所引起的塑性应变变化范围来的大。本文认为，这主要是由于材料的屈服极限相对来说改变的幅值不如材料的摩擦系数大，但综合来说，结合面的摩擦系数以及材料的屈服极限的增大，都可以有效的增强结合面抵抗塑性应变的能力。

Figure 13. Von mises stress contour plot-plastic deformation (σ_max = 200 MPa, μ_int = 0.1, U_y = 0.25 mm)

图13. 塑性变形mises应力云图(σ_max = 200 MPa, μ_int = 0.1, U_y = 0.25 mm)

Figure 14. Von mises stress contour plot-plastic deformation (σ_max = 200 MPa, μ_int = 15, U_y = 0.25 mm)

图14. 塑性变形mises应力云图(σ_max = 200 MPa, μ_int = 15, U_y = 0.25 mm)

如图15、图16所示为此模型的位移云图，在图中可以看出，在接触面的中心，即模型突起部分的位移最大，为0.269 mm，且呈椭圆形向外扩展，且基本沿中心线呈对称形状，稍稍向右侧有所偏移。在凸起部分的外部区域变形为零。从其云图的对称性可以看出，模型的位移主要是由法向载荷引起的，而切向载荷的影响较小，这也与模型的形状与加载方式有所关联。

Figure 15. Element displacement contour plot (σ_max = 200 MPa, μ_int = 0.1, U_y = 0.25 mm)

图15. 单元位移云图(σ_max = 200 MPa, μ_int = 0.1, U_y = 0.25 mm)

Figure 16. Element displacement contour plot (σ_max = 600 MPa, μ_int = 0.1, U_y = 0.25 mm)

图16. 单元位移云图(σ_max = 600 MPa, μ_int = 0.1, U_y = 0.25 mm)

4. 结论

机床制造凭借其灵活性，高精度性以及成熟的体系，广泛应用于航空航天、汽车工业、能源工程、风力发电及船舶工业等领域。而机床结构在工作的过程中，其各个部件的受力情况难免会收到复杂的外力作用。由于高精度的制造对于其机床本身的精度要求很高，因此对于机床结合面的研究对于机床发展有深远的意义。本文主要利用大型商用有限元分析软件ABAQUS/Explicit，建立了简化情况下以及真实情况下的结合面模型，并对其进行了研究，得出以下结论。

(1) 纯弹性变形的模型中可以看出当压入深度与材料弹性模量一致时，在到达峰值前，纯弹性简化结合面的切向力与切向位移成正比，并且其比值随着摩擦系数的增加而增加。这个增加的幅度与摩擦系数的大小成反比，即摩擦系数越大，其增长对于这个比值的影响会越小。

(2) 塑性变形的模型方面，本文研究的结合面切向力在切向力比切向位移到达顶峰前始终处于一个上升的趋势，同时各个模型在压入深度、屈服极限越大时，切向力越大。在摩擦系数较小时切向力在达到最大值，即法向力 * 摩擦系数时趋于平稳，即标志着结合面开始相对移动。而摩擦系数较大时，切向力随着切向位移的增长速率也会逐渐趋于平稳。

(3) 接触面积方面，其随着法向位移的增大，结合面的接触面积呈现上升的趋势，并在上升过程中有一定的曲线震荡的情况。在模型其余设置相同的情况下，屈服极限为300 Mpa的模型结合面接触面积增长更快且到达峰值更快。

综上所述，本文的结论可以对机床的制造以及运用到机床的零件加工进行减少摩擦力和降低损耗上的理论指导。