基于模型修正的螺栓连接板结构动力学建模方法

郭 宁，黄 庆,2，曹善成，徐 超

(1. 西北工业大学航天学院，西安 710072；2. 湖北航天飞行器研究所，武汉 430040)

一般情况下，影响飞行器结构动力学建模准确性的一个主要因素是其连接部位的仿真预示[1]。螺栓连接板结构是飞行器结构中最常用最重要的一种可拆卸式连接结构形式。因此，准确建立螺栓连接板结构的动力学模型，是准确获取飞行器结构动力学特性和响应的前提之一。

目前，针对螺栓连接部件动力学建模的主要方法有：线性建模方法、非线性建模方法及薄层单元方法。线性方法是将实际的非线性连接结构近似成线性对象进行动力学分析，主要有三种方法：全固连建模方法[2]、节点约束方法[3]和用杆单元、弹簧单元和梁单元等模拟结构连接形式[4−5]。线性建模方法具有工程操作简单等优点，但是不能准确地表征结构连接部位的刚度和阻尼特性，导致动响应分析结果偏差较大。非线性方法是通过假设或采用试验途径来建立连接界面的唯象力学模型。目前应用较多的界面力学模型有：Iwan模型[6−7]、LuGre毛刷模型[8]和Valanis模型[9]等。上述非线性方法存在非线性机理描述有限、非线性模型参数获取难和计算耗费大等问题，还未有效地应用于工程结构分析中。薄层单元法最初是用于建立岩石与岩石界面之间的接触模型。目前已被广泛应用于连接结构建模中，Mayer和Gaul[10]的研究表明薄层单元能够有效地模拟螺栓连接接触面的力学特性。翟学等[11]研究了基于薄层单元法的螺栓连接结构建模和分析方法，并应用于涡扇发动机的建模。Iranzad和Ahmadian[12]使用实验数据(固有频率、频响函数)来确定薄层单元的材料本构参数。姚星宇等[13]给出了薄层单元法在航空发动机机匣螺栓连接结构的建模方法，并研究了薄层单元材料参数对螺栓连接刚度的影响规律。王攀和臧朝平[14]基于螺栓连接超模型的刚度理论、赫兹接触理论以及M-B分形模型，考虑了螺栓的数量和法兰边的接触因素，建立了基于分区域薄层单元的螺栓连接简化建模方法。姜东等[15− 16]建立了基于各向同性本构关系的薄层单元理论，提出了螺栓连接结构接触面不确定性参数的识别方法。从上述的研究可以看出，采用薄层单元法建立连接结构动力学模型时薄层单元的本构参数的确定至关重要。

本文结合薄层单元理论，针对螺栓连接板结构提出了基于模型修正的动力学建模方法，其中采用薄层单元的力学行为等效模拟螺栓连接接触界面的力学行为，建立用于模型修正的初始动力学有限元模型，将不确定的薄层单元材料本构参数作为修正变量，通过模态参数构建目标函数，提出采用水循环算法(Water Cycle Algorithm, WCA)来修正薄层单元的不确定本构材料参数。

1 基于薄层单元法的螺栓连接板结构动力学建模

图1所示为工程中典型的螺栓连接板结构，首先采用板(壳)单元对被连接板部件进行动力学建模；继而采用一定厚度的8节点六面体实体单元(即薄层单元)来等效模拟螺栓连接力学行为；最后在薄层单元与对应接触的壳单元之间建立绑定约束关系，建立绑定约束时，需要设置绑定容差，容差应略大于两个被连接板中面之间的距离，以保证薄层单元与板(壳)单元之间能够有效的约束；螺栓的质量采用集中质量单元进行模拟。具体如图1所示。

图 1 螺栓连接板结构动力学建模Fig.1 Dynamic modeling of bolted joint plate

根据上述建模过程可知，螺栓连接板结构的连接界面刚度是通过薄层单元的刚度来等效的，建模过程中需要确定薄层单元的厚度以及其材料本构。薄层单元的厚度对分析结果影响较大，厚度过大则单元会有6个应变分量，难以准确体现接触界面的力学特征；厚度过小则会导致单元雅可比矩阵行列式值趋向于零，造成矩阵病态，求逆困难，无法准确计算其位移-应变关系[15]。本文采用文献[17]中给出的薄层单元厚度计算方法，即：

式中：l1和l2分别为薄层单元的长度和宽度；t为薄层单元厚度；R的取值范围为10～100。

由于薄层单元的厚度远小于其长度和宽度，因此，可忽略薄层单元的面内应变分量( εx,εy,γxy)。并且假设连接界面的法向和切向接触性能相互独立，且2个切向的接触性能相同，则表征界面接触性能的薄层单元本构方程可表示为[15]：

