一种有效的螺纹固定结合部有限元建模方法

1 引言

机床动态特性对其加工精度、性能以及寿命存在重要影响，在设计阶段建立高精度整机动力学模型是评价、优化其动态特性的关键手段。机床是由众多零部件通过各类结合部装配在一起的组合体，其中螺纹固定结合部最为常见，它指的是两零部件通过螺纹紧固件连接成整体时其相互接触的部分。结合部动态特性对整机动态特性影响显著，研究表明结合部刚度占整机刚度的（60～80）%，整机阻尼的90%以上源自结合部[1]。因此，建立合理的结合部动力学模型将有助于提高机床整机动力学建模的精度。然而结合部的形状、材料属性、结合面粗糙度及所受的外部载荷对其传力特性影响显著，如何建立能准确描述结合部传力特性的动力学模型，既是机床整机动力学建模研究的热点也是当前的一个技术难题。国内外众多学者开展了大量关于结合部动力学建模的研究，并且取得了许多有价值的成果[2-3]。文献[4]将结合面等效成一个节点，建立了结合面间六自由度的等效动力学模型，并计算了每个自由度的等效刚度与阻尼；文献[5]通过在结合面间建立一系列弹簧阻尼单元来描述结合部的动态特性，单元的刚度与阻尼可以通过解析计算和参数辨识的方法获得，由于该方法简单实用，所以被众多学者和工程师采纳；文献[6]认为各种条件下的结合面（不同材料、表面粗糙度、介质等）动态特性均可用一个理想结合面的动态特性来描述，理想结合面上任何位置处的力学性质均相同，据此提出了用于建立理想结合面动力学模型的有限元单元，并基于有限元的相关理论推导出了单元的刚度矩阵和阻尼矩阵，矩阵中的参数是通过参数辨识的方法获得。

上述研究虽然在很大程度上解决了螺纹固定结合部动力学建模问题，但是忽略了结合面压力不均匀分布对结合部建模精度的影响。螺纹固定结合部在螺纹紧固件的预紧力作用下会出现结合面压力分布不均匀的现象，而结合面接触压力是决定其接触刚度与阻尼的关键因素之一[7]。将结合面面压不均匀分布纳入到螺纹固定结合部动力学建模过程中可以有效改善其建模精度。因此，给出一种基于结合面压力不均匀分布的螺纹固定结合部有限元建模方法，建模流程简述如下：首先，在Ansys中对研究对象进行网格划分，在螺纹连接处施加实际预紧力并进行静力学分析；其次，提取结合面上的单元信息与节点信息，基于节点间距离最短原则对上、下结合面中的节点进行匹配，通过Ansys的APDL语言编写程序，利用该程序在匹配好的两节点间建立Matrix27刚度单元与Matrix27阻尼单元来模拟结合部的动力学特性。基于Yoshimura法与结合面面压分布计算出结合面刚度、阻尼分布，将Matrix27单元位置处的结合部刚度、阻尼系数带入Matrix27刚度、阻尼单元中；最后，将结合部有限元模型与各零部件有限元模型综合成整体有限元模型。

2 静力分析

通过调查各种类型的数控机床，最终确定以装配体为研究对象并进行螺纹固定结合部动力学建模研究，如图1所示。图中的装配体由上块体与下块体通过4个均布的性能等级为8.8级的M8螺栓连接而成，每个螺栓施加的预紧力矩为25N·m，其结合部为典型的螺纹固定结合部。组成装配体的上块体与下块体均由45号钢制成，其材料属性，如表1所示。

图1 装配体结构与尺寸
Fig.1 Assembly Structure and Size

表1 材料属性
Tab.1 Material Properties

所给螺纹固定结合部动力学建模方法由于引入了结合面压力分布不均匀因素，因此对施加预紧力的装配体进行静力学分析并获得结合面的压力分布是结合部动力学建模的先决条件，是获得结合面刚度、阻尼分布的依据。装配体结构中存在诸多曲面元素（例如孔面和圆柱面），在保证力学性能被准确模拟的前提下，利用Ansys中对几何形状具有较好拟合特性的solid187单元进行网格划分。solid187为10节点四面体单元，每个节点具有3个沿坐标轴方向的自由度，与传统4节点四面体单元不同，solid187可以模拟单元内部的应力变化，更贴近于实际结构的传力特性，并且相比于同样能模拟应力变化的六面体单元具有计算效率高，对几何域拟合性好的优点，因此选用solid187来模拟零部件本体的力学性能。在有限元模型中每个螺栓位置处施加25N·m的预紧力矩，每个螺栓所受到的预紧力为[8]：

