Nonlinear Vibration of Functionally Graded Graphene Nanoplatelets Polymer Nanocomposite Sandwich Beams

Abstract

We provide an analytical investigation of the nonlinear vibration behavior of thick sandwich nanocomposite beams reinforced by functionally graded (FG) graphene nanoplatelet (GPL) sheets, with a power-law-based distribution throughout the thickness. We assume the total amount of the reinforcement phase to remain constant in the beam, while defining a relationship between the GPL maximum weight fraction, the power-law parameter, and the thickness of the face sheets. The shear and rotation effects are here considered using a higher-order laminated beam model. The nonlinear partial differential equations (PDEs) of motion are derived from the Von Kármán strain-displacement relationships, here solved by applying an expansion of free vibration modes. The numerical results demonstrate the key role of the amplitudes on the vibration response of GPL-reinforced sandwich beams, whose nonlinear oscillation behavior is very important in the physical science, mechanical structures and other mathematical analyses. The sensitivity of the response to the total amount of GPLs is explored herein, along with the possible effects related to the power-law parameter, the structural geometry, and the environmental conditions. The results indicate that changing the nanofiller distribution patterns with the proposed model can remarkably increase or decrease the effective stiffness of laminated composite beams.

1. Introduction

Sandwich structures, generally made of a soft core and two hard face sheets, are largely used in the aerospace, oil, gas, and petrochemical industries, due to their enhanced mechanical properties, namely, a high strength-to-weight ratio and a high resistance to heat, humidity, and noise [1,2,3,4,5]. Hence, in recent decades, much attention has been paid to the mechanical behavior of these structures [6,7,8,9,10,11,12]. Based on the available literature, it seems that the geometry of the layers, the mechanical properties of the constituents, and the geometrical properties of the whole structure can have a meaningful effect on the static and dynamic behavior of sandwich structures [13,14,15,16,17,18,19,20]. The presence of some reinforcing layers in sandwich structures represents one important issue to consider for a general improvement of their mechanical properties [21,22,23,24]. Nowadays, with the advancement of nanotechnology, carbon nanotubes (CNTs) and graphene sheets (GSs) are two alternative options for the reinforcement of structures, due to their extraordinary properties. This has led to an extensive research on the behavior of sandwich structures reinforced with nanocomposites [25,26,27,28]. Among different reinforcement possibilities, graphene nanoplatelets (GPLs) provide a uniform reinforced assembly, as well as the easiest manufacturing process, as discussed in [29,30,31,32]. Graphene is a monolayer structure of carbon atoms with extraordinary electrical, mechanical, thermal, and optical properties [33,34,35], which make it very attractive for high-tech device applications, such as micro/nano-electromechanical systems [36]. Graphene and its derivatives—namely, GPLs—are increasingly applied as a reinforcement material in many nanocomposite structures [37,38,39,40,41,42,43]. This justifies the large attention paid in the literature to the mechanical behavior of structures reinforced with graphenic materials [44,45,46,47,48,49,50,51,52,53].

Despite the extensive literature available on the behavior of composite structures reinforced with nanostructures, there is a general lack of works focusing on the nonlinear dynamic and vibration behavior of sandwich beams reinforced by GPLs. This is here investigated for thick polymer sandwich beams with face sheets reinforced by GPLs, in a context where the reduced weight of polymers and the high strength of GPLs can provide remarkable properties in the equivalent composite structure. A novel reinforcement model is proposed herein, which considers the functionality of the GPLs distribution throughout the thickness of the face sheets, and a constant total amount of the reinforced material. A higher-order laminated beam theory is applied to include the shear and rotation effects on the thick GPL-reinforced sandwich beam, where the nonlinear governing equations of the problem are solved in a straightforward manner by means of the multiple timescales method. The main advantage of the present method is that it can cover weak or strong nonlinearities with possible damping effects. The method is demonstrated to be very simple and accurate with respect to other existing predictions and theories from the literature.

The reinforcement phase varies along the thickness according to a power-law distribution, whereby the effective material properties of the nanocomposite beam are determined by means of the Halpin–Tsai micromechanics model and the rule of mixtures. The nonlinear partial equations of motion are derived by the Hamilton’s principle, in accordance with the third-order shear deformation theory and the Von Kármán strain-displacement relationships. We then apply Galerkin’s approach to discretize the nonlinear differential equations of motion, while determining the frequency equations by means of the multiple timescales method. Various numerical examples indicate the accuracy of the proposed model and check for the sensitivity of the vibration response of GPL-reinforced sandwich beams, of great interest for design and practical purposes.

The paper is organized as follows. In Section 2 the mechanical and geometrical properties of materials and their structure are briefly described. Section 3 presents the theoretical formulation of the problem, along with the numerical procedure. A number of illustrative applications and comparative evaluations with the available literature are proposed in Section 4. Finally, in Section 5 some concluding remarks are reported.

2. Material Properties and Geometry

A nanocomposite sandwich beam with length $L$, thickness $h_t$ and width $b$ is considered, as shown in Figure 1. The Cartesian coordinate system $\left(x,z\right)$ is here used to derive the equations of motion, where the structural mid-plane is parallel to the $x$-axis. The beam is made of a homogeneous core and two face sheets with a symmetric GPL-based reinforcement, whose weight fraction satisfies the following power-law:

$$\Gamma \left(z\right)={\Gamma }_{\max }{\left(\frac{2\left|z\right|-h_c}{h_t-h_c}\right)}^{\kappa },$$

(1)

where ${\Gamma }_{\max }$ is the maximum value of the distribution function and $\kappa$ is a power-law parameter, which defines the GPLs dispersion throughout the thickness of the face sheets.

For a proper analysis, the total amount of GPLs in the beam remains constant independently of the distribution pattern (see Figure 2). This means that, if we keep constant the total amount of GPLs in the beam, ${\Gamma }_b$, the maximum value of the GPL weight fraction, ${\Gamma }_{\max }$, increases by increasing the power-law parameter. Note that the total amount of reinforced GPLs decreases by increasing the power-law parameter if the maximum value of the GPL weight fraction is kept constant.

