*Article* **Critical Temperatures for Vibrations and Buckling of Magneto-Electro-Elastic Nonlocal Strain Gradient Plates**

**Giovanni Tocci Monaco 1,2, Nicholas Fantuzzi 1,***∗***, Francesco Fabbrocino <sup>3</sup> and Raimondo Luciano <sup>2</sup>**


**Abstract:** An analytical method is presented in this work for the linear vibrations and buckling of nano-plates in a hygro-thermal environment. Nonlinear von Kármán terms are included in the plate kinematics in order to consider the instability phenomena. Strain gradient nonlocal theory is considered for its simplicity and applicability with respect to other nonlocal formulations which require more parameters in their analysis. Present nano-plates have a coupled magneto-electro-elastic constitutive equation in a hygro-thermal environment. Nano-scale effects on the vibrations and buckling behavior of magneto-electro-elastic plates is presented and hygro-thermal load outcomes are considered as well. In addition, critical temperatures for vibrations and buckling problems are analyzed and given for several nano-plate configurations.

**Keywords:** smart nano-plates; semi-analytical solution; critical temperatures; buckling; vibrations

#### **1. Introduction**

In recent years research has focused heavily on MEMS (Micro-Electro-Mechanical-System) and NEMS (Nano-Electro-Mechanical-System). This interest is mainly due to the wide variety of applications in which these devices could be used [1–4]. These structures, such as nanoplates, nanorods, and nanobeams [5], can be used in medicine [6], electronics [7], aerospace [8] and in civil construction [9], where linear and nonlinear theories are generally needed [10]. The behavior of this type of structures cannot be well described through the classical theories of continuous mechanics, as they are based on the principle of location of stresses. Due to the size of these devices, the effects induced by nanoscales must be taken into account [11,12]. Then to improve the ability of new devices and systems made with these smart materials, it is necessary to accurately investigate the mechanical behavior of these advanced structures [13,14]. Non-local theories have been widely used for the study of nanostructures since Eringen developed his theory of non-local elasticity [15]. These theories consider the nano-scale effects thanks to the introduction of one or more length scale parameters in addition to the well known linear elastic Lamé parameters [16–19]. The classification of nonlocal theories is generally presented as: strain gradient [20–23], stress gradient [24], modified strain gradient [25–27], couple stress [28], modified couple stress [29,30], integral type [31,32] and micropolar [33–35]. Article [36] offers a overview on unified continuous/reduced-order modeling and non-linear dynamic theories for thermomechanical plates. Kim in [37] developed a matrix method for evaluating effective elastic constants of generally anisotropic multilayer composites with various coupled physical effects including piezoelectricity, piezomagnetism and thermoelasticity. In [38–40] a nonlocal nonlinear first-order shear theory is used for investigating the buckling and free vibration of magneto-electro-thermo elastic (METE) nanoplates under magneto-electrothermo-mechanical loadings. Mota in [41] investigated the influence of shear factor used

**Citation:** Tocci Monaco, G.; Fantuzzi, N.; Fabbrocino, F.; Luciano, L. Critical Temperatures for Vibrations and Buckling of Magneto-Electro-Elastic Nonlocal Strain Gradient Plates. *Nanomaterials* **2021**, *11*, 274. https://doi.org/ 10.3390/nano11010087

Received: 3 December 2020 Accepted: 29 December 2020 Published: 3 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional clai-ms in published maps and institutio-nal affiliations.

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

in the context of the first-order shear deformation theory on functionally graded porous materials. In [42] free and forced vibration of a functionally graded piezoelectric plate with the properties of the material varying along the thickness are investigated. Combined asymptotic-tolerance modelling of dynamic and stability problems for functionally graded shells was given in [43,44]. In [45] are studied the pyroelectric and pyromagnetic effects on multiphase MEE cylindrical shells subjected to a uniform axisymmetric temperature using semi-analytical finite element procedures. Eremeyev et al. in [46,47] investigate the effect of flexoelectricity and flexomagnetic on nanobeams. Using FEM formulation [48,49] analyzed free vibration of orthotropic cross-ply nanoplates and nanowires. Ebrahimi in [50] studied buckling behaviour of magneto-electro-elastic functionally graded nanobeams using higer-order beam theory and Eringen's non-local elasticity. Also in [51] the behaviour of MEE nanobeams is investigated using Eulero-Bernoulli beam theory and including surface effects. In [52] the focus is on free vibration of laminated circular piezoelectric plates and discs using the weak form of the equations of motion. In [53] using the third order shear deformation plate theory, the bending, buckling, free and forced vibration behavior of a nonlocal composite microplate is analyzed. Functionally graded microplates using Kirchhoff plate's theory and straing gradient theory with only one length scale parameter are studied in [54]. In [55] theory of elasticity including surface stresses is used to study behaviour of shells with nano-scaled thickness. In [56] Kirchhoff plate's theory and the modified flexoelectric theory are used to study the nonlinear free vibration of Functionally Graded (FG) flexoelectric nanoplate by taking into account size-dependent effects. In [57] a finite element model based on a higher-order plate's theory is developed to study static and free vibration problem of magneto-electro-elastic plates. In [58] the dynamic problem of thin elastic plates resting on elastic foundation is studied using tolerance averaging method. Vibrations on periodical structures as well as bad gaps problems in dynamics are still an open topic thoroughly discussed by several researchers [59–62]. In [63] the flexural vibration band gaps in periodic beams is investigated using differential quadrature method, moreover the influence of shear deformation on the gaps is analyzed. Similarly, natural frequencies of structures made of period cells was presented for beams in [64,65] and for plates in [66,67]. In [68] the aim is the problem of vibration of band gaps in periodic Mindlin plates and it is solved using spectral element method. Also in [69] is studied the problem of vibration of band gaps in periodic plates, but is solved using differential quadrature element method. Vibration band gap of stiffened thin plate is studied also in [70] and is solved using center finite difference method. In [71,72] the problem of flexural wave propagation of a periodic beam is investigated.

The aim of this paper is the study of buckling and free vibrations of functionally graded nano-plates in hygro-thermal environment. For buckling analysis it will investigate the influence of external applied electric and magnetic potentials on critical load that leads to the instability of the plate, while in the case of dynamic analysis will be studied the behavior of natural vibration frequencies and how they are influenced by external potentials and temperature. Through the graphs will also identify the critical temperature, which corresponds to the temperature at which the vibration frequencies become zero. This paper is structured as described below. After the introduction section, the theoretical background for functionally graded (FG) thin plates in hygro-thermal environment is developed introducing also the non-linear terms of von Kármán that allow to perform the linear analysis of buckling. Using second order strain gradient theory non local effect are take into account. The following is a small paragraph showing how the electrical and magnetic potentials are approximated. Using Hamilton's principle, for the case of METE (magneto-electro-thermo-elastic) materials, the equations of motion are obtained. The analytical solution is obtained using Navier developments in double trigonometric series. Then the results obtained through calculation code implemented in MATLAB for buckling and free vibration are provided. Finally, a conclusion section is reported at the end of this paper.

#### **2. Theoretical Background**

As show in Figure 1, consider a METE thin nanoplate with length *a*, width *b* and thickness *h*, in a Cartesian reference system (*x*, *y*, *z*).

**Figure 1.** General laminate layout.

The METE nanoplate is in a hygro-thermal environment and is subjected to an electric potential *V*<sup>0</sup> and to a magnetic potential Ω<sup>0</sup> between the upper and lower surfaces. In this study, classic laminate plate theory is considered. We can define the displacement field of a generic point of the solid by means of the triad of displacement components *U*, *V*, *W*, which are functions of the coordinates (*x*, *y*, *z*) [10].

$$\begin{aligned} M(x, y, z, t) &= u(x, y, t) - z \frac{\partial w}{\partial x} \\ V(x, y, z, t) &= v(x, y, t) - z \frac{\partial w}{\partial y} \\ W(x, y, z, t) &= w(x, y, t) \end{aligned} \tag{1}$$

where *u*,*v* and *w* are the displacements along the *x*, *y* and *z* axis of the point on the middle surface and *∂w*/*∂x* and *∂w*/*∂y* are the corresponding rotation. The constitutive equations for a METE material are:

$$\begin{aligned} \sigma &= \mathbf{C}\boldsymbol{\varepsilon} - \mathbf{e}\mathbf{E} - \mathbf{q}\mathbf{H} - \mathbf{C}\boldsymbol{\alpha}\Delta T - \mathbf{C}\boldsymbol{\beta}\Delta H \\ \mathbf{D}\_{E} &= \mathbf{e}^{\top}\boldsymbol{\varepsilon} + \mathbf{\xi}\mathbf{E} + \mathbf{\xi}\mathbf{H} - \mathbf{p}\Delta T - \mathbf{h}\Delta H \\ \mathbf{B}\_{M} &= \mathbf{q}^{\top}\boldsymbol{\varepsilon} + \mathbf{\xi}\mathbf{E} + \mathbf{\chi}\mathbf{H} - \lambda\Delta T - \eta\Delta H \end{aligned} \tag{2}$$

in which *σ* is the classical stress vector, **D***<sup>E</sup>* = [*Dx*, *Dy*, *Dz*] and **B***<sup>M</sup>* = [*Bx*, *By*, *Bz*] are respectively the vector of stresses, electrical displacement and magnetic flux. *ε*, **E** and **H** are the vector of strain, electric field and magnetic field. **C**, *ξ* and *χ* represent the rigidity matrix, the electrical permittivity matrix and the magnetic permittivity matrix. Finally, **e**, **q**, *ζ*, **p**, *λ*, **h** and *η* are respectively the piezo-electric, piezo-magnetic, magneto-electro-elastic (MEE), pyro-electric, pyro-magnetic, hygro-electric and hygro-magnetic coefficients. For the stress plane state (*σ*<sup>3</sup> = 0) the matrices can be reduced by carrying out

$$\varepsilon\_{3} = -\frac{\mathsf{C}\_{13}}{\mathsf{C}\_{33}}\varepsilon\_{1} - \frac{\mathsf{C}\_{23}}{\mathsf{C}\_{33}}\varepsilon\_{2} + \frac{\varepsilon\_{33}}{\mathsf{C}\_{33}}E\_{3} + \frac{q\_{33}}{\mathsf{C}\_{33}}H\_{3} + \frac{\mathsf{C}\_{13}}{\mathsf{C}\_{33}}a\_{1}\Delta T + \frac{\mathsf{C}\_{23}}{\mathsf{C}\_{33}}a\_{2}\Delta T + \frac{\mathsf{C}\_{13}}{\mathsf{C}\_{33}}\beta\_{1}\Delta H + \frac{\mathsf{C}\_{23}}{\mathsf{C}\_{33}}\beta\_{2}\Delta H\tag{3}$$

Therefore the constitutive equations can be rewritten as follows

$$\begin{aligned} \sigma\_{1} &= \left(\mathbb{C}\_{11} - \frac{\mathbb{C}\_{13}^{2}}{\mathbb{C}\_{33}}\right) \varepsilon\_{1} + \left(\mathbb{C}\_{12} - \mathbb{C}\_{13}\frac{\mathbb{C}\_{23}}{\mathbb{C}\_{33}}\right) \varepsilon\_{2} - \left(\mathbb{c}\_{31} - \mathbb{C}\_{13}\frac{\mathbb{c}\_{33}}{\mathbb{C}\_{33}}\right) \mathbb{E}\_{3} - \left(q\_{31} - \mathbb{C}\_{13}\frac{q\_{33}}{\mathbb{C}\_{33}}\right) H\_{3} + \\ &- \left(\mathbb{C}\_{11} - \frac{\mathbb{C}\_{13}^{2}}{\mathbb{C}\_{33}}\right) a\_{1} \Delta T - \left(\mathbb{C}\_{12} - \mathbb{C}\_{13}\frac{\mathbb{C}\_{23}}{\mathbb{C}\_{33}}\right) a\_{2} \Delta T - \left(\mathbb{C}\_{11} - \frac{\mathbb{C}\_{13}^{2}}{\mathbb{C}\_{33}}\right) \beta\_{1} \Delta H - \left(\mathbb{C}\_{12} - \mathbb{C}\_{13}\frac{\mathbb{C}\_{23}}{\mathbb{C}\_{33}}\right) \beta\_{2} \Delta H \\ &= Q\_{11} \varepsilon\_{1} + Q\_{12} \varepsilon\_{2} - \overline{\varepsilon}\_{31} \mathbb{E}\_{3} - \bar{q}\_{31} H\_{3} - Q\_{11} \alpha\_{1} \Delta T - Q\_{12} a\_{2} \Delta T - Q\_{11} \beta\_{1} \Delta H - Q\_{12} \beta\_{2} \Delta H \end{aligned} \tag{4}$$

similarly for *σ*<sup>2</sup> it will be

$$\sigma\_2 = Q\_{12}\varepsilon\_1 + Q\_{22}\varepsilon\_2 - \overline{\varepsilon}\_{32}E\_3 - \overline{q}\_{32}H\_3 - Q\_{12}\mathfrak{a}\_1\Delta T - Q\_{22}\mathfrak{a}\_2\Delta T - Q\_{12}\beta\_1\Delta H - Q\_{22}\beta\_2\Delta H \tag{5}$$

finally, *Dz* can written

$$\begin{aligned} D\_z &= \left( c\_{31} - c\_{33} \frac{\mathcal{C}\_{13}}{\mathcal{C}\_{33}} \right) \varepsilon\_1 + \left( c\_{32} - c\_{33} \frac{\mathcal{C}\_{23}}{\mathcal{C}\_{33}} \right) \varepsilon\_2 + \left( \tilde{\zeta}\_{33} + \frac{\mathcal{C}\_{33}^2}{\mathcal{C}\_{33}} \right) E\_3 + \left( \tilde{\zeta}\_{33} + \mathfrak{c}\_{33} \frac{q\_{33}}{\mathcal{C}\_{33}} \right) H\_3 + \\ &- \left( p\_3 - \frac{\mathcal{C}\_{13}}{\mathcal{C}\_{33}} a\_1 - \frac{\mathcal{C}\_{23}}{\mathcal{C}\_{33}} a\_2 \right) \Delta T - \left( h\_3 - \frac{\mathcal{C}\_{13}}{\mathcal{C}\_{33}} \beta\_1 - \frac{\mathcal{C}\_{23}}{\mathcal{C}\_{33}} \beta\_2 \right) \Delta H \\ &= \tilde{\varepsilon}\_{31} \varepsilon\_1 + \tilde{\varepsilon}\_{32} \varepsilon\_2 + \tilde{\zeta}\_{33} \mathcal{E}\_3 + \tilde{\zeta}\_{33} H\_3 - \tilde{\mu}\_3 \Delta T - \tilde{h}\_3 \Delta H \end{aligned} \tag{6}$$

and similarly for *BM*,3 it will be

$$B\_z = \vec{q}\_{31}\varepsilon\_1 + \vec{q}\_{32}\varepsilon\_2 + \vec{\zeta}\_{33}E\_3 + \vec{\chi}\_{33}H\_3 - \vec{\lambda}\_3\Delta T - \vec{\eta}\_3\Delta H \tag{7}$$

So it is possible to write

$$
\begin{aligned}
\bar{\mathbf{e}} &= \begin{bmatrix} 0 & 0 & \bar{\varepsilon}\_{31} \\ 0 & 0 & \bar{\varepsilon}\_{32} \\ 0 & 0 & 0 \end{bmatrix}, \quad \bar{\mathbf{q}} = \begin{bmatrix} 0 & 0 & \bar{q}\_{31} \\ 0 & 0 & \bar{q}\_{32} \\ 0 & 0 & 0 \end{bmatrix}, \quad \bar{\mathbf{f}} = \begin{bmatrix} \bar{\varepsilon}\_{1} & 0 & 0 \\ 0 & \bar{\varepsilon}\_{2} & 0 \\ 0 & 0 & \frac{\mathcal{I}}{\mathcal{J}\_{3}} \end{bmatrix}, \quad \bar{\mathbf{x}} = \begin{bmatrix} \chi\_{1} & 0 & 0 \\ 0 & \chi\_{2} & 0 \\ 0 & 0 & \bar{\chi}\_{3} \end{bmatrix}, \\\ \bar{\boldsymbol{\xi}} &= \begin{bmatrix} \check{\varepsilon}\_{1} & 0 & 0 \\ 0 & \check{\varepsilon}\_{2} & 0 \\ 0 & 0 & \check{\varepsilon}\_{3} \end{bmatrix}, \quad \mathbf{p} = \begin{Bmatrix} p\_{1} \\ 0 & p\_{2} \\ \bar{p}\_{3} \end{Bmatrix}, \quad \mathbf{h} = \begin{Bmatrix} h\_{1} \\ h\_{2} \\ h\_{3} \end{Bmatrix}, \quad \boldsymbol{\eta} = \begin{Bmatrix} \eta\_{1} \\ \eta\_{2} \\ \bar{\eta}\_{3} \end{Bmatrix}
\end{aligned} \tag{8}$$