式中：F0—螺栓的预紧力；T—螺栓的预紧力矩；K—拧紧力矩系数。基于螺栓的实际安装工况，查询机械设计手册，取K值为0.2；d—螺栓的公称直径。

在螺栓预紧力作用下，结合面中螺纹孔周围会出现应力集中现象，根据弹性力学的相关理论可知，在此受力状态下结合面将呈现出波浪形的变形曲面，结合面在螺纹孔周围发生接触，远离螺纹孔的位置处发生分离，只有接触部分存在接触面压，未接触部分的接触面压为零。基于上述分析，将有限元模型中的接触面设置成Frictional型接触，使用该接触模型后，在计算过程中接触面之间既可以发生法向分离，也可以产生切向滑动，符合结合面在静力作用下的受力状态。调用Ansys中静力学分析模块对装配体有限元模型进行静力学分析，提取结合面中的压力分布，如图2所示。由此可见，压力在结合面中的分布是不均匀的，螺纹孔周围出现应力集中，远离螺纹孔的结合面面压逐渐减小并趋近于零。由此说明，分析得到的受力状态符合结合面实际受力状态。

图2 结合面面压分布
Fig.2 Pressure Distribution in Joint Surface

3 螺纹固定结合部有限元建模

螺纹固定结合部的结合面在微观上不是绝对光滑的平面，而是凹凸不平的曲面。螺纹紧固件将两个零部件连接成为一整体，结合面间的凸起彼此接触形成了一个个微接触对，微接触对在零部件间起到了力及运动的传递作用，并且微接触对在结合面上的分布具有随机、均匀、离散的特性。在有限元模型中，节点同样起着力及运动的传递作用，并且结合面上的网格节点也具有随机，均匀，离散的分布特性，节点在结合面间的传力、分布特性与微接触对在结合部中的传力、分布特性相似，因此采用有限元模型来描述结合部的传力特性是合理的（随着网格密度的提高，结合部有限元模型的传力特性将逐渐趋近其真实传力特性）。研究表明结合部具有存储能量与消耗能量的特性[9]，因此可以在上、下结合面的节点间建立刚度、阻尼单元以此来模拟结合部的动态特性。对于结合部而言，微接触对的质量可以忽略不计，只考虑其刚度与阻尼特性。上块体结合面中的网格节点与下块体结合面中的网格节点可能会出现部分重合与分离的情况，而Matrix27刚度、阻尼单元可以描述任意距离的两节点间的传力特性，因此所给方法采用Matrix27刚度、阻尼单元来描述结合部中两节点间的传力特性。

3.1 节点配对

Matrix27刚度、阻尼单元是两节点单元，每个节点有6个自由度（其中3个沿坐标轴方向的水平自由度，3个绕坐标轴运动的转动自由度），单元的刚度、阻尼矩阵为用户自定义。由于Matrix27刚度、阻尼单元的作用是模拟结合面间的传力特性，因此单元节点应分别位于上、下结合面中。对上、下结合面中的节点进行配对是在结合部中建立Matrix27刚度、阻尼单元的前提，节点配对的基本原则为：（1）节点对要均匀、离散的分布在结合面中。该原则能够保证结合部有限元模型中起运动、力传递作用的Matrix27刚度、阻尼单元在结合部中的分布特性与真实结合部中起运动、力传递作用的微接触对的分布特性吻合，并且节点对的均匀、离散分布有助于实现结合部压力不均匀分布对结合部动态特性的影响；（2）每个节点对的空间长度为所有配对方案中的最小值且不高于设定值ε。其中，ε是保证Matrix27刚度、阻尼单元在结合部中呈均匀、离散分布。结合面间节点配对具体流程，如图3所示。

图3 结合面节点配对流程图
Fig.3 Node Pairing Flow Chart of Joint Surface

首先分别提取上、下结合面的节点信息（包括节点编号与节点坐标），其次在Matlab中编写算法，分别计算出某结合面上的每个节点相对另一结合面上所有节点间的距离，计算公式，如式（1）所示。

