一种基于谱元法的层合板SEA参数计算方法

袁金龙，陈海波，王 昆，操小龙

(1.中国科学院 材料力学行为和设计重点实验室，中国科学技术大学近代力学系，合肥 230027；2.北京机电工程研究所，北京 100074)

复合材料具有强度高、重量轻、抗疲劳等特点，已经广泛应用于航空航天等领域。飞行器在飞行过程中承受气动热、宽频噪声、振动等复杂恶劣环境，这些环境因素对其性能和可靠性有着重要的影响。而针对构件高频域的动态响应，有限元和边界元需要用很小的网格尺寸来描述结构特征长度，使得数值分析的计算量和存储量变得非常巨大[1]，甚至不能实现。统计能量分析方法将复杂系统划分为不同的模态群，并从统计意义上把大系统分解为若干个便于分析的独立子系统，对共振模态响应进行平均，从而得到光滑准确的能量值[2]。

统计能量分析方法的预示精度在很大程度上取决于统计能量分析参数的准确性。对于三明治板这类夹芯层合板结构的研究，国内外学者已经做了很多工作。 Renji等[3]推导了考虑芯层横向剪切刚度影响的夹层板四阶振动控制方程，并由控制方程在波数空间积分得到了正交各向异性复合材料夹层板的模态密度公式，最后通过点导纳方法测得的模态密度试验值对理论公式进行了验证。Zhou等[4]由夹层板六阶振动控制方程推导了模态密度的计算公式，计算了蜂窝夹层板的传声损失。然而该研究只针对各向同性面板的夹层板，忽略了芯层面内刚度对夹层板模态密度的影响。Han等[5]考虑面板弯曲刚度和芯层面内刚度, 提出了一种改进的经典夹层板理论,能更准确地计算三明治板的模态密度。

对于一般的层合板结构，三明治板理论不再适用，这时需要采用层合板理论或数值计算方法。经典层合板理论(CLPT)是基于Kirchhoff薄板假设,认为板变形后的中面法线依然垂直于板中面，适用于薄板结构低频分析[6]；一阶剪切板理论(FSDT)考虑横向剪切和转动惯量影响，认为板变形后其中面法线仍为直线但不与中面垂直，适用于更高频率的分析[7]；而更高阶的剪切板理论控制方程复杂，波数求解困难。Chronopoulos等[8]利用波有限元方法计算了正交各向异性复合材料层合板的频散特性，并利用波数和群速度数值求解了层合板的模态密度，但是计算量巨大。

谱元法(SFEM)是一种典型的半解析方法，是在有限元方法的基础上提出的，它综合了有限元方法处理复杂结构问题的灵活性以及收敛速度快和计算精度高的特点，对结构模拟波传播特性十分有效，是提取色散曲线的一种重要方法。谱元法假设平面简谐波沿着波导轴方向传播，并基于有限元方法对波导截面内的位移场进行离散，利用类似有限元的方法进行单元刚度矩阵和单元质量矩阵的组装，从而获得一个标准的波数特征值方程。Zak等[9]将裂纹近似为刚度线性变化的弹簧，采用二维谱元法研究了含裂纹的各向同性板中的弹性波传播特性，并进而提出了一种损伤识别的方法。Peng等[10]采用三维谱元法模拟弹性波在各向同性板中的传播特性，并与采用二维谱元法的模拟结果进行对比，通过控制加载方式得到了兰姆波的不同模态。Kudela 等[11]基于Mindlin 板理论采用二维谱单元法研究了弹性波在复合材料层合板中的传播，并研究了铺层角以及纤维体积含量对层合板弹性波传播特性的影响。Datta等[12]采用谱元法研究了弹性波在对称交错层叠层合板结构的传播特性。Shorter[13]采用一维线性谱单元对线弹性各向同性层合板进行波导分析,推导了层合板总内损耗因子的计算公式。徐超等[14]推导了一种任意形状的时域谱单元，研究了高频导波在二维功能梯度材料结构中的传播特性。