By introducing second order strain gradient theory in the constitutive equations and by considering the mechanical properties variable with respect to the thickness direction we have (the dependency on the time *t* is omitted for the sake of simplicity)

$$\sigma(x, y, z) = \left(1 - \ell^2 \nabla^2\right) \left[\mathbf{Q}(z)\boldsymbol{\varepsilon} - \tilde{\mathbf{e}}(z)\mathbf{E} - \tilde{\mathbf{q}}(z)\mathbf{H}\right] - \mathbf{Q}(z)\boldsymbol{a}(z)\Delta T - \mathbf{Q}(z)\boldsymbol{\mathcal{R}}(z)\Delta H$$

$$\mathbf{D}\_E(x, y, z) = \left(1 - \ell^2 \nabla^2\right) \left[\tilde{\mathbf{e}}^\top(z)\boldsymbol{\varepsilon} + \tilde{\mathbf{f}}(z)\mathbf{E} + \tilde{\mathbf{f}}(z)\mathbf{H}\right] - \mathbf{p}(z)\Delta T - \mathbf{h}(z)\Delta H \tag{9}$$

$$\mathbf{B}\_M(x, y, z) = \left(1 - \ell^2 \nabla^2\right) \left[\tilde{\mathbf{q}}^\top(z)\boldsymbol{\varepsilon} + \tilde{\mathbf{f}}(z)\mathbf{E} + \tilde{\mathbf{y}}(z)\mathbf{H}\right] - \lambda(z)\Delta T - \eta(z)\Delta H$$

