复合材料层合板自由振动的谱切比雪夫解法

郭琛琛， 刘 涛， 王青山， 秦 斌

(1. 中南大学 轻合金研究院， 长沙 410083； 2. 中南大学 高性能复杂制造国家重点实验室， 长沙 410083；3. 中南大学 交通运输工程学院， 长沙 410075)

复合材料层合板由于其较高的比强度和比刚度等力学特性，被广泛应用于土木工程、海洋工程和航空航天等诸多工业领域。这类构件通常面临各种极端恶劣的工况，并且可能受到形式各异的边界载荷约束。因此，不同边界条件下复合材料层合板的振动特性预测一直是研究人员关注的焦点。

以往，层合板的振动问题求解主要是基于理论解[1-2]。然而，由于数学问题的复杂性以及对实际问题进行建模的高要求，这些常规方法难以同时满足计算效率和计算精度要求。随着研究的不断深入，人们提出了一系列的数值求解方法，如Galerkin法[3]、里兹法[4-6]、改进傅里叶级数法[7]、无网格法[8-15]等。Shi等用Galerkin法对四边固定的任意层合板进行了自由振动分析。Ye等用Raleigh-Ritz法得到了具有一般边界约束和内线支承的复合材料中厚层板自由振动问题的修正傅立叶解。Shi等用改进的傅里叶级数方法(IFSM)研究了弹性地基上均匀支承和多点支承的中厚复合材料层合矩形板的自由和强迫振动特性。Xiang等提出了一种基于样条径向基函数和高阶剪切变形理论的无网格法来分析固支层合板的自由振动问题。张科使用基于移动最小二乘(MLS-L)近似函数构造的三维无网格法，对正交层合板的振动问题进行了分析。陈莘莘等提出了无网格自然邻接点Pertov-Galerkin法分析了复合材料层合板的自由振动。使用插值型重构核粒子无网格法，李情等分析了复合材料层合板的自由振动特性。通过对上述文献分析可知，目前对于分析层合板振动特性的数值方法种类丰富，不同方法之间的区别主要是基/试函数和求解方法的不同，但是目前绝大多数的数值方法缺陷是在计算层合板振动问题建模过程复杂，并且求解过程一般采用传统的数值积分方法，影响计算效率。

谱切比雪夫法(ST)采用具有指数收敛特性、优越的数值稳定性和高计算精度[16-17]的切比雪夫多项式作为基函数，来分析结构的动力学问题。在文献[18]中，Yagci等首次系统阐述了谱切比雪夫法，其核心思想在于：一是采用准确的导数和内积矩阵[19]来求解方程中导数和积分运算，以提高算法的计算效率；二是采用Gauss-Lobatto节点[20]对函数进行离散化处理以保证计算结果的准确性。基于以上优点，谱切比雪夫法已成功应用于各种结构[21-26]的动力学特性研究中，例如Bediz使用二维谱切比雪夫(2D-ST)法预测了任意几何形状厚板在不同边界条件下的动力特性，然而对于复合材料层合板的自由振动问题暂无涉及。

本文将谱切比雪夫法应用于复合材料层合板的自由振动特性研究。基于一阶剪切变形理论，引入边界弹簧技术实现不同边界条件的模拟，得到层合板的能量方程，应用二维谱切比雪夫法求解能量方程获得结构的振动分析模型。通过数值分析算例，与其他方法的计算结果进行比较，证明了此方法的有效性，还在此基础上研究了弹性模量比和铺设角对复合材料层合板自由振动特性的影响。

1 理论公式

1.1 模型描述

