江 希 匡传树 帅 涛 耿培帅

(中国直升机设计研究所，江西景德镇300300)

梁作为各种复杂结构系统的基本组成单元，被普遍应用于航空航天、土木和机械等各种工程实践领域。引导梁理论发展的主要因素一方面为新材料的应用，另一方面为新的工程实践需求。

本文的理论及有限元方法主要研究目标为功能梯度材料 (functionally graded material，FGM) 梁，因此先对 FGM 做简单介绍。FGM 指一类组成结构和性能在材料某个 (或多个) 方向上连续或准连续变化的非均质复合材料[1]。与传统层合复合材料相比，FGM 具有更加优异可设计性[2-3]，可有效减小层间残余应力和热应力，并且可以作为传统层合结构的粘合材料来提高层间粘合强度。余莲英等[4]采用弹性力学逆解法，求得了功能梯度曲梁在端部受弯矩作用的解析解。李世荣等[5]基于Euler 梁理论给出了FGM 梁和均质梁静动态解之间的转换关系。Zhang[6]定义了 FGM 梁物理中面的明确表示式，将其引入Euler 梁理论，实现了非均质梁的拉弯解耦，给出了静动态解析解。徐华等[7]将物理中面的概念引入Timoshenko 梁理论，得到了FGM 梁与均匀梁静力解的转换关系。陈淑萍等[8]基于Timoshenko 梁理论分析了材料分布梯度对梁固有频率的影响。

除材料应用外，工程中需要更精细化的梁理论模型用以满足更高精度的需求。因此学者构造出了更多精细化的梁理论模型[9-13]。Levinson[12]最早在1981 年提出的轴向位移模式为三次函数的高次剪切变形理论。Reddy[13]变分自冾的概念引入到三次位移模式中，得到了变分自洽的高阶梁板理论，该理论也是目前被研究和应用最为广泛的高阶理论。Pei等[14-15]基于虚功和势能等价的思想，推导出了适用于复合材料的变分自洽的高阶梁理论和变分渐进的低阶梁理论。相较于Reddy 梁理论，该理论继承了传统均质材料梁理论(Euler 梁和Timoshenko 梁)中物理量 (广义位移、广义力、刚度等) 意义清晰、本构关系解耦等特点，并实现了向复合材料领域的延拓。

基于新型的高阶梁理论，Pei 等[14]仅作了静力学方面的分析。考虑将其理论应用于梁的模态分析，因此本文基于Pei 等[14]的高阶梁理论，构造了新型高阶单元，通过有限元方法研究了功能梯度材料梁的自由振动问题。

1 基本方程

考虑等高度矩形截面梁，材料的杨氏模量E、剪切模量G及泊松比ν沿z方向变化。待需求解的位移场为轴向位移ux和横向位移uz，两者均为坐标(x,z) 的函数，后续为简写不再注明自变量。

考虑梁结构上下表面剪应力自由，并忽略横向正应变，Pei 等[14]将梁的面内位移场表示为

式中，zc为梁截面中性面，为高阶翘曲函数，可由式(4)计算得到。u0为复合材料梁的拉伸，ϕ为复合材料梁的转角，w为复合材料梁的挠度，三者均仅为轴向坐标x的函数，可由梁的位移场(ux,uz)定义为

其中B0为复合材料梁的拉伸刚度，B2为复合材料梁的剪切刚度，两者与传统均质梁(Euler 梁和Timoshenko 梁) 中对应物理量的物理意义相当，它们的具体数学定义为

另外，式 (1) 的 1 式中，位移因函数 “1”、“z −zc”和均仅为坐标z的函数，且它们具有如下性质

式中gz为梁截面翘曲函数，较为常用函数包括：三次型z3、指数型(z−ze−2(z/t)3)、双曲型(z−tsinh(z/t))和正弦型 (−z+(π2t/24) sin(zπ/t))，本文中算例均选用三次型翘曲函数。c 和λ为常数，它们分别为

由式 (1) 可得，梁截面上的正应力和剪应力分别为

定义复合材料梁的轴力N，弯矩M，高阶弯矩P，剪力Q分别为

将式(6) 代入式(7) 中，并考虑式(4) 中位移因函数与材料分布的关系，可得

式中Bs为剪切刚度，为剪切刚度系数，η为高阶弯矩对应的刚度系数，它们具体数学表示为

为后续简写，将梁的广义位移、广义应力和广义应变表示为

其中L为微分算子矩阵，具体形式为

至此，已表示完新型高阶功能梯度梁理论中广义应力、广义应变及广义位移之间的关系。为后续推导有限元方程，这里给出新型梁理论的应变能、外力势能及动能表示式。

梁的应变能Πε可表示为

以载荷向量q={N∗ M∗ q∗}T为例，梁的外力势能可表示为

根据式(1)，梁的动能Πv可表示为

