纤维增强复合材料是由纤维增强材料与基体构成^[1].与金属材料相比，具有比强度和比刚度高，力学性能可设计等优点，广泛应用在航空航天、汽车和海洋工程等领域^[2].目前关于复合材料的联接方式主要包括胶接联接、机械联接及两者结合的方式，其中，螺栓联接具有便于拆装和维护的优点而得到最为广泛的应用^[3].复合材料联接结构在使用过程中容易产生振动而诱发不稳定性，包括由共振引起的剧烈振动^[4]、界面磨损、纤维和树脂的断裂及分层等^[5].因此，有必要对螺栓联接复合材料结构的振动特性进行相关的建模分析.

国内外对螺栓联接复合材料结构的振动特性建模进行了少量研究.Zhang等^[6]使用虚拟材料模拟结合部的接触界面进行建模，利用分形接触理论推导了虚拟材料的复模量，最后在有限元软件上采用多尺度法研究接触界面粗糙度对结构的共振频率和阻尼的影响.郝秉磊等^[7]用接触单元模拟C/SiC陶瓷基复合材料螺栓结合部的力学行为，进而计算该联接件的固有频率和振动响应.Abdus等^[8]考虑不同预紧力和温度对固有频率的影响，使用有限元软件模拟结合部.而关于螺栓联接复合材料结构的解析分析主要针对静力学，例如，Olmedo等^[9]提出了一种解析方法用于预测单螺栓搭接复合材料板的二次弯曲，并用实验进行了验证.Xiang等^[10]开发了一种考虑孔间隙和摩擦效应的改进弹簧方法，研究了复合板宽、板厚、螺栓间距和间隙对载荷分布的影响.

由以上可知，对螺栓联接复合材料的固有特性解析建模研究较少，而固有特性包括固有频率和振型，是预测结构松动和局部损伤的关键振动特性参数.因此，本文基于虚拟材料的相关理论建立了螺栓联接复合材料板的解析模型.这里虚拟材料建模法^[11-12]基本思想是:用一层虚拟材料代替真实的螺栓结合部，通过对该材料的弹性模量、泊松比、损耗因子、密度等赋予不同的参数来模拟螺栓结合部的刚度、质量及阻尼等性能.提出使用基于分解的多目标进化算法^[13]来辨识该虚拟材料模型的参数.最后，以螺栓联接TC500碳纤维/树脂基复合材料板为例，创建了分析模型并用实验验证了该建模方法的有效性.

1 基于虚拟材料的螺栓联接复合板解析建模 1.1 本构方程的建立

螺栓结合面的接触实际上是由无数个微凸体的接触构成，如图 1所示.将微凸体的接触等效成一层具有一定厚度的各向同性虚拟材料，通过改变虚拟材料的弹性模量、泊松比、密度和厚度来模拟螺栓结合部的法向和切向动态特性.

基于虚拟材料的螺栓联接复合材料板结构如图 2所示，其中上下两块复合材料板的材料是碳纤维/树脂基复合材料反对称正交铺设，长度为a，宽度为b，厚度为H_c.其纤维方向与整体坐标系x方向夹角为θ，纤维纵向弹性模量和泊松比分别为E₁, v₁₂，垂直于纤维方向的弹性模量和泊松比分别为E₁, v₂₁，面内剪切模量为G₁₂，密度为ρ_c.

将螺栓结合部假定为虚拟材料，位于两个复合板之间，该虚拟材料的长度为h，宽度b，厚度为H_s，弹性模量和泊松比分别为E_s, v_s，密度为ρ_s.其中复合材料板采用经典的正交各向异性复合材料板理论进行分析，而虚拟材料采用各向同性板理论进行分析.如图 3所示，复合材料层合板任意层的纤维方向与x轴角度为θ，从下到上层外表面的厚度方向坐标用Z₁, Z₂, …，Z_N+1表示，这里N是层数.

对于正交各向异性的单层板平面应力问题(σ_z=τ_xz=τ_yz=0)，取图 3所示xOy自然坐标系和x₁Oy₁弹性主轴坐标系.在x₁Oy₁弹性主轴坐标系中，每一层本构方程为

(1)

式中:σ₁, σ₂, τ₁₂, ε₁, ε₂, γ₁₂分别为纤维方向、纤维横向的应力和应变；Q_jl, j, l=1, 2, 6表示与弹性常数相关的系数.

利用转轴公式，在xOy参考坐标系中，纤维每层中应力和应变的本构关系为