式中，En和Gt分别是薄层单元的法向弹性常数和切向剪切常数。若接触面的法向与切向接触性能相互耦合，可在本构关系式(2)中通过引入耦合项来实现。

类比方程式(2)，在有限元建模中，可采用各向同性材料的本构关系近似替代薄层单元本构关系，具体为[15]：

式中：G为材料剪切模量；λ为拉梅常数。

式中：E为材料弹性模量；μ为材料泊松比。

在实际情况中，由于螺栓预应力的作用，螺栓连接界面的刚度分布并非一致，螺栓附近区域的界面连接刚度高于远离螺栓区域，因此，将薄层单元分为两部分，采用不同的材料本构参数进行建模。分区大小由螺栓预应力的作用区域大小决定，如图2所示，图中 α角的取值范围为25°～33°[18]，本文选择 α=30◦。为了进一步简化有限元建模，将圆形的螺栓预应力作用区域等效为正方形区域来进行有限元网格离散，但需要保证作用区域面积相同。

图 2 螺栓连接界面应力分布示意图Fig.2 Stress distribution at bolted joint interface

2 基于水循环算法的螺栓连接板结构动力学模型修正

根据第1节的分区薄层单元法建立的螺栓连接板结构动力学模型可知，需要确定的建模参数包括：螺栓预应力作用区域的薄层单元材料参数E1、μ1以及非作用区域的薄层单元材料参数E2、μ2。

2.1 目标函数的构建

本文利用固有频率的误差和模态置信准则残差的加权和来构建模型修正的目标函数：

式中：N表示实验测试模态阶数；p={p1,p2,···,pm}T为待修正变量向量，其中，下标m为修正变量个数；ωiE和ωiA分别为实验和仿真的第i阶固有频率；α为加权系数； MACi为第i阶模态的模态置信度。

2.2 水循环优化算法

从式(6)可以看出，当目标函数 Π(p)取得最小值时，即可得到模型修正问题的解，即薄层单元的待确定参数。因此需要将式(6)转化为如下的最小值问题：

本文选择水循环优化算法来执行式(9)所示的最小值优化问题。水循环优化算法是Eskandar等[19]于2012年提出的一种新颖的优化算法，其理论受启发于大自然，是通过观察大自然水循环过程中水的流动和转换而提出的。目前，水循环算法已经被广泛地应用于各类实际工程优化问题中，如Cho和Kim[20]将水循环算法应用于混合电磁铁的优化设计方法中，通过实验验证了该算法的性能优于常规优化算法；Ghaffarzadeh[21]提出采用水循环算法来开发电力系统稳定器，结果表明，采用水循环算法可以在较短时间内检测出最优参数，进一步提高了电力系统的稳定性。然而，据本文作者所知，水循环算法尚未被应用到解决结构动力学模型修正的问题中。

执行水循环算法，首先要随机生成一个初始降雨层(类似于遗传算法的种群或粒子群算法中的粒子群)，其包含的个体称为雨滴(Raindrop)，代表优化问题中的一组解(类似于遗传算法的染色体或者粒子群算法中的粒子)。初始降雨层中的每个雨滴Ri是在给定范围内随机产生的向量，对于本文中的m维修正变量问题，其表达式为：

式中： δi∈U(0,1) 为服从均匀分布的随机数；n为降雨层中的雨滴个数。

其次，将所有雨滴根据其目标函数值(即式(6))的大小划分为三个等级，其中最好的雨滴(即当前目标函数值最小的雨滴)形成海洋，较好的一些形成河流，其余的雨滴形成流向河流和海洋的溪流：

式中，下标NSr为河流数和海洋数之和，且海洋数为1。

假定所有溪流都会汇入河流或者海洋，由于河流和海洋的流量不同，流入每条河流和海洋的溪流数量NS就不同，具体为：

式中，NStream为溪流的数目。

然后，水循环算法进入汇流阶段。其中：一部分溪流先汇入河流，再汇入海洋；另一部分溪流则直接汇入海洋。具体汇流的表达式如下：

为了防止算法陷入局部最优，通过蒸发和降雨来扩大算法的搜索空间。降雨过程可根据降雨位置的不同分为两种：一种是在溪流或河流附近降雨，即随机产生溪流或河流的新个体，增大算法的搜索域；另一种是在海洋附近降雨，即在海洋附近随机产生新个体，使算法在当前最优解附近寻找更优解。可以根据以下两个条件来判断是否满足降雨条件：

式中，dmax为人为设置的阈值，通常在算法初始阶段取dmax=10−16，且随着迭代地进行自适应减小，具体为：

