Vibration and Wave Analyses in the Functionally Graded Graphene-Reinforced Composite Plates Based on the First-Order Shear Deformation Plate Theory

Abstract

Graphene platelets (GPLs) can be used to enhance the mechanical and electrical properties of the matrix material, which efficiently determines and improves the dynamic behavior in composite structures. Based on the first-order shear deformation theory, this paper investigates the vibration and wave problems in a functionally graded graphene-reinforced composite plate. The composite plate is composed of the polymer matrix reinforced with GPLs that are dispersed along the thickness direction, following four kinds of functionally graded patterns. The governing equation of dynamic problems in the composite plate can be described in the state space formulation, and be solved using the method of reverberation-ray matrix (MRRM). Unlike the traditional state space method, this method is unconditionally stable due to introducing the dual coordinates, which can inherently avoid the numerical instability. After a validation study to verify the present analysis, a parametric study is conducted to analyze the effect of weight fraction, size and distribution patterns of the reinforments, as well as the boundary conditions and aspect ratios on the dynamic behaviors of the composite plate, hence providing a better way to achieve improved dynamic resistances of the GPLs composite plates.

1. Introduction

Since its discovery in 2004 [1], graphene has gained huge amounts of attention from both research and engineering communities due to its superior mechanical, electrical and thermal properties [2,3,4]. Compared to carbon nanotubes (CNTs), graphene has a much larger surface area and comparable properties, and soon becomes an emerging reinforcement in composites [5]. Rafiee et al. [6] found that the graphene/epoxy composite exhibits significantly higer tensile strength and Young’s modulus than the epoxy composites reinforced with the same amount of CNTs. Later, other studies also revealed that the buckling behavior of the graphenen/epoxy composite is vastly superior to the composites with CNTs reinforcements [7,8].

Recently, Luong et al. [9] developed a pioneering new method of producing graphene in bulk—flash Joule heating, which can zap any inexpensive source of solid carbon, from food scraps to car tires, and turn it into graphene in less than one second, and will open up various new applications in engineering. Its outstanding tensile strength, Young’s modulus, buckling behavior and a much larger surface area of graphene, together with the avent of bulk producing method [9], enable it to become a most promising reinforcement material. As a consequence, the polymer composites reinforced with graphene nanofillers have attracted extensive research interest [10].

Functionally graded materials (FGMs) consist of two or more materials, in addition to the continuous variation in both composition and properties along certain directions. The material properties of FGMs can be tailored to meet the requirements in various extreme working conditions. Recently, the functionally graded graphene platelets (FG-GPLs) reinforced composites have attracted academic interest. Wu et al. [11] analyzed the thermal buckling and post-buckling of FG-GPLs reinforced plates. Based on the Chebyshev–Ritz method, Yang et al. [12] studied the buckling and free vibration problem in FG-GPLs reinforced porous nanocomposite plates. Based on the high-order shear deformation theory, Shen et al. [13] performed analysis on the nonlinear vibration of FG-GPLs reinforced composite beams in thermal environments, Safarpour et al. [14] studied the bending and vibration response of FG-GPLs reinforced rectangular plates for different substrates and thermal conditions both theoretically and numerically, and Wang and Ma [15] focused on the effect of thickness stretching on static and dynamic responses of FG graphene-reinforced plates. However, the higher-order theory is more exact and complicated; there is no experimental evidence that it is necessary to adopt the higher-order theory in FG-GPLs composites [16]. That is, it is sufficient to use the classical plate theory (CPT) or first-order deformation theory (FSDT) to analyze the graphene-reinforced composites.

This paper analyzes the dynamic behaviors in the functionally graded plates reinforced with GPLs using the first-order shear deformation theory. The improved Halpin-Tsai micromechanical model and traditional rule of mixture are utilized to determine the effective material properties. The governing equations are constructed in the state space, and then solved using method of reverberation-ray matrix (MRRM) [17,18]. With special formulation, MRRM can avoid numerical instability inherently even for high frequencies or big wavenumbers [19,20]. Then a parametric study is conducted, the effects of the weight fraction, nanofiller size, distribution patterns of GPLs, the boundary conditions and aspect ratios on the dynamic behaviors of the FG-GPLs plates are studied in detail in the numerical examples.

2. FG Composite Plates Reinforced with GPLs

Consider a functionally graded graphene platelets (FG-GPLs)-reinforced composite plate in a coordinate system ($0\le x\le a$, $-\infty \le y\le \infty$, $-h/2\le z\le h/2$,), which is shown in Figure 1. For the wave propogation problem, the plate is finite in the x and z-axis, infinite in the y-axis, and made from a composite material comprising an isotropic matrix with GPLs reinforcement whose mixing ratio may vary along the z-axis (functionally graded materia, FGM).

