石墨烯增强多孔功能梯度圆板非线性自由振动

Nonlinear Free Vibration of Graphene-Reinforced Porous Functionally Graded Circular Plates

摘要: 基于Halpin-Tsai微观力学模型、闭单元高斯随机场理论和混合定律得到了石墨烯纳米片(Graphene nanoplates, GPL)增强多孔功能梯度材料(functionally graded materials, FGM)板的有效物性参数，其中考虑了GPL沿厚度方向均匀、X形和O形分布和孔隙沿厚度方向两种梯度分布。然后运用一阶剪切变形理论、Von-Karman变形理论和Hamilton原理推导出多孔GPL增强FGM材料圆板的运动方程组并进行无量纲化。最后通过Kartorovich时间平均法和微分求积法(differential quadrature method, DQM)将运动方程组离散得到一个特征方程组。通过直接迭代法求解出非线性自由振动频率，并通过模型退化和现有文献对比验证，验证了计算的正确性和有效性。讨论了边界条件、GPL重量分数、孔隙率、孔隙分布、GPL分布等参数对多孔GPL圆板非线性固有频率的影响。

Abstract: Based on the Halpin-Tsai micromechanical model, the closed-cell Gaussian random field theory, and the rule of mixtures, the effective material properties of graphene nanoplate (GPL)-reinforced porous functionally graded material (FGM) plates were obtained, considering the uniform, X-shaped, and O-shaped distributions of GPLs along the thickness direction and two types of gradient distributions of porosity along the thickness direction. Subsequently, the equations of motion for porous GPL-reinforced FGM circular plates were derived using the first-order shear deformation theory, the Von-Karman deformation theory, and Hamilton’s principle, and these equations were nondimensionalized. Finally, the equations of motion were discretized into a set of eigenvalue equations using the Karman time-averaging method and the differential quadrature method (DQM). The nonlinear free vibration frequencies were obtained through a direct iteration method. The accuracy and validity of the calculations were verified by model degeneration and comparison with existing literature. The effects of boundary conditions, GPL weight fraction, porosity, porosity distribution, and GPL distribution on the nonlinear natural frequencies of porous GPL circular plates were discussed.

1. 引言

本纳米材料是近年来新兴且快速发展的研究领域之一，如石墨烯及其衍生物石墨烯纳米片[1]，因其突出的性能已经被广泛地研究和应用到各个领域，包括电子器件、生物工程和微/纳米机电系统[2]。因此，开展关于石墨烯增强功能梯度材料结构的研究有着重要的意义，可以揭示石墨烯增强功能梯度材料结构的内部特征并为结构分析提供基础。

许多学者对其进行了广泛的研究。Arefi等[3]运用非局部弹性变形理论与正弦剪切变形理论分析了石墨烯纳米片增强功能梯度材料板在Pasternak上的自由振动行为，分析了非局部参数、层数、重量分数和Pasternak对自由振动的影响。Song等[4]分析了低含量石墨烯纳米片均匀分布在聚合物的多层复合材料板的屈曲与后屈曲，详细研究了GPL重量分数、分布形式、几何尺寸和总层数对功能梯度GPLRC板屈曲和后屈曲性能的影响。Feng [5]等运用Hamilton原理和Timoshenko梁理论研究了沿厚度方向非均匀分布的石墨烯片增强的多层聚合物复合梁的非线性自由振动。其[6]还研究了石墨烯片增强的功能梯度梁的非线性弯曲行为。研究发现，通过在聚合物基体中添加极少量的GPL作为增强体，弯曲性能显著提高。Thai等[7]还应用修正应力梯度理论和高阶剪切变形理论对多层功能梯度的石墨烯增强复合材料微板进行自由振动分析，对比了基于修正耦合应力模型、修正应力梯度模型、高阶剪切变形模型对于微板固有频率的预测差异。Gholami和Ansari [8]研究了不同边界条件下三阶剪切变形功能梯度石墨烯片增强复合材料矩形板的几何非线性简谐激励振动。分析各种参数的影响，例如GPL分布模式、重量分数、GPL纳米填料的几何形状和FG-GPLRC板的边界约束。

而在功能梯度材料的制备过程中，常常会因为制备工艺的不完善而出现孔隙。孔隙的出现会削弱功能梯度材料结构整体的强度，不少学者也对含孔隙的GPL增强功能梯度材料结构进行过研究，其中黄小林等[9]基于三参数粘弹性地基模型，建立了粘弹性地基上含孔隙石墨烯增强功能梯度板的运动方程，用伽辽金法求解其固有频率和动力响应，并通过数值算例分析了粘弹性地基参数、孔隙率、孔隙类型及石墨烯纳米片分布模式、含量等因素对自由振动和动力响应的影响。Mingnan [10]等考虑了温度依赖特性，将GPL增强功能梯度多孔复合材料板的热屈曲。还考虑了沿厚度方向的均匀、线性、非线性热载荷。讨论了孔隙率、GPL质量分数、几何构型和边界条件等参数对板热屈曲的影响。Chen等[11]基于Timoshenko梁理论、von Karman几何非线性假设和修正偶应力理论研究了具有均匀与不均匀孔隙分布下的功能梯度微梁的无量纲频率和非线性频率的比值。Anirudh [12]等对GPL增强功能梯度材料多孔梁进行非线性分析，其考虑了von Karman几何非线性，并结合有限元和直接迭代技术求解得到各种设计参数下梁的弯曲情况。