式中—下结合面中编号为n的网格节点与上结合面中编号为N的网格节点的x，y，z坐标值。取距离最短的两节点编号并判断其空间距离L是否小于ε，如果满足条件，则确定该两节点对为配对节点。基于Ansys的APDL语言编写二次开发程序，在配对好的每个节点对间分别建立一个Matrix27刚度单元与Matrix27阻尼单元。经过上述工作，结合面间的所有Matrix27刚度、阻尼单元建立完成。经统计，所给研究对象两结合面共配对676对节点对，且节点间距离最大为2.1mm。由此说明任一节点对中的两节点在空间内近乎吻合，满足节点配对要求。结合部中Matrix27刚度、阻尼单元总数均为676个，每个Matrix27刚度、阻尼单元模拟了节点对在法向、切向的刚度、阻尼特性，所有配对节点在结合面上的位置分布，其中实心点为节点对中属于上块体结合面的节点分布，圆圈为节点对中属于下块体结合面的节点分布，如图4所示。

图4 结合面节点分布图
Fig.4 Node Distribution Map of Joint Surface

将结合面均匀分块，计算各子区域中节点对的密度，计算方法，如式（2）所示。

式中：Nn—编号为n的子区域中节点对的个数；Sn—该子区域的面积；ρn—区域的节点对密度。

接着计算结合面所有区域节点对密度的均值及相对于均值的偏差，如式（3）、式（4）所示。

式中：n—子区域的编号；N—子区域的个数；E—节点对密度分布的均值；δn—编号为n的子区域节点对密度分布相对于均值的偏差，选取方差为最小值时所对应的ε，此处取ε=3mm。

3.2 结合面面元划分

结合面面元是指零部件有限元模型中的单元与结合面重合时的轮廓形状。对结合面进行面元划分是在已有的有限元网格基础上，提取单元与结合面重合的所有网格节点构成一个集合，该集合存放着一个面元中的所有节点编号，是求解结合面刚度、阻尼分布的关键。具体流程如下：分别提取上块体结合面与下块体结合面中的所有节点编号及坐标并以集合的形式表示，记为A1，A2。分别将集合A1，A2中属于同一单元的节点编号进行归类，以集合Bn表示，其中n代表集合的序号，Bn中的所有节点属于单元与结合面重合的轮廓，则将Bn称为编号为n的结合面面元。分别将上块体结合面与下块体结合面中每个面元的轮廓在Matlab中绘出，如图5（a）、图5（c）所示。由于对装配体进行有限元网格化分时采用的是Solid187单元，因此结合面上的面元形状均为三角形。对比Ansys在结合面处的有限元网格，如图5（b）、图5（d）所示。由图5可见，Matlab绘制的网格与Ansys网格完全吻合，说明上述关于结合面面元的划分工作是正确的。

图5 结合面网格分布图
Fig.5 Mesh Grid Distribution Map of Joint Surface

3.3 结合面刚度与阻尼分布计算

Matrix27刚度单元的刚度矩阵，如式（5）所示。式中kxx、kyy、kzz—两节点间沿x、y、z方向的刚度系数。Matrix27阻尼单元的表达式与刚度单元相似。螺纹固定结合部中所有Matrix27刚度单元的刚度系数计算流程：首先分别计算下块体结合面内每个面元的每条轮廓函数表达式fn1、fn2、fn3，其中n表示面元的编号，1、2、3分别表示面元轮廓线的序号，如图6所示。采用式（6）来描述每个结合面面元内的压力场分布，如式（6）所示。

图6 结合面元斜边表达式
Fig.6 Oblique Edge expression of Element

式中：n—面元的编号；x，y—面元所在平面的横、纵坐标；Pn(x，y)—面元内任意一点的压力值。提取下块体结合面内所有节点的法向压力值，基于每个面元节点的坐标及面压值，辨识出参数P00-n，P01-n，P11-n。最后在Matlab中绘制出结合面的压力分布，如图7所示。

图7 结合面拟合压力分布图
Fig.7 Pressure Distribution of Joint Surface

基于式（7）和式（8）计算出每个结合面面元的平均面压，式中Pmy-n表示结合面上第n个面元的平均面压，Pn(x，y)为拟合获得的结合面中编号为n面元的压力场分布函数，Sn为此面元的面积。结合面中的某个三角形面元，如图6所示。