Due to possible difficulties during the manufacturing process of a functional reinforced lamina, each face sheet is considered to be made of N layers with equal thickness and the GPL reinforcement is assumed to be uniform within each layer (see Figure 1).

Therefore, the GPL weight fraction in the $k$th lamina can be defined as

$${\Gamma }^{\left(k\right)}={\Gamma }_{\max }{\left(\frac{2\left|h_k\right|-h_c}{h_t-h_c}\right)}^{\kappa },$$

(2)

where $h_k$ is the distance between the mid-plane of the beam and the mid-plane of the $k$th layer.

The volume fraction of the reinforced GPLs can be related to their weight fraction as

$$V_{\mathrm{GPL}}^{\left(k\right)}=\frac{{\Gamma }^{\left(k\right)}}{{\Gamma }^{\left(k\right)}+\left({\rho }_{\mathrm{GPL}}/{\rho }_M\right)\left(1-{\Gamma }^{\left(k\right)}\right)},$$

(3)

${\rho }_M$ and ${\rho }_{\mathrm{GPL}}$ being the mass density of the matrix and GPLs, respectively. Here, the Halpin–Tsai micromechanical model is adopted to define the effective elastic modulus of the reinforced face sheets [54]. Moreover, the GPL reinforcement is assumed to be randomly oriented in each lamina [27].

Therefore, the elastic modulus for the $k$th layer can be expressed as

$$E_C^{\left(k\right)}=\frac{3}{8}\left(\frac{1+{\xi }_L{\eta }_L{V}_{\mathrm{GPL}}^{\left(k\right)}}{1-{\eta }_L{V}_{\mathrm{GPL}}^{\left(k\right)}}\right)E_M+\frac{5}{8}\left(\frac{1+{\xi }_W{\eta }_W{V}_{\mathrm{GPL}}^{\left(k\right)}}{1-{\eta }_W{V}_{\mathrm{GPL}}^{\left(k\right)}}\right)E_M,$$

(4)

where

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

(5)

$${\xi }_W=\frac{2w_{\mathrm{GPL}}}{h_{\mathrm{GPL}}},{\xi }_L=\frac{2l_{\mathrm{GPL}}}{h_{\mathrm{GPL}}},$$

(6)

and $h_{\mathrm{GPL}}$, $l_{\mathrm{GPL}}$, $w_{\mathrm{GPL}}$ stand for the average thickness, length, and width of GPLs, respectively; $E_M$ and $E_{\mathrm{GPL}}$ denote the Young modulus of the matrix and GPLs, respectively.

Using the rule of mixtures, the effective mass density and Poisson’s ratio for the $k$th layer can be defined as

$${\rho }_C^{\left(k\right)}={\rho }_{\mathrm{GPL}}{V}_{\mathrm{GPL}}^{\left(k\right)}+{\rho }_M\left(1-V_{\mathrm{GPL}}^{\left(k\right)}\right),$$

(7)

$${\nu }_C^{\left(k\right)}={\nu }_{\mathrm{GPL}}{V}_{\mathrm{GPL}}^{\left(k\right)}+{\nu }_M\left(1-V_{\mathrm{GPL}}^{\left(k\right)}\right),$$

(8)

${\nu }_M$ and ${\nu }_{\mathrm{GPL}}$ being the Poisson’s ratio of the matrix and GPLs, respectively.

3. Theoretical Formulations

In this section, the nonlinear governing equations of the problem for functionally graded (FG) GPL-reinforced sandwich beams are derived by Hamilton’s principle, while using a higher-order shear deformation approach.

3.1. Displacement Field and Strains

In agreement with the third-order shear deformation theory [55,56], the displacement components $u_1\left(x,t\right)$ and $u_3\left(x,t\right)$ of an arbitrary point in the $x$ and $z$ directions for shear deformable sandwich beams can be expressed as

$$u_1\left(x,z,t\right)=u\left(x,t\right)+z\varphi \left(x,t\right)-\frac{4z^3}{3h_t^2}\left(\varphi +\frac{\partial w}{\partial x}\right),$$

(9)

$$u_3\left(x,z,t\right)=w\left(x,t\right),$$

(10)

where $u\left(x,t\right)$ and $w\left(x,t\right)$ are the displacement components of a point at the mid-plane of the beam in the $x$ and $z$ directions, respectively. Moreover, $\varphi \left(x,t\right)$ denotes the slope of a transverse normal at $z=0$. Based on the Von Kármán strain-displacement relationships, the nonlinear strain components associated with the displacement field (9)–(10) can be written as

$${\epsilon }_{xx}={\epsilon }_{xx}^{\left(0\right)}+z{\epsilon }_{xx}^{\left(1\right)}+z^3{\epsilon }_{xx}^{\left(3\right)},$$

(11)

$${\gamma }_{xz}={\gamma }_{xz}^{\left(0\right)}+z^2{\gamma }_{xz}^{\left(2\right)},$$

(12)

where

$${\epsilon }_{xx}^{\left(0\right)}=\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2,{\epsilon }_{xx}^{\left(1\right)}=\frac{\partial \varphi }{\partial x},{\epsilon }_{xx}^{\left(3\right)}=-c_1\left(\frac{\partial \varphi }{\partial x}+\frac{{\partial }^{2}w}{\partial x^2}\right),$$

(13)

$${\gamma }_{xz}^{\left(0\right)}=\varphi +\frac{\partial w}{\partial x},{\gamma }_{xz}^{\left(2\right)}=-c_2\left(\varphi +\frac{\partial w}{\partial x}\right),$$

(14)

and

$$c_1=\frac{4}{3h_t^2},c_2=3c_1=\frac{4}{h_t^2}.$$

(15)

3.2. Equations of Motion

The equations of motion of FG-GPL reinforced sandwich beams are derived from Hamilton’s principle. Accordingly, we have