根据已有的GPL圆板的文献表明，对多孔的GPL圆板的非线性振动研究较少，本文通过Hamilton原理推导控制方程组并无量纲化，采用Kantrorovich法消去时间变量后进行DQM变换，采用直接迭代求解出圆板线性与非线性频率。最后将GPL增强功能梯度材料退化和已有文献对比，验证了计算结果的准确性，并分析讨论了GPL尺寸、GPL分布形式、层数、振幅等参数对系统频率的影响。

2. GPL圆板模型建立

考虑一个厚度为h，半径为b的圆板，其用柱坐标$\left(r,\theta ,z\right)$ 分布表示径向、环向、横向坐标，圆点O放置于圆板的几何中心处，如图1所示为本文的GPL增强圆板结构及坐标系示意图。

为了实现结构沿厚度方向梯度变化，如图1所示，在厚度方向考虑三种分布模式。其中GPL在同一层中均匀分布，在不同层中GPL重量分数按照所规定的形式变化。UD分布是指每一层的重量分数相同。FG-O是指中间层的GPL重量分数最高，并向边缘层逐级降低GPL重量分数。FG-X是指中间层的GPL重量分数最低，并向边缘层逐级增加GPL重量分数。

Figure 1. Graphene-enhanced gradient plate model and distribution model

图1. 石墨烯增强梯度板模型及分布模式

三种石墨烯分布下的体积分数$V_{GPL}$ 表达式如下