如图1所示为复合材料层合板的几何形状和坐标系。直角坐标系o-xyz位于结构的几何中面上，层合板沿x,y和z方向的长度、宽度和厚度分别是a、b和h；α是层纤维方向和x轴之间的铺设角度；U、V、W分别表示层合板中任意一点在x、y、z方向上的位移。在本文中，采用边界弹簧技术实现不同类型的边界条件，引入了三组线性约束弹簧(ku、kv、kw)和两组转角约束弹簧(kφ、kθ)，kt(t=u,v,w,φ,θ)的值可以从零变化到无穷大，它们的不同组合可以表示各种边界约束条件。例如，kt= 0(t=u,v,w,φ,θ)表示自由边界条件(F)，而kt= ∞(t=u,v,w,φ,θ)表示固支边界条件(C)。

图1 复合材料层合板的几何模型与坐标系Fig.1 Geometric model and coordinate system of composite laminated plates

1.2 能量方程

应用一阶剪切变形理论，复合材料层合板中任意一点的位移场的表达式为

(1)

式中：广义位移向量u= (uvwφθ)T，u,v和w表示在层合板几何中面上沿着x,y和z方向的位移分量；φ和θ分别表示分别为x-z平面和y-z平面的转角位移。

由式(1)得应变-位移线性关系(几何方程)为

(2)

(3)

根据广义胡克定律，层合板第k层的应力-应变关系(物理方程)为

(4)

力和力矩的合力是通过积分板厚度上的应力获得的，通过引入剪切校正因子κ，最终可以得出广义力与应变之间的关系。其关系式表示如下

(5)

式中：N={NxNyNxy}T,M={MxMyMxy}T,QS={QxQy}T分别为合力，合力矩和等效剪力向量；A为拉伸刚度矩阵，B为耦合刚度矩阵，D为弯曲刚度矩阵，AS为剪切刚度矩阵，它们的表达式如下所示

(6)

因此，存储在板中的应变能US为

(7)

式中，R为5×5线性对称微分算子矩阵。其表达式为

(8)

动能T的计算表达式为

(9)

利用方程(1)定义的位移场，沿z方向求积分，动能项可改写为

(10)

G为5×5对称质量矩阵，其表达式为

(11)

式中：(I0,I1,I2)为惯性矩；ρk为第k层材料密度。Usprings是具有弹性均匀边界条件的板的四边的势能，它是用五种弹簧均匀分布在四个边来模拟。

(12)

式中,k=diag(kukvkwkφkθ)。

因此层合板的拉格朗日能量方程可写为

L=T-US-Usprings

(13)

1.3 谱切比雪夫近似原理

切比雪夫多项式为递归和正交多项式，在区间x∈[-1,1]上，其表达式为

Cl(x)=cos(larc cos(x)),l=0,1,2,…,

(14)

应用区间变换，对于任意区间[l1,l2]的函数y(x)，可以使用有限项的切比雪夫多项式展开表示为

(15)

在离散化系统(即确定采样点)时，采用Gauss-Lobatto采样方法[20]

l=0,1,…,M-1

(16)

扩展系数al与采样点yl=y(xl) 之间存在一对一的映射关系，采样函数y={y0,y1, …,yM-1}的向量与(截断的)切比雪夫展开系数向量a={a0,a1, …,aM-1}之间的关系可以写成

y=ΓBa

a=ΓFy

(17)

(18)

b=a

(19)

(20)

其中s= (l2-l1)/2为区间缩放系数。

应用切比雪夫多项式的积分性质，y(x)的积分表达式为

(21)

为了得到内积矩阵，令h(x) =f(x)g(x)，假设fd=diag[f1f2…fM]，跟据式(17)和(21)有

(22)

令vTΓF= [V1…VM]，故式(22)可改写为

(23)

对于二维问题，函数f(x,y)将产生第二阶张量，其组成元素表达式为

flm=f(xl,ym)

(24)

式中，l= 1, …,Mx和m= 1, …,My是张量的指数。fkl为组成Mx×My二维采样函数矩阵F的元素值。使用以下映射算法获得采样二维函数的矢量形式，为

fc=flm,c=(l-1)My+m

(25)

应用二维谱切比雪夫法，在坐标(x,y)中的每个位移点都用双切比雪夫多项式展开表示为

