ABSTRACT

In this paper, the free vibration characteristics of functionally graded plates in arbitrary shape were studied based on the improved Rayleigh Ritz method. Assuming that the material of the functionally graded plate changes exponentially along the thickness direction, the vibration of the plate is described by Mindlin plate theory and the improved Fourier series is used as the displacement tolerance function. The energy functional expressions of the functionally graded plate are derived and the natural frequency is obtained by solving them. The convergence and accuracy of the proposed method are verified by compared with the existing literature. Then the proposed method was applied to the triangular and circular plates respectively, in which the versatility of the proposed method was showing. Finally the influences of gradient index, boundary condition and thickness of the functionally graded plate on the free vibration are discussed.

Keywords:Arbitrary Shape, Functionally Graded, Typical Boundary Conditions, Rayleigh-Ritz Method, Improved Fourier Series Method

任意形状功能梯度板的自由振动分析

玉云艳^1,2,3*，朱翔^1,2,3，李天匀^1,2,3，郭文杰^1,4#

¹华中科技大学船舶与海洋工程学院，湖北 武汉

²船舶与海洋水动力湖北省重点实验室，湖北 武汉

³高新船舶与深海开发装备协同创新中心，上海

⁴华东交通大学，铁路环境振动与噪声教育部工程研究中心，江西 南昌

收稿日期：2019年2月16日；录用日期：2019年3月6日；发布日期：2019年3月13日

摘 要

本文基于改进的Rayleigh-Ritz法对任意形状功能梯度材料板的自由振动特性进行了分析。假设功能梯度板的材料沿厚度方向指数变化，用Mindlin板理论描述板的振动。引入改进的Rayleigh-Ritz法，利用改进的傅里叶级数作为位移容许函数，结合Mindlin板理论推导得到功能梯度板的能量泛函表达式，根据能量泛函变分原理得到了振动系统的特征方程，求解得到功能梯度板的振动固有频率。通过与文献中功能梯度板的固有频率值对比，验证了本文方法的收敛性与准确性。然后分别对三角形板和圆形板的固有频率进行了求解，表明本文方法的通用性。然后通过算例分析讨论了功能梯度板的梯度指数、边界条件及板厚等参数对固有频率的影响。

关键词 :任意形状，功能梯度，典型边界条件，Rayleigh-Ritz法，改进的傅里叶级数法

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

日本材料学家新野正之，平井敏雄和渡边龙三等人在20世纪80年代中后期提出了功能梯度材料的概念 [1] [2] 。功能梯度材料选择两种不同性能的材料，根据使用要求采用材料复合技术，使中间的组成呈现连续的梯度变化而内部不存在明显的界面，材料的性质和功能沿厚度方向呈梯度变化。这样就可以从材料的本身改善材料的各种性能。由于功能梯度材料具有良好的力学性能，日益成为研究者研究的重点 [3] 。

针对功能梯度材料板壳的振动问题，已经有大量文献开展了研究。王明禄 [4] 考虑了由材料的非均匀性引起的中面应变，求解了固有频率。曹志远 [5] 利用特殊类型各项异性及薄膜弯曲耦联的二维常系数方程研究功能梯度材料圆柱壳的固有频率的解析解。Henry Khov [6] 拓展了选用的位移函数，使其能够精准的满足边界点的控制微分方程。徐坤 [7] 假设剪应力沿厚度方向呈抛物线分布，利用剪切余能与剪切应变能相等原理与虚位移原理及Navier方法研究平板的振动。梁斌 [8] 利用Rayleigh-Ritz方法建立了功能梯度圆柱壳自由振动固有频率的特征方程，其中材料参数为指数型体积分数，推导了一端固定一端自由和两端简支两种基本边界情况下的圆柱壳的自振特性。李伟 [9] 基于哈密顿原理，采用了动力变分方程，并把任意四边边界条件转化为统一的振型参数，对正交各向异性功能梯度开口圆柱壳的固有频率与模态进行求解。陈淑萍 [10] 基于一阶剪切变形理论，求解了两端简支Timoshenko梁自由振动的固有频率，并研究了材料梯度指数对梁动力响应的影响。尹硕辉 [9] 基于一阶剪切变形理论，应用等几何有限元法对功能梯度板进行研究。此外还有大量的数值方法已经得到应用 [12] [13] [14] 。以上研究多是应用解析方法对单一形状的功能梯度板壳的研究，或者是应用数值方法，对于应用解析方法对任意形状功能梯度板壳问题的研究成果还非常少。

