基于能量有限元法的功能梯度梁振动分析

王 迪, 朱 翔，2, 李天匀，2, 衡 星, 高 双

(1.华中科技大学 船舶与海洋工程学院，武汉 430074；2.船舶与海洋水动力湖北省重点实验室，武汉 430074)

功能梯度材料(Functionally Graded Material，FGM)[1]通常是由两种或两种以上不同性能的材料组成，材料的基本要素沿结构的某一方向连续变化。功能梯度材料最早是由日本学者提出，由于其优良的结构性能和重要的应用价值，引起人们广泛的研究和关注。

梁结构是工程中最为广泛使用的基本结构，自功能梯度材料提出以来，许多学者对功能梯度梁结构的振动特性进行了研究。邓昊等[2]基于动力学的相关理论，通过对状态空间变量的替换，研究了功能梯度Timoshenko梁的理论模型。Yang等[3]利用经典的欧拉梁理论，对不同边界条件的功能梯度梁进行了自由振动分析，研究了功能梯度梁的振动特性。Simsek[4]将梁的经典理论、一阶剪切和高阶剪切理论与拉格朗日方程结合起来，分析了在不同边界条件下功能梯度材料梁的固有频率和振动特性。Alshorbagy等[5]应用虚功原理，推导了材料性质沿厚度方向幂函数变化的功能梯度材料梁自由振动的平衡方程并采用有限元数值方法进行求解对比。Ravikiran等[6]基于高阶剪切理论，采用有限元法得到陶瓷-金属功能梯度材料梁的平衡方程，应用刚度矩阵法求得了功能梯度材料梁的挠度解。

在研究结构振动的分析方法中，能量有限元法(Energy Finite Element Method, EFEM)是近些年发展起来的一种新的数值方法，其主要目的是解决结构中高频的振动问题[7]。该方法将能量密度作为动力学方程的基本变量，通过有限元的思路将不同结构离散成单元来表示结构的特性，通过求解能量密度进而得到结构的振动特性。孙丽萍等[8]以能量密度为变量,推导了薄板结构受激励产生弯曲振动时的能量密度控制方程,并得到了控制方程的有限元解。蔡忠云等[9]以能量有限元法建立了复合材料层合梁受激励时的控制方程, 并研究了复合材料梁的横向振动问题。解妙霞等[10]利用能量有限元方法研究圆柱壳体结构，并根据薄壳理论和能量平衡关系推导了圆柱壳弯曲振动时能量密度的控制方程。而目前关于FGM结构的能量有限元研究尚未见到报道。

本文基于能量有限元法对功能梯度梁以及耦合功能梯度梁的振动特性进行研究。通过FGM梁的基本动力学方程推导得到FGM梁的能量密度控制方程和有限元矩阵方程。在FGM梁的算例分析中，将本文的能量有限元法和传统有限元法得到的能量密度解进行对比和分析，验证能量有限元法的准确性。在此基础上对耦合功能梯度梁结构的能量密度和能量流进行求解。

1 FGM梁的能量有限元方程

1.1 一维FGM梁结构的能量密度控制方程

FGM梁在定义的时候，通常假定材料的各项属性沿着梁截面厚度方向以指数形式变化。坐标系如图1所示，x为梁的长度方向，z为梁的高度方向。材料的属性变化主要是体现在杨氏模量E(z)和材料密度ρ(z)这两个方面

E(z)=E0eβz

(1)

ρ(z)=ρ0eβz

(2)

式中：E0和ρ0为z=0处的弹性模量和密度值；β=ln(E2/E1)/h为FGM梯度变化的常数。

图1 功能梯度梁模型坐标系

考虑一个受横向简谐力激励下的含阻尼的功能梯度梁，如图2所示，采用简单梁理论，其运动控制微分方程可以表示为[11]

(3)

图2 受到简谐激励作用下的功能梯度梁模型

假设方程的通解为

w(x,t)=(A1e-ikfx+A2eikfx+A3e-kfx+A4ekfx)eiωt

(4)

式中：kf为复波数；A1、A2、A3和A4为方程待求解的系数，由边界条件来确定。

方程的通解式(4)代入运动方程式(3)得到FGM梁的弯曲波波速的表达式为

(5)

由结构内部的能量平衡关系可知，横向振动梁的能量密度和功率远场解局部均值有如下的关系式[12]

(6)

通过加权余量法将上述偏微分方程转化为代数方程组，将权函数作为形函数，则可通过Galerkin加权残值法将梁的能量密度控制微分方程表示为[13]