$V_{GPL}=\{\begin{array}{lll}{U}_1\hfill & \hfill & \text{UD}\hfill \\ U_2\cos \left(\frac{\pi z}{h}\right)\hfill & \hfill & \text{FG-O}\hfill \\ U_3\left(1-\cos \left(\frac{\pi z}{h}\right)\right)\hfill & \hfill & \text{FG-X}\hfill \end{array}$ (1)

式中$U_1$ 、$U_2$ 和$U_3$ 分别为均匀分布、O型分布和X型分布的最大体积分数。式中$h$ 为厚度，$z$ 为厚度方向的坐标，其中最大体积分数可由下式确定。

${\displaystyle {\int }_{-0.5h}^{0.5h}{V}_{GPL}\left[1-{\theta }_m\phi \left(z\right)\right]\text{d}z}={\displaystyle {\int }_{-0.5h}^{0.5h}\frac{g_{GPL}\left[1-{\theta }_m\phi \left(z\right)\right]}{g_{GPL}+\left({\rho }_{GPL}/{\rho }_M\right)\left(1-g_{GPL}\right)}\text{d}z}$ (2)

式中：${\rho }_{GPL}$ 和${\rho }_M$ 分别为石墨烯和基体材料的密度，$g_{GPL}^{}$ 为GPL重量分数。上式中涉及到了孔隙相关参数，本文考虑两种孔隙分布，如图2所示。

Figure 2. Two modes of pore distribution

图2. 两种孔隙的分布模式

不同孔隙分布下的杨氏模量和密度$E\left(z\right),\rho \left(z\right)$ 随密度变化的关系可由下式表示

$\{\begin{array}{l}E\left(z\right)=E_{\max }\left[1-\theta \phi \left(z\right)\right]\\ \rho \left(z\right)={\rho }_{\max }\left[1-{\theta }_m\phi \left(z\right)\right]\end{array}$ (3)

式中：

$\theta =1-\frac{E_{\max }}{E_{\min }},{\theta }_m=\frac{1.121\left(1-\sqrt[2.3]{1-\theta \phi \left(z\right)}\right)}{\phi \left(z\right)}$ (4)

$\phi \left(z\right)=\{\begin{array}{l}\cos \left(\frac{\pi z}{h}\right)V_1分布\\ \cos \left(\frac{\pi z}{2h}+\frac{\pi }{4}\right)V_2分布\end{array}$ (5)

式中：$\theta$ 为孔隙率系数，$E_{\max },E_{\min }$ 为杨氏模量的最大值与最小值，${\rho }_{\max }$ 为密度的最大值。对于其详细的计算方法可由下面给出。其中最大杨氏模量可以通过Halpin-Tsai模型计算给出[13]：

$E_{\max }=\left(\frac{3}{8}\frac{1+{\xi }_L{\eta }_L{V}_{GPL}}{1-{\eta }_L{V}_{GPL}}+\frac{5}{8}\frac{1+{\xi }_W{\eta }_W{V}_{GPL}}{1-{\eta }_W{V}_{GPL}}\right)E_M$ (6)

式中：$V_{GPL}^{\left(k\right)}$ 为k层的体积分数，$E_M$ 、$E_{GPL}$ 分别为基体和石墨烯的弹性模量[14]，${\eta }_W$ 、${\eta }_L$ 、${\xi }_M$ 和${\xi }_L$ 的表达式分别如下：

${\eta }_W=\frac{\left(E_{GPL}/E_M\right)-1}{\left(E_{GPL}/E_M\right)+{\xi }_M}\text{}\text{}\text{}\text{}\text{}\text{}\text{}$ ，${\eta }_L=\frac{\left(E_{GPL}/E_M\right)-1}{\left(E_{GPL}/E_M\right)+{\xi }_L}\text{}\text{}\text{}\text{}\text{}\text{}\text{}$

${\xi }_M=2\left(w_{GPL}/h_{GPL}\right)$ ，${\xi }_L=2\left(L_{GPL}/h_{GPL}\right)$

式中：$l_{GPL},w_{GPL},h_{GPL}$ 分别为GPL的长度、宽度和高度。最大密度与泊松比可通过混合率模型由下式给出：

${\rho }_{\max }={\rho }_{GPL}{V}_{GPL}+{\rho }_m\left(1-V_{GPL}\right)$ (7)

$v_{\max }=v_{GPL}{V}_{GPL}+v_m\left(1-V_{GPL}\right)$ (8)

其中泊松比随厚度方向变化的关系可由下式给出：

$v\left(z\right)=0.221\left(1-\frac{\rho \left(z\right)}{{\rho }_{\max }}\right)+v_{\max }\left(1-1.121\left(1-\frac{\rho \left(z\right)}{{\rho }_{\max }}\right)+0.342{\left(1-\frac{\rho \left(z\right)}{{\rho }_{\max }}\right)}^2\right)$ (9)

3. 基本方程

考虑圆板变形为轴对称的，并根据一阶剪切变形理论，其位移场可表示为

$u_r\left(r,z,t\right)=u\left(r,t\right)+z\varphi \left(r,t\right)$ (10)

$u_z\left(r,z,t\right)=w\left(r,t\right)$ (11)

式中：$u\left(r,t\right)$ 和$w\left(r,t\right)$ 分别为圆板几何中面上的径向位移和横向位移，$u_r\left(r,z,t\right)$ ，$u_z\left(r,z,t\right)$ 分别表示为圆板任意一点的径向位移和横向位移，$\varphi \left(r,t\right)$ 表示为中间平面法线的转角。

基于薄板的Von-Karman变形理论，圆板的几何非线性应变位移关系为

${\epsilon }_r=\frac{\text{d}u}{\text{d}r}+\frac{1}{2}{\left(\frac{\text{d}w}{\text{d}r}\right)}^2-z\frac{{\text{d}}^{2}w}{\text{d}{r}^2}$ (12)

${\epsilon }_{\theta }=\frac{u}{r}+\frac{z}{r}\frac{\text{d}w}{\text{d}r}$ (13)

式中：${\epsilon }_r$ 、${\epsilon }_{\theta }$ 分别为径向、环向应变。

物理方程为：

${\sigma }_r=\frac{E^{\left(k\right)}}{1-v^{\left(k\right)}{}^2}\left({\epsilon }_r+v^{\left(k\right)}{\epsilon }_{\theta }\right)$ (14)

${\sigma }_{\theta }=\frac{E^{\left(k\right)}}{1-v^{\left(k\right)}{}^2}\left({\epsilon }_{\theta }+v^{\left(k\right)}{\epsilon }_r\right)$ (15)

式中：${\sigma }_r$ 、${\sigma }_{\theta }$ 分别为径向和环向应力。

圆板的轴力与弯矩可由物理方程积分得到

$\left[\begin{array}{c}{N}_r\\ N_{\theta }\end{array}\right]={\displaystyle \sum _{k=1}^N{\displaystyle {\int }_{z_k}^{z_{k+1}}{\left[\begin{array}{c}{\sigma }_r\\ {\sigma }_{\theta }\end{array}\right]}_k\text{d}z}}=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_r^0\\ {\epsilon }_{\theta }^0\end{array}\right\}$ (16)

$\left[\begin{array}{c}{M}_r\\ M_{\theta }\end{array}\right]={\displaystyle \sum _{k=1}^N{\displaystyle {\int }_{z_k}^{z_{k+1}}{\left[\begin{array}{c}{\sigma }_r\\ {\sigma }_{\theta }\end{array}\right]}_k\cdot z^2\text{d}z}}=\left[\begin{array}{cc}{D}_{11}& D_{12}\\ D_{21}& D_{22}\end{array}\right]\left\{\begin{array}{c}-\frac{{\partial }^{2}w}{\partial r^2}\\ -\frac{1}{r}\frac{\partial w}{\partial r}\end{array}\right\}$ (17)

式中：$A_{11}=A_{22},A_{12}=A_{21},D_{11}=D_{22},D_{12}=D_{21}$ 且：

$\left(A_{11},D_{11}\right)={\displaystyle \sum _{k=1}^N{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E^{\left(k\right)}}{1-v^{{\left(k\right)}^2}}\left(1,z^2\right)\text{d}z}}$ (18)

$\left(A_{12},D_{12}\right)={\displaystyle \sum _{k=1}^N{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E^{\left(k\right)}{v}^{\left(k\right)}}{1-v^{{\left(k\right)}^2}}\left(1,z^2\right)\text{d}z}}$ (19)