Yoshimura法认为任意固定结合部在相同表面粗糙度及面压条件下，单位面积的刚度与阻尼系数相同。基于此思想，首先根据结合面面元平均面压计算得到每个面元区域内的单位面积法向、切向刚度系数，将单位面积法向，切向刚度系数乘以该面元面积可以得到面元的总法向、切向刚度；其次将面元的总刚度均分至该面元的三个角点，由于同一个节点会被多个面元共用，因此，利用式（11）计算出各节点的总刚度：

式中：KM—编号为M的节点总刚度；n—面元编号；N—拥有共同节点M的面元个数。

文献［10］通过实验获得了不同接触条件下（结合面材料、有无介质、表面粗糙度）的结合部单位面积法向、切向刚度、阻尼系数与结合面面压之间的数学表达式。则所给研究对象结合面压力与单位面积法向、切向刚度的数值关系，如式（9）、式（10）所示。

式中：Pn—结合面接触面压；Kn—单位面积法向刚度；Kt—单位面积切向刚度，结合面粗糙度Rɑ=1.6µm。

基于Matlab编写可实现上述操作的程序，计算获得的结合面节点刚度分布，如图8、图9所示。由此可以发现，螺拴孔附近的刚度相对其周围明显要大，这是由于在螺拴预紧力作用下，螺拴孔周围压力集中，因此刚度均较周围要大。

图8 结合面法向刚度分布
Fig.8 Normal Stiffness Distribution of Joint Surface

图9 结合面切向刚度分布
Fig.9 Tangential Stiffness Distribution of Joint Surface

下面通过式（12）～式（15）求解每个节点阻尼。

式中：cn、ct—结合面上每个配对节点间的法向阻尼与切向阻尼系数；C′n、C′t—结合面的法向、切向总阻尼值；N—结合面间配对的节点对总个数；S—结合面总面积；C、C′—某一面压值下单位面积法向、切向阻尼值。

最后将所求出的节度、阻尼带入到与之相对应的Matrix27刚度、阻尼单元中，至此，螺纹固定结合部有限元模型建模工作结束。

4 模型验证

将螺纹固定结合部有限元模型与装配体中的上、下块体有限元模型综合成整体有限元模型，调用Ansys模态分析模块对整体有限元模型进行模态分析，提取前四阶模态频率的计算结果，如表2所示。为验证该建模方法的有效性，基于LMS振动噪声测试系统对装配体进行锤击法模态测试，试件采用尼龙绳进行悬挂，通过力锤激励，ICP压电式三向加速度传感器拾取测试点的振动信号，运用LMStest.Lab软件提取装配体的前四阶试验模态频率，如表2所示。

表2 理论模态值与实验模态值对比
Tab.2 Comparison between Theoretical Modal and Experimental Modal

试验获得的所有激励点与测试点间的频响函数曲线及其相位图，如图10所示。试验采集设备包括PCB公司的086C03型力锤和356A16型三向加速度传感器，LMS公司的SCM202型16通道数据采集仪。

图10 频响函数幅值、相位图
Fig.10 The Graph of Frequency Response Function Amplitude，Phase

通过对比装配体前四阶理论模态频率与试验模态频率，最大误差出现在第一阶模态处，为7.2%，其余均在5%以内。将装配体的前四阶理论与试验模态振型对比于图11，通过观察发现，前四阶理论模态振型与前四阶试验模态振型相似且一一对应。由此说明，所给螺纹固定结合部动力学建模方法是合理可行的。

图11 理论模态振型与实验模态振型对比
Fig.11 Comparison between Theoretical Modal and Experimental Modal

5 结论

（1）给出了一种构建螺纹固定结合部动力学模型的新方法。通过在结合面节点间建立Matrix27刚度单元和Matrix27阻尼单元来描述结合部的动态特性，每个刚度、阻尼单元中的刚度系数、阻尼系数通过计算获得，并且考虑了结合面面压的分布对刚度数值大小的影响，即结合面面压越大，该处的刚度值就越大；（2）对所研究的装配体试件进行了模态测试，通过对比试件前四阶的理论模态与试验模态发现，理论模态的前四阶振型与试验模态的前四阶振型相似且一一对应，前四阶理论固有频率与试验固有频率间的误差在（0.7～7.2）%范围内，从而证明了所给建模方法的有效性；（3）所给建模方法由于在建模过程中不需要通过实验来获得相关建模所需参数，因此具有成本低、建模效率高、通用性好的优点。