基于非均匀分布复弹簧单元的螺栓连接薄板结构动力学有限元建模

刘晓峰， 孙 伟， 方自文

(1. 东北大学机械 工程与自动化学院， 沈阳 110819；2. 东北大学 航空动力装备振动及控制教育部重点实验室， 沈阳 110819)

螺栓连接几乎应用在所有机械装备中，其连接界面(结合面)的力学特性对整个机械装备的静力学及动力学性能均有重要的影响。因而，在机械装备设计研发过程中必须重视对螺栓结合部的性能分析与设计。而创建一个简单、高效并且具有较高精度的结合部模型则是分析和设计螺栓连接结构的重要基础。

长期以来已有大量学者致力于螺栓结合部的力学建模与分析，真实的螺栓结合部处存在螺栓、螺母等构件，在建模时若将其体现出来，建模过程会较复杂。这类高保真建模方式通常适用于静力学分析，例如，Al-Nassar等[1]建立了完整的螺栓连接板三维有限元模型并对其在剪切和拉伸载荷作用下的受力性能进行了分析。Gray等[2]提出了一种高效的螺栓连接整体三维模型建立和验证的方法。为了方便进行螺栓结合部的力学研究，一些建模方法在建模时不体现螺栓、螺母等构件，而是用一些如弹簧、阻尼类的具有集中参数的单元，或是利用一种虚拟介质等来模拟结合面的力学特性。这类简化模型的建模思路又可分为线性建模和非线性建模。

对于线性建模，通常将模拟结合部力学特性的集中参数或虚拟介质的材料参数等视为常量，例如，Ye等[3]将虚拟材料参数导入有限元分析软件，建立了包含结合部虚拟材料模拟的螺栓连接结构有限元模型。Zhou等[4]提出了单自由度黏滑摩擦模型，对螺栓连接复合板进行了多尺度建模和有限元分析，并与试验结果对照，验证了模型的正确性。McCarthy等[5]用弹簧单元模拟结合部建立了三螺栓连接复合材料梁的分析模型，研究螺栓孔间隙对多螺栓连接结合部载荷分布的影响。

实际螺栓连接结构，由于受激振力变化等因素的作用，其结合部处通常会产生非线性的力学现象，此时线性建模的思想不再适用，需进行非线性建模。如，Zhao等[6]提出了一种用非线性虚拟材料来描述螺栓连接力学特性的方法。Taheri-Behrooz等[7]建立了考虑节点非线性材料特性的单搭接多节点非线性弹簧-质量模型来模拟螺栓结合部。 Kong等[8]提出了一种包含非线性接触力、阻尼和连接界面几何特性的螺栓连接结构降阶模型，降低了求解非线性动力学方程的计算成本。Armand等[9]提出了一种新的多尺度方法，将螺栓结合部界面粗糙度、接触压力和接触刚度联系起来，并与多谐波平衡求解器相结合，可以计算不同界面粗糙度下结构的非线性动态响应。

螺栓结合部的非线性建模往往更能再现连接结构真实的振动行为，但是非线性建模的研究通常源于对线性模型的改进。当前，在用弹簧单元模拟螺栓结合部的研究中，大多是用一根弹簧单元或者刚度相等的分布弹簧单元来模拟螺栓及连接区域的刚度。这种建模方法通常无法考虑螺栓影响区刚度非均匀分布对结构的影响。还有其弹簧单元刚度通常是以结构、螺栓材料和螺栓载荷为参数通过经验计算公式得到，其计算结果的准确性也有待商榷。

螺栓影响区对整个结构系统振动特性的贡献主要体现为刚度及阻尼。因而，在考虑影响区面积的前提下，本文提出一种用刚度非均匀分布的复弹簧单元来模拟螺栓结合部力学特性的建模方法，并在螺栓连接薄板振动特性分析上进行了实践。文中详细描述了上述建模方法的原理，以及利用基于遗传算法的反推辨识技术确定复弹簧单元参数的流程。最后，以用两个螺栓连接在一起的薄板结构为例进行了实例研究，通过与试验的比对证明了上述基于非均匀分布复弹簧单元模拟螺栓结合部力学特性的合理性。

1 基于非均匀分布复弹簧单元的螺栓结合部建模

以下具体描述利用非均匀分布复弹簧单元模拟螺栓结合部力学特性的建模方法，包括建模理念以及复弹簧单元刚度及阻尼参数的确定方法。

1.1 建模理念