where is the nonlocal parameter and the operator <sup>∇</sup><sup>2</sup> <sup>=</sup> *<sup>∂</sup>*2/*∂x*<sup>2</sup> <sup>+</sup> *<sup>∂</sup>*2/*∂y*<sup>2</sup> is the second order gradient operator. For the hygro-thermal loads a linear variation is considered along the thickness as

$$
\Delta T = T\_0 + zT\_1/h, \quad \Delta H = H\_0 + zH\_1/h \tag{10}
$$

#### **3. Electric and Magnetic Potentials**

To satisfy Maxwell's equations [73] the electrical and magnetic potential are approximated along the thickness with a linear and cosinusoidal combination. The first amends for the open-circuit condition and the latter for the closed-circuit one

$$\begin{aligned} \Phi(x, y, z, t) &= -\cos\frac{\pi z}{h} \phi(x, y, t) + \frac{2z}{h} V\_0 \\ \Upsilon(x, y, z, t) &= -\cos\frac{\pi z}{h} \gamma(x, y, t) + \frac{2z}{h} \Omega\_0 \end{aligned} \tag{11}$$

in which *V*<sup>0</sup> represents the difference in electrical potential between the two faces of the plate and Ω<sup>0</sup> represents the difference in magnetic potential. The relationships between electric field and electric potential can be written in accordance with the above

$$\begin{aligned} E\_x &= -\frac{\partial \Phi}{\partial x} = \cos\frac{\pi z}{h} \frac{\partial \phi}{\partial x} \\ E\_y &= -\frac{\partial \Phi}{\partial y} = \cos\frac{\pi z}{h} \frac{\partial \phi}{\partial y} \\ E\_z &= -\frac{\partial \Phi}{\partial z} = -\frac{\pi}{h} \sin\frac{\pi z}{h} \phi - \frac{2}{h} V\_0 \end{aligned} \tag{12}$$

that in matrix notation can be rewritten in this form

$$\mathbf{E} = \mathbf{f}\_E \mathbb{D}\_E \boldsymbol{\Phi} + \mathbf{E}\_0 \tag{13}$$

with

$$\mathbf{f}\_E = \begin{bmatrix} \cos\frac{\pi z}{h} & 0 & 0\\ 0 & \cos\frac{\pi z}{h} & 0\\ 0 & 0 & -\frac{\pi}{h}\sin\frac{\pi z}{h} \end{bmatrix}, \quad \mathbb{D}\_E = \begin{Bmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ 1 \end{Bmatrix}, \quad \mathbf{E}\_0 = \begin{Bmatrix} 0\\ 0\\ -\frac{2}{h}V\_0 \end{Bmatrix} \tag{14}$$

Similarly for the magnetic field we can write as

$$\begin{aligned} H\_x &= -\frac{\partial Y}{\partial x} = \cos\frac{\pi z}{h} \frac{\partial \gamma}{\partial x} \\ H\_y &= -\frac{\partial Y}{\partial y} = \cos\frac{\pi z}{h} \frac{\partial \gamma}{\partial y} \\ H\_z &= -\frac{\partial Y}{\partial z} = -\frac{\pi}{h} \sin\frac{\pi z}{h} \gamma - \frac{2}{h} \Omega\_0 \end{aligned} \tag{15}$$

which in matrix notation becomes

$$\mathbf{H} = \mathbf{f}\_H \mathbb{D}\_H \gamma + \mathbf{H}\_0 \tag{16}$$

with

$$\mathbf{f}\_{H} = \begin{bmatrix} \cos\frac{\pi z}{h} & 0 & 0\\ 0 & \cos\frac{\pi z}{h} & 0\\ 0 & 0 & -\frac{\pi}{h}\sin\frac{\pi z}{h} \end{bmatrix}, \quad \mathbb{D}\_{H} = \begin{Bmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ 1 \end{Bmatrix}, \quad \mathbf{H}\_{0} = \begin{Bmatrix} 0\\ 0\\ -\frac{2}{h}\Omega\_{0} \end{Bmatrix} \tag{17}$$

#### **4. Equations of Motion**

The equations of motion are derived through Hamilton's principle

$$\int\_{t\_1}^{t\_2} (\delta H\_{\text{cut}} + \delta V - \delta \mathcal{K}) dt = 0 \tag{18}$$

Writing the variation of enthalpy *δHent*

$$\begin{aligned} \delta H\_{\text{cnt}} &= \int\_{\mathcal{A}} \int\_{-\frac{h}{2}}^{\frac{h}{2}} \left\{ \sigma\_{\text{xx}} \delta \varepsilon\_{\text{xx}} + \sigma\_{yy} \delta \varepsilon\_{yy} + \sigma\_{\text{xy}} \delta \gamma\_{\text{xy}} \right. \\\\ &- D\_{\text{x}} \delta E\_{\text{x}} - D\_{\text{y}} \delta E\_{\text{y}} - D\_{\text{z}} \delta E\_{\text{z}} - B\_{\text{x}} \delta H\_{\text{x}} - B\_{\text{y}} \delta H\_{\text{y}} - B\_{\text{z}} \delta H\_{\text{z}} \right\} d\boldsymbol{z} d\boldsymbol{\mathcal{A}} \end{aligned} \tag{19}$$

by introducing the classical stress resultants *Nxx*, *Nyy*, *Nxy* and *Mxx*, *Myy*, *Mxy* and the piezo and magneto resultants as

$$
\begin{Bmatrix} \mathcal{D}\_x \\ \mathcal{D}\_y \\ \mathcal{D}\_z \end{Bmatrix} = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \mathbf{f}\_E \mathbf{D}\_E \, dz, \quad \begin{Bmatrix} \mathcal{B}\_x \\ \mathcal{B}\_y \\ \mathcal{B}\_z \end{Bmatrix} = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \mathbf{f}\_B \mathbf{B}\_M \, dz \tag{20}
$$

The definition of the integrated elastic properties are given in the Appendix A in Equations (A1)–(A3). Thus, it is obtained

$$\begin{split} \delta H\_{\text{crit}} &= \int\_{\mathcal{A}} \left\{ N\_{\text{xx}} \left( \frac{\partial \delta u}{\partial \mathbf{x}} + \frac{\partial w}{\partial \mathbf{x}} \frac{\partial \delta w}{\partial \mathbf{x}} \right) + N\_{\text{yy}} \left( \frac{\partial \delta v}{\partial \mathbf{y}} + \frac{\partial w}{\partial \mathbf{y}} \frac{\partial \delta w}{\partial \mathbf{y}} \right) \right. \\ &\left. + N\_{\text{xy}} \left( \frac{\partial \delta u}{\partial \mathbf{y}} + \frac{\partial \delta v}{\partial \mathbf{x}} + \frac{\partial \delta w}{\partial \mathbf{y}} \frac{\partial w}{\partial \mathbf{x}} + \frac{\partial \delta w}{\partial \mathbf{x}} \frac{\partial w}{\partial \mathbf{y}} \right) + \\ &\quad + M\_{\text{xx}} \left( -\frac{\partial^{2} \delta w}{\partial \mathbf{x}^{2}} \right) + M\_{\text{yy}} \left( -\frac{\partial^{2} \delta w}{\partial y^{2}} \right) + M\_{\text{xy}} \left( -2 \frac{\partial^{2} \delta w}{\partial \mathbf{x} \partial y} \right) + \\ &\quad - \left( \mathcal{D}\_{\text{x}} \frac{\partial \delta \boldsymbol{\phi}}{\partial \mathbf{x}} + \mathcal{D}\_{\text{y}} \frac{\partial \delta \boldsymbol{\phi}}{\partial y} + \mathcal{D}\_{\text{x}} \delta \boldsymbol{\phi} + \mathcal{B}\_{\text{x}} \frac{\partial \delta \gamma}{\partial \mathbf{x}} + \mathcal{B}\_{\text{y}} \frac{\partial \delta \gamma}{\partial y} + \mathcal{B}\_{\text{z}} \delta \gamma \right) \right) d\mathcal{A} \end{split} \tag{21}$$

Integrating by parts Equation (21) is obtained

*δHent* = - A #*∂Nxx <sup>∂</sup><sup>x</sup>* <sup>+</sup> *∂Nxy ∂y δu* + *∂Nxy <sup>∂</sup><sup>x</sup>* <sup>+</sup> *∂Nyy ∂y δv* + *∂ ∂x Nxx ∂w <sup>∂</sup><sup>x</sup>* <sup>+</sup> *Nxy ∂w ∂y* + *∂ ∂y Nxy ∂w <sup>∂</sup><sup>x</sup>* <sup>+</sup> *Nyy ∂w ∂y* + *∂*2*Mxx <sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *∂*2*Myy <sup>∂</sup>y*<sup>2</sup> <sup>+</sup> <sup>2</sup> *∂*2*Mxy ∂x∂y δw* − *∂*D*<sup>x</sup> <sup>∂</sup><sup>x</sup>* <sup>+</sup> *∂*D*<sup>y</sup> ∂y* + D*<sup>z</sup> δφ* − *∂*B*<sup>x</sup> <sup>∂</sup><sup>x</sup>* <sup>+</sup> *∂*B*<sup>y</sup> ∂y* + B*<sup>z</sup> δγ*\$ *d*A + − - Γ #% *Nxx nx* + *Nxy ny* & *δu* + % *Nxy nx* + *Nyy ny* & *δv* + % *Nxxnx* + *Nxyny* & *<sup>∂</sup><sup>w</sup> ∂x* +% *Nxynx* + *Nyyny* & *<sup>∂</sup><sup>w</sup> ∂y* + *∂Mxx <sup>∂</sup><sup>x</sup>* <sup>+</sup> *∂Mxy ∂y nx* + *∂Myy ∂y* + *∂Mxy ∂x ny δw* −% *Mxxnx* + *Mxyny* & *∂δ<sup>w</sup> <sup>∂</sup><sup>x</sup>* <sup>−</sup> % *Mxynx* + *Myyny* & *∂δ<sup>w</sup> ∂y* −% D*<sup>x</sup> nx* + D*<sup>y</sup> ny* & *δφ* − % B*<sup>x</sup> nx* + B*<sup>y</sup> ny* & *δγ* \$ *d*Γ (22)

The external work due to the external boundary loads (where mechanical, electrical and magnetic loads are neglected) can be written as

$$\begin{split} \delta V &= \int\_{\Gamma} \left\{ \left( \hat{\mathsf{N}}\_{\text{xx}} n\_{\text{x}} + \hat{\mathsf{N}}\_{\text{xy}} n\_{\text{y}} \right) \delta u + \left( \hat{\mathsf{N}}\_{\text{xy}} n\_{\text{x}} + \hat{\mathsf{N}}\_{\text{yy}} n\_{\text{y}} \right) \delta v + \\ & \quad - \left( \hat{\mathsf{M}}\_{\text{xx}} n\_{\text{x}} + \hat{\mathsf{M}}\_{\text{xy}} n\_{\text{y}} \right) \frac{\partial \delta w}{\partial x} - \left( \hat{\mathsf{M}}\_{\text{xy}} n\_{\text{x}} + \hat{\mathsf{M}}\_{\text{yy}} n\_{\text{y}} \right) \frac{\partial \delta w}{\partial y} + \left( \hat{\mathsf{Q}}\_{\text{x}} n\_{\text{x}} + \hat{\mathsf{Q}}\_{\text{y}} n\_{\text{y}} \right) \delta w \right\} d\Gamma \end{split} \tag{23}$$

Variation of the Kinetic energy can be written as

$$\begin{split} \delta K &= \int\_{\mathcal{A}} \left\{ \left( -I\_{0}\ddot{u} + I\_{1}\frac{\partial \vec{w}}{\partial x} \right) \delta u + \left( -I\_{0}\ddot{v} + I\_{1}\frac{\partial \vec{w}}{\partial y} \right) \delta v \\ &+ \left( -I\_{0}\ddot{w} - I\_{1}\frac{\partial \vec{u}}{\partial x} - I\_{1}\frac{\partial \vec{v}}{\partial y} + I\_{2}\frac{\partial \vec{w}}{\partial x} + I\_{2}\frac{\partial \vec{w}}{\partial y} \right) \delta w \right\} d\mathcal{A} + \\ &+ \int\_{\Gamma} \left\{ I\_{1}\ddot{u}n\_{x} + I\_{1}\ddot{v}n\_{y} - I\_{2}\frac{\partial \vec{w}}{\partial x}n\_{x} - I\_{2}\frac{\partial \vec{w}}{\partial y}n\_{y} \right\} \delta w \, d\Gamma \end{split} \tag{24}$$

Introducing N (*w*) and P(*w*) as defined below

$$\begin{array}{rcl}\mathcal{N}(w) &=& \frac{\partial}{\partial x} \Big( N\_{\text{xx}} \frac{\partial w}{\partial x} \Big) \delta w + \frac{\partial}{\partial y} \Big( N\_{\text{yy}} \frac{\partial w}{\partial y} \Big) \delta w + \frac{\partial}{\partial y} \Big( N\_{\text{xy}} \frac{\partial w}{\partial x} \Big) \delta w + \frac{\partial}{\partial x} \Big( N\_{\text{xy}} \frac{\partial w}{\partial y} \Big) \delta w \\\mathcal{P}(w) &=& \Big( N\_{\text{xx}} \frac{\partial w}{\partial x} + N\_{\text{xy}} \frac{\partial w}{\partial y} \Big) n\_x + \Big( N\_{\text{xy}} \frac{\partial w}{\partial x} + N\_{\text{yy}} \frac{\partial w}{\partial y} \Big) n\_y \end{array} \tag{25}$$

the motion equations can be written as follow

$$\begin{aligned} \frac{\partial N\_{xx}}{\partial x} + \frac{\partial N\_{xy}}{\partial y} &= I\_0 \ddot{u} - I\_1 \frac{\partial \ddot{w}}{\partial x} \\ \frac{\partial N\_{yy}}{\partial y} + \frac{\partial N\_{xy}}{\partial x} &= I\_0 \ddot{v} - I\_1 \frac{\partial \ddot{w}}{\partial y} \\ \frac{\partial^2 M\_{xx}}{\partial x^2} + 2 \frac{\partial^2 M\_{xy}}{\partial x \partial y} + \frac{\partial^2 M\_{yy}}{\partial y^2} + \mathcal{N}(w) &= I\_0 \ddot{w} + I\_1 \left(\frac{\partial \ddot{u}}{\partial x} + \frac{\partial \ddot{v}}{\partial y}\right) - I\_2 \left(\frac{\partial^2 \ddot{w}}{\partial x^2} + \frac{\partial^2 \ddot{w}}{\partial y^2}\right) \\ \frac{\partial \mathcal{D}\_x}{\partial x} + \frac{\partial \mathcal{D}\_y}{\partial y} + \mathcal{D}\_z &= 0 \\ \frac{\partial \mathcal{B}\_x}{\partial x} + \frac{\partial \mathcal{B}\_y}{\partial y} + \mathcal{B}\_z &= 0 \end{aligned} \tag{26}$$

and relative boundary conditions become

*δu* = 0 or % *Nxx* <sup>−</sup> *<sup>N</sup>*<sup>ˆ</sup> *xx*& *nx* + % *Nxy* <sup>−</sup> *<sup>N</sup>*<sup>ˆ</sup> *xy*& *ny* = 0 *δv* = 0 or % *Nyy* <sup>−</sup> *<sup>N</sup>*<sup>ˆ</sup> *yy*& *ny* + % *Nxy* <sup>−</sup> *<sup>N</sup>*<sup>ˆ</sup> *xy*& *nx* = 0 *<sup>δ</sup><sup>w</sup>* <sup>=</sup> 0 or *∂Mxx <sup>∂</sup><sup>x</sup>* <sup>+</sup> *∂Mxy <sup>∂</sup><sup>y</sup>* <sup>−</sup> *<sup>I</sup>*1*u*¨ <sup>+</sup> *<sup>I</sup>*<sup>2</sup> *∂w*¨ *∂x nx*+ + *∂Myy ∂y* + *∂Mxy <sup>∂</sup><sup>x</sup>* <sup>−</sup> *<sup>I</sup>*1*v*¨ <sup>+</sup> *<sup>I</sup>*<sup>2</sup> *∂w*¨ *∂y ny* + P(*w*) − % *Q*ˆ *<sup>x</sup>* + *Q*ˆ *<sup>y</sup>* & = 0 *∂δw <sup>∂</sup><sup>x</sup>* <sup>=</sup> 0 or % *Mxx* <sup>−</sup> *<sup>M</sup>*<sup>ˆ</sup> *xx*& *nx* + % *Mxy* <sup>−</sup> *<sup>M</sup>*<sup>ˆ</sup> *xy*& *ny* = 0 *∂δw <sup>∂</sup><sup>y</sup>* <sup>=</sup> 0 or % *Myy* <sup>−</sup> *<sup>M</sup>*<sup>ˆ</sup> *yy*& *ny* + % *Mxy* <sup>−</sup> *<sup>M</sup>*<sup>ˆ</sup> *xy*& *nx* = 0 *δφ* = 0 or D*<sup>x</sup> nx* + D*<sup>y</sup> ny* = 0 *δγ* = 0 or B*<sup>x</sup> nx* + B*<sup>y</sup> ny* = 0 (27)

#### **5. Navier Solution**

Analytical solution is obtained using Navier's expansion. This type of solution allow to solve simply supported plate case. Navier expansion for the displacements take the form

$$
\begin{Bmatrix} \mu \\ v \\ w \end{Bmatrix} = \sum\_{m=1}^{M} \sum\_{n=1}^{N} \begin{bmatrix} \cos ax \sin \beta y & 0 & 0 \\ 0 & \sin ax \cos \beta y & 0 \\ 0 & 0 & \sin ax \sin \beta y \end{bmatrix} \begin{Bmatrix} \mathcal{U}\_{mn} \\ \mathcal{V}\_{mn} \\ \mathcal{W}\_{mn} \end{Bmatrix} \tag{28}
$$

whereas, the electric and magnetic potentials are both approximated with a double sinusoidal trigonometric expansion.

$$\phi = \sum\_{m=1}^{M} \sum\_{n=1}^{N} \sin ax \sin \beta y \,\,\Phi\_{mn}, \quad \gamma = \sum\_{m=1}^{M} \sum\_{n=1}^{N} \sin ax \sin \beta y \,\,\Gamma\_{mn} \tag{29}$$

#### *5.1. Buckling Analysis*

Replacing the displacement field in the motion equations and performing the derivates the algebraic system is obtained

$$
\begin{bmatrix}
\hat{\varepsilon}\_{11} & \hat{\varepsilon}\_{12} & \hat{\varepsilon}\_{14} & \hat{\varepsilon}\_{15} & \hat{\varepsilon}\_{13} \\
\hat{\varepsilon}\_{12} & \hat{\varepsilon}\_{22} & \hat{\varepsilon}\_{24} & \hat{\varepsilon}\_{25} & \hat{\varepsilon}\_{23} \\
\hat{\varepsilon}\_{14} & \hat{\varepsilon}\_{24} & \hat{\varepsilon}\_{44} & \hat{\varepsilon}\_{45} & \hat{\varepsilon}\_{34} \\
\hat{\varepsilon}\_{15} & \hat{\varepsilon}\_{25} & \hat{\varepsilon}\_{45} & \hat{\varepsilon}\_{55} & \hat{\varepsilon}\_{35} \\
\hat{\varepsilon}\_{13} & \hat{\varepsilon}\_{23} & \hat{\varepsilon}\_{34} & \hat{\varepsilon}\_{35} & \hat{\varepsilon}\_{33} + \tilde{s}\_{33}
\end{bmatrix}
\begin{Bmatrix}
\mathcal{U}\_{mn} \\
\mathcal{V}\_{mn} \\
\Phi\_{mn} \\
\Gamma\_{mn} \\
\mathcal{W}\_{mn}
\end{Bmatrix} = \begin{Bmatrix}
0 \\
0 \\
0 \\
0 \\
0 \\
\end{Bmatrix} \tag{30}
$$

The coefficients *c*ˆ*ij* and *s*˜33 are defined in Appendix A at Equation (A4). By introducing the quantities *<sup>N</sup>*<sup>0</sup> <sup>=</sup> <sup>−</sup>*N*<sup>ˆ</sup> *xx*, *<sup>κ</sup>* <sup>=</sup> *<sup>N</sup>*<sup>ˆ</sup> *yy*/*N*<sup>ˆ</sup> *xx* and *amn* as