According to the modified Halpin-Tsai micromechanics model [6], the effective Young’s modulus of FG-GPLs reinforced composite can be well approximated by

$$E=\frac{3}{8}\frac{1+{\xi }_l{\eta }_l{V}_{\mathrm{G}}}{1-{\eta }_l{V}_{\mathrm{G}}}{E}_{\mathrm{M}}+\frac{5}{8}\frac{1+{\xi }_w{\eta }_w{V}_{\mathrm{G}}}{1-{\eta }_w{V}_{\mathrm{G}}}{E}_{\mathrm{M}}$$

(1)

where ${\eta }_l=\frac{E_{\mathrm{G}}/E_{\mathrm{M}}-1}{E_{\mathrm{G}}/E_{\mathrm{M}}+{\xi }_l}$, ${\eta }_w=\frac{E_{\mathrm{G}}/E_{\mathrm{M}}-1}{E_{\mathrm{G}}/E_{\mathrm{M}}+{\xi }_w}$, $E_{\mathrm{G}}$ and $E_{\mathrm{M}}$ represent the Young’s moduli of GPLs and polymer matrix, respectively, and $V_{\mathrm{G}}$ represents the GPL volume fraction. The scaled parameters ${\xi }_l$ and ${\xi }_w$ denote the effects of both the geometry and size of GPL nanofillers, and can be expressed as

$${\xi }_l=2a_{\mathrm{G}}/t_{\mathrm{G}},\hspace{0.17em}\hspace{0.17em}{\xi }_w=2b_{\mathrm{G}}/t_{\mathrm{G}}$$

(2)

in which $a_{\mathrm{G}},b_{\mathrm{G}},t_{\mathrm{G}}$ are the length, width and thickness of GPL nanofillers, respectively.

However the modified Halpin–Tsai model cannot deterimine the mass density and Poisson’s ratio of composite materials. By using the rules of mixture, the Poisson’s ratio and density of the nanocomposites can be obtained as

$$\nu ={\nu }_{\mathrm{G}}{V}_{\mathrm{G}}+{\nu }_{\mathrm{M}}\left(1-V_{\mathrm{G}}\right)$$

(3)

$$\rho ={\rho }_{\mathrm{G}}{V}_{\mathrm{G}}+{\rho }_{\mathrm{M}}\left(1-V_{\mathrm{G}}\right)$$

(4)

where ${\nu }_{\mathrm{G}},{\rho }_{\mathrm{G}}$ and ${\nu }_{\mathrm{M}},{\rho }_{\mathrm{M}}$ represent the Poisson’s ratios, densities of GPLs and polymer matrix, respectively.

In this paper, four different GPLs distribution patterns are considered, named FG-V, FG-X, FG-O and UD [11,12], as schematically depicted in Figure 2, and can be expressed as

$$\begin{gathered} \mathrm{FG}-\mathrm{V}: V_{\mathrm{G}}\left(z\right)=V_{\mathrm{G}}^{\ast }\left(1+2z/h\right) \\ \mathrm{FG}-\mathrm{X}: V_{\mathrm{G}}\left(z\right)=V_{\mathrm{G}}^{\ast }4\left|z\right|/h \\ \mathrm{FG}-\mathrm{O}: V_{\mathrm{G}}\left(z\right)=V_{\mathrm{G}}^{\ast }\left(2-4\left|z\right|/h\right) \\ \mathrm{UD}: V_{\mathrm{G}}\left(z\right)=V_{\mathrm{G}}^{\ast } \end{gathered}$$

(5)

where $z\in [-0.5h,\hspace{0.17em}0.5h]$, and $V_{\mathrm{G}}^{\ast }$ is the total volume fraction of GPLs, which can be determined by the total weight fraction of GPLs $W_{\mathrm{G}}$,

$$V_{\mathrm{G}}^{\ast }=\frac{W_{\mathrm{G}}}{W_{\mathrm{G}}+\left({\rho }_{\mathrm{G}}/{\rho }_{\mathrm{M}}\right)\left(1-W_{\mathrm{G}}\right)}$$

(6)

3. Theoretical Formulations of Composite Plates

3.1. Governing Equations in the State Space

Based on the first-order shear deformation plate theory (FSDT), the displacements can be described as

$$u_1(x,y,z,t)=u(x,y,t)+z{\phi }_x(x,y,t)$$

(7a)

$$u_2(x,y,z,t)=v(x,y,t)+z{\phi }_y(x,y,t)$$

(7b)

$$u_3(x,y,z,t)=w(x,y,t)$$

(7c)