本文采用改进的Rayleigh-Ritz法，基于一阶剪切理论对功能梯度材料板的自由振动进行研究。考虑材料沿厚度方向梯度指数变化的功能梯度板，引入改进的Rayleigh-Ritz法，利用改进的傅里叶级数和Mindlin板理论推导得到功能梯度板的能量泛函表达式，求解得到功能梯度板的振动固有频率。通过数值算例的计算并和文献对比证明本方法的收敛性与准确性，然后对其他形状板的自由振动进行了求解。

2. 理论模型

2.1. 功能梯度材料

对于功能梯度材料板结构，通常研究其材料属性沿厚度方向的变化。采用空间连续函数模型，则由两种材料组成的功能梯度材料沿厚度方向材料组分的分布规律：

$X\left(z\right)=X_m{V}_m\left(z\right)+X_c{V}_c\left(z\right)$ (1)

其中， $X\left(z\right)$ 为功能梯度材料的物理属性如密度、杨氏模量等， $$ 和 $X_c$ 为两种组分材料的物理属性， $V_m\left(z\right)$ 和 $$ 为两种材料的体积分数，且满足：

$V_m\left(z\right)+V_c\left(z\right)=1$ (2)

采用幂函数形式，则体积分数表示为：

$V_m={\left(\frac{2z+h}{2h}\right)}^k\left(k\ge 0\right)$ (3)

式中k为梯度指数。代入式(1)可得

$X\left(z\right)=\left(X_m-X_c\right){\left(\frac{2z+h}{2h}\right)}^k+X_c$ (4)

2.2. 改进的Rayleigh-Ritz法

改进的Rayleigh-Ritz法与传统的Rayleigh-Ritz法相比的改进之处在于对边界条件的处理以及位移容许函数的选取方法。

2.2.1. 函数域的拓展

对Rayleigh-Ritz法中的函数域进行扩展，模型如图1。

图1. 任意形状板的计算模型

S为薄板的求解域，S_a为试函数的积分域，长度分别为L_x和L_y。当求解域S为任意的不规则形状时，可以假想将结构域扩展为一个略大于真实域的矩形域S_a，矩形域要涵盖待求取结构的真实求解域。一般来说，取求解域距x、y轴最大距离作为矩形域的边界。当采用改进的傅里叶级数作为结构的试函数时，使用矩形域能更便捷地构造出合理的试函数。而在求板的应变能、动能时，只需对真实求解域进行离散并积分，然后求取泛函极值得到振动方程。

2.2.2. 曲边边界的虚拟弹簧处理方法

经典的边界条件，比如固支和简支等往往不能准确模拟结构的实际边界，对于实际边界条件，可采用弹簧模型，通过在结构域的边界上采用平动弹簧和转动弹簧来模拟弹性边界。以矩形板为例施加弹簧边界，如图2所示。

图2. 矩形板物理模型

假设平动弹簧和转动弹簧刚度值分别为k_ij(N/m)和K_ij(Nm/rad)，各种复杂边界条件可通过设置两种弹簧的刚度值模拟。经典边界条件下的刚度值如表1所示。模拟弹性边界条件时则将刚度系数取为相应的值即可。

表1. 经典边界条件下的边界弹簧取值

通过弹性边界设置，这将任意边界对结构振动的影响转化为边界弹簧的弹性势能对结构总刚度的影响。

3. 公式推导

功能梯度板的振型函数可以表示为：

$\begin{array}{l}w\left(x,y,t\right)={\displaystyle \sum _{m=0}^M{\displaystyle \sum _{n=0}^N{A}_{mn}}}{f}_m\left(x\right)g_n\left(y\right){\text{e}}^{i\omega t}\\ {\beta }_x\left(x,y,t\right)={\displaystyle \sum _{m=0}^M{\displaystyle \sum _{n=0}^N{B}_{mn}}}{f}_m\left(x\right)g_n\left(y\right){\text{e}}^{i\omega t}\\ {\beta }_y\left(x,y,t\right)={\displaystyle \sum _{m=0}^M{\displaystyle \sum _{n=0}^N{C}_{mn}}}{f}_m\left(x\right)g_n\left(y\right){\text{e}}^{i\omega t}\end{array}$ (4)

其中A_mn，B_mn，C_mn为未知系数， $f_m\left(x\right)$ ， $g_n\left(y\right)$ 分别为与x、y有关的函数， ${\text{e}}^{i\omega t}$ 为时间因子，在下文中忽略。

本文方法选用改进的傅里叶级数作为试函数 [11] ，具体形式为：

$\begin{array}{l}{f}_m\left(x\right)=\sin \left(\frac{m\pi }{Lx}x\right)0<m<5\\ f_m\left(x\right)=\cos \left[\frac{\left(m-5\right)\text{π}}{Lx}x\right]m\ge 5\end{array}\}$ (5a)