$${{\displaystyle \int }}_{t_1}^{t_2}\left(\delta U-\delta W-\delta K\right)dt=0,$$

(16)

where $U$ is the strain energy, $W$ is the work done by external forces, and $K$ is the kinetic energy. The virtual strain energy $\delta U$ for the third-order shear deformable sandwich beams reads as follows

$$\begin{array}{lll}\delta U\hfill & =\hfill & {{\displaystyle \int }}_A{{\displaystyle \int }}_0^L\left({\sigma }_{xx}\delta {\epsilon }_{xx}+{\sigma }_{xz}\delta {\gamma }_{xz}\right)dxdA\\ & =& {{\displaystyle \int }}_A{{\displaystyle \int }}_0^L\left[{\sigma }_{xx}\left(\delta {\epsilon }_{xx}^{\left(0\right)}+z\delta {\epsilon }_{xx}^{\left(1\right)}+z^3\delta {\epsilon }_{xx}^{\left(3\right)}\right)+{\sigma }_{xz}\left(\delta {\gamma }_{xz}^{\left(0\right)}+z^2\delta {\gamma }_{xz}^{\left(2\right)}\right)\right]dxdA\\ & =& {{\displaystyle \int }}_0^L[-\frac{\partial N_{xx}}{\partial x}\delta u-\frac{\partial }{\partial x}\left(N_{xx}\frac{\partial w}{\partial x}\right)\delta w-\frac{\partial M_{xx}}{\partial x}\delta \varphi +c_1\frac{\partial P_{xx}}{\partial x}\delta \varphi -c_1\frac{{\partial }^2{P}_{xx}}{\partial x^2}\delta w\\ & +& Q_x\delta \varphi -\frac{\partial Q_x}{\partial x}\delta w-c_2{R}_x\delta \varphi +c_2\frac{\partial R_x}{\partial x}\delta w]dx\\ & +& {\left[N_{xx}\delta u+N_{xx}\frac{\partial w}{\partial x}\delta w+M_{xx}\delta \varphi -c_1{P}_{xx}\delta \varphi -c_1{P}_{xx}\frac{\partial }{\partial x}\delta w+c_1\frac{\partial P_{xx}}{\partial x}\delta w+Q_x\delta w-c_2{R}_x\delta w\right]}_0^L,\end{array}$$

(17)

where

$$\left[\begin{array}{c}{N}_{xx}\\ M_{xx}\\ P_{xx}\end{array}\right]={{\displaystyle \int }}_A{\sigma }_{xx}\left[\begin{array}{c}1\\ z\\ z^3\end{array}\right]dA,\left[\begin{array}{c}{Q}_x\\ R_x\end{array}\right]={{\displaystyle \int }}_A{\sigma }_{xz}\left[\begin{array}{c}1\\ z^2\end{array}\right]dA.$$

(18)

The virtual kinetic energy $\delta K$ is defined as

$$\begin{array}{ll}\delta K& ={{\displaystyle \int }}_A{{\displaystyle \int }}_0^L\rho \left(z\right)\left({\dot{u}}_1\delta {\dot{u}}_1+{\dot{u}}_3\delta {\dot{u}}_3\right)dAdx\\ & ={{\displaystyle \int }}_0^L[\left(m_0\stackrel{..}{u}+m_1\stackrel{..}{\varphi }-c_1{m}_3\left(\stackrel{..}{\varphi }+\frac{\partial \stackrel{..}{w}}{\partial x}\right)\right)\delta u\\ & +\left(m_1\stackrel{..}{u}+m_2\stackrel{..}{\varphi }-c_1{m}_3\stackrel{..}{u}-c_1{m}_4\left(2\stackrel{..}{\varphi }+\frac{\partial \stackrel{..}{w}}{\partial x}\right)+c_1^2{m}_6\left(\stackrel{..}{\varphi }+\frac{\partial \stackrel{..}{w}}{\partial x}\right)\right)\delta \varphi \\ & +\left(m_0\stackrel{..}{w}+c_1{m}_3\frac{\partial \stackrel{..}{u}}{\partial x}+c_1{m}_4\frac{\partial \stackrel{..}{\varphi }}{\partial x}-c_1^2{m}_6\left(\frac{\partial \stackrel{..}{\varphi }}{\partial x}+\frac{{\partial }^2\stackrel{..}{w}}{\partial x^2}\right)\right)\delta w]dx,\end{array}$$

(19)

with

$$m_i={{\displaystyle \int }}_A\rho \left(z\right)z^{i}dA={{\displaystyle \sum }}_{k=1}^{2N+1}b{{\displaystyle \int }}_{z_k}^{z_{k+1}}{\rho }_c^{\left(k\right)}{z}^{i}dz,i=0,1,2,3,4,6.$$

(20)

In the total absence of external forces on the structure, it is $\delta W=0$. Therefore, by substitution of Equations (17) and (19) into Equation (16), by integrating the result by parts, and equating the coefficients of $\delta u$, $\delta w$, and $\delta \varphi$ to zero separately, we get the following nonlinear equations of motion:

$$\frac{\partial N_{xx}}{\partial x}=m_0\stackrel{..}{u}+m_1\stackrel{..}{\varphi }-c_1{m}_3\left(\stackrel{..}{\varphi }+\frac{\partial \stackrel{..}{w}}{\partial x}\right),$$

(21)

$$\frac{\partial M_{xx}}{\partial x}-c_1\frac{\partial P_{xx}}{\partial x}-Q_x+c_2{R}_x=m_1\stackrel{..}{u}+m_2\stackrel{..}{\varphi }-c_1{m}_3\stackrel{..}{u}-c_1{m}_4\left(2\stackrel{..}{\varphi }+\frac{\partial \stackrel{..}{w}}{\partial x}\right)+c_1^2{m}_6\left(\stackrel{..}{\varphi }+\frac{\partial \stackrel{..}{w}}{\partial x}\right),$$

(22)