系统动能T的变分为

$\delta T=\frac{1}{2}\cdot \delta \left[{\displaystyle {\int }_{\Omega }\rho }\left(z\right)\cdot \left(u_r^2+u_z^2\right)\text{d}\Omega \right]={\displaystyle \int {\displaystyle {\int }_A\left[\left(I_{0}u-I_1\phi \right)\cdot \delta u+I_{0}w\delta w+\left(-I_{1}u+I_{2}w\right)\delta \phi \right]\text{d}A}}$ (20)

式中：$\left(I_0,I_1,I_2\right)={\displaystyle \sum _{k=1}^N{\displaystyle {\int }_{z_k}^{z_{k+1}}{\rho }^{\left(k\right)}\left(z\right)\left(1,z,z^2\right)\text{d}z}}$ 。

系统应变能U的变分为

$\delta U={\displaystyle {\int }_V\left({\sigma }_r\delta {\epsilon }_r+{\sigma }_{\theta }\delta {\epsilon }_{\theta }\right)\text{d}V}$ (21)

根据薄板理论略去影响很小的面内惯性和耦合惯性项，并由Hamilton原理推导出石墨烯增强圆板的运动方程组为：

$A_{11}\left(r\frac{{\partial }^{2}u}{\partial r^2}+\frac{\partial u}{\partial r}-\frac{u}{r}+r\frac{{\partial }^{2}w}{\partial r^2}\frac{\partial w}{\partial r}+\frac{1}{2}{\left(\frac{\partial w}{\partial r}\right)}^2\right)-A_{12}\frac{1}{2}{\left(\frac{\partial w}{\partial r}\right)}^2=0$ (22)

$\begin{array}{l}{D}_{11}\left(\frac{{\partial }^{4}w}{\partial r^4}+\frac{2}{r}\frac{{\partial }^{3}w}{\partial r^3}-\frac{1}{r^2}\frac{{\partial }^{2}w}{\partial r^2}+\frac{1}{r^3}\frac{\partial w}{\partial r}\right)-A_{11}\left[\frac{\partial u}{\partial r}\frac{{\partial }^{2}w}{\partial r^2}+\frac{1}{2}{\left(\frac{\partial w}{\partial r}\right)}^2\frac{{\partial }^{2}w}{\partial r^2}+\frac{u}{r^2}\frac{\partial w}{\partial r}\right]\\ -A_{12}\left[\frac{u}{r}\frac{{\partial }^{2}w}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}\frac{\partial w}{\partial r}+\frac{1}{2r}{\left(\frac{\partial w}{\partial r}\right)}^2\frac{\partial w}{\partial r}\right]+I_0\frac{{\partial }^{2}w}{\partial t^2}-P\left(\frac{{\partial }^{2}w}{\partial r^2}+\frac{1}{r}\frac{\partial w}{\partial r}\right)=0\end{array}$ (23)

运动控制方程可以通过以下式子进行无量纲化

$\begin{array}{l}x=\frac{r}{b},\overline{w}=\frac{w}{h},\overline{u}=\frac{b}{h^2}u,p=\frac{P{b}^2}{D_{11}},F_1=\frac{A_{11}{h}^2}{D_{11}},\\ \overline{t}=t{\left(\frac{{\rho }_{m}h{b}^4}{D_m}\right)}^{-1/2},F_2=\frac{A_{12}{h}^2}{D_{11}},F_3=\frac{D_m{I}_0}{D_{11}{\rho }_{m}h},F_4=\frac{A_{12}}{A_{11}}\end{array}$

式中：$D_m$ 和${\rho }_m$ 分别为纯基底材料时的弯曲刚度与密度。将上述无量纲量代入运动方程组中可以得到：

$\frac{{\partial }^{2}f\left(x\right)}{\partial x^2} + \frac{1}{x}\frac{\partial f\left(x\right)}{\partial x} -\frac{f\left(x\right)}{x^2}+ \frac{{\partial }^{2}g\left(x\right)}{\partial x^2} \frac{\partial g\left(x\right)}{\partial x} +\frac{1}{2x}{\left(\frac{\partial g\left(x\right)}{\partial x}\right)}^2-F_4\frac{1}{2x}{\left(\frac{\partial g\left(x\right)}{\partial x}\right)}^2=0$ (24)