(7)

同时，梁的能量流远场解与能量密度的远场解存在如下关系[14]

(8)

可以将式(7)用矩阵的形式表示出来

[Ke]{ee}={Pe}+{Qe}

(9)

式中：[Ke]为每个梁单元有关质量和刚度信息的系数矩阵；{ee}则是需要求解的各个节点的能量密度矩阵；{Pe}为在节点处输入的激励能量；{Qe}则为每个梁单元两端进出的能量流。

因此，根据能量有限元法求解式(9)可得到节点的能量密度，即可进一步得到每个梁单元的能量流。

1.2 耦合梁结构的能量有限元方程

对于更加复杂的耦合结构，由于耦合结构会在耦合处存在材料或几何的不连续，从而导致结构振动波在耦合处产生反射和透射。通常能量密度在耦合节点处是不连续的，而能量流是连续的。因此在建立能量有限元系数矩阵的同时，需要单独考虑耦合节点处的能量传递和连接，可以通过连接矩阵[15]的形式来解决耦合结构的连接问题。

图3 功能梯度梁模型中的能量传递过程

当仅有弯曲波由第一段梁向第二段梁入射时，两段梁的横向位移可以分别表示为

(10)

(11)

梁中弯矩M和剪力N可以表示为

(12)

(12)

式中：I为梁的截面惯性矩。

根据耦合节点处弯矩、剪力、位移以及转角的连续，可以得到四个连续条件

M1=M2

(14)

N1=N2

(15)

W1=W2

(16)

(17)

通过求解上述四个方程可得

(18)

(19)

(20)

(21)

(22)

忽略近场效应，可以得出梁耦合节点处的能量透射系数τ12和反射系数γ11分别为

(23)

(24)

在耦合节点处，能量流有正负两个方向，+代表正方向、-代表负方向，如图4所示，并且有如下的关系式[16]。

图4 耦合节点处的能量流

(25)

(26)

波在耦合节点处会同时产生反射和透射，主要有如下的关系式

(27)

(28)

式中：γ和τ分别为耦合梁处的反射系数和透射系数，其满足关系式τ+γ=1。将式(27)和式(28)联立，再利用耦合节点处功率连续性q1m+q2m=0的条件可以得到

(29)

同时，在耦合节点处，能量密度是由正负两个方向的能量密度值之和组成，有如下关系式

(30)

(31)

将式(25)～式(28)代入式(30)、式(31)可得

(32)

联立式(29)和式(32)可以得到

(33)

将关于能量流的项进行统一的整理，可以将系统的整体矩阵表示为

[Ke+J]{ee}={Pe}+{Qe}

(34)

由式(34)即可对耦合FGM梁的能量密度求解，从而得到振动能量。

2 FGM梁结构的能量有限元振动分析

2.1 FGM梁的能量密度计算验证

本节按照上节推导的公式，对某FGM梁的能量密度进行计算，得到结构中能量密度值的分布。同时也通过传统有限元模型对该梁的能量密度进行计算并对比，从而验证EFEA方法的正确性。

考虑图1所示的简支梁模型，其长度为4 m，截面是边长为0.02 m的正方形，梁的上表面材料为陶瓷，密度ρ为3 800 kg/m3，弹性模量E为3.8×1011Pa，μ为0.3，梯度变化常数β取85，因此可得到下表面的材料参数。结构的阻尼损耗因子η为0.02。在梁的中心处作用一输入功率为1.12×10-3W的横向激励力，激励频率为3 000 Hz。

按照上一节的推导，对整个梁划分能量有限元单元，单元数为10，形函数取n=3的Lagrange插值函数，计算得到能量密度分布如图5所示，其中参考能量密度值为1×10-12J/m。

为了验证本文能量有限元计算的准确性，同时也采用有限元法计算梁中的能量密度。有限元模型中的FGM梁采用8节点实体单元来模拟，梁的高度方向分成10份，每层梁单元按照梯度公式赋予不同的E和ρ来模拟功能梯度材料的变化，长度方向划分100份，通过提取中性面节点的速度来计算能量密度，如图5。从图中可见由于能量有限元法计算得到的是时间和空间平均后的解，因此能量密度分布比较均匀，而有限元解反映出局部特征，其得到的能量密度随着EFEM的解上下波动。但是两种方法的局部误差均不超过1 dB。

图5 两种方法沿梁长度方向能量密度分布的结果对比