$$\frac{\partial }{\partial x}\left(N_{xx}\frac{\partial w}{\partial x}\right)+c_1\frac{{\partial }^2{P}_{xx}}{\partial x^2}+\frac{\partial Q_x}{\partial x}-c_2\frac{\partial R_x}{\partial x}=m_0\stackrel{..}{w}+c_1{m}_3\frac{\partial \stackrel{..}{u}}{\partial x}+c_1{m}_4\frac{\partial \stackrel{..}{\varphi }}{\partial x}-c_1^2{m}_6\left(\frac{\partial \stackrel{..}{\varphi }}{\partial x}+\frac{{\partial }^2\stackrel{..}{w}}{\partial x^2}\right).$$

(23)

Thus, we define the stress resultants in terms of the displacement and rotation components of the sandwich beam as

$$N_{xx}=A_{11}\left(\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2\right)+B_{11}\left(\frac{\partial \varphi }{\partial x}\right)-c_1{E}_{11}\left(\frac{\partial \varphi }{\partial x}+\frac{{\partial }^{2}w}{\partial x^2}\right),$$

(24)

$$M_{xx}=B_{11}\left(\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2\right)+D_{11}\left(\frac{\partial \varphi }{\partial x}\right)-c_1{F}_{11}\left(\frac{\partial \varphi }{\partial x}+\frac{{\partial }^{2}w}{\partial x^2}\right),$$

(25)

$$P_{xx}=E_{11}\left(\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2\right)+F_{11}\left(\frac{\partial \varphi }{\partial x}\right)-c_1{H}_{11}\left(\frac{\partial \varphi }{\partial x}+\frac{{\partial }^{2}w}{\partial x^2}\right),$$

(26)

$$Q_x=\left(A_{55}-c_2{D}_{55}\right)\left(\varphi +\frac{\partial w}{\partial x}\right),$$

(27)

$$R_x=\left(D_{55}-c_2{F}_{55}\right)\left(\varphi +\frac{\partial w}{\partial x}\right),$$

(28)

where

$$\left(A_{11},B_{11},D_{11},E_{11},F_{11},H_{11}\right)={{\displaystyle \sum }}_{k=1}^{N}b{{\displaystyle \int }}_{h_k}^{h_{k+1}}{Q}_{11}^{\left(k\right)}\left(1,z,z^2,z^3,z^4,z^6\right)dz,$$

(29)

$$\left(A_{55},D_{55},F_{55}\right)={{\displaystyle \sum }}_{k=1}^{N}b{{\displaystyle \int }}_{h_k}^{h_{k+1}}{Q}_{55}^{\left(k\right)}\left(1,z^2,z^4\right)dz,$$

(30)

and

$$Q_{11}^{\left(k\right)}=\frac{E_C^{\left(k\right)}}{1-{\left({\nu }_C^{\left(k\right)}\right)}^2},Q_{55}^{\left(k\right)}=G_C^{\left(k\right)}=\frac{E_C^{\left(k\right)}}{2\left(1+{\nu }_C^{\left(k\right)}\right)},$$

(31)

In view of Equations (21)–(28), the nonlinear partial differential equations of motion of FG GPL-reinforced sandwich beams can be written as

$$\begin{array}{ll}{A}_{11}\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}w}{\partial x^2}\frac{\partial w}{\partial x}\right)& +\left(B_{11}-c_1{E}_{11}\right)\frac{{\partial }^2\varphi }{\partial x^2}-c_1{E}_{11}\frac{{\partial }^{3}w}{\partial x^3}\\ & =m_0\frac{{\partial }^{2}u}{\partial t^2}+\left(m_1-c_1{m}_3\right)\frac{{\partial }^2\varphi }{\partial t^2}-c_1{m}_3\frac{{\partial }^{3}w}{\partial x\partial t^2},\end{array}$$

(32)

$$\begin{array}{ll}\left(B_{11}-c_1{E}_{11}\right)\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}w}{\partial x^2}\frac{\partial w}{\partial x}\right)& +\left(D_{11}-2c_1{F}_{11}+c_1^2{H}_{11}\right)\frac{{\partial }^2\varphi }{\partial x^2}+\left(-c_1{F}_{11}+c_1^2{H}_{11}\right)\frac{{\partial }^{3}w}{\partial x^3}\\ & +\left(-A_{55}+2c_2{D}_{55}-c_2^2{F}_{55}\right)\left(\varphi +\frac{\partial w}{\partial x}\right)\\ & =\left(m_1-c_1{m}_3\right)\frac{{\partial }^{2}u}{\partial t^2}+\left(m_2-2c_1{m}_4+c_1^2{m}_6\right)\frac{{\partial }^2\varphi }{\partial t^2}+\left(-c_1{m}_4+c_1^2{m}_6\right)\frac{{\partial }^{3}w}{\partial x\partial t^2},\end{array}$$

(33)

$$\begin{array}{ll}{c}_1{E}_{11}\left(\frac{{\partial }^{3}u}{\partial x^3}-{\left(\frac{{\partial }^{2}w}{\partial x^2}\right)}^2\right)\hfill & +\left(c_1{F}_{11}-c_1^2{H}_{11}\right)\frac{{\partial }^3\varphi }{\partial x^3}-c_1^2{H}_{11}\frac{{\partial }^{4}w}{\partial x^4}\hfill \\ \hfill & +\left(A_{55}-2c_2{D}_{55}+c_2^2{F}_{55}\right)\left(\frac{\partial \varphi }{\partial x}+\frac{{\partial }^{2}w}{\partial x^2}\right)\hfill \\ \hfill & +A_{11}\left(\frac{{\partial }^{2}u}{\partial x^2}\frac{\partial w}{\partial x}+\frac{\partial u}{\partial x}\frac{{\partial }^{2}w}{\partial x^2}+\frac{3}{2}\left(\frac{{\partial }^{2}w}{\partial x^2}{\left(\frac{\partial w}{\partial x}\right)}^2\right)\right)\hfill \\ \hfill & +\left(B_{11}-c_1{E}_{11}\right)\left(\frac{{\partial }^2\varphi }{\partial x^2}\frac{\partial w}{\partial x}+\frac{\partial \varphi }{\partial x}\frac{{\partial }^{2}w}{\partial x^2}\right)\hfill \\ \hfill & =m_0\frac{{\partial }^{2}w}{\partial t^2}+c_1{m}_3\frac{{\partial }^{3}u}{\partial x\partial t^2}+\left(c_1{m}_4-c_1^2{m}_6\right)\frac{{\partial }^3\varphi }{\partial x\partial t^2}-c_1^2{m}_6\frac{{\partial }^{4}w}{\partial x^2\partial t^2}.\hfill \end{array}$$