式中，M为最大迭代次数。

如果满足式(16)，则在整个搜索空间内进行降雨，具体为：

如果满足公式(17)，则在海洋附近区域进行降雨，具体为：

式中：λ为服从均匀分布的随机数向量；ν为在海洋附近的搜索范围的大小，通常取ν=0.1。

同其他启发式智能优化算法一样，本文选用最大迭代次数作为水循环优化算法的终止条件。

综上，结合薄层单元法和模型修正方法的螺栓连接板结构动力学建模流程如图3所示，且具体实施步骤如下：

第一步，建立部件结构的动力学模型，并对其建模参数进行修正；

第二步，采用薄层单元法建立螺栓连接板结构的初始有限元模型，并进行模态分析，提取仿真模态参数；进行螺栓连接板结构的实验模态分析，提取相应的实验模态参数；

第三步，利用实验模态参数与对应仿真模态参数，构造模型修正目标函数；

第四步，选择修正变量，给定修正范围；

第五步，设置水循环算法的初始参数，如初始雨滴数目n，海洋加河流的数目NSr，溪流数目参数NStream，蒸发条件判断数dmax，最大迭代次数M等；

第六步，产生初始种群，根据目标函数值划分海洋、河流和溪流，计算汇流数目，进行相应的汇流；

第七步，判断是否满足蒸发条件，如果满足，则执行降雨过程；

第八步，判断是否满足算法终止条件，如果满足，则结束算法并输出最优解；如果不满足，则返回到汇流过程，继续执行算法，直到满足终止条件为止；

第九步，输出算法最优解，即为建模所需的薄层单元的材料参数。

图 3 基于模型修正的螺栓连接板结构动力学建模流程图Fig.3 Flow chart of dynamic modeling of bolted plate structure based on model updating

3 实验验证

为了验证本文所提的螺栓连接板结构动力学建模方法的可行性和实用性，设计了如图4所示的两端固支五螺栓搭接板模型，其中，单板的尺寸为240 mm×120 mm×3 mm(去除边界加持区域)；搭接区域长度为40 mm。搭接板材料为铝合金，采用5个M6不锈钢螺栓进行连接，且拧紧力矩为5.2 N∙m。

图 4 五螺栓搭接板结构Fig.4 Five-bolt lap plate

实验结构与测试设备如图5所示，包括信号采集及模态识别系统(LMS SCADAS Mobile)、PCB力锤(灵敏度为2.267 mV/N)、PCB加速度传感器(灵敏度为100.5 mV/g)和STANLEY扭矩扳手。实验中采用“移动激励点法”进行模态测试。为了抑制由于信号截断引起的泄漏，实验测试中对输入信号(即激励信号)添加力窗，并对输出信号(即响应信号)添加指数窗。为了一定程度上消除测试噪声的影响，实验中在每个测点进行3次测量，并对3次测试响应数据进行平均。本实验所关心的结构的频带范围为1 Hz～620 Hz，包含了结构前6阶模态，实验中设定振动信号的采样频率为1280 Hz(大于最大测试频率的2倍)，频率分辨率为0.31 Hz。实验获得的结构前六阶模态结果见表1。

图 5 模态试验现场布置图Fig.5 Layout of modal test

首先对部件结构(即板1和板2)的模型参数进行修正。通过对实验件进行测量，得到板1和板2的尺寸分别为240 mm×120 mm×2.96 mm和240 mm×120 mm×2.95 mm，重量分别为284 g和282 g，可计算得到结构的材料密度为2717.7×10−12t/mm3。所以仅需要对其材料弹性模量E和泊松比μ进行修正。选取板1为修正结构，进行模态实验，边界条件为自由-自由，且采用4边形壳单元对单板结构进行有限元建模，修正后的材料参数见表2。

表 1 五螺栓搭接板前六阶实验模态Table 1 First six experimental modes of five-bolt lap plate

根据第2节中所述的螺栓连接结构的动力学建模方法，建立了如图6(a)所示的有限元模型，薄层单元采用六面体8节点实体单元，其中白色单元区域为螺栓预应力作用区域，红色单元区域为非螺栓预应力作用区域，如图6(b)所示。对结构进行有限元离散时选择的最大单元边长为4.5 mm，根据式(1)，我们选取得薄层单元厚度为0.1 mm。

图 6 五螺栓搭接板有限元模型Fig.6 Finite element model of five-bolt lap plate

定义设计变量集p={E1,μ1,E2,μ2}，每个设计变量相对于固有频率的灵敏度如图7所示。

从图7可以看出，泊松比μ1和μ2对低阶模态不敏感，因此，选择敏感性较强的设计变量E1和E2作为模型修正的修正变量，并分别给定修正区间为E1∈[0.00110] 和E2∈[0.00110]；建模中设定μ1=μ2=0.3。