螺栓在固定两个连接结构时会在结合面产生一个以螺栓杆轴心为中心的一个非均匀分布近似圆形的压力区域，称为螺栓影响区[10]。在振动分析中，可以说对螺栓结合部的模拟，就是要对这个螺栓影响区的有效模拟，具体涉及确定影响区的面积以及影响区刚度及阻尼特性的分布等。

1.1.1 螺栓影响区的面积

如图1所示，dw为螺栓承载面直径，t为连接板的厚度，D为螺栓影响区域直径。可见，获得了螺栓影响区直径就确定了整个影响区的面积。

(a) 螺栓连接板结构

当螺栓型号及连接结构已定时，可以确定dw及t的值，而半顶角α取值范围一般在25°≤α≤33°[11]，基于上述参数，螺栓影响区直径的计算式可表达为

D=dw+2t·tanα

(1)

按上述理论计算的螺栓影响区面积为圆形，但在有限元建模时(尤其是在自行研发的有限元程序中)，对圆形区域划分网格不够方便，为了简化建模，这里采用一个与圆形区域面积相接近的正方形区域来描述螺栓影响区的力学特性。因而在本研究中，螺栓影响区的面积仅作为建模时添加复弹簧单元区域的参考，并没有将实际螺栓影响区面积引入分析模型。

1.1.2 螺栓影响区的刚度及阻尼分布假定

研究表明[12]，在螺栓影响区中靠近螺栓区域的刚度值较大而远离螺栓区域的刚度值逐渐降低，直至影响区边界刚度降至0。因而，采用均匀分布的弹簧模拟螺栓影响区的力学特性是不恰当的。本文提出用三种刚度非均匀分布弹簧单元来模拟螺栓影响区的力学特性，包括正弦、抛物线和线性分布，如图2所示。从图中可以较为直观的看出，三种分布均反映了螺栓影响区中心处刚度值最大，中心区向外围区刚度值逐渐递减的趋势，因而可以较为真实地描述结合面的力学行为。

(a) 正弦分布

需要说明的是在本文的建模方法中，弹簧单元仅用于模拟螺栓影响区的面外刚度，且每根弹簧的刚度定义为复刚度k*，该复刚度同时包含了结合面的刚度及阻尼，表达为

k*=k+ic

(2)

这里k表征了螺栓结合面的刚度，而c描述了结合面的阻尼，且这是一种结构阻尼。

首先讨论影响区的刚度值分布。基于图2所描述的螺栓影响区刚度分布假定，在影响区内任一弹簧刚度k可由影响区中心处刚度值kmax以及弹簧单元距离影响区中心距离的2倍(即图2中的d)来综合确定。对于正弦、抛物线和线性分布，在影响区内任一弹簧刚度的求解式可分别描述如下

k=kmaxsin[90°(1-d/D)]

(3)

k=kmax[1-(d/D)2]

(4)

k=kmax(1-d/D)

(5)

接着讨论螺栓影响区的阻尼值分布。结合面的阻尼分布与刚度分布通常是不一样的，通常在压力大的地方刚度大而阻尼小[13]。考虑到阻尼机理的复杂性，为了简化建模，将每根弹簧的结构阻尼取同一值(即将阻尼在影响区内视为均匀分布)。

综上，基于非均匀分布复弹簧单元的螺栓结合部建模过程可简要描述如下：① 根据螺栓的类型确定影响区的面积；② 基于此面积完成螺栓影响区网格划分，在相应节点上建立同时包含刚度及阻尼的复弹簧单元：③ 指定一种非均匀分布来描述螺栓影响区的弹簧单元刚度，用均匀分布来描述螺栓影响区的阻尼；④ 利用参数辨识技术确定复弹簧单元的刚度及阻尼值。

1.2 复弹簧单元刚度及阻尼参数确定方法

参见式(3)～式(5)，由于假定刚度值满足一种特定的分布，因而在实际研究中只需要确定中心区最大刚度值kmax，而阻尼值对于每根弹簧单元是一样，即需要确定参数c的值。

这里采用反推辨识法[14-15]确定上述待定参数，反推法基本思想可概括为：将需辨识的量值定义为设计变量(这里为kmax和c)，通过一个匹配算法不断修正理论模型中的设计变量值，使计算获得的结构振动特性参数与试验测得的相一致，进而反推出待辨识的参数。用于模拟螺栓结合面力学特性的复弹簧单元刚度参数直接影响螺栓连接板结构的固有频率，而复弹簧单元的阻尼参数直接影响螺栓连接板结构的频响函数幅值。因而，这里提出通过实测螺栓连接薄板结构的固有频率和频响函数幅值来反推复弹簧单元的刚度及阻尼参数。