where t is the time, $u_1,\hspace{0.17em}{u}_2,\hspace{0.17em}{u}_3$ are the displacements at any point $(x,y,z)$, $u,\hspace{0.17em}v,\hspace{0.17em}w$ are the displacements of the point $(x,y,0)$ on the mid-plane, ${\phi }_x$ and ${\phi }_y$ are the rotations of the normal about the y- and x-axis, respectively.

The stress–strain relationship is given by the following expressions,

$$\left[\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\tau }_{xy}\end{array}\right]=\left[\begin{array}{ccc}{Q}_{11}& Q_{12}& 0\\ Q_{21}& Q_{22}& 0\\ 0& 0& Q_{66}\end{array}\right]\left[\begin{array}{c}{u}_{,x}+z{\phi }_{x,x}\\ v_{,y}+z{\phi }_{y,y}\\ u_{,y}+v_{,x}+z({\phi }_{x,y},+{\phi }_{y,x})\end{array}\right]$$

(8a)

$$\left[\begin{array}{c}{\tau }_{xz}\\ {\tau }_{yz}\end{array}\right]=\kappa \left[\begin{array}{cc}{Q}_{44}& 0\\ 0& Q_{55}\end{array}\right]\left[\begin{array}{c}{w}_{,x}+{\phi }_x\\ w_{,y}+{\phi }_y\end{array}\right]$$

(8b)

in which the shear correction coefficient $\kappa$ is assumed to be 5/6 for FSDT by Whitney and Pagano [21]. The position-dependent elastic constants $Q_{ij}$ are described as