$\begin{array}{l}{g}_n\left(y\right)=\sin \left(\frac{n\text{π}}{Ly}y\right)0<n<5\\ g_n\left(y\right)=\cos \left[\frac{\left(n-5\right)\text{π}}{Ly}y\right]n\ge 5\end{array}\}$ (5b)

其中 $m=1,2,3,\cdots ,M$ ； $n=1,2,3,\cdots ,N$ 。由于改进的傅里叶级数可以满足任意边界条件，使板的试函数在整个求解域内三阶导数连续且四阶导数各点均存在，可以有效克服边界处可能出现的不连续现象。

考虑到材料沿厚度不均匀分布会导致物理中面的改变，引入 $z_0$ 表示新物理中面的位置。

$$ (6)

其中 $$ 为弹性模量，随厚度变化。

考虑物理中面位置改变，由中厚板的Mindlin理论可以得到位移的形式：

$\left[\begin{array}{c}u\\ v\\ w\end{array}\right]=\left[\begin{array}{c}\left(z-z_0\right){\beta }_x\\ \left(z-z_0\right){\beta }_y\\ w\left(x,y\right)\end{array}\right]$ (7)

应力与应变的关系可以表示为：

$\left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}=\frac{E\left(z\right)}{1-{\mu }^2}\left[\begin{array}{ccc}1& \mu & 0\\ \mu & 1& 0\\ 0& 0& \frac{1-\mu }{2}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right\}$ (8)

$\left\{\begin{array}{c}{\tau }_{yz}\\ {\tau }_{xz}\end{array}\right\}=G\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left\{\begin{array}{c}{\gamma }_{yz}\\ {\gamma }_{xz}\end{array}\right\}$ (9)

其中 $G=\frac{E\left(z\right)}{2\left(1+\mu \right)}$ 为沿厚度方向变化的剪切模量。

外力与应力的关系为：

$\begin{array}{l}{M}_x=-{\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_x\left(z-z_0\right)\text{d}z}\\ M_y=-{\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_y\left(z-z_0\right)\text{d}z}\\ M_{xy}=-{\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\tau }_{xy}\left(z-z_0\right)\text{d}z}\\ Q_x={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\tau }_{xz}\left(z-z_0\right)\text{d}z}\\ Q_y={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\tau }_{yz}\left(z-z_0\right)\text{d}z}\end{array}\}$ (10)

将式(7)~(9)代入式(10)中可得：

$\left\{\begin{array}{c}{M}_x\\ M_y\\ M_{xy}\\ Q_x\\ Q_y\end{array}\right\}=\left[\begin{array}{cc}{G}_b& \\ & G_s\end{array}\right]\left\{\begin{array}{c}\frac{\partial {\beta }_x}{\partial x}\\ \frac{\partial {\beta }_y}{\partial y}\\ \frac{\partial {\beta }_x}{\partial y}+\frac{\partial {\beta }_y}{\partial x}\\ \frac{\partial w}{\partial x}-{\beta }_x\\ \frac{\partial w}{\partial y}-{\beta }_y\end{array}\right\}$ (11)

式中G_b为弯曲刚度，G_s为扭转刚度：

$G_b=-{\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E\left(z\right){\left(z-z_0\right)}^2}{1-{\mu }^2}\left[\begin{array}{ccc}1& \mu & 0\\ \mu & 1& 0\\ 0& 0& \frac{1-\mu }{2}\end{array}\right]\text{d}z}$ (12)

$G_s=k_0{\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}\left[\begin{array}{cc}G& 0\\ 0& G\end{array}\right]\text{d}z}$ (13)

式中k₀为截面校正因子，本文中取k₀ = 5/6。

设

$d^{\text{T}}=\left[\begin{array}{ccc}{A}_{mn}& B_{mn}& C_{mn}\end{array}\right]$ (14)

${\lambda }^{\text{T}}=\left[f_1\left(x\right)g_1\left(y\right),f_1\left(x\right)g_2\left(y\right),\cdots ,f_1\left(x\right)g_N\left(y\right),\cdots ,f_M\left(x\right)g_N\left(y\right)\right]$ (15)

将式(14)和(15)代入式(7)中可得：

$\left[\begin{array}{l}w\\ u\\ v\end{array}\right]=\left[\begin{array}{ccc}{\lambda }^{\text{T}}& 0& 0\\ 0& \left(z-z_0\right){\lambda }^{\text{T}}& 0\\ 0& 0& \left(z-z_0\right){\lambda }^{\text{T}}\end{array}\right]d=Pd$ (16)

则中厚板的动能为：