图3为利用反推辨识法确定复弹簧单元刚度及阻尼参数的流程，主要包括三项内容，分别是试验测试、理论计算和匹配计算。对于试验测试，这里通过锤击试验获得螺栓连接薄板结构的频响函数，进而由频响函数峰值处对应的频率值确定各阶固有频率，这就为反推辨识提供了完整的试验数据。在理论计算中，这里主要利用自编的有限元程序来完成螺栓连接薄板的有限元建模，详见第2章。

图3 复弹簧单元刚度(阻尼)辨识流程

匹配计算的目标是使理论计算结果与实测值尽快匹配，这里采用遗传算法[16]来实现参数寻优。模型参数的确定是通过将实测的固有频率(频响函数值)与相应的有限元模型预测值之间的差异最小化来实现的，故反推辨识弹簧刚度的目标函数可表达为

(6)

(7)

在基于遗传算法反推复弹簧刚度及阻尼的匹配计算中，首先，需设定好种群数量，即迭代过程中任意一代给出的个体数量(这里个体指复弹簧单元刚度或阻尼参数)、变异概率(种群里个体发生变异行为的概率)、交叉概率(种群里个体发生交叉行为的概率)、迭代次数等参数；接着，实施迭代计算，遗传算法将通过选择、交叉和变异等遗传操作进行概率化的迭代寻优；最后，当达到所设定收敛条件(例如达到最大的迭代次数)后，即可获得待辨识的复弹簧刚度或阻尼值。

2 螺栓连接薄板结构动力学有限元建模及分析

这里用非均匀分布复弹簧单元模拟螺栓影响区的力学特性，完成螺栓连接薄板结构动力学有限元建模与分析。为了更好地实施所研发的建模理念以及反推辨识确认复弹簧单元的刚度及阻尼参数，这里利用Matlab软件自行研发有限元程序。图4为处于悬臂状态的螺栓连接薄板结构，两块薄板由若干个螺栓搭接在一起。以下从板单元的分析、复弹簧单元引入以及螺栓连接板结构振动特性求解描述整个分析的过程。

图4 螺栓连接薄板结构

2.1 薄板单元分析

选用如图5所示的4节点薄板单元，每个节点具有3个自由度，即一个位移自由度w(挠度)、两个转动自由度θx和θy，单元的长度及宽度分别为a和b。

图5 4节点薄板单元

单元节点位移矢量为

δe=[w1θx1θy1…w4θx4θy4]T

(8)

每节点的形函数可表示为

(9)

i=1,2,3,4

式中：ξ，η为介于-1和+1间的局部坐标；ξi和ηi为第i个节点的局部坐标。

根据单元刚度矩阵定义可求出单元刚度矩阵，求解式为

Ke=t∬BTDBdxdy

(10)

式中：B是由几何方程推得的单元应变矩阵;D为由材料参数确定的弹性矩阵

相应地，薄板单元质量矩阵求解式可表示为

Me=t∬NTρNdxdy

(11)

式中：ρ为板的密度；N为由式(9)描述的节点形函数组成的形函数矩阵。

2.2 复弹簧单元的引入与系统组集

在螺栓影响区的每个节点上都有一个同时包含刚度及阻尼的复弹簧单元，其刚度方向对应于上述板单元的位移自由度为w(挠度)，在有限元分析时，需要将这些单元与板单元进行有效的组集。

复弹簧单元的刚度及阻尼矩阵可分别表达为

(12)

(13)

则总的复弹簧单元复刚度矩阵表达为

(14)

引入边界条件后，将弹簧的复刚度矩阵与两块板的刚度矩阵进行组集，整个组建过程可用式(15)进行描述。式(15)中，左上角框圈起来的部分为第一块板的刚度组集，m表示第一块板的节点自由度数，左下角为第二块板的刚度组集，n表示第二块板的节点自由度数，中间部分为搭接部分与复弹簧单元复刚度矩阵的组集，l表示复弹簧单元个数。

(15)

通过这种刚度组集方式，得到螺栓连接结构的复刚度矩阵K*，该复刚度矩阵同时包含了结构系统的刚度及螺栓结合面产生的结构阻尼。螺栓连接薄板结构的总质量矩阵M不涉及复弹簧单元，因而总质量矩阵可表达为

(16)

2.3 螺栓连接薄板结构振动特性求解

悬臂状态下的螺栓连接薄板结构的频域运动方程为

(K*-ω2M)X=F

(17)

式中：ω为激振频率;X为响应向量;F为激振力向量。