(26)

式中:alm是双展开系数;Mx和My是在每个方向上使用的切比雪夫多项式数。将式(25)张量到矢量的映射应用于展开系数，式(26)可改写为

(27)

二维问题的展开系数与采样函数值的关系为

f=ΨBa

a=ΨFf

(28)

(29)

对于微分运算，需要获得二维问题中关于x和y方向的导数矩阵Dx和Dy。按照文献[25]的方法来扩展导数矩阵，就可以分别得到相对于x和y的扩展导数矩阵Dx和Dy。使用向前和向后变换矩阵以及扩展的导数矩阵，有如下定义：

fx=ΨBbx=ΨBDxa

fy=ΨBby=ΨBDya

(30)

对于面积积分，使用内积矩阵方法

∬Af(x,y)g(x,y)dxdy=fTVg

(31)

V为内积矩阵，其表达式为

V=x⊗y

(32)

1.4 振动求解

应用二维谱切比雪夫法求解式(13)可得到

(33)

式中，M,KU和KS分别代表质量矩阵。应变能刚度矩阵和边界弹簧刚度矩阵，其具体表达式如下

(34)

t=u,v,w,φ,θ

Vx=x⊗INy×Ny,Vy=INx×Nx⊗y

(35)

I为单位矩阵。

对于式(33)，采用弱形式的求解方法得到对未知变量求得偏导数方程

(36)

将式(36)以矩阵的形式表达，则结构的振动特性方程可以写为以下形式

(KU+KS-Mω2)A=0

(37)

式中:A代表未知系数的矢量;ω表示结构的固有频率。通过直接求解式(37)，可以得到板的频率特性和相应的特征向量。

2 数值分析与算例

基于前面得到的自由振动特征求解方程，本部分将对复合材料层合板的自由振动特性进行分析。为简化研究，如无特别说明，在接下来的计算分析中，层合板的几何参数为：a= 1 m,b= 1 m,h= 0.1 m；材料参数为：E2= 10 GPa、E1= 25E2、G12= 0.5E2、G13=G12、G23=0.2E2、μ12= 0.25、μ12=μ21、ρ= 1 500 kg/m3，κ= 5/6。如前所述，对于不同的边界可以设置与之对应的弹簧约束刚度实现，符号C、S及F分别用来表示固支、简支、自由边界。特别地，本文所涉及的边界组合按逆时针方向设定。此外，本文方法具有简便性和通用性，当结构的几何参数和边界条件发生变化时，只需要修改相关参数即可，并不用重新推导，从而可以大大简化计算过程。此外，本文中采用式(38)对层合板的固有频率进行无量纲化。

(38)

2.1 数值验证

在本节中，通过与文献中的结果进行比较，对二维谱切比雪夫法进行了数值验证。首先，表1为改变x和y方向上所需的多项式数Mx×My，四端固支的方形层合板的前五阶无量纲频率。层合材料的纤维铺设角为[30°/-30°]，板厚h=0.05 m。从表1可以看出，随着多项式数的递增，各阶频率依次从低阶到高阶收敛至稳定值, 当Mx×My=16×16时，所有无量纲频率均已达到收敛值。与Jacobi-Ritz法进行对比，两种方法均是基于一阶剪切变形理论，所以两者前五阶无量纲频率计算结果高度吻合，相对误差在0.001%以内。此外，本文的谱切比雪夫法达到收敛值的矩阵维数为1 280×1 280 (5MxMy× 5MxMy)，而文献[6]的矩阵维数为5 120×5 120。本文只需要少量的多项式数即可获得足够准确的结果，充分说明了本文所提出的谱切比雪夫法具有收敛速度快，同时计算准确性好的特点。

表1 采用不同多项式数 Mx×My的复合材料层合板的前五阶无量纲频率对比Tab.1 Comparison of the first five-order dimensionless frequency of composite laminated plates with different Mx×My