由于能量有限元模型和传统有限元模型的单元数量不同，且能量有限元方法求得的能量密度是对时间和空间平均后得到的相对平滑的解，因此也可将有限元法得到的能量密度值按照能量有限元法设定的单元长度选出对应包含的节点，将这些节点上的能量密度进行一次平均，然后来比较两种方法的结果，如表1所示。从均值的量化结果来分析，两种方法也吻合较好，从而说明本文针对FGM梁的能量有限元分析是准确的。

表1 EFEA和FEA(均值)结果量化对比分析

与传统的有限元方法相比，能量有限元法的单元数量少，计算速度快，因此在计算效率上有了很大程度的提升，但是和有限元相比，其求解的能量反映的是时间和空间上的均值。与统计能量法相比，由于统计能量分析的基本变量是在一个子系统内对各组相似模态的能量密度在指定频宽内频率平均，而能量有限元则是同时对能量密度在一个波长的空间平均和在某个频宽内的频率平均，因此统计能量法只能求解子系统的能量，而能量有限元法能够求解出每个单元的能量，且对结构模态数没有要求，分析频段可扩展到中高频段。

2.2 耦合FGM梁的振动分析

如图6所示，考虑两段同轴的耦合FGM梁，每段梁长度为2 m，两段梁截面面积不同。第一段梁的截面是边长为0.01 m的正方形，第二段梁的截面是边长为0.02 m的正方形。两段梁的属性和3.1节相同，在第二段梁上施加输入功率为2.24×10-3W的横向激励力，激励频率为3 000 Hz。按照能量有限元法，将每段FGM梁划分成十个单元，然后计算梁单元节点处的能量密度值和梁单元的能量流。将计算得到的能量密度和能量流沿梁的分布表示出来，如图7和图8所示。

图6 耦合梁的EFEA分析模型

图7 耦合梁的能量密度分布

从图7和图8可见，耦合梁在耦合位置处能量密度发生了突变，这是由于耦合梁的截面产生突变，造成能量密度在耦合节点处不连续，同时可见在能量输入的位置附近能量密度较为集中。在图8中由于能量平衡，因此能量流在耦合位置处是连续的，没有发生突变，并且在激励输入处能量流达到峰值。同时，由于梁中阻尼的存在，远离激励位置的能量流逐渐衰减。

2.3 截面尺寸变化对耦合FGM梁的振动影响

考虑到耦合梁的截面参数变化会改变耦合面处的透射系数和反射系数，对连接矩阵产生影响，因此本节讨论截面尺寸变化对耦合FGM梁的振动影响。本节算例中除了梁的截面边长发生变化，其他参数均和上节相同。保持第一段梁的截面边长不变，改变第二段梁的截面边长，对比不同边长比下的能量密度和能量流结果，如图9和图10所示。

图8 耦合梁的能量流分布

图9 不同边长比下耦合梁的能量密度分布

图10 不同边长比下耦合梁的能量流分布

从图9和图10可见，两段梁的边长比变化对耦合梁的能量密度和能量流都会产生一定的影响。随着边长比的增大，两段梁的能量密度和能量流的差别也逐渐增大。当边长比达到3时，能量流在激励点的位置达到最大，此时连接处截面突变最明显，因此耦合处能量密度的突变也最明显。这是由于边长比变化会通过能量透射系数和反射系数引起连接矩阵的改变，随着边长比增大，能量在两段梁的分布差别也会逐渐增加。

为了进一步研究不同边长比下梁的振动规律，本节还计算了耦合梁单元的速度分布，如图11所示。从图中可以看出，由于耦合处截面的不连续，不同单元速度在耦合处依然是不连续的。相比第一段梁，由于第二段梁的截面发生变化，弯曲刚度也随之变化，因此第二段梁单元的速度变化更加明显。随着边长比的不断增大，单元速度值也随之减小。

图11 不同边长比下耦合梁的单元速度分布

3 结 论

本文基于能量有限元法对功能梯度梁和耦合功能梯度梁的振动特性进行了研究，推导了梁和耦合梁的能量有限元法基本方程，并推导了耦合FGM梁的连接矩阵和能量有限元矩阵方程。通过将EFEA计算得到的单元能量密度与传统有限元FEM的结果进行对比，验证了本文推导和求解的正确性。通过对耦合梁的能量密度和能量流计算可见，耦合FGM梁在耦合位置处能量密度发生了突变。能量流在耦合节点处保持连续，并且在能量的输入位置能量流达到峰值。由于阻尼的存在，远离激励位置的能量流逐渐衰减。耦合梁的截面变化越明显，耦合梁之间的能量密度变化也越明显。本文的研究为能量有限元法应用到复杂功能梯度材料结构的振动研究提供了借鉴。