用于求解螺栓连接薄板结构模态振型的特征方程可表达为

(18)

式中：KR为复刚矩阵K*的实部，ωr和φr分别为第r阶固有频率和实模态振型。

参见式(17)，螺栓连接薄板结构复频响函数矩阵为

(19)

本文中螺栓连接板结构模型可认为是小阻尼系统，可认为引入的阻尼并不影响系统的特征向量，因而可用实模态理论分析系统振动行为。

利用求得的模态振型矩阵φ，对上述方程做正交变换有

(20)

式中，ηr为由复弹簧单元的阻尼部分产生的第r阶模态损耗因子。

在实测时，通常无法也没有必要获得整个系统的频响函数矩阵。只需获得在i点激励和j点拾振获得的频响函数，这对应频响函数矩阵的一个元素，可表示为

(21)

实际研究中一般只需要考虑前若干阶次即可，即考虑前n0个阶次。

3 实例研究

3.1 问题描述

这部分以螺栓连接钢板为例描述本文提出的用非均匀分布复弹簧单元模拟螺栓结合部力学特性进而完成组合结构动力学建模与分析的方法。图6(a)为所研究的螺栓连接薄板结构，其由两块钢制带孔薄板经螺栓连接紧固构成，相关的几何及材料参数如表1所示。采用锤击法对该螺栓连接薄板结构进行固有特性测试，对螺栓施加3 N·m的预紧力，测试中通过PCB SN 30272力锤对薄板连接结构施加宽频激励，使用Polytec PDV-100激光多普勒测振仪拾振，LMS SCSDAS数据采集分析仪用于获取激励及响应信号，在获取模态振型时，锤击点为图6(a)中标注的所有节点，拾振点位置如图不变，利用获取的频响函数可得到各阶固有频率。而文中分析所选取的频响函数在试验获取时的锤击点位置在图6(a)中做了突出标注，相关测试结果列在后续的与理论分析相对照的各表格中。

(a) 螺栓连接薄板

表1 螺栓连接薄板结构中板的相关材料及几何参数

3.2 螺栓连接薄板动力学建模

用于连接两块薄板的螺栓为M6外六角头螺栓，因而参照式(1)可确定每个螺栓影响区的面积为171.37 mm2。所创建模的有限元模型见图6(b)，每个板共划分了400个单元，882个节点，两块板在螺栓影响区对应节点处用复弹簧单元进行连接。图6(b)中的节点包含了图6(a)中的所有节点，以保证分析获得的频响函数能与试验结果比较。考虑到圆形影响区在有限元建模时难以模拟，因而这里用面积相近的矩形区替换实际的圆形影响区来对结合部进行有限元建模。每个影响区共设置了16个板单元，其总的面积与前述圆形影响区面积计算值大体一致。用25个复弹簧单元模拟螺栓结合面的刚度及阻尼特性，且其中复弹簧单元的刚度值共划分为6个梯度，分别按正弦、抛物线和线性分布赋值。参照式(3)～式(5)，每个复弹簧单元的刚度值都可以通过螺栓影响区中心处最大刚度值kmax计算获得。

为了明确用哪种复弹簧单元分布形式模拟螺栓影响区力学特性更能提升分析模型的精度，这里进行了对比研究。同时，为了说明这种用分布弹簧单元模拟螺栓影响区建模理念的先进性，还对比了用单根以及均布弹簧单元模拟螺栓结合部力学特性的结果。分别采用反推辨识技术，利用实测的前5阶固有频率，针对不同的螺栓影响区建模方式，确定弹簧单元的刚度值，其中对前5阶分配的权重依次为0.1、0.1、0.3、0.3和0.2，遗传算法中种群数量、变异概率、交叉概率、迭代次数依次设置为50、0.05、0.9和50，辨识得到的螺栓影响区中心弹簧刚度值结果如表2所示。

表2 反推法辨识获得的模拟螺栓影响区中心处的弹簧单元刚度

需要说明的是表2中的弹簧刚度值是使对应的螺栓影响区模拟方法达到最高模拟精度的刚度值。用此刚度值获得的螺栓连接板的前5阶固有频率与试验值的比对分别如表3和表4所示。

表3 单根弹簧和弹簧刚度值均匀分布时模型固有频率与试验固有频率对比

表4 弹簧刚度值非均匀分布时模型固有频率与试验固有频率对比