表2所示为在固支边界条件下，[0/α]和[α/-α]两种铺设方式下的复合材料层合板的前三阶无量纲频率，密度为ρ=1 000 kg/m3。将计算结果分别和Galerkin法和局部移动Kriging无网格法的计算结果进行了比较,从表2可以看出，三种方法的计算结果可以很好地匹配，微小的误差来源于主要是因为不同数值方法，不同方法之间的计算原理在本质上存在一定差异。

表2 [0/α]和[α/-α]两种铺设方式下复合材料层合板的前三阶无量纲频率对比Tab.2 Comparison of the first three-order dimensionless frequencies of composite laminated plates with the two layer types [0/α] and [α/-α]

最后，由于本文方法是为了求解层合板在任意边界条件下的自由振动问题，因此有必要验证其在各种边界条件下的有效性和计算精度，相关对比数据结果如表3所示，参考值为修正傅里叶解计算结果。材料密度ρ=1 000 kg/m3，纤维铺设角为[45°/-45°]。在该算例中，考虑了四种边界条件，分别是C-C-C-C、S-S-S-S、F-C-F-C、F-S-F-S，和两种宽厚比，即h/a=0.05和0.2。两种方法同时基于一阶剪切变形理论，只是基函数和求解形式不同，所以计算结果非常接近。结果表明，本文方法适用于任意边界条件。

表3 不同边界条件下[45°/-45°]复合材料层合板的前三阶无量纲频率对比Tab.3 Comparison of the first three-order dimensionless frequencies of [45°/-45°] composite laminated plates with different boundary conditions

2.2 参数化研究

根据引言调研的文献可知，在工程实践中，板结构的材料属性往往是多种多样的，而纤维的敷设角度和纤维自身的弹性模量取值对于层合板的振动特性而言至关重要，这也直接反应在复合材料的纤维类型选择和加工工艺上。因此，详细研究纤维的敷设角度和纤维自身的弹性模量取值对复合材料层合板的自由振动特性的影响对指导实际工程设计和应用有着重要意义。本节以四端固支，铺设方式为[α/-α]的层合板为例，研究了弹性模量比E1/E2和铺设角α对复合材料层合板无量纲频率的影响。层合板的其余结构参数参考默认。E1/E2的变化区间为[1,60]，α从0°递增至180°，提取层合板的1阶、10阶和100阶无量纲频率，最终计算结果如图2所示。需要注意的是，当计算100阶固有频率时，所采用的多项式数为Mx×My=40×40。从图2可以看出，随着E1/E2的增大，1阶、10阶和100阶的无量纲频率也随之递增，E1/E2取最大值60时，频率也达到了最大值；对于铺设角α，1阶、10阶和100阶的无量纲频率变化规律关于α=90°对称，并且模态阶次不同，频率变化规律也会略有不同，1阶频率在0°、90°和180°附近达到最大值，10阶和100阶频率则在区间[30°,60°]和[120°,150°]内取最大值。

1阶

10阶

100阶图2 弹性模量比E1/ E2和铺设角α对层合板无量纲频率的影响Fig.2 Effects of elastic modulus ratio E1/ E2 and laying angle α on non-dimensional natural frequencies of laminated plates

3 结 论

本文在一阶剪切变形理论的基础上，应用二维谱切比雪夫法建立了任意边界条件下复合材料层合板的自由振动分析模型。通过数值对比分析和参数化研究，得出以下结论：

(1) 所采用的谱切比雪夫法具有收敛速度快，计算准确性好的特点，可以有效分析任意边界条件下复合材料层合板的自由振动问题。

(2) 弹性模量比E1/E2对复合材料层合板的自由振动频率有重要影响。随着弹性模量比E1/E2的增加，结构的无量纲固有频率也会显著增加。

(3) 对于铺设角α，复合材料层合板的无量纲频率关于α=90°对称分布，并且与模态阶次有一定的关联性。