与传统的有限元方法相比，谱元方法只需要对波导截面位移场进行离散，降低了模型的自由度，从而极大地提高了计算效率。而且由于在波的传播方向采用精确解，从而避免了有限元离散时产生的数值误差，进一步提高了计算结果的准确性。但将该法应用于统计能量参数获取的研究还未见报道。本文采用二次谱单元对一般层合板厚度方向进行离散，建立波数的特征根方程，并在波数求解时采用皮尔逊相关系数对其对应的模态进行分类，进而计算一般层合板的统计能量分析参数。通过对薄板、三明治夹芯板以及厚板的数值模拟分析，验证了该方法的可行性和有效性，这对一般层合板的统计能量分析具有工程指导意义。

1 谱元法(SFEM)

1.1 单元位移场

平面波在平板中传播时，高频段波传播对板面内尺寸如长、宽不敏感，因此在厚度方向可采用一维单元离散,由波的传播特性，单元位移场可表示为：

(1)

式中:ω为圆频率，k为波数，φ为方位角，u(z)、v(z)和w(z)为厚度方向z坐标处波传播的振幅。长度为Li的单元中，位移振幅通过一维形函数插值可表示为：

(2)

图1 一维网格

p=u,v,w

(3)

N(z)为单元的形函数，文献[13]采用的是二节点单元，为了进一步提高计算精度和效率，本文采用三节点二次单元：

(4)

1.2 单元质量矩阵与刚度矩阵

单元的应变矩阵为：

(5)

其中：

(6)

将式(1)和式(2)代入式(5)中得

(7)

其中：

B1=LzN,z,B2=cos(φ)LxN+cos(φ)LyN

(8)

N,z为形函数矩阵对坐标z的偏导，单元内的应力为：

(9)

式中:C为模量矩阵，S为柔量矩阵，对于正交各向异性材料的独立弹性常数为9个，即三个拉压弹性模量、三个剪切模量和三个主泊松比：

在层合板坐标系下，模量矩阵为：

Cx-y=T(θ)-1CRT(θ)R-1

(10)

式中：m=cos(θ)，n=sin(θ)，θ为材料主轴方向与参考方向的夹角。层合板动能为：

(11)

层合板势能为：

(12)

H代表共轭转置。将式(7)代入式(12)可得：

U=

(13)

由哈密顿变分原理可得

(14)

其中：

将所有单元的刚度矩阵和质量矩阵组装起来即可获得总体的刚度矩阵和质量矩阵：

[K1(φ)+K2(φ)k+K3(φ)k2-ω2M]q0=0

(15)

2 统计能量分析

2.1 模态密度求解

当方位角和频率给定时，方程式(15)就变成k的二次特征根问题，二次特征根问题又可以进一步转化成下列线性特征值问题[13]：

(16)

其中：

q1=kq0;C2=K3;

C1=K2;C0=K1-ω2M

(17)

特征值的实部和虚部的符号代表给定波的传播方向，特征向量代表给定频率和方位角的横截面波的形状。在低频段，振动波主要为弯曲波、剪切波和拉伸波，当频率继续升高时将出现高阶振动波，因此需要对波的传播类型进行分类。层合板理论和三明治板理论都是通过求解高次方程来获得波数，存在波数分类困难的问题。谱元法在求解特征值的同时得到了特征向量，从而可用它来对波数进行分类。文献[8]采用特征向量余弦相似度的平方进行波数分类。余弦相似度指通过计算两个向量的夹角余弦值来评估它们的相似度，它适用于数据稀疏的情况，但易受向量平移的影响。而皮尔逊相关系数是利用两个变量之间的协方差和标准差的商来衡量两个变量的相关性，相比于余弦相似度，它对数据进行了中心化，具有平移不变性和尺度不变性的特点，且不受维度缺失的影响[15]。

根据模态相似原则，本文采用皮尔逊相关系数的绝对值来进行振型分类：

(18)

式中:Qm和Qn为相邻方位角或者相邻频率点的特征向量，E表示取均值。对振型相似的模态，P的绝对值接近1；在实际分类时，先取较大值，再取较小值，直到筛选出最相似的模态。利用方程(16)和(18)即可获得各阶模态对应的波数。

根据统计能量分析方法中模态密度的定义，在波数空间中，圆频率ω下的模态数N(ω)可由圆频率ω曲线包含的面积获得：

(19)

式中:a和b为矩形板的长和宽。由统计能量法模态密度的定义有：

(20)

2.2 结构与声腔的耦合