$$\begin{aligned} a\_{mn} &= \pounds\_{33} + \boldsymbol{\alpha}^2 \left(\boldsymbol{\hat{N}}\_{xx}^T + \boldsymbol{\hat{N}}\_{xx}^E + \boldsymbol{\hat{N}}\_{xx}^H\right) + \boldsymbol{\beta}^2 \left(\boldsymbol{\hat{N}}\_{yy}^T + \boldsymbol{\hat{N}}\_{yy}^E + \boldsymbol{\hat{N}}\_{yy}^H\right) \\ &- \left\{\boldsymbol{\varepsilon}\_{13} \quad \boldsymbol{\varepsilon}\_{23} \quad \boldsymbol{\varepsilon}\_{34} \quad \boldsymbol{\varepsilon}\_{35}\right\} \begin{bmatrix} \boldsymbol{\varepsilon}\_{11} & \boldsymbol{\varepsilon}\_{12} & \boldsymbol{\varepsilon}\_{14} & \boldsymbol{\varepsilon}\_{15} \\ \boldsymbol{\varepsilon}\_{12} & \boldsymbol{\varepsilon}\_{22} & \boldsymbol{\varepsilon}\_{24} & \boldsymbol{\varepsilon}\_{25} \\ \boldsymbol{\varepsilon}\_{14} & \boldsymbol{\varepsilon}\_{24} & \boldsymbol{\varepsilon}\_{44} & \boldsymbol{\varepsilon}\_{45} \\ \boldsymbol{\varepsilon}\_{15} & \boldsymbol{\varepsilon}\_{25} & \boldsymbol{\varepsilon}\_{45} & \boldsymbol{\varepsilon}\_{55} \end{bmatrix}^{-1} \begin{Bmatrix} \boldsymbol{\varepsilon}\_{13} \\ \boldsymbol{\varepsilon}\_{23} \\ \boldsymbol{\varepsilon}\_{34} \\ \boldsymbol{\varepsilon}\_{35} \end{Bmatrix} \tag{31}$$

we can write the solution of the eigenvalue problem as

$$\text{N}\_{\text{0}} = \frac{a\_{mn}}{\mathfrak{a}^2 + \kappa \mathfrak{F}^2} \tag{32}$$

The load that buckles the plate depends on *m* and *n* and in particular the critical load is the lowest of the buckling loads. The terms *N*ˆ *<sup>E</sup> xx* , *N*ˆ *<sup>E</sup> yy* , *N*ˆ *<sup>H</sup> xx* , *N*ˆ *<sup>H</sup> yy* are defined below

