复合材料封闭变截面薄壁梁自由振动分析

马静敏， 任勇生，姚文莉

(1.山东科技大学 理学院，山东 266510;2.山东科技大学 机械电子工程学院，山东 266510;3.青岛理工大学 理学院，青岛 266510)

复合材料由于轻质、高强度、高刚度以及抗疲劳性能突出等优点［1－2］，在工程领域中得到了广泛的应用。封闭截面复合材料薄壁梁是航空、宇航结构及风力机叶片的主要结构形式，研究其振动动态特性对其结构动力学设计具有十分重要的意义。

Mansfield等［3－5］提出复合材料封闭截面薄壁梁的运动模型，由于该模型忽略了梁的横向剪切变形和横截面的翘曲影响，仅适用于解决具有特定铺层顺序的复合材料梁的振动问题。Libove［6］建立了用于计算封闭截面薄壁梁剪切流和正应力的简单理论，该理论虽然考虑了弯曲剪切变形，但是没有考虑到由于不均匀翘曲引起的剪切变形。Librescu等［7－11］建立了广泛适合于工程领域问题的模型，该模型引入了弯曲中的剪切弹性变形的影响，但是翘曲扭转变形被忽略了。

上述这些研究都是针对等截面直梁的复合材料薄壁结构模型进行静力或动力建模与分析，而实际的风力机叶片、机翼截面尺寸和形状是变化的，所以研究变截面复合材料薄壁梁的动力学性能对于解决工程实际问题具有重要意义。Librescu 等［8，12－14］对复合材料封闭截面薄壁梁进行了一系列的研究，建立了一套理论体系，用于解决等直梁、变截面梁和变截面旋转梁的振动问题。这些研究考虑经典的由于扭转产生的翘曲和横向剪切变形的影响，使用哈密顿原理获得振动微分方程。采用这种方法由于需要涉及分部积分运算，推导过程较为繁琐。

本文综合考虑了均匀轴向拉伸、绕y轴和z轴的弯曲所引起的截面翘曲，但是暂时未考虑横向剪切变形，基于拉格朗日方程推导出复合材料封闭变截面薄壁梁的自由振动方程。可以看到，相对于哈密顿原理，本文所采用的方法使公式推导变得十分简洁。在建立复合材料封闭变截面薄壁梁振动微分方程的基础上，通过应用迦辽金法，通过数值计算得到周向均匀刚度(CUS)和周向反对称刚度(CAS)配置的矩形变截面薄壁悬臂梁的固有频率，并且与ANSYS有限单元法的计算结果进行了比较。分析了复合材料的各向异性以及纤维铺层角和截面变化对复合材料封闭变截面薄壁梁固有振动的影响。

1 自由振动方程推导

1.1 薄壁梁动能和势能

截面形状任意且沿轴向任意变化的细长变截面薄壁弹性梁及参考坐标系如图1所示。其中，L表示杆件的长度，h表示截面任意位置处的壁厚，R表示截面任意位置处的中线半径，d表示截面的特征尺寸，且满足下列条件:

图1 薄壁梁几何图及参考坐标系Fig.1 Beam geometry and coordinate systems

为了分析该梁的自由振动，定义x位置所在的截面沿x，y，z坐标轴方向的平均位移分别为 U1(x)，U2(x)，U3(x)，该截面扭转变形的扭转角为φ(x)。

悬臂梁上任意一点沿坐标系(x，y，z)的三个坐标轴方向的位移分量为［16］:

其中:g(s，x)表示由于均匀轴向拉伸、绕y轴和z轴的弯曲所引起的截面翘曲。

薄壁梁的动能和势能为［16］:

其中:(·)，(··)表示对时间t求1阶和2阶导数，依次类推;()'，()″表示对位置坐标x求1阶和2阶导数，依次类推。

该势能表达式是根据复合材料薄壁结构的特点，依据二维各向异性壳理论，综合考虑了截面薄壁特点和翘曲影响建立的。此外，由文献［16］可知，［C］4×4中的元素都是由封闭曲线积分计算得到的，因此充分反映了封闭截面的特点。

1.2 薄壁梁拉格朗日方程

对于三维连续系统 κ =κ(x，y，z，t)，相应的拉格朗日函数为 L= ∫∫∫ζd x d y d z，其中 ζ为拉格朗日密度［18］。根据薄壁梁的动能和势能，拉格朗日密度为:

对该薄壁杆件系统应用哈密顿原理有:

交换积分与变分符号可得:

其中:

将上式适当变换，代入式(6)，并注意到:

可得:

因为在(t1，t2)内:δU1，δU2，δU3，δφ 都不等于零，所以可得薄壁梁的拉格朗日方程为:

1.3 变截面薄壁梁自由运动方程

对于变截面直杆件而言，［C］4×4随截面位置的变化而变化。将式(5)代入式(7)中，可得复合材料薄壁梁的自由运动方程为:

其中:

等直梁的情况下，式(8)与文献［16］直接应用哈密顿原理建立的振动微分方程是一致的。

2 应用

其中:σa=aT/aR，σb=bT/bR表示锥度，η=x/l表示横截面位置的的无量纲量且η∈［0，1］。下标R和T分别表示根截面和端截面。若规定由根部到端部截面尺寸逐渐减小，则 σa∈［0，1］和 σb∈［0，1］，且锥度取0值时对应着端节面尺寸变为0，而取1时对应着截面尺寸不变化。

本文将式(8)的运动微分方程用于求解矩形变截面薄壁层压复合材料悬臂梁的固有频率。层压板采用两种铺设方式，一种是CUS，一种是CAS。矩形截面梁的各项参数，材料的力学性能，铺层方式见图2、图3和表1。矩形截面的宽度，高度随着截面位置呈线性变化。任一位置x处的截面宽度和高度可以表示为:

表1 薄壁梁尺寸和材料性能参数Tab.1 Cantilever geometry and properties

图2 矩形梁几何形状参数Fig.2 The geometry of box beam

图3 两种铺层方式Fig.3 Two cross-section layup

2.1 两种构型梁的自由振动微分方程

对于CUS结构来说，由于轴向刚度、耦合刚度、剪切刚度、壁厚及材料密度在全截面上是常数 ，所以C13=C14=C23=C24=C34=0，且静矩Sy=Sz=0。式(8)可以简化为:

从运动方程可以得到，对于CUS构型的薄壁梁来说，式(11)和式(12)描述了拉伸与扭转的耦合振动，式(13)描述了绕z轴的弯曲振动，式(14)描述了绕y轴的弯曲振动。

对于CAS结构来说，轴向刚度、壁厚及材料密度在全截面上是常数，而矩形截面的耦合刚度B对立边符号相反，所以C12=C13=C14=C24=C34=0，且静矩Sy=Sz=0。式(8)可以简化为:

从运动方程可以得到，对于CAS构型的薄壁梁来说，式(16)和式(17)描述了绕y轴弯曲与扭转的耦合振动，式(15)描述了轴向振动，式(18)描述了绕z轴的弯曲振动。

显然，由于复合材料的力学性能随着铺层方式发生变化，不同的铺层方式，将引起不同的耦合振动，CUS构型存在拉扭耦合振动，而CAS构型将引起弯扭耦合振动。因此可以通过对复合材料进行弹性剪裁来获得所需要的力学性能。

2.2 固有频率求解方法

采用连续系统振动和稳定性分析通用的Galerkin近似求解方法［19－20］，利用假设振型消除叶片空间位置变量，将自由振动方程转化为关于时间的常微分方程。假定梁的轴向变形、扭转变形和弯曲变形可以表示如下:

拉伸、弯曲和扭转选取标准的非旋转、非耦合、均匀悬臂梁的振型函数［21］:

将式(19)，式(2

0)代入CUS构型梁的振动微分方程式(11)～式(14)，由Galerkin法，通过振型函数加权积分，得:

其中:D1(U1，φ，U2，U3)，D2(U1，φ，U2，U3)，D3(U1，φ，U2，U3)，D4(U1，φ，U2，U3)分别表示振动方程式(11)～式(14)等号左边的表达式。

式(21)经化简可得到下列的CUS构型薄壁变截面梁的4N个常微分方程组:

其中:

上述矩阵中的元素为:

同理可得CAS构型梁的4 N个振动微分方程也为式(22)，M阵和K阵与CUS构型相同。CAS构型薄壁变截面梁的K阵为:

其中:

求解式(22)可得不同构型下矩形变截面薄壁梁各种振动解析解的固有频率。

2.3 两种构型梁振动固有频率数值分析

图4 ω1随锥度变化规律Fig.4 ω1 vs taper ratio

图5 ω2随锥度变化规律Fig.5 ω2 vs taper ratio

图6 ω3随锥度变化规律Fig.6 ω3 vs taper ratio