$T={\displaystyle {\iint }_S{P}^{\text{T}}{m}_{0}P\text{d}x\text{d}y}$ (17)

式中 $m_0$ 为单位面积密度。

$m_0={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}\rho \left(z\right)}\left[\begin{array}{ccc}z& 0& 0\\ 0& z^2& 0\\ 0& 0& z^2\end{array}\right]\text{d}z$ (18)

则中厚板的弯曲应变能可以相应地表示为：

$V_p=\frac{1}{2}{\displaystyle {\iint }_s{B}^{\text{T}}\left[\begin{array}{cc}{G}_b& 0\\ 0& G_s\end{array}\right]B\text{d}x\text{d}y}$ (19)

式中

$B=\left[\begin{array}{ccc}0& \frac{\partial \lambda }{\partial x}& 0\\ 0& 0& \frac{\partial \lambda }{\partial y}\\ 0& \frac{\partial \lambda }{\partial y}& \frac{\partial \lambda }{\partial x}\\ \frac{\partial \lambda }{\partial x}& \lambda & 0\\ \frac{\partial \lambda }{\partial y}& 0& \lambda \end{array}\right]$ (20)

采用弹簧模型处理边界条件，则三种弹簧组成的弹簧刚度系数矩阵为

$K_s=\left[\begin{array}{ccc}{k}_l& 0& 0\\ 0& K_{xl}& 0\\ 0& 0& K_{yl}\end{array}\right]$ (21)

其中k_l为平动弹簧的刚度值，k_xl为的转动弹簧刚度值，k_yl为的转动弹簧刚度值。

功能梯度板的边界弹簧势能表示为：

$V_s=\frac{1}{2}{\displaystyle {\int }_l{\lambda }^{\text{T}}{K}_s\lambda \text{d}l}$ (22)

于是系统的能量泛函为：

$\prod =V_p+V_s-T$ (23)

对未知系数求极值：

$\begin{array}{l}\frac{\partial \prod }{\partial B_{mn}}=0\\ \frac{\partial \prod }{\partial A_{mn}}=0\\ \frac{\partial \prod }{\partial C_{mn}}=0\end{array}\}$ (24)

求解式(24)并将所得结果整理化简，于是可将振动问题转化为求解特征值的问题：

$\left(K-{\omega }^{2}M\right)A=0$ (25)

其中 $\omega$ 为圆频率。

4. 数值计算

4.1. 材料属性

本节算例中考虑如下几种材料，材料属性如表2所示。考虑两种功能梯度材料，分别为Al/Al₂O₃和Al/ZrO₂。

4.2. 收敛性分析

由前文分析可知，计算过程中存在针对级数的求和，因此计算结果的准确性受到截断项数M和N的影响。此外在用弹簧模拟诸如刚性固定等边界时，实际弹簧刚度取值也不可能取无穷大，因此需要对截断项数及弹簧刚度系数的取值进行收敛性分析。

表2. 材料属性

在算例中选取了功能梯度矩形板作为计算模型，如图2所示。分别选用四边自由和四边刚性固定的条件进行收敛性分析。模型参数如下：材料为Al/Al₂O₃，长为a = 1 m，宽为b = 1 m，厚h = 0.1 m，材料梯度指数k = 1。

首先分析截断项数M和N的收敛性，考虑四边自由板，取k_l和K为0，M和N从10取到15，固有频率大小随M和N的变化如表3所示。由表3可知，当M = N = 14时，结构的各阶固有频率不再随着M和N的增加而变化，此时计算结果已收敛。

表3. 四边自由(F-F-F-F)边界条件下功能梯度矩形板固有频率随截断项数的变化(Hz)

然后将板的四边设置为固支约束，取M = N = 14，且K和k_l的值从10¹²变化到10²⁰，以验证K和k的收敛性，得到的固有频率如表4所示。表4中数据表明，当弹簧的刚度系数不断增大时，固支功能梯度矩形板的弯曲自由振动的频率逐渐增大后趋于稳定。当弹簧的刚度系数达到k_l = 10¹⁴时结果已经收敛。故可以取上述数值模拟弹簧刚度系数为无穷大时的取值。

表4. 四边固支(C-C-C-C)功能梯度矩形板固有频率随弹簧刚度值变化(Hz)

4.3. 准确性验证

本小节计算不同边界条件下，不同厚度的功能梯度矩形板的无量纲频率，将无量纲频率与已有的文献结果对比，验证本方法在求解功能梯度板的振动的准确性。

选用边长为a = 1 m，b = 1 m，厚度分别为h = 0.05、0.1、0.2 m的功能梯度矩形板，材料为Al/Al₂O₃。计算简支功能梯度矩形板的无量纲频率 $$ ，与文献 [15] 进行比对，如表5所示。