耦合损耗因子是子系统间耦合作用大小的一种度量。结构与声腔相互作用时结构的损耗因子ηrad与辐射阻抗Rrad有关：

(21)

式中:Ms为板结构的总质量,对于四边简支板的辐射阻抗可由下式获得[15]:

(22)

式中:As为板结构的面积，fc为板结构的临界频率，c为空气中声音的传播速度，ρa为空气密度。P为板边缘总长度。

2.3 内损耗因子

子系统的动力学特性与系统的刚度、质量和阻尼有关，阻尼在频响曲线共振区起主导作用，因此阻尼特性对研究结构的动力学响应和声传递特性具有重要意义。在统计能量分析方法中，一般用内损耗因子来表示系统的阻尼损耗特性。内损因子可定义为：

(23)

式中:Pdiss为时间平均的功率耗散。总功率损耗因子可以由各层的应变能获得：

(24)

式中:n是层合板的总层数，将式(24)代入(23)得：

(25)

式中:Ki是层合板第i层的组装刚度，即第i层对总体刚度矩阵的贡献量。考虑到高频振动时结构动能与势能的互等关系，将第k阶传播波的振型向量代入式(5)，该阶波的内损耗因子可表示为：

(26)

式中:Pk是第k阶传播波的特征向量。由式(26)可知，总内损耗因子与各层的势能有关，单层的势能越大，其对总内损耗因子的贡献越大。

2.4 传声损失

(27)

图2 传声模型

其中

(28)

式中：Vi表示声源室和接收室的体积，ηi为内损耗因子，声场的内损耗因子可以由混响时间T3获得，ni为模态密度，声场的模态密度可由下式获得：

(29)

式中:Ai为声场的总表面积，Si声场的边缘的总长度。

3 数值算例

3.1 各向同性薄板

对于厚度为h的各向同性薄板，其弯曲、剪切和拉伸的波数理论值可由下式给出：

(30)

以厚度为1 mm的均质铝板A为研究对象，铝板的密度为2 700 kg/m3,弹性模量为71 GPa，泊松比为0.329 6。图3为采用不同谱单元SFEM计算结果与式(30)理论值相比的误差。由图3可知，二次单元只需一个单元就能获得较准确的结果，而线性单元用了4个还未达到相应的精度，由此可见二次单元比线性单元的效率更高。图4为分别采用式(30)与采用二次单元的SFEM方法的波数计算结果对比，虽然SFEM只用了1个单元，二者几乎一致。本文后面的算例均采用二次单元计算。由图5可知，对于各向同性材料，波数的分布与方位角无关。

图3 板A弯曲波数SFEM与理论值误差(方位角为0°)

图4 板A波数(方位角为0°)

图5 板A波数与方位角的关系(SFEM)

3.2 正交各向异性薄板

正交各向异性薄板B为碳-碳(Carbon)复合材料板，材料参数如表1所示。由图6可知，当厚度比较小时，SFEM的波数计算结果与层合板理论CLPT和FSDT几乎一致。对于正交各向异性材料，波数的分布与方位角有关，如图7所示，此时板B在靠近x轴时波数较大。

表1 材料力学参数

图6 板B波数(方位角为0°)

图7 板B波数与方位角的关系(SFEM)

至此，通过各向同性板和正交各向异性板的波数计算结果对比，验证了SFEM计算波数的有效性。

3.4 三明治板

三明治板是一类特殊的层合板结构，一般由两块较薄、模量较大的面板和轻质、模量较小的厚夹芯层构成。以文献[4]中的蜂窝板C为研究对象，板C为夹芯层为6.35 mm的蜂窝板结构，面板和夹芯材料的力学性质，如表1所示。

采用SFEM分析板C时，面板划分1个单元，夹芯层划分5个单元，其他参数与文献[4]一致。计算出波数后，模态密度可由式(20)计算。如图8所示，在较低频段,SFEM与文献[4]和文献[5]的三明治板理论计算结果差距较小，都比实验值略小；在较高频段，SFEM计算结果相比于文献[4]和文献[5]与实验值吻合更好。由SFEM获得的临界频率为1 680 Hz，通过式(22)计算辐射阻抗如图9所示，相比于文献[4]中的数值解，由SFEM计算的辐射阻抗与实验值一致性更好。如图10所示，由SFEM计算的传声损失在低频段和高频段与实验值都比较接近，在临界频率处二者有一些偏差，总体结果要比文献[4]中的好。通过三明治板的模态密、辐射效率和传声损失的计算结果对比，进一步验证了本文所提SFEM计算SEA参数的准确性。