$$
\hat{N}\_{\text{xx}}^{E} = \hat{N}\_{yy}^{E} = \int\_{-h/2}^{h/2} \mathbb{I}\_{31}(z) \frac{2V\_0}{h} \, dz \quad , \quad \hat{N}\_{\text{xx}}^{H} = \hat{N}\_{yy}^{H} = \int\_{-h/2}^{h/2} \vec{q}\_{31}(z) \frac{2\Omega\_0}{h} \, dz \tag{33}
$$

Note that the electric and magnetic in-plane loads have the same intensity since in the following applications the material is isotropic in-plane and anisotropic out-of-plane.

#### *5.2. Thermal Free Vibration*

In this paragraph it will be treated the problem of free vibrations of the FG plate simply supported. For this problem it is necessary to rewrite the solving system in a homogeneous form, and the rotational inertia *I*<sup>1</sup> are neglected. The system becomes then

$$
\begin{pmatrix}
\begin{bmatrix}
\hat{\varepsilon}\_{11} & \hat{\varepsilon}\_{12} & \hat{\varepsilon}\_{13} & \hat{\varepsilon}\_{14} & \hat{\varepsilon}\_{15} \\
\hat{\varepsilon}\_{12} & \hat{\varepsilon}\_{22} & \hat{\varepsilon}\_{23} & \hat{\varepsilon}\_{24} & \hat{\varepsilon}\_{25} \\
\hat{\varepsilon}\_{13} & \hat{\varepsilon}\_{23} & \hat{\varepsilon}\_{33} & \hat{\varepsilon}\_{34} & \hat{\varepsilon}\_{35} \\
\hat{\varepsilon}\_{14} & \hat{\varepsilon}\_{24} & \hat{\varepsilon}\_{34} & \hat{\varepsilon}\_{44} & \hat{\varepsilon}\_{45} \\
\hat{\varepsilon}\_{15} & \hat{\varepsilon}\_{25} & \hat{\varepsilon}\_{35} & \hat{\varepsilon}\_{45} & \hat{\varepsilon}\_{55}
\end{bmatrix} - \omega^{2} \\
\begin{bmatrix}
\hat{m}\_{11} & 0 & 0 & 0 & 0 \\
0 & \hat{m}\_{22} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0
\end{bmatrix}
\end{pmatrix}
\begin{cases}
\begin{bmatrix}
U\_{mu} \\
V\_{mn} \\
W\_{mn} \\
\Phi\_{mn} \\
\Phi\_{mn} \\
\Gamma\_{mn}
\end{bmatrix} = \begin{Bmatrix}
0 \\
0 \\
0 \\
0 \\
0
\end{Bmatrix}
\end{aligned}
\tag{34}$$

with *m*ˆ <sup>11</sup> = *m*ˆ <sup>22</sup> = *I*<sup>0</sup> e *m*ˆ <sup>33</sup> = *I*<sup>0</sup> + *I*2(*α*<sup>2</sup> + *β*2). This system can then be rewritten in a more compact matrix form as follows

$$
\omega \begin{pmatrix} \begin{bmatrix} \mathbf{K}\_{\mu u} & \mathbf{K}\_{u\phi} \\ \mathbf{K}\_{\phi u} & \mathbf{K}\_{\phi\phi} \end{bmatrix} - \omega^2 \begin{bmatrix} \mathbf{M}\_{uu} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \end{pmatrix} \begin{Bmatrix} \mathbf{U} \\ \mathbf{0} \end{Bmatrix} = \begin{Bmatrix} \mathbf{0} \\ \mathbf{0} \end{Bmatrix} \tag{35}$$

Rewriting the matrix system by applying static condensation we get

$$\left(\vec{\mathcal{K}} - \omega^2 \mathbf{M}\_{\mu u}\right) \mathcal{U} = \begin{array}{c} \mathbf{0} \\ \end{array} \tag{36}$$

where *K***¯** is

$$\mathcal{K} = (\mathbf{K}\_{\mu u} - \mathbf{K}\_{\nu \phi} (\mathbf{K}\_{\phi \phi})^{-1} \mathbf{K}\_{\phi u}) \tag{37}$$

and *U* represents the ways of vibrating and *ω* natural frequencies.

#### **6. Numerical Results**

In this paper for the numerical solution has been considered a FG nanoplate composed of CoFe2O4 and BaTiO3, the properties of the materials are shown in Table 1. Since it has not been possible to find in literature the hygrometric coefficients of the materials, the applications presented foresee only the case of thermal environment.

**Table 1.** Piezo-electro-magnetic-thermal properties of BaTiO3 and CoFe2O4.


The variation of the material properties along the thickness is regulated by the following relationship

$$P(z) = (P\_l - P\_b) \left(\frac{z}{h} + \frac{1}{2}\right)^{p\_n} + P\_b \tag{38}$$

where *Pt* and *Pb* represent the properties of the material placed on the top and bottom of the plate, respectively.

Note that if *pn* = 0 the plate will be composed entirely of the material of the top side while if *pn* → ∞ the plate will be composed entirely of the material of the bottom side.

#### *6.1. Buckling*

In the following applications the values of the critical load will be presented in dimensionless form through the following relationship

$$
\bar{N}\_{\rm cr} = \frac{N\_{\rm cr} a^2}{C\_{11, \rm m} h^3} \tag{39}
$$

where *C*11,*<sup>m</sup>* is the average stiffness value of the two materials. As a first result the comparison with Li [74], Park and Han [75] is reported. The plate considered is isotropic with *a*/*h* = 1000 and the properties of the material used are as reported in the cited works (which are obtained as average values of two isotropic constituents, not as a functionally graded composite).

Figure 2 shows the critical buckling load by varying the electric and magnetic potentials. The present results agree well with the ones presented in the mentioned papers [74,75]. It is emphasized that the slight difference in the results is due to the fact that in the cited studies the potential is approximated using three contributions: one parabolic, one linear and one constant unlike the present study in which the potentials are approximated with a cosine function and a linear part. In addition, in the cited studies, the in-plane components *Ex*, *Ey*, *Hx* and *Hy* of electric and magnetic fields are null.

**Figure 2.** Critical load *N*¯ *cr* of a square Functionally Graded (FG) nanoplate for different values of electric potential *V*<sup>0</sup> (**a**) and magnetic potential Ω<sup>0</sup> (**b**).

In the applications below *a*/*h* = 1000 and *np* = 1 are always considered. Table 2 shows the dimensionless critical load of a square plate FG for different values of externally applied potentials and non-local parameter. It can be seen that as the non-local parameter increases, the value of the critical load increases. It can also be seen that by increasing the magnetic potential the critical load increases while the electric potential has the opposite effect. This last phenomenon is also clearly visible in Figure 3, where the value of the critical load is reported as the external potentials applied vary and for different values of the non-local parameter. Finally it is remarked that for same values of the electric and magnetic potentials the critical load takes a negative value, thus, buckling occurs for traction loads instead of compression. Figure 4 shows the dimensional critical load when the aspect ratio varies and for different values of the non-local parameter. As expected the critical load increases as the plane becomes of rectangular shape and as the nonlocal parameter increases.