(2)

式中:σ_x, σ_y, τ_xy, ε_x, ε_y, γ_xy为层合板第i层在xy平面的应力和应变；Q_jl, j, l=1, 2, 6具体表达如下:

(3)

1.2 螺栓联接板结构的能量分析

整体结构的应变能为

(4)

式中:U_c1，U_c2和U_s分别为上下两块复合材料板及虚拟材料的应变能:

(5)

(6)

这里假设w₀(x, y, t)=W(x, y)sinω₀t为复合材料板中性面自由振动的挠度，W(x, y)为位移振型，ω₀为固有频率.

由于上下两块复合材料板的材料和尺寸完全相同，因此其弯曲刚度系数D_jl相同，可表示为

(7)

虚拟材料的应变能U_s的表达式为

(8)

式中:

整体结构的动能T表达式为

(9)

式中:T_c1，T_c2和T_s分别表示上下两块复合材料板的动能和虚拟材料的动能.

1.3 位移振型的假设

假设螺栓联接复合材料板的横向振动幅值是由正交多项式表达的二重级数P_r(α)和P_s(β)的连续函数的形式:

(10)

式中:α=x/a；β=x/b；R和S为多项式的数量，通常取R=S；c_rs为待定参数；P_r(α)和P_s(β)为特征正交多项式，对于悬臂边界状态，其具体表达式为

(11)

式中:

(12)

这样，便得到了关于α和β的无量纲的W(α, β).将α=x/a，β=x/b再代回式(10)，则横向振动幅值W重新变成关于x，y的W(x, y)函数.

1.4 基于拉格朗日方程的固有特性求解

将动能和应变能代入拉格朗日方程中，设L=T-U，拉格朗日方程的具体表达式为

(13)

式中：p=1, 2, 3, …, R; k=1, 2, 3, …, S.

经过一系列求导等计算，最终得到系统的自由振动方程:

(14)

式中: K, M, c分别为螺栓联接系统刚度矩阵、质量矩阵和特征向量.计算结构的固有频率需要式(14)有解，因此其特征向量c对应的系数矩阵的行列式为0，即

(15)

对K和M求广义特征值就能得到固有频率的平方及特征向量，然后将特征向量代入W(x, y)中，再将求解的W(x, y)代入w₀(x, y, t)中，可得到各阶次模态振型.

2 虚拟材料参数的确定

在初步计算时可仅靠复合材料构件的参数来确定虚拟材料的参数，例如设虚拟材料弹性模量E_s=E₂/2(1-v₁₂²)，泊松比v_s=v₁₂，厚度H_s=H_c/N.为了粗略保证系统总体质量不变，设定虚拟材料的密度ρ_s与螺栓密度相同.

如果要用虚拟材料模型精确模拟复合材料螺栓结合部，可通过模型修正技术反推确定虚拟材料参数.使用基于分解的多目标进化算法(MOEA/D)来辨识虚拟材料的参数，设定修正的多个目标函数为

(16)

式中:f_i(x)(i=1, 2, …, m)表示第i个目标函数，x为决策变量；Ω表示n维决策空间；f_i^S和f_i^E分别表示仿真和实验获得的第i阶固有频率.

由式(16)可知，需要同时优化多个目标函数，采用MOEA/D算法能够根据相邻子目标函数的信息产生相应的权重向量，权重向量在迭代过程中发生变化，从而完成优化过程.多目标优化问题的优化过程中会产生一个Pareto最优解的集合PS.所有PS中的解对应的目标函数值所形成的集合PF称为Pareto前沿.利用MOEA/D算法确定虚拟材料参数的过程描述如下:

1) 输入:虚拟材料的参数种群大小为P，算法的进化代数为Maxgen，Tchebycheff分解的邻居大小为T，P个均匀分布的权重向量[λ¹，λ²，…，λ^P]，求解式(16)的运动方程并形成的目标函数.

2) 将虚拟材料的参数E_s, v_s, H_s, ρ_s种群初始化，之后将进化过程中的非支配个体储存在外部群体(EP)中.

3) 计算任意两个权向量之间的欧式距离，然后计算每个与欧式距离最近的T个权向量.

4) 随机产生一个虚拟材料参数的初始种群x¹, …, x^N，并且调用求解式(16)的运动方程计算目标函数值FVⁱ=F(xⁱ).

5) 繁殖:对种群个体使用遗传算子操作，之后计算目标函数值产生一个新的解y.