图8 板C的模态密度

图9 板C的辐射阻抗

图10 板C的传声损失

3.3 一般层合板

层合板D、E和F为五层厚度相同的层合板结构，长和宽均为1 m，每层的材料均为Carbon，铺层角度为[0/90/0/90/0]。单层的厚度分别为0.1 mm、0.5 mm和1 mm。分别采用CLPT、FSDT和SFEM计算层合板的模态密度，如图11～13所示，CLPT和FSDT相对于SFEM的偏差如图14所示。CLPT随着厚度和频率升高，模态密度计算偏差越大，FSDT在厚度小时与SFEM吻合得较好，当层合板厚度变大时，偏差增大。CLPT、FSDT和SFEM的计算时间对比如表2所示，FSDT的计算时间比CLPT略大，而SFEM的计算时间比CLPT多9倍左右。影响SFEM计算效率的主要因素是总体刚度矩阵和质量矩阵的维度以及波形分类引起的额外计算量，较细的单元划分会增加矩阵的维度但可以提高计算精度，因此与普通有限元一样，实际计算中要寻求符合计算精度的较稀疏网格进行计算。另外，采用SFEM计算SEA参数的原理与层合板理论计算SEA参数是一样的，都是基于波数的方法，因此其计算时间在求解过程中并不会频率升高显著增加。下一步我们将对如何提高SFEM的计算效率做专题研究。

图11 板D的模态密度

图12 板E的模态密度

图13 板F的模态密度

图14 层合板理论相对于SFEM的偏差

表2 计算时间对比

三明治理论虽然也能较为准确的计算模态密度，但是它仅适用于特殊夹芯结构；层合板理论在厚度和频率增大时，计算偏差会变大；三明治板理论和层合板理论都是求解高次方程，难于进行模态分类，而SFEM在求解特征值问题时能利用皮尔逊相关系数准确地进行模态分类，适用于一般层合板结构的分析，并且对网格要求低，计算量小，精确度高。对于厚板，还可以适当增加单元数目来提高SFEM计算结果的准确性。

3.4 汽车玻璃层合板

板G是一种典型汽车玻璃(Glass)层合板，由两块厚度为4 mm的玻璃板和0.75 mm厚的PVB材料组成，其如图15所示[13]。

图15 板G结构示意图

Glass的材料为各向同性材料，杨氏模量为62 GPa，泊松比为0.24，密度为2 300 kg/m3；PVB为黏弹性材料，泊松比为0.499，其剪切模量与频率有关，其关系如图2所示，PVB材料的内损耗因子也是与频率有关，其关系如图17所示，Glass的内损耗因子取0[13]。

图16 PVB材料剪切模量

板G的三层各划分为一个单元，利用SFEM计算获得各模态的波数和特征向量，根据式(26)获得板G的总内损耗因子如图18所示，基于SFEM的总内损耗因子计算结果与文献[13]和VA one软件的结果几乎是一致的，验证了SFEM计算总内损耗因子的准确性。弯曲模态的内损耗因子随频率的升高先增大后减小，量级为10-1左右；而剪切模态和拉伸模态的损耗因子持续增大，量级为10-5～10-4左右。黏弹性层PVB对弯曲模态损耗因子影响较大，而对剪切模态和拉伸模态的损耗因子影响较小。

图17 PVB材料内损耗因子

图18 板G内损耗因子

4 总 结

本文基于二次谱单元建立层合板的特征方程，在波数求解时，采用了皮尔逊相关系数来衡量模态的相似度，根据相似度准确地进行波数分类，进而计算层合板的统计能量分析参数。通过数值模拟验证了SFEM方法的通用性与有效性。本文的SFEM在处理复杂铺层、任意厚度和材料分布的层合平板方面，比传统层合板理论精度更高、适用范围更广，且其计算量远远小于传统有限元方法。但要解决曲板、加筋板等复杂层合板结构的建模问题，需要发展高维的谱单元。本文基于SFEM的SEA参数计算方法对飞行器层合板结构的高频建模有较重要的工程实际意义。