(34)

3.3. Solution Procedure

In this section, the nonlinear equations of motion are solved numerically, in order to obtain the linear and nonlinear frequency equations. In this regard, the nonlinear partial differential equations (PDEs) of motion (32)–(34) are discretized as ordinary differential equations by employing the Galerkin method. Afterwards, the multiple timescales approach is used to obtain the nonlinear frequency equation. Here, we assume that the GPL-reinforced sandwich beam is simply supported at both ends with movable supports. Based on these assumptions, the displacement and rotation field of the beam can be defined as expansions of the free vibration mode shapes, namely,

$$u\left(x,t\right)={{\displaystyle \sum }}_{n=1}^{\infty }{U}_n\left(t\right){\chi }_n\left(x\right),{\chi }_n\left(x\right)=\mathrm{Cos}\left(\frac{n\pi x}{L}\right),$$

(35)

$$\varphi \left(x,t\right)={{\displaystyle \sum }}_{n=1}^{\infty }{\phi }_n\left(t\right){\psi }_n\left(x\right),{\psi }_n\left(x\right)=\mathrm{Cos}\left(\frac{n\pi x}{L}\right),$$

(36)

$$w\left(x,t\right)={{\displaystyle \sum }}_{n=1}^{\infty }{W}_n\left(t\right){\lambda }_n\left(x\right),{\lambda }_n\left(x\right)=\mathrm{Sin}\left(\frac{n\pi x}{L}\right),$$

(37)

$U_n$, ${\phi }_n$ and $W_n$ being the unknown generalized coordinates which stand for the amplitude of the vibration. Moreover, the following functions ${\chi }_n\left(x\right)$, ${\psi }_n\left(x\right)$, and ${\lambda }_n\left(x\right)$ are introduced to satisfy all boundary conditions of the system. By considering a single mode approximate solution and by substitution of Equations (35)–(37) into Equations (32)–(34), after multiplying the results by $\chi$, $\psi$, and $\lambda$ and after their integration over the domain of the system, the following nonlinear differential equations of motion are obtained

$$g_{11}{U}_n+g_{12}{W}_n+g_{13}{W}_n{}^2+g_{14}{\phi }_n=0,$$

(38)

$$g_{21}{U}_n+g_{22}{W}_n+g_{23}{W}_n{}^2+g_{24}{\phi }_n=0,$$

(39)

$$g_{31}{U}_n+g_{32}{W}_n+g_{33}{U}_n{W}_n+g_{34}{W}_n{}^2+g_{35}{W}_n{}^3+g_{36}{\phi }_n+g_{37}{W}_n{\phi }_n+g_{38}{U}_n{}^{″}+g_{39}{W}_n{}^{\prime \prime }+g_{310}{\phi }_n{}^{″}=0,$$

(40)

where coefficients $g_{11}$, $g_{12}$, …, $g_{310}$ are detailed in Appendix A. The ordinary differential equation of transverse motion can be obtained by solving $U_n$ and ${\phi }_n$ in terms of $W_n$ from Equations (38) and (39) and substituting the results in Equation (40). Thus, we get

$$W_n{}^{″}+{\alpha }_1{W}_n+{\alpha }_2{W}_n{}^3=0,$$

(41)

where

$${\alpha }_1={\omega }_L^2=\frac{g_{14}\left(g_{22}{g}_{31}-g_{21}{g}_{32}\right)+g_{12}\left(-g_{24}{g}_{31}+g_{21}{g}_{36}\right)+g_{11}\left(g_{24}{g}_{32}-g_{22}{g}_{36}\right)}{g_{24}\left(-g_{12}{g}_{38}+g_{11}{g}_{39}\right)+g_{14}\left(g_{22}{g}_{38}-g_{21}{g}_{39}\right)+\left(g_{12}{g}_{21}-g_{11}{g}_{22}\right)g_{310}},$$

(42)

$${\alpha }_2=\frac{g_{14}\left(-g_{23}{g}_{33}+g_{21}{g}_{35}\right)+g_{13}\left(g_{24}{g}_{33}-g_{21}{g}_{37}\right)+g_{11}\left(-g_{24}{g}_{35}+g_{23}{g}_{37}\right)}{g_{24}\left(g_{12}{g}_{38}-g_{11}{g}_{39}\right)+g_{14}\left(-g_{22}{g}_{38}+g_{21}{g}_{39}\right)+\left(-g_{12}{g}_{21}+g_{11}{g}_{22}\right)g_{310}},$$

(43)

and ${\omega }_L$ is the natural frequency of the nanocomposite sandwich beam. According to the multiple timescale approach [57,58], we approximate the solution of Equation (41) by means of the following expansion,

$$W\left(t,ϵ\right)=ϵW_1\left(T_0,T_1,T_2,\dots \right)+{ϵ}^2{W}_2\left(T_0,T_1,T_2,\dots \right)+{ϵ}^3{W}_3\left(T_0,T_1,T_2,\dots \right)+\dots$$

(44)

where $ϵ$ is a small perturbation parameter and $T_{\mathrm{n}}={ϵ}^{n}t$ refers to the independent variables for $n=0,1,2,\dots$, whose derivatives with respect to $t$ are defined as follows:

$$\begin{array}{c}\frac{d}{d\mathrm{t}}=\frac{d{T}_0}{d\mathrm{t}}\frac{\partial }{\partial T_0}+\frac{d{T}_1}{d\mathrm{t}}\frac{\partial }{\partial T_1}+\dots =D_0+ϵD_1+\dots \\ \frac{d^2}{d{\mathrm{t}}^2}=D_0^2+2ϵD_0{D}_1+{ϵ}^2\left(D_1^2+2D_0{D}_2\right)+2{ϵ}^3{D}_1{D}_2+{ϵ}^4{D}_2^2+\dots \end{array}$$

(45)

where $D_n=\partial /\partial T_n$. In our case, we apply the expansion up to $O\left({ϵ}^3\right)$, such that we need $T_0,T_1$ and $T_2$. By substitution of Equations (44) and (45) into Equation (41), expanding and equating coefficients of $ϵ$, ${ϵ}^2$, and ${ϵ}^3$ to zero, we get the following relations:

$$\mathrm{Order}ϵ:D_0^2{W}_1+{\omega }_L^2{W}_1=0$$

(46)

$$\mathrm{Order}{ϵ}^2:D_0^2{W}_2+{\omega }_L^2{W}_2=-2D_0{D}_1{W}_1$$

(47)

$$\mathrm{Order}{ϵ}^3:D_0^3{W}_3+{\omega }_L^2{W}_3=-2D_0{D}_2{W}_2-D_1^2{W}_1-2D_0{D}_2{W}_1-{\alpha }_2{W}_1^3$$

(48)

The solution of Equation (46) takes the following form:

$$W_1=A\left(T_1,T_2\right)\exp \left(i{\omega }_L{T}_0\right)+\overline{A}\exp \left(-i{\omega }_L{T}_0\right)$$

(49)

where $A$ is an unknown complex function and $\overline{A}$ is its complex conjugate. By substitution of Equation (49) into Equation (46) we obtain the following relation:

$$D_0^2{W}_2+{\omega }_L^2{W}_2=-2i{\omega }_L{D}_{1}A\exp \left(i{\omega }_L{T}_0\right)+cc$$

(50)

$cc$ being the complex conjugate of the previous term. Any particular solution of Equation (50) has a secular term containing the factor $T_0\exp \left(i{\omega }_L{T}_0\right)$ unless $D_{1}A=0$. This means that $A$ is independent of $T_1$, whereby the solution of Equation (50) is verified to be identically null.

By substitution of $W_2=0$, together with Equation (49), into Equation (48) we get the following expression

$$D_0^2{W}_3+{\omega }_L^2{W}_3=-\left[2i{\omega }_L{D}_{2}A-3{\alpha }_2{A}^2\overline{A}\right]\exp \left(i{\omega }_L{T}_0\right)-{\alpha }_2{A}^3\exp \left(i{\omega }_L{T}_0\right)$$

(51)

In this last relation the secular terms containing $\exp \left(i{\omega }_L{T}_0\right)$ must be equal to zero to have a periodic solution, which corresponds to enforce the following relation:

$$2i{\omega }_L{D}_{2}A-3{\alpha }_2{A}^2\overline{A}=0$$

(52)

whose solution can be found by defining $A$ as

$$A=\frac{1}{2}a\exp \left(i\beta \right)$$

(53)

where $a$ and $\beta$ are real functions of $T_2$.

By substituting Equation (53) into Equation (52) and by equating the real and imaginary parts to zero, we obtain

$${\omega }_{L}a\prime =0\mathrm{and}{\omega }_{L}a{\beta }^{\prime }-3/8{\alpha }_2{a}^3=0$$

(54)

where the prime denotes the derivative with respect to $T_2$. Solving both relations in Equation (54), it follows that $a$ is a constant and

$$\beta =3/8\frac{{\alpha }_2}{{\omega }_L^{}}{a}^2{T}_2+{\beta }_0$$

(55)

where ${\beta }_0$ is a constant. By combination of Equations (53) and (55) with Equation (49), we obtain the following closed-form solution for the nonlinear frequency of the transverse vibration of GPL-reinforced sandwich beams based on third-order shear deformation theory:

$${\omega }_{NL}={\omega }_L\left(1+\frac{3}{8}\frac{{\alpha }_2}{{\omega }_L^2}{\epsilon }^2{a}^2\right),$$

(56)

where $a$ is the amplitude of the vibration.

4. Numerical Results

In this section, we present the numerical results from a large parametric investigation into the nonlinear vibration behavior of thick sandwich beams reinforced with GPLs. To check for the accuracy of the proposed model, our numerical results obtained for FG GPL-reinforced beams are compared with those available from the literature. In this regard, the thickness of the core layer in the present model is assumed to be zero ($h_c=0$). In

Table 1. First three dimensionless natural frequencies (${\overline{\omega }}_L={\omega }_L\times L\times \sqrt{m_{10}/A_{10}}$) of simply supported laminated beams reinforced with GPLs ($L/h_t=20$).

Pattern	Reference	Mode
		1	2	3
UD ($\kappa =0$)	Feng et al. [30]	0.21542	0.85226	1.88292
	Present	0.23482	0.92903	2.05352
FG-X ($\kappa =1$)	Feng et al. [30]	0.25853	1.01309	2.20666
	Present	0.26759	1.05199	2.30258

and

Table 2. First three dimensionless nonlinear frequencies (${\overline{\omega }}_{NL}={\omega }_{NL}\times L\times \sqrt{m_{10}/A_{10}}$) of simply supported laminated beams reinforced with GPLs ($L/h_t=20$).

Pattern	Reference	Mode
		1	2	3
UD ($\kappa =0$)	Feng et al. [30]	0.27259	1.07270	2.33122
	Present	0.29020	1.08015	2.41741
FG-X ($\kappa =1$)	Feng et al. [30]	0.31973	1.20509	2.54097
	Present	0.31619	1.18544	2.62712