$$Q_{11}=\frac{E_{11}(z)}{1-{\nu }_{12}{\nu }_{21}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{Q}_{22}=\frac{E_{22}(z)}{1-{\nu }_{12}{\nu }_{21}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{Q}_{12}=\frac{{\nu }_{12}{E}_{22}(z)}{1-{\nu }_{12}{\nu }_{21}}$$

(9a)

$$Q_{44}=G_{13}(z),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{Q}_{55}=G_{23}(z),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{Q}_{66}=G_{12}(z)$$

(9b)

The governing equations of motion for the FG-GPLs reinforced plate [21] can be achieved

$$N_{x,x}+N_{xy,y}=I_1\ddot{u}+I_2{\ddot{\phi }}_x$$

(10a)

$$N_{xy,x}+N_{y,y}=I_1\ddot{v}+I_2{\ddot{\phi }}_y$$

(10b)

$$Q_{x,x}+Q_{y,y}=I_1\ddot{w}$$

(10c)

$$M_{x,x}+M_{xy,y}-Q_x=I_2\ddot{u}+I_3{\ddot{\phi }}_x$$

(10d)

$$M_{xy,x}+M_{y,y}-Q_y=I_2\ddot{v}+I_3{\ddot{\phi }}_y$$

(10e)

where the over dot represents partial derivative of time t. The stress or moment resultants and the inertia-related terms are defined by

$$(N_x,N_y,N_{xy},Q_x,Q_y)={\displaystyle {\int }_{-h/2}^{h/2}({\sigma }_{xx},{\sigma }_{yy},{\tau }_{xy},{\tau }_{xz},{\tau }_{yz})}\mathrm{d}z$$

(11a)

$$(M_x,M_y,M_{xy})={\displaystyle {\int }_{-h/2}^{h/2}({\sigma }_{xx},{\sigma }_{yy},{\tau }_{xy})}z\mathrm{d}z$$

(11b)

$$(I_1,\hspace{0.17em}{I}_2,\hspace{0.17em}{I}_3)={\displaystyle {\int }_{-h/2}^{h/2}\rho (1,\hspace{0.17em}z,\hspace{0.17em}{z}^2)}\mathrm{d}z$$

(11c)

Using Equations (7)–(9) and (11), we can reconstruct Equation (10) in the state space form as

$$\frac{\partial }{\partial x}\mathsf{\Theta }=\mathsf{\Sigma }\mathsf{\Theta }=\left[\begin{array}{cc}{\mathsf{\Sigma }}_1& {\mathsf{\Sigma }}_2\\ {\mathsf{\Sigma }}_3& {\mathsf{\Sigma }}_4\end{array}\right]\mathsf{\Theta }$$

(12)

where $\mathsf{\Theta }={[u,v,w,{\phi }_x,{\phi }_y,M_x,M_{xy},N_x,N_{xy},Q_x]}^{\mathrm{T}}$ is the vector, and the coefficient matrices ${\mathsf{\Sigma }}_i$ (i = 1, …, 4) are given in Appendix A.

For the wave propagating along the y direction, the vector in the governing equations can be denoted in the following dimensionless form

$$\mathsf{\Theta }=\left[\begin{array}{c}u\\ v\\ w\\ {\phi }_x\\ {\phi }_y\\ M_x\\ M_{xy}\\ N_x\\ N_{xy}\\ Q_x\end{array}\right]=\left\{\begin{array}{c}h\overline{u}(\xi )\\ h\overline{v}(\xi )\\ h\overline{w}(\xi )\\ {\phi }_{\xi }(\xi )\\ {\phi }_{\eta }(\xi )\\ E_{\mathrm{M}}{h}^2{\overline{M}}_{\xi }(\xi )\\ E_{\mathrm{M}}{h}^2{\overline{M}}_{\xi \eta }(\xi )\\ E_{\mathrm{M}}h{\overline{N}}_{\xi }(\xi )\\ E_{\mathrm{M}}h{\overline{N}}_{\xi \eta }(\xi )\\ E_{\mathrm{M}}h{\overline{Q}}_{\xi }(\xi )\end{array}\right\}{\mathrm{e}}^{\mathrm{i}\hspace{0.17em}(\chi \eta -\omega t)}$$

(13)

in which $\xi =x/a$, $\eta =y/a$, $\chi$ the wavenumber of the $\eta$ direction and $\omega$ the circular frequency. Plugging Equation (13) into Equation (12), one obtains the dimensionless state space equation,

$$\frac{\mathrm{d}}{\mathrm{d}\xi }\overline{\mathsf{\Theta }}=\overline{\mathsf{\Sigma }}\overline{\mathsf{\Theta }}=\left[\begin{array}{cc}{\overline{\mathsf{\Sigma }}}_1& {\overline{\mathsf{\Sigma }}}_2\\ {\overline{\mathsf{\Sigma }}}_3& {\overline{\mathsf{\Sigma }}}_4\end{array}\right]\overline{\mathsf{\Theta }}$$

(14)

where $\overline{\mathsf{\Theta }}={[\overline{u},\overline{v},\overline{w},{\phi }_{\xi },{\phi }_{\eta },{\overline{M}}_{\xi },{\overline{M}}_{\xi \eta },{\overline{N}}_{\xi },{\overline{N}}_{\xi \eta },{\overline{Q}}_{\xi }]}^{\mathrm{T}},$ and the submatrices of the coefficient matrix $\overline{\mathsf{\Sigma }}$ are listed in Appendix B.

On the edges of $\xi =0,\hspace{0.17em}1$ ($x=0,\hspace{0.17em}a$), we take three cases for example: simply supported (S), clamped (C), or traction-free (F), which can be described as,

$$\overline{w}={\phi }_{\eta }={\overline{M}}_{\xi }={\overline{N}}_{\xi }={\overline{N}}_{\xi \eta }=0\hspace{0.17em}\left(\mathrm{S}\right)$$

(15a)

$$\overline{u}=\overline{v}=\overline{w}=\hspace{0.17em}{\phi }_{\xi }={\phi }_{\eta }=0\hspace{0.17em}\left(\mathrm{C}\right)$$

(15b)

$${\overline{M}}_{\xi }={\overline{M}}_{\xi \eta }={\overline{N}}_{\xi }={\overline{N}}_{\xi \eta }={\overline{Q}}_{\xi }=0\hspace{0.17em}\left(\mathrm{F}\right)$$

(15c)

The transfer matrix method, as one of the most efficient approaches dealing with the matrix equation, is usually adopted to seek the solution of state space Equation (14) [22,23,24]. However, when the product of the frequency and length is large enough [25,26,27], it will encounter the inherent numerical instability. So, in this work, the method of reverberation-ray matrix (MRRM) is utilized in what follows to overcome this problem.

3.2. MRRM Formulation

Owing to the constant state matrix $\overline{\mathsf{\Sigma }}$, we can describe the general solution of Equation (14) as

$$\overline{\mathsf{\Theta }}(\xi )=\mathsf{\Psi }\exp (\mathsf{\Lambda }\xi )\mathsf{\gamma }$$

(16)

AS stated in the eigenproblem in alignment of molecules [28], $\mathsf{\Lambda }=\mathrm{diag}({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_{10})$ is a diagonal matrix with ${\lambda }_i$ (i = 1,2,…,10) being the eigenvalues of the coefficient matrix $\overline{\mathsf{\Sigma }}$, $\mathsf{\Psi }=[{\mathsf{\Psi }}_1,{\mathsf{\Psi }}_2,{\mathsf{\Psi }}_3,\cdots ,{\mathsf{\Psi }}_{10}]$ are the eigenvectors, and $\mathsf{\gamma }$ are the wave amplitude vectors. Due to the physical meaning of wave propagation, the coefficient matrix $\overline{\mathsf{\Sigma }}$ is a Hamiltonian matrix, so the eigenvalue matrix can always be separated as follows,

$$\mathsf{\Lambda }=\left[\begin{array}{cc}{\mathsf{\Lambda }}_{+}& \mathbf{0}\\ \mathbf{0}& {\mathsf{\Lambda }}_{-}\end{array}\right]$$

(17)

where ${\mathsf{\Lambda }}_{+}$ comprises five eigenvalues with positive real part or positive imaginary part for pure imaginary, while ${\mathsf{\Lambda }}_{-}$ contains the remaining five eigenvalues, and ${\mathsf{\Lambda }}_{-}=-{\mathsf{\Lambda }}_{+}$. Similarly, we can divide the wave amplitudes $\mathsf{\gamma }$ to $\mathsf{\gamma }={[{({\mathsf{\gamma }}_{+})}^{\mathrm{T}},{({\mathsf{\gamma }}_{-})}^{\mathrm{T}}]}^{\mathrm{T}}$. So, the solution can be rewritten as

$$\overline{\mathsf{\Theta }}(\xi )=[{\mathsf{\Psi }}_{+}\hspace{0.17em}{\mathsf{\Psi }}_{-}]\left[\begin{array}{cc}\exp ({\mathsf{\Lambda }}_{+}\xi )& \mathbf{0}\\ \mathbf{0}& \exp ({\mathsf{\Lambda }}_{-}\xi )\end{array}\right]\left\{\begin{array}{c}{\mathsf{\gamma }}_{+}\\ {\mathsf{\gamma }}_{-}\end{array}\right\}$$

(18)

Before using MRRM, we should introduce the dual local coordinates first. For one element problem shown in Figure 3, the two local axes ${\xi }_1,\hspace{0.17em}{\xi }_2$ are in opposite directions.

In the two local coordinates ${\xi }_1$ and ${\xi }_2$, the solution to Equation (18) can be rewritten as

$${\overline{\mathsf{\Theta }}}^1({\xi }_1)=[{\mathsf{\Psi }}_{+}^1\hspace{0.17em}{\mathsf{\Psi }}_{-}^1]\left[\begin{array}{cc}\exp ({\mathsf{\Lambda }}_{+}^1{\xi }_1)& \mathbf{0}\\ \mathbf{0}& \exp ({\mathsf{\Lambda }}_{-}^1{\xi }_1)\end{array}\right]\left\{\begin{array}{c}{\mathsf{\gamma }}_{+}^1\\ {\mathsf{\gamma }}_{-}^1\end{array}\right\}$$

(19a)

$${\overline{\mathsf{\Theta }}}^2({\xi }_2)=[{\mathsf{\Psi }}_{+}^2\hspace{0.17em}{\mathsf{\Psi }}_{-}^2]\left[\begin{array}{cc}\exp ({\mathsf{\Lambda }}_{+}^2{\xi }_2)& \mathbf{0}\\ \mathbf{0}& \exp ({\mathsf{\Lambda }}_{-}^2{\xi }_2)\end{array}\right]\left\{\begin{array}{c}{\mathsf{\gamma }}_{+}^2\\ {\mathsf{\gamma }}_{-}^2\end{array}\right\}$$

(19b)

According to $1-{\xi }_1={\xi }_2$ in the dual local coordinate system in Figure 3, the relationship between ${\overline{\mathsf{\Theta }}}^1$ and ${\overline{\mathsf{\Theta }}}^2$ can be derived as

$${\overline{\mathsf{\Theta }}}^1(1-{\xi }_1)=T_f{\overline{\mathsf{\Theta }}}^2({\xi }_1)$$

(20)

where $T_f$ is the transformation matrix. Plugging Equation (20) into Equation (14) yields the following relationships [29]

$${\overline{\mathsf{\Sigma }}}^2=-T_f{\overline{\mathsf{\Sigma }}}^1{T}_f^{-1}, {\mathsf{\Lambda }}^2=-{\mathsf{\Lambda }}^1, {\mathsf{\Psi }}^2=T_f{\mathsf{\Psi }}^1$$

(21)

Combining Equations (19)–(21) yields

$$\left\{\begin{array}{c}{\mathsf{\gamma }}_{+}^1\\ {\mathsf{\gamma }}_{+}^2\end{array}\right\}=\left[\begin{array}{cc}\mathbf{0}& \exp (-{\mathsf{\Lambda }}_{+}^1)\\ \exp ({\mathsf{\Lambda }}_{-}^1)& \mathbf{0}\end{array}\right]\left\{\begin{array}{c}{\mathsf{\gamma }}_{-}^1\\ {\mathsf{\gamma }}_{-}^2\end{array}\right\}=P\left\{\begin{array}{c}{\mathsf{\gamma }}_{-}^1\\ {\mathsf{\gamma }}_{-}^2\end{array}\right\} \mathrm{or} a=Pd$$

(22)

where $P$ is the phase matrix. By dividing the unknown constants into two groups, exponential operation of large numbers is avoided in the phase matrix, which assures the unconditional numerical stability. Using Equation (19), the boundary conditions on the edges $\xi =0,\hspace{0.17em}1$ are

$$D^1{\mathsf{\gamma }}_{-}^1=S^1{\mathsf{\gamma }}_{+}^1\hspace{0.17em}(\mathrm{at}\hspace{0.17em}{\xi }_1=0), D^2{\mathsf{\gamma }}_{-}^2=S^2{\mathsf{\gamma }}_{+}^2\hspace{0.17em}(\mathrm{at}\hspace{0.17em}{\xi }_2=0\hspace{0.17em}/\hspace{0.17em}{\xi }_1=1)$$

(23)

or reformulated in a combined form

$$\left[\begin{array}{cc}{D}^1& \mathbf{0}\\ \mathbf{0}& D^2\end{array}\right]\left\{\begin{array}{c}{\mathsf{\gamma }}_{-}^1\\ {\mathsf{\gamma }}_{-}^2\end{array}\right\}=\left[\begin{array}{cc}{S}^1& \mathbf{0}\\ \mathbf{0}& S^2\end{array}\right]\left\{\begin{array}{c}{\mathsf{\gamma }}_{+}^1\\ {\mathsf{\gamma }}_{+}^2\end{array}\right\} \mathrm{or} Dd=Sa$$

(24)

where the scattering matrix $D,\hspace{0.17em}S$ are determined by the boundary conditions [30]. For simply supported, clamped or free edges, their expressions are summarized in Appendix C. Combining Equations (22) and (24), we can derive the characteristic equation of the wave propagating in the FG-GPLs reinforced composite plate,

$$\left|D-SP\right|=0$$

(25)

from which the dispersion relationship between the frequency and wavenumber can be derived. Unlike the solving process in Reference [31], we use $\left|D-SP\right|=0$ instead of $\left|I-D^{-1}SP\right|=0$ to avoid the inversion of matrix D, which may introduce numerical problem when $\left|D\right|=0$.

4. Results and Discussion

First, the free vibration of FGM square plates is analyzed to validate the present research. The square plates are made of aluminum oxide (ceramic) and Ti-6Al-4V (metal), with two opposite edges simply supported. The material properties, the ceramic variation and the plate size can be found in References [31,32,33,34]. The results of the classical plate theory [32], semi-analytical [33] and analytical results [34] based on the higher-order shear deformation plate theory are also included for comparison. Since the validation analysis has been performed in our former work [31], the numerical examples are omitted here for the sake of simplicity. As

Table 1. Material properties of the polymer [35] and GPLs [6].

Material	$\mathit{E}$	$\mathit{\nu }$	$\mathit{\rho }$
PMMA	2.5 GPa	0.34	1190 kg/m3
GPLs	1.01 TPa	0.186	1062.5 kg/m3

in reference [31] shows, the results are in good agreement with those of other theories, which validates our formulation and programing.

Then we study the free vibration problem in the FG-GPLs-reinforced plates. In what follows, the FG-GPLs reinforced plates are made from a mixture of the polymer PMMA and GPLs with a length of $a_{\mathrm{G}}=2.5 \mathsf{\mu }\mathrm{m}$, width of $b_{\mathrm{G}}=1.5 \mathsf{\mu }\mathrm{m}$, and thickness of $t_{\mathrm{G}}=1.5 \mathrm{nm}$. The material properties of PMMA and GPLs are shown in Table 1. The effective material properties of the composite plates are determined by Equations (1)–(4). Unless declared otherwise, $s=0.1$, $W_{\mathrm{G}}=1\%$, $a_{\mathrm{G}}=2.5 \mathsf{\mu }\mathrm{m}$, $b_{\mathrm{G}}=1.5 \mathsf{\mu }\mathrm{m}$, $t_{\mathrm{G}}=1.5 \mathsf{n}\mathrm{m}$, the plate is simply supported at two opposite ends with FGM pattern UD used in

Table 2. Natural frequencies of the SSSS GPLs reinforced PMMA square plates for different GPLs weight fractions.

	WG = 0
Mode	m = 1	m = 2	m = 3	m = 4	m = 5
1	0.5766	1.3758	2.5724	4.0427	5.6884
2	1.3757	2.1112	2.5726	4.6178	6.1920
3	1.3758	2.1113	2.5727	5.5164	6.1933
4	1.4760	3.2274	3.2270	5.5175	6.9880
5	2.5725	3.4071	3.2273	6.6755	6.9885
	WG = 0.1%
Mode	m = 1	m = 2	m = 3	m = 4	m = 5
1	0.5876	1.4019	2.6214	4.1198	5.7980
2	1.4020	2.1514	3.2892	4.7049	6.3113
3	1.5044	3.2890	4.3181	4.7054	7.1215
4	2.6216	3.2892	5.3838	5.6221	8.1849
5	2.8107	3.4728	5.5734	5.6229	9.4425
	WG = 0.5%
Mode	m = 1	m = 2	m = 3	m = 4	m = 5
1	0.6295	1.5022	2.8094	4.4151	6.2143
2	1.5021	2.3056	3.5249	4.4152	6.7636
3	1.6133	3.5247	3.5252	4.4153	7.6356
4	2.8093	3.7245	4.6279	5.0432	7.6358
5	3.0151	4.2794	4.6281	6.0257	8.7724
	WG = 1%
Mode	m = 1	m = 2	m = 3	m = 4	m = 5
1	0.6782	1.6186	3.0270	4.7585	6.6983
2	1.6185	2.4842	3.7984	5.4362	7.2900
3	1.7399	2.4843	3.7988	6.4957	7.2902
4	3.0271	2.4844	4.9881	7.8604	8.2288
5	3.2532	3.7983	6.2285	7.8609	8.2301

,

Table 3. Natural frequencies of the SSSS GPLs reinforced PMMA square plates for different GPLs aspect ratios (b_G = 1.5 μm).

	aG/bG = 1
Mode	m = 1	m = 2	m = 3	m = 4	m = 5
1	0.6782	1.6183	3.0267	4.7581	6.6961
2	1.6183	1.6184	3.0270	4.7584	6.6971
3	1.7398	2.4839	3.7981	5.4348	6.6973
4	3.0269	2.4840	4.9879	6.4941	7.2899
5	3.2529	2.4841	6.2280	6.4945	8.2277
	aG/bG = 4
Mode	m = 1	m = 2	m = 3	m = 4	m = 5
1	0.6783	1.6186	3.0272	4.7588	6.6978
2	1.6186	2.4845	3.7991	4.7591	6.6985
3	1.7400	3.7987	4.9885	5.4351	7.2911
4	3.0273	4.0178	6.2289	5.4353	8.2296
5	3.2534	4.6164	8.1301	6.4955	9.4578
	aG/bG = 10
Mode	m = 1	m = 2	m = 3	m = 4	m = 5
1	0.6783	1.6187	3.0273	4.7592	6.6987
2	1.6187	2.4846	3.0275	5.4354	7.2903
3	1.7401	3.7987	3.7990	6.4980	8.2312
4	3.0275	3.7988	4.9884	7.8609	8.2315
5	3.2535	4.0179	6.2291	7.8618	9.4575

and

Table 4. Natural frequencies of the SSSS GPLs reinforced PMMA square plates for different GPLs width-to-thickness ratios (t_G = 1.5 nm).

	bG/tG = 10
Mode	m = 1	m = 2	m = 3	m = 4	m = 5
1	0.6629	1.5822	2.9590	4.6511	6.5472
2	1.5820	2.4285	2.9591	4.6515	7.1267
3	1.7007	3.7128	3.7133	5.3141	8.0444
4	2.9590	3.9271	4.8759	6.3494	9.2438
5	3.1799	4.5122	4.8763	7.6830	9.2441
	bG/tG = 100
Mode	m = 1	m = 2	m = 3	m = 4	m = 5
1	0.6763	1.6140	3.0187	4.7449	6.6789
2	1.6139	2.4773	3.7873	5.4216	6.6806
3	1.7350	2.4774	3.7881	6.4760	7.2672
4	3.0186	3.7878	4.9743	6.4773	7.2707
5	3.2440	4.0061	6.2109	7.8382	7.2708
	bG/tG = 1000
Mode	m = 1	m = 2	m = 3	m = 4	m = 5
1	0.6782	1.6186	3.0270	4.7585	6.6983
2	1.6185	2.4842	3.7984	5.4362	7.2900
3	1.7399	2.4843	3.7988	6.4957	7.2902
4	3.0271	2.4844	4.9881	7.8604	8.2288
5	3.2532	3.7983	6.2285	7.8609	8.2301

. The dimensionless frequency is

$$\mathsf{\Omega }=\sqrt{{\rho }_M{\omega }^2{a}^2/E_M}$$

(26)

Table 2 compares the first five natural frequencies of SSSS square plates for m = 1~5 with different GPLs weight fractions (W_G = 0, 0.1, 0.5, 1%). For every m, the frequencies increase as the weight fraction rises. This can be expected as the increased GPLs content improves the stiffness of the plates substantially.

Table 3 and Table 4 are focused on the effect of GPLs shape on the natural frequencies of the SSSS square plates for m = 1~5 under free vibration. The influence of different shape ratios is investigated, including GPLs aspect ratio a_G/b_G and width-to-thickness ratio b_G/t_G. It can be observed from Table 3 and Table 4 that increasing a_G/b_G with a fixed b_G/t_G or raising b_G/t_G under a constant a_G/b_G, grows the natural frequencies of the square composite plates. A higher value of a_G/b_G represents a larger surface area, and a larger b_G/t_G indicates a wider graphene sheet. Hence, the bigger or wider GPLs can provide the better stiffening effect for the composite plates.

To better illustate the influence of different parameters on the wave spectra, Figure 4, Figure 5, Figure 6 and Figure 7 are presented. Unless stated otherwise, $s=0.1$, $W_{\mathrm{G}}=1\%$, the plate is simply supported at $\xi =0.1$ with FGM pattern UD used in all examples. The dimensionless frequency and wavenumber are $\mathsf{\Omega }=\sqrt{{\rho }_M{\omega }^2{a}^2/E_M}$ and $\chi =ka$, respectively.

Figure 4 investigates the effect of weight fractions on the dispersion spectra of GPLs reinforced plates. The frequencies always increase with the increasing weight fraction, so the dispersion curves will move upward as shown in Figure 4. For every mode, the blue solid branch for W_G = 0 is the lowest, the red dot dash line for W_G = 0.5% is middle, while the black dash line for W_G = 1% is the highest.

Figure 5 shows the influences of different GPLs distribution patterns on the dispersion spectra of the FG plates. As demonstrated by the figure, the GPLs distribution patterns have little influence on the in-plane modes, while they have a considerable impact on the flexural modes [31]. For the same branch of the dispersion spectra, the red dot dash is the lowest (FG-O), the green dash (FG-X) is the highest, yet the blue solid (UD) and the black dash (FG-V) are quite close to each other. Compared with the other patterns, the FG-X plate has the highest frequencies. This is due to the symmetric configuration in the thickness direction. Meanwhile, the top and bottom surfaces with the maximum normal stress are GPLs rich and the neutral layer with the minimum normal stress is pure PMMA. This distribution pattern is the most effective since it makes better use of GPLs reinforcements and determines the higher bending stiffness. On the contrary, FG-O patterns performs the lowest frequency due to the weakest bending stiffness.

Figure 6 analyzes the effect of the boundary conditions on the dispersion spectra of the composite plates with different boundary conditions at $\xi =0,1$ in which SS, CC and FF represent both edges, simply supported, clamped and free. As we can see, the boundary condition has little effect on the in-plane modes of the dispersion spectra, and the flexural modes are quite sensitive to the boundary conditons. Especially the first two black dash branches are unique for the FF plate, the corresponding branches for the SS and CC conditions are absent.

Figure 7 depicts the first four dispersion spectra of the composite plates with different aspect ratios (s = 0.1, 0.15, 0.2), which include two in-plane and two flexural modes. Similar to the phenomenon shown in Figure 6, the dispersion relations are insensitive to the aspect ratio for the two in-plane modes, while the aspect ratio noticeably affects the flexural modes. The frequency grows with the apsect ratio increases for the flexural modes. The in-plane mode branches go steeper, while the flexural mode branches become milder. Concluded from Figure 4, Figure 5, Figure 6 and Figure 7, the GPLs distribution patterns, boundary conditions and aspect ratios influence the flexural modes more obvious than the in-plane modes, while the GPLs weight fraction affects both the flexural and in-plane modes remarkably.

5. Conclusions

In this paper, the dynamic behaviors of FG-GPLs-reinforced plates are analyzed based on the first-order shear deformation theory within the framework of the state space. The method of reverberation-ray matrix (MRRM) is utilized to obtain the natural frequencies and dispersion spectra. By using the dual coordinates to avoid big number computation and eliminating the inversion of the zero matrix, MRRM can determine a solution which is inherently numerically stable. The free vibration and dispersion relation are revealed with a detailed examination of the influence of the GPLs weight fractions, GPLs size, distribution patterns, boundary conditions and aspect ratios. Numerical results show that the weight fraction and size of GPLs can improve the natural frequencies, and all the dispersion curves can be considerably affected by changing the graphene weight fraction, while only the flexural modes of the dispersion curves are sensitive to the graphene distribution, the boundary conditions and aspect ratios. It is envisaged that the present analysis can provide guidelines to design novel, better-performing composite structures.