将CAS和CUS构型梁在30°铺层角时前五阶固有频率随σa和σb的变化绘制成曲面图，如图4～图8所示。由图4～图8可知，无论锥度如何变化，CAS构型的前五阶固有频率要大于CUS构型的固有频率且两种构型固有频率随锥度的变化规律基本是一致的。

令 σb=1只有 σa变化和令 σa=1只有 σb变化，CUS构型前五阶固有频率随σa和σb单独变化规律如图9和图10所示，CAS构型如图11和图12所示。由图9～图12可以说明，梁的固有频率随σa的变化规律和随σb的变化规律是不同的，而两种构型梁随σa的变化规律基本相同，随σb的变化规律也基本相同。这种变化规律的差异是因为前五阶振动分别属于不同的振动模态，而锥度的变化，对于不同的振动模态的刚度影响不同。此外，图10和图12中的曲线存在着尖点，是因为固有频率由小到大顺序排列的缘故，并非固有频率随锥度不光滑地变化。

图7 ω4随锥度变化规律Fig.7 ω4 vs taper ratio

图8 ω5随锥度变化规律Fig.8 ω5 vs taper ratio

图9 CUS 构型梁 ω1，ω2随锥度变化规律Fig.9 ω1 and ω2 vs taper ratio of CUS beam

图10 CUS 构型梁 ω3，ω4，ω5随锥度变化规律Fig.10 ω3，ω4 and ω5 vs taper ratio of CUS beam

图11 CAS构型梁 ω1，ω2随锥度变化规律Fig.11 ω1 and ω2 vs taper ratio of CAS beam

图12 CAS 构型梁 ω3，ω4，ω5随锥度变化规律Fig.12 ω3，ω4 and ω5 vs taper ratio of CAS beam

图13 ω1和 ω2随铺层角度变化规律Fig.13 ω1 and ω2 vs ply angel

图14 ω3，ω4和 ω5随铺层角度变化规律Fig.14 ω3，ω4 and ω5 vs ply angel

图13和图14说明了在锥度σa=σb=0.5的条件下，CAS和CUS构型梁前五阶固有频率随铺层角度的变化规律。两种构型梁的前五阶固有频率随铺层角度的变化规律都是下降曲线，由文献［17］可知，CUS构型的前五阶振动为弯曲振动模态和拉扭耦合振动模态中的一种，而CAS构型梁为弯曲为主的振动模态。由铺设方式可知，采用0°和90°铺层时，CUS和CAS配置方式实质是相同的，因此固有频率应该相等。图12和图13中，CUS和CAS构型的曲面在0°和90°时相交与实际情况相符。而其它角度铺层时，CAS构型的固有频率要大于CUS构型的固有频率，在20°左右时，差别最显著。因此将20°时，锥梁和等直梁的固有频率值列于表2中。可见锥度对于固有频率值存在着显著的影响，表3则从具体数值上说明铺层角度和铺层方式对于固有频率的影响。

需要说明的是，Gakerkin法求解的固有频率的精度与保留模态个数有关。固有频率随保留模态数的变化情况列于表4中。可见，CUS与CAS构型梁的各阶固有频率求解精度与保留模态数的关系不完全相同。且不同阶的固有频率随模态数增加，收敛情况也不相同。对于低阶的固有频率值保留较少的模态数就达到较高的精度，而对于高阶固有频率值则需要保留较多的模态以达到足够的精度。

因为薄壁结构梁在ANSYS中采用线性壳单元SHELL99，设定铺层方式和材料性能参数及约束，将图2中的锥度σa=σb=0.5，铺层角度为20°的梁划分为1416个四边形单元进行分析。ANSYS有限单元法和Galerkin法保留四阶模态的求解结果列于表5中。由表5可知，CUS和CAS构型梁的低阶固有频率值与ANSYS求解结果非常接近。充分说明本文建立的运动微分方程和采用Galerkin法求解的正确性。此外表中CAS构型梁的第三阶和第五阶固有频率值与有限单元法的求解结果差值较大，原因为保留模态数较少。若增加保留模态数，则差值会明显降低。

表2 20°铺层时CAS与CUS梁前五阶固有频率Tab.2 The first five natural frequency of CUS and CAS beams with 20°ply angel

表3 σa=σb=0.5时CAS与CUS梁前五阶固有频率Tab.3 The first five natural frequency of CUS and CAS beams with σa=σb=0.5

表4 模态个数对于固有频率的影响(σa=σb=0.5，θ=20°)Tab.4 Effect of the number on natural frequencies(σa=σb=0.5，θ=20°)