6) 新解的改进:通过具体的问题对新解y进行罚函数越界判断后变成y′.

7) 更新外部群体集合(EP):将所有F(y′)支配的向量从EP中移除，如果F(y′)没有被EP中其他向量支配，就将F(y′)加入到EP中.

8) 停止条件:当达到算法的最大进化代数时，停止并输出EP.否则，转向步骤5).

3 实例研究

这里以螺栓联接两块TC500碳纤维/树脂基复合材料板为例进行分析，复合材料板的长度a=0.15 m，宽度b=0.06 m，厚度H_c=0.002 36 m，该类型纤维增强复合板为正交铺设，整个复合板层数为20，每个铺层具有相同的厚度和纤维体积分数，纤维方向角为[-45°/+45°]₁₀，复合材料TC500的纤维纵向弹性模量E₁=136 GPa，横向弹性模E₂=7.92 GPa，泊松比为0.32，剪切模量G₁₂=3.39 GPa，密度ρ_c=1 780 kg/m³，结合部区域的长度l_c=0.03 m.

螺栓组合复合板的一端固定在夹具上(夹持区为30 mm)，通过力矩扳手施加2 N · m预紧力，利用锤击法对组合板的固有特性进行测试，所组建的实验系统见图 4.这里用Polytec PDV-100激光多普勒测振仪拾振，用PCB SN 30272力锤激励试件，用LMS SCADAS进行数据采集.

3.1 螺栓联接复合材料板的固有特性求解

通过MATLAB软件编写相应的程序，取多项式的项数R=S=12，求解得到的结构前6阶固有频率和模态阵型见表 1和表 2.

由表 1和表 2可知，解析模型中求得的固有频率与实验对比最大偏差为10.50 %，最小偏差为0.99 %, 偶数阶次的频率偏差与奇数阶次相比较要大.推断偏差较大的原因:复合材料的力学性质是通过制造商机械测试提取出来的，其同一批次在生产加工过程中可能存在一定偏差；模型中虚拟材料参数不准确导致的结构偏差.因此有必要采用优化算法对复合材料及模拟螺栓联接的虚拟材料的参数进行辨识从而达到模型修正的目的.

3.2 模型修正

首先对单个悬臂板进行模态分析以确定较为准确的复合材料参数.采用解析法得到其对应的固有频率，与实验进行对比，其中决策变量选取E₁，E₂和G₁₂，变量上下界为0.8~1.1倍初始值；目标函数选取前三阶模拟值与实验值偏差的平方.采用基于分解的多目标进化算法(MOEA/D)修正复合材料参数得到的计算结果，如表 3所示.

对复合材料参数进行修正后，使用基于分解的多目标进化算法(MOEA/D)对虚拟材料模型进行修正，目标函数选取偏差较大的第1，2，4阶固有频率与实验值偏差的平方.采用切比雪夫方法对多个目标函数进行分解，设置决策变量种群的大小P=200，邻居大小T=20，最大进化代数Maxgen=200，目标函数的Pareto前沿如图 5所示.

由图 5可知，模型修正后的Pareto前沿是凹的并且解集分散的均匀性较好，当目标函数f₁较小时，f₂与f₃中总有一个偏大，即一个目标函数的优化性能会牺牲其他目标函数的优化性能.对PF_known进行分析并选取最满意的一组非支配解作为虚拟材料的参数，得到多目标算法修正后的固有频率，如表 4所示.

通过对比表 1可以看出，修正后模型的1, 2, 4阶次固有频率与实验值的偏差明显减少，并且其他阶次的偏差也相应减少.表 5为最终确定的模拟复合材料螺栓结合部的虚拟材料参数.

4 结论

1) 将螺栓结合部的力学特性用虚拟材料来模拟，进一步采用能量法创建了螺栓联接复合材料板的解析分析模型.实践表明这种不考虑螺栓结合部细节创建组合结构分析模型的方法可显著提高螺栓联接结构的分析效率.初步计算获得的前6阶固有频率与实测值最大偏差为10.5 %，而模态振型与实测值基本一致.

2) 为了提高所创建的螺栓联接复合材料板动力学模型的分析精度，可以采用基于分解的多目标进化算法(MOEA/D)来修正用于模拟螺栓结合部的虚拟材料参数，该优化方法具有收敛快、精度高的特点，计算获得的螺栓联接板前6阶固有频率与实测值最大偏差仅为2.83 %.