[1] 李信，刘海昌，滕元成，等. 功能梯度材料的研究现状及展望[J]. 材料导报,2012(增刊1):370-373.

LI Xin, LIU Haichang, TENG Yuancheng，et al. Research status and future directions on functionally gradient materials[J]. Materials Review,2012(Sup1):370-373.

[2] 邓昊，程伟. 一种功能梯度Timoshenko梁的损伤识别方法[J]. 北京航空航天大学学报,2016,42(10):2214-2221.

DENG Hao, CHENG Wei. A damage identification method for functionally graded Timoshenko beams[J]. Journal of Beijing University of Aeronautics and Astronautics, 2016,42(10):2214-2221.

[3] YANG J, CHEN Y. Free vibration and buckling analyses of functionally graded beams with edge cracks[J]. Bulletin of the Polish Academy of Sciences Mathematics, 2008, 83(1):48-60.

[4] SIMSEK M. Fundamental frequency analysis of functionally graded beams by using different higher-order beam theories[J]. Nuclear Engineering and Design, 2010,240(4): 697-705.

[5] ALSHORBAGY A E, ELTAHER M A, MAHMOUD F F. Free vibration characteristics of a functionally graded beam by finite element method[J]. Applied Mathematical Modeling, 2011,35(1):412-425.

[6] RAVIKIRAN K, KASHIF A, GANESAN N. Static analysis of functionally graded beams using higher order shear deformation theory[J]. Applied Mathematical Modeling, 2008,32(12): 2509-2525.

[7] 原凯，王建民，韩丽，等. 能量有限元在振动与噪声预示中的研究进展[J]. 强度与环境,2015,42(3):10-19.

YUAN Kai, WANG Jianmin, HAN Li, et al. Review on prediction of vibration and noise using energy finite element analysis[J]. Structure & Environment Engineering,2015,42(3):10-19.

[8] 孙丽萍，聂武. 能量有限元法在船舶结构中的应用[J]. 哈尔滨工业大学学报,2008,40(9):1491-1494.

SUN Liping，NIE Wu. Application of energy finite element method in ship structures[J]. Journal of Harbin Institute of Technology,2008,40(9):1491-1494.

[9] 蔡忠云，唐文勇，张圣坤. 能量有限元方法在复合材料层合梁耦合结构振动分析中的应用[J]. 振动与冲击,2010,29(10):23-27.

CAI Zhongyun, TANG Wenyong, ZHANG Shengkun. Application of energy finite element method in vibration analysis of coupled composite laminated beam structures[J]. Journal of Vibration and Shock,2010,29(10):23-27.

[10] 解妙霞，陈花玲，吴九汇. 圆柱壳高频弯曲振动的能量有限元分析[J]. 西安交通大学学报,2008,42(9):1113-1116.

XIE Miaoxia, CHEN Hualing, WU Jiuhui. Energy finite elenment analysis to high-frequency bending vibration in cylindrical shell[J]. Journal of Xian Jiaotong University,2008,42(9):1113-1116.

[11] AYDIN K. Free vibration of functionally graded beams with arbitrary number of surface cracks[J]. European Journal of Mechanics - A/Solids, 2013, 42(42):112-124.

[12] 朱继清. 基于能量流有限元的耦合结构振动传递特性分析[D].济南：山东大学,2011.

[13] 杨旸. 能量有限元分析法的研究及其在复合材料层合板中的应用[D]. 宁波：宁波大学,2013.

[14] 孙丽萍. 能量有限元法研究及其应用[D]. 哈尔滨：哈尔滨工程大学,2004.

[15] 游进,李鸿光,孟光. 耦合板结构随机能量有限元分析[J]. 振动与冲击,2009,28(11):43-46.

YOU Jin，LI Hongguang，MENG Guang. Random energy finite element analysis of coupled plate structures[J]. Journal of Vibration and Shock, 2009,28(11):43-46.

[16] 徐福慧,郑世杰.基于能量有限元法的耦合梁分析[J]. 江苏航空,2012(增刊1):108-110.

XU Fuhui，ZHENG Shijie. Analysis of coupled beam based on energy finite element method[J]. Jiangsu Aviation, 2012(Sup1):108-110.

[17] 郑国涛. 能量有限元方法若干问题的研究[D]. 上海：上海交通大学,2013.