表5 Galerkin法和有限单元法计算结果(σa=σb=0.5，θ=20°)Tab.5 Solution of Galerkin method and finite element method(σa=σb=0.5，θ=20°)

3 结论

采用拉格朗日方程推导了任意封闭变截面薄壁复合材料梁的自由振动方程。该方程适用于沿截面周线刚度任意变化的变截面直梁。应用该方程使用Galerkin法求解了两种铺层方式的复合材料梁的固有频率。有限单元法的求解结果和理论分析解析结果一致，说明了所建立公式的正确性。此外求解结果也反映了薄壁复合材料锥形梁自由振动的弹性耦合机制。说明不同的铺设方法将引起不同的耦合振动，即使是同样的铺设方法，采取不同铺设角度和不同的锥度，对于固有频率也存在影响。这些结果对于风力机叶片等工程结构的设计，都具有指导意义。

［1］Barbero E J.Introduction to composite material design［M］.Taylor＆ Francis Inc.，Philadelphia，1999.

［2］Jonesm P M.Mechanics of composite materials［M］.Taylor＆ Francis Inc.，Philadelphia，1999.

［3］Mansfield E H，Sobey A J.The fiber composite helicopter blade.Part I:stiffness properties.Part II:prospects of aeroelastic tailoring［J］.Aeronaut Quart，1979，30:413－449.

［4］Mansfield E H.The stiffness of a two-cell anisotropic tube［J］.Aeronaut.Quart，1981，32:338 － 353.

［5］Kapania R K，Raciti S.Recent advances in analysis of laminated beams and plates 1.shear effects and buckling［J］.AIAA，1989，27(7):923 －934.

［6］Libove C.Stresses and rate of twist in single cell thin-walled beams with anisotropic walls［J］.AIAA ，1988，26(9):1107－1118.

［7］Librescu L，Song O.On the aeroelastic tailoring of composite aircraft swept wings modeled as thin-walled beam structures［J］.Compos Eng，1992，2(5 －7):497 －512.

［8］Song O，Librescu L.Free vibration of anisotropic composite thin-walled beams of closed cross-section contour［J］.Sound Vibr，1993，167(1):129 －147.

［9］Bhaskar K，Librescu L.A geometrically non-linear theory for laminated anisotropic thin-walled beams［J］.Eng Sci，1995，33(9):1331－1344.

［10］Qin Z， Librescu L. On a shear-deformable theory of anisotropic thin-walled beams:further contribution and validations［J］.Compos Struct，2002，56(4):345 － 358.

［11］Na S S，Librescu L.Dynamic response of elastically tailored adaptive cantilevers of non-uniform cross section exposed to blast pressure pulses［J］.Impact Eng，2001，25(9):847－867.

［12］Librescu L，Na S S.Active vibration control of doubly tapered thin-walled beams using piezoelectric actuation［J］.Thin-Walled structures，2001，39:65 －82.

［13］Na S S，Librescu L.Dynamic behavior of aircraft wings modeled as doubly-tapered composite thin-walled beams［C］.Recent Advances in Solids and Structures，New York:ASME，1999(398):59 －68.

［14］Chandiramani N K，Librescu L，Shete C D.On the freevibration of rotating composite beams using a higher-order shear formulation［J］.Aerospace Science and technology，2006(6):545－561.

［15］Zhao Y H，Hu H Y.Structural modeling and aeroelastic analysis of high-aspect-ratio composite wings［J］.Chinese journal of aeronautics，2005，18(1):25－30.

［16］Armanios E A，Badir A M.Free vibration analysis of anisotropic thin-walled closed-section beams［J］.AIAA J，1995，33(10):1905－1910.

［17］Dancila D S，Armanion E A.The influence ofcoupling on the free vibration of anisotropic thin-walled closed-section beams［J］.Int.J.Solids Structures，1998，35(23):3105 － 3119.

［18］易 中，周丽珍.分析力学初步［M］.北京:冶金工业出版社，2006:38－42.

［19］Thomsen J J.Vibrations and stability:advanced theory，analysis，and tools［M］.Springer-Verlag:Berlin-Heidelberg-New York，2003.

［20］Meirovitvh L.Computational methods in structural dynamics［M］. Sijthoff-Noordhoff:Alphen aan den Rijn， The Netherlands，1980.

［21］Hodges D H ，Dowell E H，Nonlinear equations of notion for the elastic bending and torsion of twisted nonuniform rotor blades［R］.NACA TND －7818，1974.