$\begin{array}{l}\frac{{\partial }^{4}w}{\partial x^4}+\frac{2}{x}\frac{{\partial }^{3}w}{\partial x^3}-\frac{1}{x^2}\frac{{\partial }^{2}w}{\partial x^2}+\frac{1}{x^3}\frac{\partial w}{\partial x}-F_1\left[\frac{\partial u}{\partial x}\frac{{\partial }^{2}w}{\partial x^2}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2\frac{{\partial }^{2}w}{\partial x^2}+\frac{u}{x^2}\frac{\partial w}{\partial x}\right]\\ -F_2\left[\frac{u}{x}\frac{{\partial }^{2}w}{\partial x^2}+\frac{1}{x}\frac{\partial u}{\partial x}\frac{\partial w}{\partial x}+\frac{1}{2x}{\left(\frac{\partial w}{\partial x}\right)}^2\frac{\partial w}{\partial x}\right]+F_3\frac{{\partial }^{2}w}{\partial t^2}-p\left(\frac{{\partial }^{2}w}{\partial r^2}+\frac{1}{r}\frac{\partial w}{\partial r}\right)=0\end{array}$ (25)

对于自由振动问题，我们忽略外力的影响，采用假设时间模态法求解非线性振动问题，其假设非线性振动的形式为谐振动，然后利用Kartorovich时间平均法将方程中时间变量消去。对于石墨烯增强材料圆板的自由振动，可令：

$u\left(x,t\right)=f\left(x\right)\cdot {\sin }^2\left(\lambda \cdot t\right)$ (26)

$w\left(x,t\right)=g\left(x\right)\cdot \sin \left(\lambda \cdot t\right)$ (27)

式中：$f\left(x\right)$ 和$g\left(x\right)$ 分别为关于u，w的模态函数，$\lambda$ 为圆板的固有频率。

将上式代入无量纲控制方程，并运用Kartorovich时间平均法：

${\displaystyle {\int }_0^{2\pi /\lambda }\left\{R_1\left(x,t\right),R_2\left(x,t\right)\right\}}\sin \lambda t\text{d}t=0$ (28)

则得到石墨烯梯度增强材料圆板简化后无量纲非线性自由振动的控制方程组为

$\frac{{\text{d}}^{2}f\left(x\right)}{{\text{d}}^{}{x}^2} + \frac{1}{x}\frac{\text{d}f\left(x\right)}{{\text{d}}^{}x} -\frac{f\left(x\right)}{x^2}+ \frac{{\text{d}}^{2}g\left(x\right)}{{\text{d}}^{}{x}^2} \cdot \frac{\text{d}g\left(x\right)}{{\text{d}}^{}x} +\frac{1}{2x}{\left(\frac{\text{d}g\left(x\right)}{{\text{d}}^{}x}\right)}^2-\frac{F_4}{2x}{\left(\frac{\text{d}g\left(x\right)}{{\text{d}}^{}x}\right)}^2=0$ (29)

$\begin{array}{l}\frac{{\text{d}}^{4}g\left(x\right)}{\text{d}{x}^4}+\frac{2}{x}\frac{{\text{d}}^{3}g\left(x\right)}{\text{d}{x}^3}-\frac{1}{x^2}\frac{{\text{d}}^{2}g\left(x\right)}{\text{d}{x}^2}+\frac{1}{x^3}\frac{\text{d}g\left(x\right)}{\text{d}x}\\ -\frac{3}{4}{F}_1\left[\frac{\text{d}f\left(x\right)}{\text{d}x}\frac{{\text{d}}^{2}g\left(x\right)}{\text{d}{x}^2}+\frac{1}{2}{\left(\frac{\text{d}g\left(x\right)}{\text{d}x}\right)}^2\frac{{\text{d}}^{2}g\left(x\right)}{\text{d}{x}^2}+\frac{f\left(x\right)}{x^2}\frac{\text{d}g\left(x\right)}{\text{d}x}\right]\\ -\frac{3}{4}{F}_2\left[\frac{f\left(x\right)}{x}\frac{{\text{d}}^{2}g\left(x\right)}{\text{d}{x}^2}+\frac{1}{x}\frac{\text{d}f\left(x\right)}{\text{d}x}\frac{\text{d}g\left(x\right)}{\text{d}x}+\frac{1}{2x}{\left(\frac{\text{d}g\left(x\right)}{\text{d}x}\right)}^2\frac{\text{d}g\left(x\right)}{\text{d}x}\right]-F_3{\lambda }^{2}g\left(x\right)=0\end{array}$ (30)

4. DQM数值计算

运用微分求积法可以将连续函数的定义域分解为离散的节点，考虑在$a\le x\le b$ 上的函数$F\left(x\right)$ ，$F\left(x\right)$ 在$x_i$ 处的k阶导数可以近似为所有离散点的加权线性和。从而将微分方程转换为代数方程表示。