采用第2节建立的模型修正方法对修正变量进行修正，其中水循环算法的具体参数设定见表3。

图 7 设计变量对于固有频率灵敏度分析Fig.7 Sensitivity of natural frequencies versus design parameters

表 3 水循环算法的参数设定Table 3 Water cycle algorithm parameters

将修正后的参数输入有限元模型(即E1=0.001 MPa，E2=0.2957 MPa)，进行模态分析，并将参数修正后的模态分析结果与实验结果进行对比，详见表4。

表 4 修正后模态分析结果Table 4 Modal analysis results after updating

从表4可以看出，采用本文提出的动力学建模方法建立的螺栓连接板结构的动力学模型可以很好的预示结构的动力学特性，固有频率的最大相对误差仅为1.69%，且MAC均大于0.91，满足工程中对结构动力学建模的需求。

3.1 优化算法对比

本文分别采用遗传算法[22](Genetic Algorithm,GA)、粒子群算法[23](Particle Swarm Optimization,PSO)和水循环算法(Water Cycle Algorithm, WCA)三种全局智能优化算法对算例中的薄层单元的参数进行修正。为了使修正结果具有可比性，将算法的种群数、粒子群数和雨滴数均设置为30，迭代次数均设置为40次。图8为分别采用三种优化算法执行模型修正问题的收敛情况对比结果。从图8中可以看出，相比遗传算法和粒子群算法，水循环算法的收敛曲线下降最快，目标函数值最小，表明水循环算法在求解结构动力学模型修正问题时具有收敛快、精度高的优点。

图 8 不同智能优化算法优化结果对比图Fig.8 Comparison of optimization results of different intelligent optimization algorithms

为了进一步验证水循环算法在求解最优化问题(式(9))时具有唯一性，本文在保证相同优化参数设定的前提下，重复运行求解过程10次，均得到了相同的最优解和最佳适应度值，即，E1=0.001 MPa；E2=0.2957 MPa；∏(p)=0.006 561 7，说明了采用水循环算法在求解结构动力学优化问题时具有很好的鲁棒性。

3.2 抗噪性验证

本节对所提动力学建模方法的抗噪性进行验证。通过对实验模态数据添加人为噪声来模拟带有测试噪声的模态参数识别结果，即：

式中：XE代表固有频率向量或第i阶模态振型向量；η为噪声水平；δ为均值为0、方差为1，并满足标准正态分布的随机向量。

众所周知，测试噪声对识别结构的固有频率和模态振型的影响是不相同的，固有频率一般可以比较精确地识别，而模态振型的识别精度比固有频率的精度要差很多。因此，选取了两种噪声水平工况来考察本文所提出的螺栓连接板结构动力学建模方法的抗噪性，具体为：工况1(固有频率( η=1% )，模态振型( η=5%)；工况2(固有频率( η=2% )，模态振型( η=10%)。表5和表6分别给出了不同噪声水平下的固有频率和MAC值对比结果。图9示出了不同噪声水平下修正后获得的MAC值对比情况。

表 5 不同噪声水平下的固有频率对比Table 5 Comparison of modal frequencies under different noise levels

表 6 不同噪声水平下的MAC值对比Table 6 Comparison of MAC under different noise levels

图 9 不同误差水平下修正后的MAC值对比Fig.9 Comparison of MAC under different noise levels

综合上述对比结果可以看出，在第一种噪声水平下，固有频率的最大相对误差仅为2.49%，MAC值均大于0.9，其最小值为0.907。在第二种噪声水平下，固有频率的最大相对误差为4.61%，MAC值均大于0.9，其最小值为0.911。上述结果表明，在两种噪声工况下建立的螺栓连接板结构的动力学有限元模型的动力学特性都和真实结构的测试结果吻合得很好，表明本文所建立的动力学建模方法具有良好的抗噪性。

4 结论

本文针对螺栓连接板的动力学建模问题，提出了一种结合薄层单元法和模型修正的建模方法。采用薄层单元法建立了用于模型修正的初始动力学模型，将薄层单元的不确定材料本构参数作为修正变量，引入水循环算法来进行优化求解。从本文的研究中可得出以下结论：

(1)相对于遗传算法和粒子群算法，采用水循环算法进行动力学模型修正，极大的提高了求解效率。

(2)采用两端固支五螺栓搭接板结构作为实验验证模型，验证了提出的螺栓连接板结构动力学建模方法是合理可行的，建立的动力学预示模型可以更好地用于进行结构的动力学分析，并且具有良好的抗噪性。