从表3与表4的对比可非常直观地看出用单根弹簧或弹簧刚度值均匀分布模拟螺栓影响区时模型获得的固有频率与试验结果的接近程度明显差于弹簧刚度值按正弦、抛物线和线性分布时的结果。进一步，利用均方根误差(RMSE)方法对比上述三种非均匀分布模拟时的各阶固有频率与试验值的偏差，RMSE的求解式如下

(22)

将表4中的数据代入到式(22)，获得对应三种分布形式的仿真与试验前5阶固有频率的RMSE，结果如表5所示。

表5 各非均匀分布形式中仿真与试验前5阶固有频率的RMSE

通过表5中的RMSE对比可以看出，采用正弦分布时，模型前5阶固有频率与试验值的接近程度更好。故在后续分析时，对于本文的结构，螺栓影响区均采用刚度值按正弦分布的非均匀分布弹簧单元来模拟。

为了进一步描述清楚该建模方法的可用性，以下进一步讨论了螺栓影响区面积对分析结果的影响。将螺栓影响区增加至用25个板单元来模拟，相应的复弹簧单元数量增加到36个，并按正弦分布赋弹簧刚度值，求解该搭接板结构的固有频率。将获得的结果与上面用16个板单元模拟影响区面积的结果进行比对，如表6所示。从对比结果可以看出，变动影响区面积会使分析结果发生改变，但变动的幅度并不大。可见只要基于一定的理论确定螺栓影响区面积参考值，在实际建模中基于此参考值大致确定一个影响区面积作为添加复弹簧单元的依据则可获得具有一定精度的分析结果。

表6 影响区面积变动前后模型固有频率对比

3.3 振动特性求解

前一部分已描述了固有频率的求解，这里继续求解螺栓连接板的模态振型。仿真计算获得的模态振型与实测模态振型的对比如表7所示，可以看出两者的振型基本一致。

表7 试验与仿真前5阶振型对照

接下来进行频响函数的计算，计算之前需获得复弹簧单元的阻尼。前面已述，这里采用阻尼值均匀分布的复弹簧单元模拟螺栓影响区的阻尼特性。参照图6所描述的锤击点及拾振点，利用式(21)进行频响函数计算。首先，采用反推辨识技术，利用实测包含前5阶固有频率的频响函数曲线，依次分配的权重为0.4、0.1、0.3、0.1和0.1，遗传算法中种群数量、变异概率、交叉概率、迭代次数依次设置为50、0.05、0.9和50，不断迭代模型中弹簧阻尼单元的阻尼值，使频响函数曲线尽可能接近实测曲线，完成迭代后，螺栓影响区任一弹簧阻尼单元的阻尼值为c=3.56×107N/m。接着用此阻尼值获得最终的频响函数曲线并与实测值比对，如图7所示。从图中也可看出，仿真与实测的频响函数也有较好的接近。

图7 实测与仿真频响函数对比

综上，从固有频率、模态振型以及频响函数的比对中，可以得出结论，本文所提出的用非均匀分布弹簧单元模拟螺栓影响区进而实施有限元建模可实现较高的仿真计算精度。

4 结 论

本文采用非均匀分布复弹簧单元模拟螺栓结合部影响区，以螺栓连接薄板为对象，描述了相关的有限元建模及参数辨识方法，得出如下结论：

(1) 提出用刚度非均匀分布复弹簧单元模拟螺栓影响区压力非均匀分布的力学特性，并给出了详细的建模流程及方法。实践表明用本文假定的正弦、抛物线和线性等非均匀分布复弹簧单元模拟螺栓影响区并建立动力学模型，其建模精度远高于用单根及均布弹簧单元模拟螺栓影响区的模型。

(2) 提出的复弹簧单元同时包含了刚度及阻尼，这种建模方式既可以较为精确的模拟结合面的力学特性又可简化建模过程，通过所提出的基于遗传算法的反推辨识流程可较为精确的确定复弹簧单元的刚度及阻尼参数。

(3) 在有效引入非均匀分布复弹簧单元的基础上，建立了螺栓连接薄板结构的动力学有限元模型，计算得到的前5阶固有频率与试验测得的固有频率偏差均小于4%，求解得到的模态振型也与试验振型基本一致，仿真和实测的频响函数也有着很好的接近，从而证明了该建模方法的合理性。需要说明的是，由于模拟螺栓结合部的复弹簧单元的刚度及阻尼是常值，因而本文所创建的动力学有限元模型仅能模拟结构线性振动特性。无法再现由于激振力幅的增加而使螺栓连接结构出现的软式(或硬式)、谐波失真、跳跃等非线性振动现象。后续，可对复弹簧单元模型加以改进，以适应螺栓连接结构非线性振动研究。