**1. Introduction**

Sandwich structures, generally made of a soft core and two hard face sheets, are largely used in the aerospace, oil, gas, and petrochemicalindustries, due to their enhancedmechanical properties, namely, a high strength-to-weight ratio and a high resistance to heat, humidity, and noise [1–5]. Hence, in recent decades, much attention has been paid to the mechanical behavior of these structures [6–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–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–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–28]. Among different reinforcement possibilities, graphene nanoplatelets (GPLs) provide a uniform reinforced assembly, as well as the easiest manufacturing process, as discussed in [29–32]. Graphene is a monolayer

structure of carbon atoms with extraordinary electrical, mechanical, thermal, and optical properties [33–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–43]. This justifies the large attention paid in the literature to the mechanical behavior of structures reinforced with graphenic materials [44–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 *ht* and width *b* is considered, as shown in Figure 1. The Cartesian coordinate system (*x*, *z*) 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(z) = \Gamma\_{\text{max}} \left( \frac{2|z| - h\_{\text{f}}}{h\_{\text{f}} - h\_{\text{c}}} \right)^{\kappa},
\tag{1}
$$

where Γmax is the maximum value of the distribution function and κ is a power-law parameter, which defines the GPLs dispersion throughout the thickness of the face sheets.

**Figure 1.** General configuration and graphene nanoplatelet (GPL) dispersion description in a laminated GPL-reinforced sandwich beam.

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, Γ*b*, the maximum value of the GPL weight fraction, Γ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.

**Figure 2.** Variation of the dimensionless GPL weight fraction (Γ = Γ(*z*)/Γ*b*) through the thickness of the top face sheet with respect to various power-law parameters (*h* = (2*z* − *hc*)/2*hf*).

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^{(k)} = \Gamma\_{\text{max}} \left( \frac{2|h\_k| - h\_\emptyset}{h\_t - h\_\emptyset} \right)^\chi,\tag{2}
$$

where *hk* 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\_{\rm GPL}^{(k)} = \frac{\Gamma^{(k)}}{\Gamma^{(k)} + (\rho\_{\rm GPL}/\rho\_M) \Big(1 - \Gamma^{(k)}\Big)} \,\tag{3}$$

ρ*<sup>M</sup>* and ρ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^{(k)} = \frac{3}{8} \left( \frac{1 + \xi\_L \eta\_L V\_{\rm GPL}^{(k)}}{1 - \eta\_L V\_{\rm GPL}^{(k)}} \right) \mathbf{E}\_M + \frac{5}{8} \left( \frac{1 + \xi\_W \eta\_W V\_{\rm GPL}^{(k)}}{1 - \eta\_W V\_{\rm GPL}^{(k)}} \right) \mathbf{E}\_{M\prime} \tag{4}$$

where

$$
\eta\_W = \frac{(E\_{\rm GPL}/E\_M) - 1}{(E\_{\rm GPL}/E\_M) + \xi\_W}, \quad \eta\_L = \frac{(E\_{\rm GPL}/E\_M) - 1}{(E\_{\rm GPL}/E\_M) + \xi\_L}, \tag{5}
$$

$$
\xi\_W = \frac{2w\_{\rm GPL}}{h\_{\rm GPL}}, \ \xi\_L = \frac{2l\_{\rm GPL}}{h\_{\rm GPL}},\tag{6}
$$

and *h*GPL, *l*GPL, *w*GPL stand for the average thickness, length, and width of GPLs, respectively; *EM* and *E*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\_{\mathcal{C}}^{(k)} = \rho\_{\text{GPL}} V\_{\text{GPL}}^{(k)} + \rho\_M \Big(1 - V\_{\text{GPL}}^{(k)}\Big) \tag{7}
$$

$$\nu\_{\mathbb{C}}^{(k)} = \nu\_{\text{GPL}} V\_{\text{GPL}}^{(k)} + \nu\_{\text{M}} \Big(1 - V\_{\text{GPL}}^{(k)}\Big),\tag{8}$$

ν*<sup>M</sup>* and ν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(*x*, *t*) and *u*3(*x*, *t*) of an arbitrary point in the *x* and *z* directions for shear deformable sandwich beams can be expressed as

$$u\_1(\mathbf{x}, z, t) = u(\mathbf{x}, t) + z\phi(\mathbf{x}, t) - \frac{4z^3}{3h\_t^2} (\phi + \frac{\partial w}{\partial \mathbf{x}}) \tag{9}$$

$$w\_3(\mathbf{x}, z, \mathbf{t}) = w(\mathbf{x}, \mathbf{t}),\tag{10}$$

where *u*(*x*, *t*) and *w*(*x*, *t*) are the displacement components of a point at the mid-plane of the beam in the *x* and *z* directions, respectively. Moreover, φ(*x*, *t*) 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

$$
\varepsilon\_{\rm xx} = \varepsilon\_{\rm xx}^{(0)} + z \varepsilon\_{\rm xx}^{(1)} + z^3 \varepsilon\_{\rm xx}^{(3)}.\tag{11}
$$

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

where

$$
\varepsilon\_{\rm xx}^{(0)} = \frac{\partial u}{\partial \mathbf{x}} + \frac{1}{2} \left(\frac{\partial w}{\partial \mathbf{x}}\right)^2,\\ \varepsilon\_{\rm xx}^{(1)} = \frac{\partial \phi}{\partial \mathbf{x}},\\ \varepsilon\_{\rm xx}^{(3)} = -\mathbf{c}\_1 \left(\frac{\partial \phi}{\partial \mathbf{x}} + \frac{\partial^2 w}{\partial \mathbf{x}^2}\right). \tag{13}$$

$$\gamma\_{xx}^{(0)} = \phi + \frac{\partial w}{\partial \mathbf{x}} \; \; \; \gamma\_{xx}^{(2)} = -\mathbf{c}\_2 \Big(\phi + \frac{\partial w}{\partial \mathbf{x}}\Big) \tag{14}$$

and

$$c\_1 = \frac{4}{3h\_t^{2}}, \ c\_2 = 3c\_1 = \frac{4}{h\_t^2}. \tag{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

$$\int\_{t\_1}^{t\_2} (\delta \mathcal{U} - \delta \mathcal{W} - \delta \mathcal{K}) dt = 0,\tag{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 δ*U* for the third-order shear deformable sandwich beams reads as follows

$$\begin{array}{rcl}\delta I I &=&\int\_{A} \int\_{0}^{L} (\sigma\_{xx}\delta \varsigma\_{xx} + \sigma\_{x}\delta \gamma\_{xx}) dx dA\\ &=&\int\_{A} \int\_{0}^{L} \Big[\sigma\_{x}\Big(\delta \varsigma\_{xx}^{(0)} + z\delta \epsilon\_{xx}^{(1)} + z^{3}\delta \epsilon\_{xx}^{(3)}\Big) + \sigma\_{x}\Big[\delta \gamma\_{xx}^{(0)} + z^{2}\delta \gamma\_{xx}^{(2)}\Big]\Big] dx dA\\ &=&\int\_{0}^{L} \Big[-\frac{\partial \mathcal{N}\_{\text{Tx}}}{\partial x} \delta w - \frac{\partial}{\partial x} \Big(\mathcal{N}\_{\text{xx}}\frac{\partial \mathcal{w}}{\partial x}\Big) \delta w - \frac{\partial \mathcal{M}\_{\text{Tx}}}{\partial x} \delta \phi + c\_{1}\frac{\partial \mathcal{P}\_{\text{xx}}}{\partial x} \delta \phi - c\_{1}\frac{\partial^{2} \mathcal{P}\_{\text{xx}}}{\partial x^{2}} \delta w\\ &+&Q\_{2}\delta \phi - \frac{\partial Q\_{x}}{\partial x} \delta w - c\_{2}\mathcal{R}\_{\text{x}} \delta \phi + c\_{2}\frac{\partial \mathcal{R}\_{\text{x}}}{\partial x} \delta w\Big] dx\end{array} \tag{17}$$

where

$$
\begin{bmatrix} N\_{\text{xx}} \\ M\_{\text{xx}} \\ P\_{\text{xx}} \end{bmatrix} = \int\_{A} \sigma\_{\text{xx}} \left| \begin{array}{c} 1 \\ z \\ z^{3} \end{array} \right| dA, \quad \left[ \begin{array}{c} Q\_{\text{x}} \\ R\_{\text{x}} \end{array} \right] = \int\_{A} \sigma\_{\text{xz}} \left[ \begin{array}{c} 1 \\ z^{2} \end{array} \right] dA. \tag{18}
$$

The virtual kinetic energy δ*K* is defined as

$$\begin{array}{lcl} \delta K &= \int\_{\mathcal{A}} \int\_{0}^{L} \rho(z) \Big( \bar{u}\_{1} \delta \bar{u}\_{1} + \bar{u}\_{3} \delta \bar{u}\_{3} \Big) dA \, d\mathbf{x} \\ &= \int\_{0}^{L} \Big[ \Big( m\_{0} \bar{u} + m\_{1} \bar{\phi} - c\_{1} m\_{3} (\bar{\phi} + \frac{\partial \bar{w}}{\partial x}) \Big) \delta u \\ &+ \Big( m\_{1} \bar{u} + m\_{2} \bar{\phi} - c\_{1} m\_{3} \bar{u} - c\_{1} m\_{4} (2 \bar{\phi} + \frac{\partial \bar{w}}{\partial x}) + c\_{1}^{2} m\_{6} (\bar{\phi} + \frac{\partial \bar{w}}{\partial x}) \Big) \delta \phi \\ &+ \Big( m\_{0} \bar{w} + c\_{1} m\_{3} \frac{\partial \bar{u}}{\partial x} + c\_{1} m\_{4} \frac{\partial \bar{\phi}}{\partial x} - c\_{1}^{2} m\_{6} \Big( \frac{\partial \bar{\phi}}{\partial x} + \frac{\partial^{2} \bar{w}}{\partial x^{2}} \Big) \delta w \Big] dx, \end{array} \tag{19}$$

with

$$m\_i = \int\_A \rho(z) z^i dA = \sum\_{k=1}^{2N+1} b \int\_{z\_k}^{z\_{k+1}} \rho\_c^{(k)} z^i dz, \quad i = 0, 1, 2, 3, 4, 6. \tag{20}$$

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

$$\frac{\partial N\_{\text{xx}}}{\partial \mathbf{x}} = m\_0 \ddot{\boldsymbol{\mu}} + m\_1 \ddot{\boldsymbol{\phi}} - \mathbf{c}\_1 m\_3 \left(\ddot{\boldsymbol{\phi}} + \frac{\partial \ddot{\boldsymbol{w}}}{\partial \mathbf{x}}\right) \tag{21}$$

$$\frac{\partial M\_{\text{xx}}}{\partial \mathbf{x}} - c\_1 \frac{\partial P\_{\text{xx}}}{\partial \mathbf{x}} - Q\_{\text{x}} + c\_2 R\_{\text{x}} = m\_1 \ddot{u} + m\_2 \ddot{\phi} - c\_1 m\_3 \ddot{u} - c\_1 m\_4 \left(2 \ddot{\phi} + \frac{\partial \ddot{w}}{\partial \mathbf{x}}\right) + c\_1^2 m\_6 \left(\ddot{\phi} + \frac{\partial \ddot{w}}{\partial \mathbf{x}}\right) \tag{22}$$

$$\frac{\partial}{\partial \mathbf{x}} \left( \mathbf{N}\_{\mathbf{x}\mathbf{x}} \frac{\partial \mathbf{w}}{\partial \mathbf{x}} \right) + c\_1 \frac{\partial^2 P\_{\mathbf{x}\mathbf{x}}}{\partial \mathbf{x}^2} + \frac{\partial Q\_{\mathbf{x}}}{\partial \mathbf{x}} - c\_2 \frac{\partial \mathbf{R}\_{\mathbf{x}}}{\partial \mathbf{x}} = m\_0 \ddot{w} + c\_1 m\_3 \frac{\partial \ddot{u}}{\partial \mathbf{x}} + c\_1 m\_4 \frac{\partial \ddot{\phi}}{\partial \mathbf{x}} - c\_1^2 m\_6 \left( \frac{\partial \ddot{\phi}}{\partial \mathbf{x}} + \frac{\partial^2 \ddot{w}}{\partial \mathbf{x}^2} \right). \tag{23}$$

*Appl. Sci.* **2020**, *10*, 5669

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

$$N\_{\rm xx} = A\_{11} \left( \frac{\partial u}{\partial \mathbf{x}} + \frac{1}{2} \left( \frac{\partial w}{\partial \mathbf{x}} \right)^2 \right) + B\_{11} \left( \frac{\partial \phi}{\partial \mathbf{x}} \right) - c\_1 E\_{11} \left( \frac{\partial \phi}{\partial \mathbf{x}} + \frac{\partial^2 w}{\partial \mathbf{x}^2} \right) \tag{24}$$

$$M\_{\rm xx} = B\_{11} \left( \frac{\partial \mathbf{u}}{\partial \mathbf{x}} + \frac{1}{2} \left( \frac{\partial \mathbf{w}}{\partial \mathbf{x}} \right)^2 \right) + D\_{11} \left( \frac{\partial \phi}{\partial \mathbf{x}} \right) - c\_1 F\_{11} \left( \frac{\partial \phi}{\partial \mathbf{x}} + \frac{\partial^2 \mathbf{w}}{\partial \mathbf{x}^2} \right) \tag{25}$$

$$P\_{\rm xx} = E\_{11} \left( \frac{\partial \mu}{\partial \mathbf{x}} + \frac{1}{2} \left( \frac{\partial w}{\partial \mathbf{x}} \right)^2 \right) + F\_{11} \left( \frac{\partial \phi}{\partial \mathbf{x}} \right) - c\_1 H\_{11} \left( \frac{\partial \phi}{\partial \mathbf{x}} + \frac{\partial^2 w}{\partial \mathbf{x}^2} \right) \tag{26}$$

$$Q\_{\rm x} = (A\_{55} - c\_2 D\_{55}) \left(\phi + \frac{\partial w}{\partial \mathbf{x}}\right) \tag{27}$$

$$R\_{\mathbf{x}} = (D\_{\mathbf{55}} - c\_2 F\_{\mathbf{55}}) \Big( \phi + \frac{\partial w}{\partial \mathbf{x}} \Big) \tag{28}$$

where

$$(A\_{11}, B\_{11}, D\_{11}, E\_{11}, F\_{11}, H\_{11}) = \sum\_{k=1}^{N} b \int\_{\eta\_k}^{\eta\_{k+1}} Q\_{11}^{(k)}(1, z, z^2, z^3, z^4, z^6) dz,\tag{29}$$

$$(A\_{55}, D\_{55}, F\_{55}) = \sum\_{k=1}^{N} b \int\_{h\_k}^{h\_{k+1}} Q\_{55}^{(k)}(1, z^2, z^4) dz,\tag{30}$$

and

$$Q\_{11}^{(k)} = \frac{E\_{\mathbb{C}}^{(k)}}{1 - \left(\nu\_{\mathbb{C}}^{(k)}\right)^2}, \quad Q\_{55}^{(k)} = G\_{\mathbb{C}}^{(k)} = \frac{E\_{\mathbb{C}}^{(k)}}{2\left(1 + \nu\_{\mathbb{C}}^{(k)}\right)},\tag{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{split} A\_{11} \Big( \frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 w}{\partial \mathbf{x}^2} \frac{\partial w}{\partial \mathbf{x}} \Big) &\quad + (B\_{11} - c\_1 E\_{11}) \frac{\partial^2 \phi}{\partial \mathbf{x}^2} - c\_1 E\_{11} \frac{\partial^3 w}{\partial \mathbf{x}^3} \\ &= m\_0 \frac{\partial^2 u}{\partial t^2} + (m\_1 - c\_1 m\_3) \frac{\partial^2 \phi}{\partial t^2} - c\_1 m\_3 \frac{\partial^3 w}{\partial x \partial t^2} \Big) \end{split} \tag{32}$$

$$\begin{array}{ll} (B\_{11} - c\_1 \mathbf{E}\_{11}) \Big( \frac{\partial^2 \mu}{\partial \mathbf{z}^2} + \frac{\partial^2 \mu}{\partial \mathbf{z}^2} \frac{\partial \mu}{\partial \mathbf{z}} \Big) & + \Big( D\_{11} - 2c\_1 \mathbf{F}\_{11} + c\_1^2 \mathbf{H}\_{11} \Big) \frac{\partial^2 \phi}{\partial \mathbf{z}^2} + \Big( -c\_1 \mathbf{F}\_{11} + c\_1^2 \mathbf{H}\_{11} \Big) \frac{\partial^3 \mu}{\partial \mathbf{z}^3} \\ & + \Big( -A\_{55} + 2c\_2 \mathbf{D}\_{55} - c\_2^2 \mathbf{F}\_{55} \Big) \Big( \phi + \frac{\partial \mu}{\partial \mathbf{z}} \Big) \\ = (m\_1 - c\_1 m\_3) \frac{\partial^2 \mu}{\partial t^2} + \Big( m\_2 - 2c\_1 m\_4 + c\_1^2 m\_6 \right) \frac{\partial^3 \phi}{\partial t^2} + \Big( -c\_1 m\_4 + c\_1^2 m\_6 \Big) \frac{\partial^3 \mu}{\partial \mathbf{z} \partial \mathbf{z}^2}, \end{array} \tag{33}$$

$$\begin{split} c\_{1}E\_{11}\Big(\frac{\partial^{3}u}{\partial x^{3}}-\left(\frac{\partial^{2}w}{\partial x^{2}}\right)^{2}\Big) &+ \Big(c\_{1}F\_{11}-c\_{1}^{2}H\_{11}\Big)\frac{\partial^{3}\phi}{\partial x^{3}}-c\_{1}^{2}H\_{11}\frac{\partial^{4}w}{\partial x^{4}} \\ &+ \Big(A\xi\overline{s}-2c\_{2}D\xi\overline{s}+c\_{2}^{2}F\_{55}\Big)\Big(\frac{\partial\phi}{\partial x}+\frac{\partial^{2}w}{\partial x^{2}}\Big) \\ &+ A\_{1}\Big(\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}\Big(\frac{\partial^{2}w}{\partial x^{2}}\Big)\Big) \\ &+ (B\_{11}-c\_{1}E\_{11})\Big(\frac{\partial^{2}\phi}{\partial x^{2}}\frac{\partial w}{\partial x}+\frac{\partial\phi}{\partial x}\frac{\partial^{2}w}{\partial x^{2}}\Big) \\ &= m\_{0}\frac{\partial^{2}w}{\partial t^{2}}+c\_{1}m\_{3}\frac{\partial^{3}u}{\partial x\partial t^{2}}+\Big(c\_{1}m\_{4}-c\_{1}^{2}m\_{6}\Big)\frac{\partial^{3}\phi}{\partial x\partial t^{2}}-c\_{1}^{2}m\_{6}\frac{\partial^{4}w}{\partial x^{2}\partial t^{2}}.\end{split} \tag{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,

$$\mu(\mathbf{x},t) = \sum\_{n=1}^{\infty} \mathcal{U}\_n(t) \chi\_n(\mathbf{x}), \quad \chi\_n(\mathbf{x}) = \text{Cost}\left(\frac{n\pi\mathbf{x}}{L}\right) \tag{35}$$

$$\phi(\mathbf{x},t) = \sum\_{n=1}^{\infty} \varphi\_n(t)\psi\_n(\mathbf{x}), \quad \psi\_n(\mathbf{x}) = \text{Cost}\left(\frac{n\pi\mathbf{x}}{L}\right) \tag{36}$$

$$w(\mathbf{x},t) = \sum\_{n=1}^{\infty} \mathcal{W}\_n(t)\lambda\_n(\mathbf{x}), \quad \lambda\_n(\mathbf{x}) = \text{Sim}\left(\frac{n\pi\mathbf{x}}{L}\right),\tag{37}$$

*Un*, ϕ*<sup>n</sup>* and *Wn* being the unknown generalized coordinates which stand for the amplitude of the vibration. Moreover, the following functions χ*n*(*x*), ψ*n*(*x*), and λ*n*(*x*) 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 χ, ψ, and λ and after their integration over the domain of the system, the following nonlinear differential equations of motion are obtained

$$
\lg\_{11} \mathcal{U}\_n + \lg\_{12} \mathcal{W}\_n + \lg\_{13} \mathcal{W}\_n^2 + \lg\_{14} \mathcal{q}\_n = 0,\tag{38}
$$

$$
\mathcal{g}\_{21}\mathcal{U}\_n + \mathcal{g}\_{22}\mathcal{W}\_n + \mathcal{g}\_{23}\mathcal{W}\_n^2 + \mathcal{g}\_{24}\mathcal{q}\_n = 0,\tag{39}
$$

$$\mathcal{G}\_{31}\mathcal{U}\_{\boldsymbol{n}} + \mathcal{g}\_{32}\mathcal{W}\_{\boldsymbol{n}} + \mathcal{g}\_{33}\mathcal{U}\_{\boldsymbol{n}}\mathcal{W}\_{\boldsymbol{n}} + \mathcal{g}\_{34}\mathcal{W}\_{\boldsymbol{n}}^{2} + \mathcal{g}\_{35}\mathcal{W}\_{\boldsymbol{n}}^{3} + \mathcal{g}\_{36}\wp\_{\boldsymbol{n}} + \mathcal{g}\_{37}\mathcal{W}\_{\boldsymbol{n}}\wp\_{\boldsymbol{n}} + \mathcal{g}\_{38}\mathcal{U}\_{\boldsymbol{n}}\mathcal{"} + \mathcal{g}\_{39}\mathcal{W}\_{\boldsymbol{n}}\mathcal{"} + \mathcal{g}\_{310}\wp\_{\boldsymbol{n}} = 0,\tag{40}$$

where coefficients *g*11, *g*12, ... , *g*<sup>310</sup> are detailed in Appendix A. The ordinary differential equation of transverse motion can be obtained by solving *Un* and ϕ*<sup>n</sup>* in terms of *Wn* from Equations (38) and (39) and substituting the results in Equation (40). Thus, we get

$$\mathcal{W}\_{n}{}^{\prime\prime} + a\_{1}\mathcal{W}\_{n} + a\_{2}\mathcal{W}\_{n}{}^{3} = 0,\tag{41}$$

where

$$a\_1 = a\_L^2 = \frac{g\_{14}(g\_{22}g\_{31} - g\_{21}g\_{32}) + g\_{12}(-g\_{24}g\_{31} + g\_{21}g\_{36}) + g\_{11}(g\_{24}g\_{32} - g\_{22}g\_{36})}{g\_{24}(-g\_{12}g\_{38} + g\_{11}g\_{39}) + g\_{14}(g\_{22}g\_{38} - g\_{21}g\_{39}) + (g\_{12}g\_{21} - g\_{11}g\_{22})g\_{310}},\tag{42}$$

$$a\_{2} = \frac{\text{g}\_{14}(-\text{g}\_{23}\text{g}\_{33} + \text{g}\_{21}\text{g}\_{35}) + \text{g}\_{13}(\text{g}\_{24}\text{g}\_{33} - \text{g}\_{21}\text{g}\_{37}) + \text{g}\_{11}(-\text{g}\_{24}\text{g}\_{35} + \text{g}\_{23}\text{g}\_{37})}{\text{g}\_{24}(\text{g}\_{12}\text{g}\_{38} - \text{g}\_{11}\text{g}\_{39}) + \text{g}\_{14}(-\text{g}\_{22}\text{g}\_{38} + \text{g}\_{21}\text{g}\_{39}) + (-\text{g}\_{12}\text{g}\_{21} + \text{g}\_{11}\text{g}\_{22})\text{g}\_{310}},\tag{43}$$

and ω*<sup>L</sup>* 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,

$$\mathcal{W}(t,\varepsilon) = \varepsilon \mathcal{W}\_1(T\_0, T\_1, T\_2, \dots) + \varepsilon^2 \mathcal{W}\_2(T\_0, T\_1, T\_2, \dots) + \varepsilon^3 \mathcal{W}\_3(T\_0, T\_1, T\_2, \dots) + \dots \tag{44}$$

where is a small perturbation parameter and *T*<sup>n</sup> = *nt* refers to the independent variables for *n* = 0, 1, 2, ..., whose derivatives with respect to *t* are defined as follows:

$$\begin{aligned} \frac{d}{dt} &= \frac{dT\_0}{dt} \frac{\partial}{\partial T\_0} + \frac{dT\_1}{dt} \frac{\partial}{\partial T\_1} + \dots = D\_0 + \epsilon D\_1 + \dots \\ \frac{d^2}{dt^2} &= D\_0^2 + 2\epsilon D\_0 D\_1 + \epsilon^2 \left( D\_1^2 + 2D\_0 D\_2 \right) + 2\epsilon^3 D\_1 D\_2 + \epsilon^4 D\_2^2 + \dots \end{aligned} \tag{45}$$

where *Dn* = ∂/∂*Tn*. In our case, we apply the expansion up to *O* 3 , such that we need *T*0, *T*<sup>1</sup> and *T*2. By substitution of Equations (44) and (45) into Equation (41), expanding and equating coefficients of , 2, and <sup>3</sup> to zero, we get the following relations:

$$\text{Order\\_er: } D\_0^2 \mathcal{W}\_1 + \omega\_L^2 \mathcal{W}\_1 = 0 \tag{46}$$

$$\text{Order } \epsilon^2 \text{: } D\_0^2 \mathcal{W}\_2 + a\_L^2 \mathcal{W}\_2 = -2D\_0 D\_1 \mathcal{W}\_1 \tag{47}$$

$$\text{Order } \epsilon^3 \text{: } D\_0^3 \mathcal{W}\_3 + \omega\_L^2 \mathcal{W}\_3 = -2D\_0 D\_2 \mathcal{W}\_2 - D\_1^2 \mathcal{W}\_1 - 2D\_0 D\_2 \mathcal{W}\_1 - a\_2 \mathcal{W}\_1^3 \tag{48}$$

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

$$\mathcal{W}\_1 = A \begin{pmatrix} T\_1 \ T\_2 \end{pmatrix} \exp(i\omega \iota\_\* T\_0) + \overline{A} \exp(-i\omega \iota\_\* T\_0) \tag{49}$$

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

$$D\_0^2 \mathcal{W}\_2 + \omega\_L^2 \mathcal{W}\_2 = -2i\omega\_L D\_1 A \exp\left(i\omega\_L T\_0\right) + c\varepsilon \tag{50}$$

*cc* being the complex conjugate of the previous term. Any particular solution of Equation (50) has a secular term containing the factor *T*<sup>0</sup> exp(*i*ω*LT*0) 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*<sup>2</sup> = 0, together with Equation (49), into Equation (48) we get the following expression

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

In this last relation the secular terms containing exp(*i*ω*LT*0) must be equal to zero to have a periodic solution, which corresponds to enforce the following relation:

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

whose solution can be found by defining *A* as

$$A = \frac{1}{2}a\exp(i\beta)\tag{53}$$

where *a* and β 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 \nu = 0 \text{ and } a \omega\_L a \beta' - 3/8a\_2 a^3 = 0 \tag{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}{\alpha\_L} a^2 T\_2 + \beta\_0 \tag{55}
$$

where β<sup>0</sup> 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\_{\rm NL} = \omega\_{\rm L} \left( 1 + \frac{3}{8} \frac{\alpha\_2}{\omega\_L^2} \varepsilon^2 a^2 \right) \tag{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 (*hc* = 0). In Tables 1 and 2, the dimensionless free linear (ω*<sup>L</sup>* = <sup>ω</sup>*<sup>L</sup>* <sup>×</sup> *<sup>L</sup>* <sup>×</sup> <sup>+</sup> *<sup>m</sup>*10/*A*10) and nonlinear (ω*NL* = <sup>ω</sup>*NL* <sup>×</sup> *<sup>L</sup>* <sup>×</sup> <sup>+</sup> *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: *EM* = 2.85 GPa, ρ*<sup>M</sup>* = 1200Kg/m3, *E*GPL = 1.01 TPa, ρGPL = 1062.5Kg/m3, *w*GPL = 1.5 μm, *l*GPL = 2.5 μm, and *h*GPL = 1.5 nm.

**Table 1.** First three dimensionless natural frequencies (ω*<sup>L</sup>* = <sup>ω</sup>*<sup>L</sup>* <sup>×</sup> *<sup>L</sup>* <sup>×</sup> <sup>+</sup> *m*10/*A*10) of simply supported laminated beams reinforced with GPLs (*L*/*ht* = 20).


**Table 2.** First three dimensionless nonlinear frequencies (ω*NL* = <sup>ω</sup>*NL* <sup>×</sup> *<sup>L</sup>* <sup>×</sup> <sup>+</sup> *m*10/*A*10) of simply supported laminated beams reinforced with GPLs (*L*/*ht* = 20).


In Table 3, a comparison has been attempted between results from the present formulation for the first four dimensionless frequencies (ω*<sup>L</sup>* = ω*<sup>L</sup>* × *L*2/*h* + ρ/*E*11) of simply supported orthotropic beams and those from the literature, based on different higher-order shear deformation theories [59–61]. The material properties are assumed to be *E*<sup>11</sup> = 144.9 GPa, *E*<sup>22</sup> = 9.65 GPa, *G*<sup>12</sup> = *G*<sup>13</sup> = 4.14 GPa, *G*<sup>23</sup> = 3.45 GPa, ρ = 1389.23Kg/m3, and ν<sup>12</sup> = ν<sup>21</sup> = 0.3.

**Table 3.** Comparison of the first four dimensionless natural frequencies (ω*<sup>L</sup>* = ω*<sup>L</sup>* × *L*2/*h* + ρ/*E*11) of orthotropic thick beams based on different higher-order shear deformable beam theories (*L*/*ht* = 10).


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–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 Tables 4 and 5, 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 *EM* = 2.85 GPa, ρ*<sup>M</sup>* = 1200Kg/m3, *E*GPL = 1.01 TPa, ρGPL = 1062.5Kg/m3, *w*GPL = 1.5 μm, *l*GPL = 2.5 μm, and *h*GPL = 1.5 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 (Γ*b*). As mentioned before, Γ*<sup>b</sup>* 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 Tables 4 and 5, an increased total weight fraction of GPLs (Γ*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.

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


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


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 Γ*b*. This aspect is illustrated in detail as follows.

#### *4.1. E*ff*ect of the Amplitude of the Vibrations*

Figures 3 and 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.

**Figure 3.** Variation of the first-order nonlinear frequency and nonlinear-to-linear frequency ratio of a GPL-reinforced sandwich beam with respect to an increasing amplitude of vibrations, and for different power-law parameters (*hc*/*ht* = 0.6, *N* = 10, Γ*<sup>b</sup>* = 1%, *L*/*ht* = 10).

**Figure 4.** Variation of the first-order nonlinear frequency and nonlinear-to-linear frequency ratio of a GPL-reinforced sandwich beam with respect to an increasing amplitude of vibrations, and for a different total amount of GPLs in the beam (*hc*/*ht* = 0.6, *N* = 10, *a*/*ht* = 1, *L*/*ht* = 10, κ = 2).

#### *4.2. E*ff*ect 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, *hc*/*ht*, and vibration amplitudes *a*/*ht*. 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.

**Figure 5.** Variation of the first-order nonlinear frequency of a GPL-reinforced sandwich beam due to increasing power-law parameter for different core-to-beam thickness ratios and amplitudes of vibration (*N* = 10, Γ*<sup>b</sup>* = 1%, *L*/*ht* = 10).

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.

**Figure 6.** Variation of the first-order nonlinear-to-linear frequency ratio of a GPL-reinforced sandwich beam due to an increased power-law parameter for different core-to-beam thickness ratios and amplitudes of vibration (*N* = 10, Γ*<sup>b</sup>* = 1%, *L*/*ht* = 10).

In Figures 7 and 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.

**Figure 7.** Variation of the first-order nonlinear frequency of a GPL-reinforced sandwich beam due to an increasing power-law parameter, for different core-to-beam thickness ratios and total amount of GPLs reinforcement in the beam (*N* = 10, *a*/*ht* = 1, *L*/*ht* = 10).

**Figure 8.** Variation of the first-order nonlinear-to-linear frequency ratio of a GPL-reinforced sandwich beam due to an increasing power-law parameter for different core-to-beam thickness ratios and total amounts of GPLs reinforcements in the beam (*N* = 10, *a*/*ht* = 1, *L*/*ht* = 10).

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. E*ff*ect of the Thickness of the Face Sheets*

Figures 9 and 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.

**Figure 9.** Variation of the first-order nonlinear frequency of a GPL-reinforced sandwich beam with respect to an increasing core-to-face thickness ratio for different power-law parameters and amplitudes of vibration (*N* = 10, Γ*<sup>b</sup>* = 1%, *L*/*ht* = 10).

**Figure 10.** Variation of the first-order nonlinear-to-linear frequency ratio of a GPL-reinforced sandwich beam with respect to an increasing core-to-beam thickness ratio for different power-law parameters and amplitudes of vibration (*N* = 10, Γ*<sup>b</sup>* = 1%, *L*/*ht* = 10).

#### *4.4. E*ff*ects 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.

**Figure 11.** Variation of the first-order nonlinear frequency and nonlinear-to-linear frequency ratio of a GPL-reinforced sandwich beam with respect to an increasing total amount of GPLs in the beam, for different power-law parameters (*hc*/*ht* = 0.6, *N* = 10, *a*/*ht* = 3, *L*/*ht* = 10).

#### **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:



**Author Contributions:** Formal analysis, M.S.N., H.M., R.D. and F.T.; investigation, M.S.N., H.M., R.D. and F.T.; methodology, H.M., R.D. and F.T.; software, H.M., F.T. and R.D.; supervision, F.T.; validation, M.S.N. and R.D.; writing—original draft, M.S.N. and H.M.; writing—review and editing, R.D. and F.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

In the following, we provide the extended relations for the coefficients *g*11, *g*12, ... in Equations (38)–(40), i.e.,

$$\mathcal{g}\_{11} = \frac{n^2 \pi^2 A\_{11}}{2L},\tag{A1}$$

$$\mathcal{g}\_{12} = -\frac{n^3 \pi^3 c\_1 E\_{11}}{2L^2},\tag{A2}$$

$$g\_{13} = \frac{\left(1 + (-1)^{n+1}\right) n^2 \pi^2 A\_{11}}{3L^2},\tag{A3}$$

$$\lg\_{14} = \frac{n^2 \pi^2 (B\_{11} - c\_1 E\_{11})}{2L},\tag{A4}$$

$$\text{g}\_{21} = \frac{n^2 \pi^2 (B\_{11} - c\_1 E\_{11})}{2L},\tag{A5}$$

$$\mathbf{g}\_{22} = \frac{n^3 \pi^3 \left(-c\_1 F\_{11} + c\_1^2 H\_{11}\right)}{2L^2} - \frac{n \pi \left(-A\_{55} + 2c\_2 D\_{55} - c\_2^2 F\_{55}\right)}{2},\tag{A6}$$

$$\mathcal{G}\varepsilon\varepsilon = \frac{\left(1 + (-1)^{n+1}\right) n^2 \pi^2 (B\_{11} - c\_1 E\_{11})}{3L^2},\tag{A7}$$

$$\mathcal{g}\_{24} = \frac{n^2 \pi^2 \left(D\_{11} - 2c\_1 F\_{11} + c\_1^2 H\_{11}\right)}{2L} - \frac{L\left(-A\_{55} + 2c\_2 D\_{55} - c\_2^2 F\_{55}\right)}{2},\tag{A8}$$

$$g\_{31} = \frac{n^3 \pi^3 c\_1 E\_{11}}{2L^2},\tag{A9}$$

$$\mathbf{g}\_{32} = \frac{n^2 \pi^2 \left(-A\_{55} + 2c\_2 D\_{55} - c\_2^2 F\_{55}\right)}{2L} - \frac{n^4 \pi^4 c\_1^2 H\_{11}}{2L^3} \tag{A10}$$

$$\mathbf{g}\_{33} = \frac{n^2 \pi^2 \Big( -1 + \left( -1 \right)^n + 4 \Big( 2 + \left( -1 \right)^n \Big) \text{Sim}^4 \left( \frac{n \pi}{2} \right) \Big) A\_{11}}{3L^2},\tag{A11}$$

$$\mathcal{g}\_{34} = \frac{4\left(-2 + (-1)^{n+1}\right) n^3 \pi^3 \text{Sin}^4\left(\frac{n\pi}{\mathcal{Z}}\right) c\_1 E\_{11}}{3L^3},\tag{A12}$$

$$\text{gæs} = -\frac{3n^4\pi^4A\_{11}}{16L^3},\tag{A13}$$

$$\mathbf{g}\_{36} = \frac{n\pi \left(-A\_{55} + 2c\_2 D\_{55} - c\_2^2 F\_{55}\right)}{2} - \frac{n^3 \pi^3 \left(-c\_1 F\_{11} + c\_1^2 H\_{11}\right)}{2L^2},\tag{A14}$$

$$\mathbf{g}\_{37} = \frac{n^2 \pi^2 \left(-1 + (-1)^n + 4\left(2 + (-1)^n\right) \text{Sim}^4\left(\frac{n\pi}{2}\right)\right) \left(B\_{11} - c\_1 E\_{11}\right)}{3L^2},\tag{A15}$$

$$g\_{38} = \frac{n\pi c\_1 m\_3}{2},\tag{A16}$$

$$\mathbf{g}\_{\mathfrak{B}} = -\frac{n^2 \pi^2 c\_1^2 m\_6}{2L} - \frac{Lm\_0}{2},\tag{A17}$$

$$\mathbf{g}\_{310} = \frac{n\pi \left(-c\_1 m\_4 + c\_1^2 m\_6\right)}{2}.\tag{A18}$$

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