表5. 四边简支(S-S-S-S)功能梯度矩形板无量纲固有频率(材料Al/Al₂O₃, a/b = 1)

更换材料为Al/ZrO₂，计算无量纲频率 $\gamma =\omega a^2\sqrt{{\rho }_m/E_m}/h$ ，与文献 [16] 进行比对。如表6所示：

表6. 简支边界下功能梯度矩形板无量纲固有频率( $$ ，S-S-S-S，h = 0.1，a/b = 1)

表5、表6中结果表明，对各种材料、厚度的功能梯度板，本文的计算结果与文献均吻合很好，说明了本方法的准确性。

4.4. 功能梯度三角形板的弯曲自由振动特性分析

本节对功能梯度三角形板进行弯曲自由振动特性分析，板为边长分别为， $b=\sqrt{2}\text{m}$ ，c = 2 m的等腰直角三角形，h = 0.1 m，材料为Al/Al₂O₃，计算其简支边界条件下的固有频率，换算为无量纲频率，如表7。

表7. 不同梯度指数下功能梯度三角形板前六阶无量纲固有频率

4.5. 功能梯度圆形板的弯曲自由振动分析

对功能梯度圆形板进行弯曲自由振动的特性分析，圆板的几何参数为半径R = 1 m，h = 0.1 m，材料为Al/Al₂O₃，计算简支边界条件下的固有频率，换算为无量纲频率，如表8，从表中可见随着梯度指数增大，相同阶的固有频率逐渐降低。

表8. 不同梯度指数下功能梯度圆形板前六阶无量纲固有频率

4.6. 梯度指数k对功能梯度板的固有频率影响分析

本节讨论功能梯度指数、边界条件和板厚等对固有频率的影响。

首先讨论梯度指数对频率的影响与不同的边界条件的关系。选用边长a = 1，b = 1，厚度h = 0.1的功能梯度板，材料为Al/Al₂O₃，边界分别为四边固支、四边简支、一边固支三边自由。求无量纲频率，作图3。

从图3可以看出，随着梯度指数增大，首阶固有频率逐渐降低，但是降低的程度随梯度指数增大趋于平缓，且从不同边界条件的三条曲线可见，梯度指数相同时，四周固支板固有频率最高，一边固支板固有频率最低。

图3. 不同边界条件下功能梯度矩形板首阶无量纲固有频率随梯度指数变化的曲线

以下为四边简支功能梯度板为例，讨论梯度指数对无量纲频率的影响以及与固有频率阶数的关系。几何和材料参数和上一算例相同。图4为前三阶固有频率随梯度指数变化的曲线。

图4. 功能梯度矩形板前三阶无量纲固有频率随梯度指数变化的曲线，S-S-S-S

从图4中看出，功能梯度板的前三阶固有频率随梯度指数的变化趋势基本相同，随梯度指数增大，固有频率逐渐降低，这是由于梯度指数增大，板的等效刚度逐渐降低。同时也可以得出功能梯度指数对板的低阶固有频率的影响与阶数基本无关。

以下讨论梯度指数对频率的影响与板的厚度的关系。选用边长a = 1 m，b = 1 m，厚度分别为h = 0.05 m、0.1 m、0.2 m的简支功能梯度板，材料为Al/Al₂O₃，得到图5在不同板厚下的无量纲频率随梯度指数变化的曲线。

从图5中可以看出，对于不同厚度的功能梯度板，材料梯度指数的变化对其自由振动首阶频率的影

图5. 不同板厚下功能梯度矩形板首阶无量纲固有频率随梯度指数变化的曲线

响大致相同。板厚越大，首阶固有频率越高。

5. 结论

本文基于改进的Rayleigh-Ritz法研究了任意形状功能梯度板的振动。引入改进的傅里叶级数作为结构的位移容许函数，基于一阶剪切理论，计算结构的应变能和动能，利用弹簧模拟边界条件的方法对任意形状边界上的支撑条件进行处理，通过求取系统控制方程的特征值，计算得到结构的固有频率。通过算例对比验证了方法的准确性，然后计算了三角形和圆形功能梯度板的固有频率，讨论了梯度指数、边界条件等对固有频率的影响。

文章引用

玉云艳,朱 翔,李天匀,郭文杰. 任意形状功能梯度板的自由振动分析
Free Vibration Analysis of Functional Gradient Mindlin Plate of Arbitrary Shape[J]. 声学与振动, 2019, 07(01): 28-40. https://doi.org/10.12677/OJAV.2019.71004

参考文献