**Table 2.** Dimensionless critical load *N*¯ *cr* of a square FG plate composed of BaTiO3/CoFe2O4 for different electric and magnetic potentials and nonlocal parameter (-/*a*)2. (*κ* = 1,(*m*, *n*)=(1, 1)).

**Figure 3.** Critical load *N*¯ *cr* of a square FG nanoplate for different values of electric potential *V*<sup>0</sup> (**a**) and magnetic potential Ω<sup>0</sup> (**b**), for different values of non local parameter (-/*a*)2.

**Figure 4.** Critical load *N*¯ *cr* of square FG nanoplate for different values of ratio *a*/*b* and different values of non local parameter (-/*a*)2.

#### *6.2. Thermal Free Vibration*

As first application for the free vibration of piezo-electro-magnetic-thermal plates a comparison with [76] is reported (Table 3). The plate is composed of BaTiO3/CoFe2O4 and is of rectangular shape with *a* = 2 and *b* = 1, ratio *h*/*a* = 0.1 and properties vary linearly along the thickness (*np* = 1). The results are written in dimensionless form through the following formula

$$
\omega = \omega \left(\frac{a^2}{h}\right) \sqrt{\frac{\rho}{\mathcal{C}\_{11}}} \tag{40}
$$

where *ρ* and *C*<sup>11</sup> respectively represent the density and the (1, 1) position element of the material stiffness matrix on the underside of the plate. It should be noted that in the study just mentioned the Mindlin's moderately thick plates theory (FSDT) is used, and so the results deviate slightly and this difference increases as the vibration mode increases as the effects of shear become more relevant.

**Table 3.** Dimensionless natural frequencies *ω*¯ of a simply supported rectangular FG plate composed of BaTiO3/CoFe2O4.


Table 4 shows a comparison with article [57] for a thick magneto-electro-elastic square plate. The thickness of the plate is constant and the laminae all have the same thickness. The properties of the material are those reported in Table 1, except the density which is assumed constant for the two materials and equal to 1600 kg/m3. The values calculated in this study differ slightly from those in the literature because in [57] a third order plate theory is used while in this study it is used thin plate theory.

**Table 4.** Natural frequencies (rad/s) of a simply supported square plate (*a* = 1 m; *h* = 0.3 m). B = BaTiO3; F = CoFe2O4.


Table 5 analyzes the influence of the non-local parameter on the natural frequencies of a square plate with *a* = *b* = 1 m and ratio *a*/*h* = 100 and *np* = 1. The results show an increase in natural frequencies, due to the stiffening of the plate, as the non-local parameter increases.

Figure 5 shows the influence of temperature and externally applied potentials on the natural frequency of a nanoplate composed by BaTiO3/CoFe2O4. The plate considered is a simply supported square plate and the thickness ratio is *a*/*h* = 100. In particular, the critical temperature of each structure can be identified when the frequency becomes zero.

**Table 5.** Natural frequencies *ω*¯ of a simply supported square FG nanoplate composed by BaTiO3/CoFe2O4.


**Figure 5.** Natural frequencies *ω*¯ of a simply supported square FG nanoplate composed by BaTiO3/CoFe2O4 to vary of temperature *T*<sup>0</sup> and for different values of magnetic and electric potentials and non-local parameter: (**a**) local configuration (- = 0); (**b**) nonlocal configuration with *V*<sup>0</sup> = 0; (**c**) nonlocal configuration with Ω<sup>0</sup> = 0; (**d**) nonlocal configuration with *V*<sup>0</sup> = Ω<sup>0</sup> = 0.

#### **7. Conclusions**

In this paper the dynamic and buckling problems of METE nanoplates have been analyzed. In particular, the interest focused on the coupling of magnet-electro-thermoelastic effects and the influence that the external potentials applied to the plate have on the critical load and natural vibration frequencies. Through Hamilton's principle motion equations for FG METE thin plates are derived and analytical solution using Navier method

is obtained. The materials that have been used in the simulations are BaTiO3 and CoFe2O4 and the properties of the materials used are reported in the article. The results show that increasing the non-local parameter increases the critical load and natural vibration frequencies. For external potentials instead it was seen that the critical load increases with the increase of the negative electrical potential and the positive magnetic potential. Finally, from the graphs of the natural frequency of vibration it can be seen that the frequencies tend to increase by subjecting the plate to a positive magnetic potential and to decrease by subjecting it to a positive electrical potential. For what concerns the temperature instead we see how an increase of the latter leads to a reduction of the natural frequencies.

**Author Contributions:** Conceptualization, N.F.; methodology, N.F.; software, N.F. and G.T.M.; validation, N.F. and G.T.M.; formal analysis, N.F. and G.T.M.; investigation, N.F. and G.T.M.; resources, N.F., F.F. and R.L.; data curation, G.T.M.; writing—original draft preparation, N.F. and G.T.M.; writing—review and editing, N.F., G.T.M., F.F. and R.L.; visualization, F.F. and R.L.; supervision, N.F., F.F. and R.L.; project administration, N.F., F.F. and R.L.; funding acquisition, F.F. and R.L. All authors have read and agreed to the published version of the manuscript.

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

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing not applicable.

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

#### **Appendix A**

Integrating last two integrals of Equation (22) along the thickness the following quantities can be defined

$$\begin{aligned} \mathbf{A}\_{E}^{f} &= \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \tilde{\mathbf{e}} \mathbf{f}\_{E} \, dz, \quad \mathbf{A}\_{E} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \tilde{\mathbf{e}} \, dz, \quad \mathbf{B}\_{E}^{f} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \tilde{\mathbf{e}} \mathbf{f}\_{E} \, z \, dz \quad , \quad \mathbf{B}\_{E} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \tilde{\mathbf{e}} \, z \, dz \\\ \mathbf{A}\_{H}^{f} &= \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \tilde{\mathbf{q}} \mathbf{f}\_{H} \, dz, \quad \mathbf{A}\_{H} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \tilde{\mathbf{q}} \, dz, \quad \mathbf{B}\_{H}^{f} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \tilde{\mathbf{q}} \mathbf{f}\_{H} \, z \, dz, \quad \mathbf{B}\_{H} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \tilde{\mathbf{q}} \, z \, dz \end{aligned} \tag{A1}$$

− *h* 2 − *h* 2 − *h* 2 − *h* 2 **A***f <sup>ξ</sup>* = *h* 2 − *h* 2 **f** *E* **˜** *ξdz*, **B***<sup>f</sup> <sup>ξ</sup>* = *h* 2 − *h* 2 **f** *E* **˜** *<sup>ξ</sup>***f***Edz*, **<sup>A</sup>***<sup>f</sup> <sup>ζ</sup>* = *h* 2 − *h* 2 **f** *E* **˜** *ζdz*, **B***<sup>f</sup> <sup>ζ</sup>* = *h* 2 − *h* 2 **f** *E* **˜** *ζ***f***Hdz* **A***<sup>p</sup>* = *h* 2 − *h* 2 **f** *<sup>E</sup>* **p***dz*, **B***<sup>p</sup>* = *h* 2 − *h* 2 **f** *<sup>E</sup>* **p***z dz*, **A***<sup>h</sup>* = *h* 2 − *h* 2 **f** *<sup>E</sup>* **h***dz*, **B***<sup>h</sup>* = *h* 2 − *h* 2 **f** *<sup>E</sup>* **h***z dz* (A2) **A***f <sup>ζ</sup>* = *h* 2 − *h* 2 **f** *H* **˜** *ζdz*, **B***<sup>f</sup> <sup>ζ</sup>* = *h* 2 − *h* 2 **f** *H* **˜** *<sup>ζ</sup>***f***Edz*, **<sup>A</sup>***<sup>f</sup> <sup>χ</sup>* = *h* 2 − *h* 2 **f** *<sup>H</sup>χ***˜***dz*, **<sup>B</sup>***<sup>f</sup> <sup>χ</sup>* = *h* 2 − *h* 2 **f** *<sup>H</sup>χ***˜f***Hdz h h h h* (A3)

$$\mathbf{A}\_{\lambda} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \mathbf{f}\_{H}^{\top} \lambda dz, \quad \mathbf{B}\_{\lambda} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \mathbf{f}\_{H}^{\top} \lambda z \, dz, \quad \mathbf{A}\_{\eta} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \mathbf{f}\_{H}^{\top} \eta dz, \quad \mathbf{B}\_{\eta} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \mathbf{f}\_{H}^{\top} \eta z \, dz$$

Considering what reported in Equations (A1)–(A3) it can be write *c*ˆ*ij* coefficients as follows

*c*ˆ11 = *α*2*A*<sup>11</sup> + *β*2*A*<sup>66</sup> + -2 *<sup>α</sup>*4*A*<sup>11</sup> <sup>+</sup> *<sup>α</sup>*2*β*2(*A*<sup>11</sup> <sup>+</sup> *<sup>A</sup>*66) <sup>+</sup> *<sup>β</sup>*4*A*66 *c*ˆ12 = *αβ*(*A*<sup>12</sup> + *A*66) + -2 *α*3*β*(*A*<sup>12</sup> + *A*66) + *αβ*3(*A*<sup>12</sup> + *A*66) *<sup>c</sup>*ˆ13 <sup>=</sup> <sup>−</sup>*α*3*B*<sup>11</sup> <sup>−</sup> *αβ*2(*B*<sup>12</sup> <sup>+</sup> <sup>2</sup>*B*66) − -2 *α*5*B*<sup>11</sup> + *α*3*β*2(*B*<sup>12</sup> + *B*<sup>11</sup> + 2*B*66) + *αβ*4(*B*<sup>12</sup> + 2*B*66) *c*ˆ22 = *β*2*A*<sup>22</sup> + *α*2*A*<sup>66</sup> + -2 *<sup>α</sup>*2*β*2(*A*<sup>22</sup> <sup>+</sup> *<sup>A</sup>*66) <sup>+</sup> *<sup>α</sup>*4*A*<sup>66</sup> <sup>+</sup> *<sup>β</sup>*4*A*22 *<sup>c</sup>*ˆ23 <sup>=</sup> <sup>−</sup>*α*2*β*(*B*<sup>12</sup> <sup>+</sup> <sup>2</sup>*B*66) <sup>−</sup> *<sup>β</sup>*3*B*<sup>22</sup> − -2 *<sup>α</sup>*4*β*(*B*<sup>12</sup> <sup>+</sup> <sup>2</sup>*B*66) <sup>+</sup> *<sup>α</sup>*2*β*3(*B*<sup>22</sup> <sup>+</sup> *<sup>B</sup>*<sup>12</sup> <sup>+</sup> <sup>2</sup>*B*66) <sup>+</sup> *<sup>β</sup>*5*B*22 *c*ˆ33 = *α*4*D*<sup>11</sup> + 2*α*2*β*2(*D*<sup>12</sup> + 2*D*66) + *β*4*D*<sup>22</sup> + -2 *α*6*D*<sup>11</sup> + *α*4*β*2(*D*<sup>11</sup> + 2*D*<sup>12</sup> + 4*D*66) + *α*2*β*4(*D*<sup>22</sup> + 2*D*<sup>12</sup> + 4*D*66) + *β*6*D*<sup>22</sup> *<sup>c</sup>*ˆ14 <sup>=</sup> <sup>−</sup>*αA<sup>f</sup> <sup>E</sup>*,13 + -2 −*α*3*A<sup>f</sup> <sup>E</sup>*,13 <sup>−</sup> *αβ*2*A<sup>f</sup> E*,13 *<sup>c</sup>*ˆ15 <sup>=</sup> <sup>−</sup>*αA<sup>f</sup> <sup>H</sup>*,13 + -2 −*α*3*A<sup>f</sup> <sup>H</sup>*,13 <sup>−</sup> *αβ*2*A<sup>f</sup> H*,13 *<sup>c</sup>*ˆ24 <sup>=</sup> <sup>−</sup>*βA<sup>f</sup> <sup>E</sup>*,23 + -2 −*β*3*A<sup>f</sup> <sup>E</sup>*,23 <sup>−</sup> *<sup>α</sup>*2*βA<sup>f</sup> E*,23 *<sup>c</sup>*ˆ25 <sup>=</sup> <sup>−</sup>*βA<sup>f</sup> <sup>H</sup>*,23 + -2 −*β*3*A<sup>f</sup> <sup>H</sup>*,23 <sup>−</sup> *<sup>α</sup>*2*βA<sup>f</sup> H*,23 *<sup>c</sup>*ˆ34 <sup>=</sup> *<sup>α</sup>*2*B<sup>f</sup> <sup>E</sup>*,13 <sup>+</sup> *<sup>β</sup>*2*B<sup>f</sup> <sup>E</sup>*,23 + -2 *α*4*B<sup>f</sup> <sup>E</sup>*,13 <sup>+</sup> *<sup>α</sup>*2*β*<sup>2</sup> *Bf <sup>E</sup>*,13 <sup>+</sup> *<sup>B</sup><sup>f</sup> E*,23 + *β*4*B<sup>f</sup> E*,23 *<sup>c</sup>*ˆ35 = *<sup>α</sup>*2*B<sup>f</sup> <sup>H</sup>*,13 <sup>+</sup> *<sup>β</sup>*2*B<sup>f</sup> <sup>H</sup>*,23 + -2 *α*4*B<sup>f</sup> <sup>H</sup>*,13 <sup>+</sup> *<sup>α</sup>*2*β*<sup>2</sup> *Bf <sup>H</sup>*,13 <sup>+</sup> *<sup>B</sup><sup>f</sup> H*,23 + *β*4*B<sup>f</sup> H*,23 *<sup>c</sup>*ˆ44 <sup>=</sup> <sup>−</sup>*B<sup>f</sup> <sup>ξ</sup>*,33 <sup>−</sup> *<sup>α</sup>*2*B<sup>f</sup> <sup>ξ</sup>*,11 <sup>−</sup> *<sup>β</sup>*2*B<sup>f</sup> ξ*,22 − -2 *α*4*B<sup>f</sup> <sup>ξ</sup>*,11 <sup>+</sup> *<sup>α</sup>*2*β*<sup>2</sup> *Bf <sup>ξ</sup>*,11 <sup>+</sup> *<sup>B</sup><sup>f</sup> ξ*,22 + *β*4*B<sup>f</sup> <sup>ξ</sup>*,22 <sup>+</sup> *<sup>α</sup>*2*B<sup>f</sup> <sup>ξ</sup>*,33 <sup>+</sup> *<sup>β</sup>*2*B<sup>f</sup> ξ*,33 *<sup>c</sup>*ˆ45 <sup>=</sup> <sup>−</sup>*Bf E <sup>ζ</sup>*,33 <sup>−</sup> *<sup>α</sup>*2*Bf E <sup>ζ</sup>*,11 <sup>−</sup> *<sup>β</sup>*2*Bf E ζ*,22 − -2 *α*4*Bf E <sup>ζ</sup>*,11 <sup>+</sup> *<sup>α</sup>*2*β*<sup>2</sup> *Bf E <sup>ζ</sup>*,11 <sup>+</sup> *<sup>B</sup>f E ζ*22 + *β*4*Bf E <sup>ζ</sup>*,22 <sup>+</sup> *<sup>α</sup>*2*Bf E <sup>ζ</sup>*,33 <sup>+</sup> *<sup>β</sup>*2*Bf E ζ*,33 *<sup>c</sup>*ˆ54 <sup>=</sup> <sup>−</sup>*Bf H <sup>ζ</sup>*,33 <sup>−</sup> *<sup>α</sup>*2*Bf H <sup>ζ</sup>*,11 <sup>−</sup> *<sup>β</sup>*2*Bf H ζ*,22 − -2 *α*4*Bf H <sup>ζ</sup>*,11 <sup>+</sup> *<sup>α</sup>*2*β*<sup>2</sup> *Bf H <sup>ζ</sup>*,11 <sup>+</sup> *<sup>B</sup>f H ζ*,22 + *β*4*Bf H <sup>ζ</sup>*,22 <sup>+</sup> *<sup>α</sup>*2*Bf H <sup>ζ</sup>*,33 <sup>+</sup> *<sup>β</sup>*2*Bf H ζ*,33 *<sup>c</sup>*ˆ55 <sup>=</sup> <sup>−</sup>*B<sup>f</sup> <sup>χ</sup>*,33 <sup>−</sup> *<sup>α</sup>*2*B<sup>f</sup> <sup>χ</sup>*,11 <sup>−</sup> *<sup>β</sup>*2*B<sup>f</sup> χ*22 − -2 *α*4*B<sup>f</sup> <sup>χ</sup>*,11 <sup>+</sup> *<sup>α</sup>*2*β*<sup>2</sup> *Bf <sup>χ</sup>*,11 <sup>+</sup> *<sup>B</sup><sup>f</sup> χ*,22 + *β*4*B<sup>f</sup> <sup>χ</sup>*,22 <sup>+</sup> *<sup>α</sup>*2*B<sup>f</sup> <sup>χ</sup>*,33 <sup>+</sup> *<sup>β</sup>*2*B<sup>f</sup> χ*,33 *s*˜33 = *α N*ˆ *xx* + *N*ˆ *<sup>T</sup> xx* + *N*ˆ *<sup>E</sup> xx* + *N*ˆ *<sup>H</sup> xx* + *β N*ˆ *yy* + *N*ˆ *<sup>T</sup> yy* + *N*ˆ *<sup>E</sup> yy* + *N*ˆ *<sup>H</sup> yy* (A4)

In case the electrical and magnetic potentials have the same through-the-thickness expansion, as in the present study, **B***f E <sup>ζ</sup>* <sup>=</sup> **<sup>B</sup>***f H <sup>ζ</sup>* <sup>=</sup> **<sup>B</sup>***<sup>f</sup> <sup>ζ</sup>* , hence, *c*ˆ45 = *c*ˆ54.

#### **References**