, the dimensionless free linear (${\overline{\omega }}_L={\omega }_L\times L\times \sqrt{m_{10}/A_{10}}$) and nonlinear (${\overline{\omega }}_{NL}={\omega }_{NL}\times L\times \sqrt{m_{10}/A_{10}}$) frequencies are provided for a simply supported, GPL-reinforced beam and compared with numerical results reported by Feng et al. [30]. The reference model is based on a Timoshenko beam theory, whereby two different patterns are considered for validation purposes. The following properties are assumed for the beam: $E_M=2.85\mathrm{GPa},$ ${\rho }_M=1200\mathrm{Kg}/{\mathrm{m}}^3,E_{\mathrm{GPL}}=1.01\mathrm{TPa},$ ${\rho }_{\mathrm{GPL}}=1062.5\mathrm{Kg}/{\mathrm{m}}^3,$ $w_{\mathrm{GPL}}=1.5\mathsf{\mu }\mathrm{m},$ $l_{\mathrm{GPL}}=2.5\mathsf{\mu }\mathrm{m}$, and $h_{\mathrm{GPL}}=1.5\mathrm{nm}$.

In

Table 3. Comparison of the first four dimensionless natural frequencies (${\overline{\omega }}_L={\omega }_L\times \left(L^2/h\right)\sqrt{\rho /E_{11}}$) of orthotropic thick beams based on different higher-order shear deformable beam theories ($L/h_t=10$).

Reference	Mode
	1	2	3	4
Shen et al. [61]	2.3100	6.9538	11.9707	17.0393
Li and Qiao [60]	2.3188	7.0204	12.0894	17.3139
Vo and Thai [59]	2.3198	7.0091	12.1250	17.2949
Present	2.4038	7.2110	12.3975	17.6279

, a comparison has been attempted between results from the present formulation for the first four dimensionless frequencies (${\overline{\omega }}_L={\omega }_L\times \left(L^2/h\right)\sqrt{\rho /E_{11}}$) of simply supported orthotropic beams and those from the literature, based on different higher-order shear deformation theories [59,60,61]. The material properties are assumed to be $E_{11}=144.9 \mathrm{GPa},$ $E_{22}=9.65 \mathrm{GPa},$ $G_{12}=G_{13}=4.14 \mathrm{GPa},$ $G_{23}=3.45 \mathrm{GPa},$ $\rho =1389.23\mathrm{Kg}/{\mathrm{m}}^3$, and ${\nu }_{12}={\nu }_{21}=0.3$.

As is clearly visible in Table 3, the numerical results based on the proposed model agree very well with predictions from the literature based on other shear deformable models. It seems that higher-order-models available in the literature [59,60,61] get more conservative results than our formulation. At the same time, the proposed multiple timescale approach proves to be an efficient analytical tool to solve nonlinear systems in a very easy and straightforward manner.

In the benchmark

Table 4. First-order nonlinear frequency of a simply supported GPL-reinforced sandwich beam ($h_c/h_t=0.6,N=10,a/h_t=1$).

$\mathit{L}/{\mathit{h}}_{\mathit{t}}$	${\mathit{\Gamma }}_{\mathit{b}}$	$\mathit{\kappa }$
		0.5	1	2	5
10	0.5	1.52551	1.53220	1.54115	1.54952
	1	1.96610	1.97452	1.98470	1.99110
	2	2.62940	2.63683	2.64207	2.63204
20	0.5	0.38485	0.38708	0.39006	0.39311
	1	0.49821	0.50165	0.50612	0.51036
	2	0.66931	0.67405	0.67978	0.68391

and

Table 5. First-order nonlinear-to-linear frequency ratio of a simply supported GPL-reinforced sandwich beam ($h_c/h_t=0.6,N=10,a/h_t=1$).

$\mathit{L}/{\mathit{h}}_{\mathit{t}}$	${\mathit{\Gamma }}_{\mathit{b}}$	$\mathit{\kappa }$
		0.5	1	2	5
10	0.5	1.59007	1.56735	1.54135	1.51106
	1	1.55652	1.53643	1.51484	1.49074
	2	1.54324	1.52830	1.51607	1.50875
20	0.5	1.55431	1.53007	1.50227	1.47047
	1	1.51103	1.48708	1.46023	1.42993
	2	1.48700	1.46431	1.43990	1.41399

, we report the numerical results in terms of the first-order nonlinear frequency and nonlinear-to-linear frequency ratio for a GPL-reinforced sandwich beam. The material and geometrical properties of the beam are considered to be $E_M=2.85\mathrm{GPa},$ ${\rho }_M=1200\mathrm{Kg}/{\mathrm{m}}^3,E_{\mathrm{GPL}}=1.01\mathrm{TPa},$ ${\rho }_{\mathrm{GPL}}=1062.5\mathrm{Kg}/{\mathrm{m}}^3,$ $w_{\mathrm{GPL}}=1.5\mathsf{\mu }\mathrm{m},$ $l_{\mathrm{GPL}}=2.5\mathsf{\mu }\mathrm{m},$ and $h_{\mathrm{GPL}}=1.5\mathrm{nm}$. These properties are kept constant for the following examples. The numerical results are obtained for a different power-law parameter, length-to-total thickness of the beam ratio, as well as for a different total weight fraction of the GPLs reinforced in the beam (${\Gamma }_b$). As mentioned before, ${\Gamma }_b$ is defined such that the total amount of GPLs remains constant with respect to any change in the power-law parameter or thickness of the face sheets. According to Table 4 and Table 5, an increased total weight fraction of GPLs (${\Gamma }_b$) yields a meaningful increase of the nonlinear frequency of the beam, while decreasing the nonlinear-to-linear frequency ratio. In addition, an increased length-to-thickness ratio provides a decreasing effect on the nonlinear frequency of the system and its associated nonlinear-to-linear ratio.

On the other hand, an increased power-law parameter gets an increased nonlinear frequency and a decreased nonlinear-to-linear frequency ratio. A non-uniform behavior can be observed, sometimes, for an increasing power-law parameter and for a large amount of ${\Gamma }_b$. This aspect is illustrated in detail as follows.

4.1. Effect of the Amplitude of the Vibrations