$\begin{array}{l}\frac{{\text{d}}^{k}F\left(x_i\right)}{\text{d}{x}^k}={\displaystyle \sum _{j=1}^N{c}_{ij}^{\left(k\right)}F\left(x_j\right)}\\ i=1,2,\cdots ,N-1,N;k=1,2,\cdots ,N-2,N-1\end{array}$ (31)

其中N为节点总数，$c_{ij}^{\left(k\right)}$ 为第k阶加权系数，其数值与节点选取有关，因为非均匀节点相较均匀节点有着更好的收敛性，本文节点选择Chebyshev-Gauss-Lobatto分布，如下所示

$x_i=\frac{1}{2}\left\{1-\cos \left(\frac{i-1}{N-1}\right)\pi \right\},i=1,2,\cdots ,N$ (32)

加权系数$c_{ij}^{\left(k\right)}$ 可由以下式子得到

$C_{ij}^{\left(1\right)}=\{\begin{array}{l}\frac{{\displaystyle \prod _{\begin{array}{l}k=1\\ k\ne i,j\end{array}}^N\left(x_i-x_k\right)}}{\left(x_i-x_j\right){\displaystyle \prod _{\begin{array}{l}k=1\\ k\ne j\end{array}}^N\left(x_j-x_k\right)}}i\ne j\\ {\displaystyle \sum _{\begin{array}{l}k=1\\ k\ne i\end{array}}^N\frac{1}{x_i-x_k}}i=j\end{array}$ (33)

$C_{ij}^{\left(k\right)}={\displaystyle \sum _{k=1}^N{C}_{ik}^{\left(1\right)}}{C}_{ik}^{\left(k-1\right)}\left(k=2,3,\cdots ,N-1\right)$ (34)

基于此原理，将控制方程组进行DQM变换。

$\begin{array}{l}{\displaystyle \sum _{j=1}^N{C}_{ij}^2\cdot f\left(x_j\right)}+\frac{1}{x_i}{\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot f\left(x_j\right)}-\frac{f\left(x_i\right)}{x_i^2}+{\displaystyle \sum _{j=1}^N{C}_{ij}^2\cdot g\left(x_j\right)}\cdot {\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot g\left(x_j\right)}\\ +\frac{1}{2x}{\left({\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot g\left(x_j\right)}\right)}^2-\frac{F_4}{2x}{\left({\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot g\left(x_j\right)}\right)}^2=0\end{array}$ (35)

$\begin{array}{l}{\displaystyle \sum _{j=1}^N{C}_{ij}^4\cdot g\left(x_j\right)}+\frac{2}{x_i}{\displaystyle \sum _{j=1}^N{C}_{ij}^3\cdot g\left(x_j\right)}-\frac{1}{x_i^2}{\displaystyle \sum _{j=1}^N{C}_{ij}^2\cdot g\left(x_j\right)}+\frac{1}{x_i{}^3}{\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot g\left(x_j\right)}-\frac{3}{4}{F}_1[{\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot g\left(x_j\right)}\cdot {\displaystyle \sum _{j=1}^N{C}_{ij}^2\cdot g\left(x_j\right)}\\ +\frac{1}{2}{\left({\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot g\left(x_j\right)}\right)}^2\cdot {\displaystyle \sum _{j=1}^N{C}_{ij}^2\cdot g\left(x_j\right)}+\frac{1}{x_i}f\left(x_i\right)\cdot {\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot g\left(x_j\right)}]-\frac{3}{4x}{F}_2+[{\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot f\left(x_j\right)}\cdot {\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot g\left(x_j\right)}\\ +f\left(x_i\right)\cdot {\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot g\left(x_j\right)}+\frac{1}{2}{\displaystyle \sum _{j=1}^N{\left(C_{ij}^1\cdot g\left(x_j\right)\right)}^2\cdot {\displaystyle \sum _{j=1}^N{C}_{ij}^1\cdot g\left(x_j\right)}}]-F_3\cdot {\lambda }^2\cdot g\left(x_j\right)=0\end{array}$ (36)

同样的，对FGM圆板的边界也需要进行DQM变换

夹紧边界：

$u_N=0,w_N=0,{\displaystyle \sum _{j=1}^N{C}_{Nj}^1{w}_j}=0$ (37)

简支边界：

$u_N=0,w_N=0,{\displaystyle \sum _{j=1}^N{C}_{Nj}^2{w}_j}+v\cdot {\displaystyle \sum _{j=1}^N{C}_{Nj}^1{w}_j}=0$ (38)

园板中心连续条件：

$u_1=0,{\displaystyle \sum _{j=1}^N{C}_{1j}^1{w}_j}=0,{\displaystyle \sum _{j=1}^N{C}_{1j}^3{w}_j}+\frac{1}{x_1}{\displaystyle \sum _{j=1}^N{C}_{1j}^3{w}_j}=0$ (39)

对于自由振动的DQM控制方程组，其可以看成以下形式的特征方程

$\left[K\right]\left[W\right]+\left[F\left\{W\right\}\right]\left[W\right]={\lambda }^2\left[W\right]$ (40)

其中${\left[W\right]}^{\text{T}}=\left\{u_1,u_2,\cdots ,u_{N-1},u_N,w_1,w_2,\cdots ,w_{N-1},w_N\right\}$ 对于上述非线性方程，可以通过以下步骤迭代求解，从而得到无量纲非线性频率解及模态。

(一) 首先将非线性项设为零，求解简化后的线性特征值问题。得到线性频率${\omega }_L$ 和其对应的特征向量。

(二) 放大特征向量，使圆板最大横向位移等于给定的振幅$w_{\max }$ ，并利用此特征向量计算(40)中的非线性项

(三) 将步骤二得到的非线性项带入(40)，求解出新的最大横向位移，如最大横向位移与$w_{\max }$ 之间的相对误差大于0.01%，则返回步骤二。并持续迭代直至达到规定误差值。

5. 数值结果与分析

通过编写MATLAB程序可以得到由微分求积法求解多孔GPL增强圆板的无量纲固有频率。本文选取环氧树脂作为基体材料，其弹性模量$E_m=3\text{Gpa}$ ，密度为${\rho }_m=1200\text{kg}/{\text{m}}^{\text{3}}$ ，泊松比$v_m=0.34$ 石墨烯作为加强材料，其弹性模量$E_g=1.01\text{Tpa}$ ，其密度${\rho }_g=1062\text{kg}/{\text{m}}^{\text{3}}$ ，泊松比$v_g=0.175$ ，石墨烯片几何尺寸$l_{GPL}=2.5\text{um}$ 、$b_{GPL}=1.5\text{um}$ 、$h_{GPL}=1.5\text{nm}$ ，石墨烯重量分数$g=0.5\%$ 。石墨烯增强板层数$N_m=50$ 层。

为了验证计算结果的正确性，本文计算了夹紧与简支情况下各向同性圆板的小振幅振动频率并与文献对比，如表1所示，其结果与文献吻合良好。

Table 1. Circular plate small amplitude vibration first three frequencies

表1. 圆板小振幅振动前三阶频率

无量纲频率	夹紧			简支
	一阶频率	二阶频率	三阶频率	一阶频率	二阶频率	三阶频率
本文	10.2158	39.7711	89.1041	4.9351	29.7200	74.1561
[15]	10.2158	39.7711	89.1041	4.9351	29.7200	74.1561
[16]	10.215	39.771	89.104	4.935	29.720	74.156

图3和图4分析了在小振幅和大振幅两种情况下重量分数对圆板第一阶频率的影响。其中图3考虑了$V_1$ 孔隙分布模式下，孔隙率$\theta =0.3$ 的情况，图4考虑了$V_2$ 孔隙分布模式下，孔隙率$\theta =0.3$ 的情况。观察圆板随着GPL重量分数变化的关系，两种孔隙分布和三种不同的石墨烯分布下，频率都表现出随着GPL重量分数增加而增加的趋势。对比不同的分布模式，可以发现相对于$V_1$ 孔隙分布，$V_2$ 孔隙分布的无量纲第一阶固有频率稍微低于$V_1$ 孔隙分布。对比不同分布模式下GPL重量分数对频率的影响，可以发现FG-O分布下的频率最低，UD分布下的频率略高于FG-O分布，FG-X分布下频率最高。仔细观察不同分布模式下频率的增长速率，可以发现FG-X分布下频率随GPL重量分数增长的速率最大，UD分布次之，FG-O分布下频率随GPL重量分数增长速率最低。对比振幅变化，可以发现三种GPL分布形式下和两种孔隙分布形式下，$w_{\max }/h=0.8$ 时的频率始终小于$w_{\max }/h=0.01$ 。观察GPL重量分数–频率曲线可以看出最大振幅对GPL重量分数–频率曲线的斜率没有影响，只为简单的频率减少。对比不同边界可以发现，简支边界相同条件下圆板的固有基频低于夹紧边界圆板。且简支圆板下频率随GPL重量分数增长的速率小于夹紧圆板。

(a) 夹紧边界 (b) 简支边界

Figure 3. Effect of GPL weight fraction on frequency at small amplitude ($V_1$ Distribution)

图3. GPL重量分数对频率的影响($V_1$ 分布)

(a) 夹紧边界 (b) 简支边界

Figure 4. Effect of GPL weight fraction on frequency at small amplitude ($V_2$ Distribution)

图4. GPL重量分数对频率的影响($V_2$ 分布)

图5表示了频率与GPL几何参数的关系，其中考虑孔隙为$V_1$ 分布，孔隙率$\theta =0.3$ 下夹紧和简支圆板。对比发现使用正方形GPL对比使用长方形有着更高的频率。随着$l_{GPL}/h_{GPL}$ 的增加，圆板的频率也随之增加，当$l_{GPL}/h_{GPL}<1000$ 时，圆板的频率增加速率较快，当$l_{GPL}/h_{GPL}>1000$ 时圆板频率增加缓慢。

(a) 夹紧边界 (b) 简支边界

Figure 5. The effect of GPL geometric parameters on frequency

图5. GPL几何参数对频率的影响

图6表示了频率与振幅的关系，其中考虑孔隙为$V_1$ 分布，孔隙率$\theta =0.3$ 下夹紧和简支圆板。对比不同的石墨烯重量分数可以发现随着石墨烯重量分数的增加，会放大振幅对频率降低的影响。对比三种不同的分布模式，可以发现三种不同的分布模式下频率都随着最大振幅的增加减少，但三种分布下频率降低的速率不一样。对比可以发现FG-O分布下频率减少速率大于均匀分布大于FG-X分布。

(a) 夹紧边界 (b) 简支边界

Figure 6. The effect of amplitude on frequency

图6. 振幅对频率的影响

图7到图8给出了孔隙率$\theta$ 对频率的影响曲线。其中考虑了不同的孔隙分布和不同的GPL分布，还考虑了两种不同的最大振幅，分别为$w_{\max }/h=0.01$ 、$w_{\max }/h=0.8$ 。通过对比可以发现，在$V_1$ 分布下固有频率随着孔隙率的增加而略微增加，这是因为孔隙率增加的时候GPL圆板的有效承载面积降低，导致弹性模量$E$ 降低，随而使材料弯曲刚度$D_{11}$ 降低。同时随着孔隙率的增加材料的总质量也随之降低，总质量的降低会导致惯性系数$I_0$ 降低。在线性问题下，当弯曲刚度$D_{11}$ 不变时，频率会随着惯性系数$I_0$ 的降低而增加；当惯性系数$I_0$ 不变时，频率会随着弯曲刚度$D_{11}$ 的降低而降低。相同的孔隙率下，$V_1$ 孔隙分布的弯曲刚度小于$V_2$ 孔隙分布的弯曲刚度，而两种孔隙分布模式下惯性系数$I_0$ 随孔隙率增加而降低的速率相当。因此在小振幅情况下，$V_1$ 孔隙分布下圆板的第一阶固有频率大于$V_2$ 孔隙分布下的圆板。但非线性问

(a) 夹紧边界 (b) 简支边界

Figure 7. The effect of $V_1$ distributed porosity on frequency

图7. $V_1$ 分布孔隙率对频率的影响

(a) 夹紧边界 (b) 简支边界

Figure 8. The effect of $V_2$ distributed porosity on frequency

图8. $V_2$ 分布孔隙率对频率的影响

题不能以相同的概念推导，通过对比非线性问题和线性问题求解出的曲线，可以发现在大振幅下圆板基频都低于小振幅下的圆板基频，但孔隙率的增加能略微降低大振幅频率和小振幅频率的差值。通过观察不同孔隙分布下的孔隙率–频率曲线可以看出，在$V_1$ 孔隙分布下频率随着孔隙率的增加而增加，其中UD分布增加速率最大，FG-O增加速率最小。在$V_2$ 孔隙分布下频率随着孔隙率的增加而降低。对比简支和夹紧圆板可以发现，孔隙率–频率变化曲线整体趋势相同，但简支无量纲频率低于夹紧，因为边界约束减弱时，结构更容易发生变形和振动，从而固有频率降低。

6. 结论

本文基于Halpin-Tsai模型和冯卡门大变形理论，对石墨烯增强梯度圆板的非线性自由振动问题进行分析，利用Hamilton原理导出了控制方程组并对控制方程无量纲化，然后通过微分求积法迭代求解方程组的非线性特征值问题，得到主要结论如下。

1) 随着石墨烯重量分数的增加，圆板的频率也随之增加。且FG-X型GPL分布模式的频率大于UD分布大于FG-O分布，且随着重量分数增长，FG-X型分布的增长速度大于均匀分布大于FG-O型分布。

2) 使用正方形GPL对比使用长方形有着更高的频率。随着$l_{GPL}/h_{GPL}$ 的增加，圆板的频率也随之增加，当$l_{GPL}/h_{GPL}<1000$ 时，圆板的频率增加速率较快，当$l_{GPL}/h_{GPL}>1000$ 时圆板频率增加缓慢。

3) 振幅与频率的关系为负相关，且随着振幅的增加，频率降低的速率也增加，且石墨烯重量分数的增加，会放大振幅对频率降低的影响。对比三种不同的分布模式，发现FG-O分布下频率减少速率大于均匀分布大于FG-X分布。

4) 在$V_1$ 孔隙分布下，孔隙率增加会提高频率；$V_2$ 孔隙分布下，孔隙率增加会降低频率，这与$V_1$ 孔隙分布下的弯曲刚度的大于$V_2$ 孔隙分布有关。

NOTES

^*通讯作者。