2 有限元方程

本小节将根据上一小节的理论构造2 节点8 自由度的C1 型高阶梁单元Beam-H2。

将梁的广义位移场u由节点位移向量插值表示为

其中N为形函数矩阵，具体形式为

梁单元局部坐标与全局坐标转换关系为

其中B为几何矩阵，可由N形函数矩阵表示为

将式(25) 代入式(15)，则梁的应变能可离散表示为

其中下标 “n” 表示第n个单元，Ne代表单元个数。Ke为单元刚度矩阵，具体形式为

其中Le为单元积分域。

根据式(17) 和式(19)，梁的动能可离散表示为

其中Λ矩阵为

式(29) 中Me表示单元质量矩阵，具体为

将式 (19) 代入式 (15)，则梁的外力势能可离散表示为

对于动力学分析，可使用哈密顿原理进行方程推导，该原理数学表示为

将式 (26)、式 (28) 及式 (32) 代入式 (34) 得

将式(35) 所表示的动力学平衡方程简写为

其中K为总体刚度矩阵，M为总体质量矩阵，F为总体外载荷矩阵，它们分别由单元刚度阵Ke、单元质量阵Me和单元载荷矩阵Fe组装得到。

对于自由振动分析，略去外载荷F，并假设模态位移为可得

此时自由振动问题即转化为求解K和M广义特征值和特征向量的问题。至此可用于功能梯度梁模态分析的Beam-H2 单元已构造完毕。

3 数值算例

3.1 均质梁模态分析

本小节将使用Beam-H2 单元以如图1 所示的悬臂梁为例做模态分析作为验证性算例。

图1 悬臂梁的几何描述

梁长度L= 0.5 m，厚度t= 0.05 m，杨氏模量E= 210 GPa，剪切模量G= 80.77 GPa, 材料密度ρ=7850 kg/m2。本算例参考文献[16-17] 做对比验证。均质悬臂梁的前五阶固有频率计算结果如表 1 所示。

表1 悬臂梁前五阶圆频率

对比段铁城文献 [17] 和文献 [16] 以及 Timoshenko 梁理论的结果，Beam-H2 与它们都非常接近，即验证本文中的理论及有限元方法是正确的。

图2 对比了悬臂梁前五阶固有振型，因不同理论所对应的振型并没有显著差异，所以这里仅给出了Beam-H2 的振型图。

图2 悬臂梁的前五阶固有振型

3.2 FGM 梁模态分析

在3.1 节中通过算例验证了本文的Beam-H2 单元，本小节将利用Beam-H2 单元进行FGM 梁模态分析，梁结构仍如图1 所示。L=1 m，t=0.1 m。

功能梯度材料为 Al-ZrO2， 其中铝的密度为2700 kg/m3，杨式模量为 70 GPa；二氧化锆的密度为 5850 kg/m3，杨氏模量为 200 GPa。材料性能分布公式[18-19]为

其中k为材料分布参数，t为梁厚度。P代表材料性能，如：杨式模量、剪切模量、泊松比或材料密度等,下标m 代表金属，下标c 代表陶瓷。

定义无量纲频率为

表 2 悬臂 FG 梁前三阶圆频率

表3 悬臂FG 梁前三阶无量纲固有圆频率

表3 悬臂FG 梁前三阶无量纲固有圆频率

k i 1 2 3 0 (陶瓷) 3.489 0 20.936 2 54.417 5 1 3.490 8 21.007 0 55.539 1 2 3.488 8 20.931 3 55.067 0 4 3.485 4 20.802 6 52.918 6∞ (铝) 3.489 0 20.936 2 54.417 5

根据式(40)，可得FGM 梁固有频率转换式

因为式 (40) 和式 (41) 仅是对一个算例的归纳总结，仅验证了材料分布参数k对的影响，结果可能存在一定偶然性，所以接下来将对式(40) 和式(41) 的适应性作进一步算例验证分析。包括不同材料组分性能对的影响和长高比对的影响。

表 4 给定不同 r，s 值时所对应的 e1 (单位：%)

由表 5 不同长高比的相对误差可得，长细比较小(L/t<10) 时，长细比的变化对无量纲频率影响越大，且频率阶次越高，影响越为显著；而长细比较大(L/t >20) 时，长细比的变化对无量纲频率影响较为微弱，若以(e2|L/t=80−e2|L/t=20) 的绝对值作为评判，则1 阶偏差在 0.3%以内，2 阶偏差2%以内，3 阶偏差5%以内。

表 5 不同长高比 (L/t) 所对应的 e2 (单位：%)

4 结论

(1)转换两者的梁结构长细比相同，可作3 阶以内频率转换。

(2) 转换两者的梁结构长细比均较大 (L/t >20)，可作3 阶以内频率转换。

(3) 若仅作基频转换，转换两者长细比L/t >6即可。