Figure 3 and Figure 4 show the effect of an increasing amplitude on the nonlinear frequency and the nonlinear-to-linear frequency ratio. According to both figures, the nonlinear frequency of a GPL-reinforced sandwich beam generally increases by increasing the amplitude of the vibrations as well as its nonlinear-to-linear frequency. However, an increased vibration amplitude significantly affects the nonlinear frequency and its rational form, for smaller values of the power-law parameter. It seems that for sandwich beams with a larger amount of GPL-reinforcement, the nonlinear frequency is increasingly affected by larger vibration amplitudes.

4.2. Effect of the Power-Law Parameter

Here, we study the effect of the power-law parameter on the nonlinear vibration behavior of the nanocomposite structure. Figure 5 depicts the variation of the nonlinear frequency as a function of the power-law parameter for different thickness ratios, $h_c/h_t$, and vibration amplitudes $a/h_t$. Note that, for a small amplitude of vibration, the nonlinear frequency of the system increases by increasing the power-law parameter. The effect of an increasing power-law parameter on the vibration response of the system is different depending on the amount of the thickness ratio. For an increased amplitude of the vibration up to a threshold value, the results become non-uniform for an increased power-law parameter. An increased amplitude of vibration will completely change the behavior of the system; namely, by increasing the power-law parameter, the nonlinear frequency of the GPL-reinforced sandwich beam decreases for a large amplitude of vibration.

In Figure 6, we plot the variation of the nonlinear-to-linear frequency ratio of the nanocomposite structure vs. the power-law parameter. It is worth noting that an increased power-law parameter has a decreasing effect on the nonlinear-to-linear frequency ratio of the system. These effects become even more pronounced for larger vibration amplitudes, while leaving the overall behavior almost unaltered.

In Figure 7 and Figure 8, the nonlinear frequency and the nonlinear-to-linear frequency ratio of a GPL-reinforced sandwich beam are plotted vs. the power-law parameter, while assuming different thickness ratios and total weight fractions of the GPL phase in the beam.

According to Figure 7, for a small amount of the GPLs weight fraction, the nonlinear frequency increases by increasing the power-law parameter. Moreover, the effect of an increasing power-law parameter becomes more pronounced for sandwich beams with thick face sheets. As the total weight fraction of GPLs increases, a different response is noticed, in terms of nonlinear frequency, by increasing the power-law parameter. More specifically, for a large GPL weight fraction, we notice a threshold value after which the nonlinear frequency decreases by increasing the power-law parameter. Of course, the value of the threshold point varies with the thickness of the face sheets. On the other hand, the results for the nonlinear-to-linear frequency ratio show some opposite effects for an increasing power-law parameter. An increased total amount of GPL reinforcement phase in the face sheets has a pronounced effect on the general behavior of the system, which has to be studied carefully.

4.3. Effect of the Thickness of the Face Sheets

Figure 9 and Figure 10 show the effect of the core-to-face thickness ratio on the nonlinear vibration response of the structure. As visible in both figures, a decreasing thickness of the face sheets generally increases the nonlinear frequency of the system for a small vibration amplitude. By increasing the amplitude vibrations, the structural response changes due to an increased thickness ratio; namely, for large amplitude vibrations, the frequency decreases by increasing the thickness ratio, and the nonlinear-to-linear frequency ratio decreases accordingly. This behavior is almost unaffected by the amplitude vibrations.

4.4. Effects of the Total Weight Fraction of the GPLs

The last parametric investigation checks for the sensitivity of the response to an increased total amount of GPL phase in the face sheets of a sandwich beam, both in terms of nonlinear frequency and nonlinear-to-linear frequency. As shown in Figure 11, an increased total weight fraction of the GPLs will increase the nonlinear frequency of the system. Different effects on the nonlinear-to-linear frequency ratio are observable, depending on the value of the selected power-law parameter.

5. Concluding Remarks

In the present paper we analyze the nonlinear free vibration response of thick sandwich beams with FG GPL-reinforced face sheets, based on a novel dispersion model for the reinforcement phase. A higher-order laminated beam model is associated with the Von Kármán strain-displacement relationships to capture the shear and rotation effects on the structural behavior of the system. The nonlinear equations of motion are determined through Hamilton’s principle, and they are discretized to ordinary differential equations by means of Galerkin’s approach. An analytical solution procedure based on the multiple timescales method is then used to obtain the nonlinear frequency equation. A large numerical investigation analyzes the effect of the vibration amplitude, the thickness of the face sheets, and the GPL dispersion on the nonlinear vibration response of the reinforced sandwich structure, where the following conclusions can be summarized as follows:

-

An increased amplitude of vibrations significantly increases the nonlinear frequency and its ratio to linear frequency. The sensitivity of the nonlinear response varies with the total amount of GPL reinforcement in the beam, as well as with the value of the power-law parameter and the thickness of the face sheets.

-

An increased power-law parameter can have different effects on the stiffness of the GPL-reinforced sandwich beam, depending on the value of the vibration amplitude and the total weight fraction of the GPLs.

-

There exists a threshold value for the vibration amplitude, after which the behavior of the system can change for an increased thickness of the face sheets and power-law parameter.

-

For low amplitude vibrations, the nonlinear frequency increases by increasing the power-law parameter and by decreasing the thickness of the faces sheets. For large amplitude vibrations, the contrary occurs by increasing both the power-law parameter and core-to-face thickness ratio.

-

The effect of an increasing power-law parameter and a decreasing thickness of the face sheets on the nonlinear-to-linear frequency ratio is independent of the vibration amplitude.

-

An increasing total weight fraction of GPLs in the beam generally increases the nonlinear frequency of the system. The sensitivity of the nonlinear-to-linear frequency ratio can be more or less pronounced, depending on the GPL dispersion pattern in the face sheets.

-

The proposed parametric study would be of great interest for optimization and design of materials and for an appropriate evaluation of stability for sandwich beams under different environmental conditions, which would prove useful in many space and aircraft applications.