*Article* **Vibration and Buckling of Shear Deformable Functionally Graded Nanoporous Metal Foam Nanoshells**

#### **Yufei Zhang 1,\* and Fei Zhang 2,\***

<sup>1</sup> College of Aerospace Engineering, Shenyang Aerospace University, Shenyang 110136, China

<sup>2</sup> College of Sciences, Northeastern University, Shenyang 110819, China

**\*** Correspondence: yufeizhang73@163.com (Y.Z.); feizhangneu@163.com (F.Z.)

Received: 22 January 2019; Accepted: 11 February 2019; Published: 15 February 2019

**Abstract:** This article aims to investigate free vibration and buckling of functionally graded (FG) nanoporous metal foam (NPMF) nanoshells. The first-order shear deformation (FSD) shell theory is adopted and the theoretical model is formulated by using Mindlin's most general strain gradient theory, which can derive several well-known simplified models. The symmetric and unsymmetric nanoporosity distributions are considered for the structural composition. Hamilton's principle is employed to deduce the governing equations as well as the boundary conditions. Then, via the Navier solution technique, an analytical solution for the free vibration and buckling of FG NPMF nanoshells is presented. Afterwards, a detailed parametric analysis is conducted to highlight the effects of the nanoporosity coefficient, nanoporosity distribution, length scale parameter, and geometrical parameters on the mechanical behaviors of FG NPMF nanoshells.

**Keywords:** nanoporous metal foam; nanoshell; buckling; free vibration; strain gradient theory; first-order shear deformation theory

#### **1. Introduction**

Functionally graded materials (FGMs) have a continuous and smooth graded distribution of material properties in the spatial field. Due to their superior properties and advantages, FGMs have been successfully extended to various engineering applications and received much attention [1–24]. Recently, a breakthrough made it possible to realize desired structural properties by adjusting the local density of structures, thereby developing novel functionally graded (FG) porous structures composed of metal foams having graded density [25–28]. The application of nanoporous metal foams (NPMFs) has been extended to some advanced engineering fields due to their extremely high specific surface area [29–32]. This kind of material has a combination of properties that is not achievable for ceramics, metals, or dense polymers.

Micro/nanostructures have been successfully used in shape memory alloys [33] and micro- and nano-electro-mechanical systems (MEMS and NEMS) [34,35]. The small-scale effects on the mechanical behaviors of micro/nanostructures have been experimentally observed in their applications [36,37]. It was revealed that the mechanical behaviors of micro/nanostructures were different from their macro counterparts due to the size effect [38,39]. Due to the lack of intrinsic material length scale parameters, the classical continuum theory has no ability to predict the mechanical characteristics of micro/nanostructures. Therefore, several size-dependent continuum theories have been proposed to compensate for the drawbacks of the classical continuum theory for micro/nanostructures. One of the size-dependent continuum theories is Mindlin's strain gradient theory (SGT) [40], which is known as the general form of the SGT containing five additional material length scale parameters compared to the classical continuum theory. Later, several special forms of Mindlin's SGT were proposed. For instance,

one of the most popular forms is the modified strain gradient theory (MSGT) [37]. In fact, this theory is a more useful form of Mindlin's SGT including three material length scale parameters related to symmetric rotation gradients, deviatoric stretch gradients, and dilatation gradients. Several successful applications of the MSGT in dynamic and static analyses of micro/nanobeams [41–44], plates [45–47], and shells [48–50] have been reported. It should be noticed that the modified couple stress theory (MCST) [51] can be achieved by ignoring two of the three material length scale parameters in the MSGT. Moreover, the MSGT can be simplified to the classical theory (CT) by neglecting all of the three material length scale parameters.

Recently, the structural performance of NPMF micro/nanobeams has been investigated by several researchers. Post-buckling analysis for nanobeams made of NPMFs is presented by Barati and Zenkour [52] via the nonlocal elasticity theory (NET). By using the nonlocal strain gradient theory together with the third-order shear deformation beam theory, nonlinear bending of FG NPMF micro/nanobeams reinforced by graphene platelets has been analyzed by Sahmani et al. [53]. Wang et al. [54] utilized the sinusoidal beam theory and the MSGT to study the vibration and bending of NPMF microbeams.

Shell-type structures have excellent mechanical properties [55–61], and thus, nanoshells are important components in various MEMS and NEMS [62–64]. Complete knowledge of the mechanical properties of nanoshells encourages researchers to use them more efficiently. Therefore, some research has been made to illustrate the buckling and vibration characteristics of nanoshells. For example, by using the NET, Hoseinzadeh and Khadem [65] investigated the thermoelastic vibration of double-walled carbon nanotubes (CNTs). By employing the classical shell theory together with the Gurtin-Murdoch elasticity theory, Sahmani et al. [66] analyzed the postbuckling and nonlinear buckling of cylindrical nanoshells subjected to radial and axial compressive loads. Implementing the NET in the first-order shear deformation (FSD) shell theory, Ansari et al. [67] explored the buckling behavior of multi-walled CNTs including the effect of the thermal environment. Wang et al. [68] studied the nonlinear vibration of nanoshells conveying fluid based on the surface stress elasticity theory as well as the classical shell theory.

In the present study, we aim to make an attempt to investigate the vibration and buckling of circular cylindrical nanoshells made from FG NPMFs. In order to accommodate the size dependency of the nanostructure, the general SGT is used to develop the size-dependent first-order shear deformable nanoporous nanoshell model. The governing equations, as well as the related boundary conditions, are obtained simultaneously by utilizing Hamilton's principle. The free vibration and axial buckling of simply supported nanoporous circular cylindrical nanoshells are solved analytically by means of the Navier solution technique. Moreover, the influence of some key parameters on the vibration and buckling properties of the system is shown.

#### **2. FG NPMF Circular Cylindrical Nanoshells**

An FG NPMF circular cylindrical nanoshell of middle-surface radius *R*, thickness *h*, and length *L* is shown in Figure 1. Two kinds of nanoporosity distribution in the thickness direction are considered, namely, nanoporosity-1 and nanoporosity-2. Additionally, the nanoshell is subjected to axial loads *N*<sup>0</sup> *xx*.

Owing to non-uniform nanoporosity distribution, mass densities *ρ*(*z*), Young's modulus *E*(*z*), and shear modulus *μ*(*z*) of the nanoshell are functions of position and can be written as [69–74]:

Nanoporosity-1:

$$E(z) = E\_1^\* \left[ 1 - \varepsilon\_0 \cos(\pi \mathcal{I}\_\*^\circ) \right] \tag{1}$$

$$\rho(z) = \rho\_1^\* [1 - e\_m \cos(\pi \zeta)] \tag{2}$$

$$
\mu(z) = \mu\_1^\* [1 - e\_0 \cos(\pi \zeta)] \tag{3}
$$

Nanoporosity-2:

$$E(z) = E\_1^\* \left[ 1 - c\_0 \cos \left( \frac{\pi \zeta}{2} + \frac{\pi}{4} \right) \right] \tag{4}$$

$$\rho(z) = \rho\_1^\* \left[ 1 - \varepsilon\_m \cos \left( \frac{\pi \zeta}{2} + \frac{\pi}{4} \right) \right] \tag{5}$$

$$
\mu(z) = \mu\_1^\* \left[ 1 - e\_0 \cos \left( \frac{\pi \zeta}{2} + \frac{\pi}{4} \right) \right] \tag{6}
$$

where *ζ* = *z*/*h*, the nanoporosity coefficients are *e*0= 1 − *E*<sup>∗</sup> 0/*E*<sup>∗</sup> <sup>1</sup> (0 ≤ *e*<sup>0</sup> < 1) and *em*= 1− *ρ*<sup>∗</sup> 0/*ρ*<sup>∗</sup> <sup>1</sup> (0 ≤ *em* < 1), *ρ*<sup>∗</sup> <sup>0</sup> and *ρ*<sup>∗</sup> <sup>1</sup> are the minimum and maximum values of the mass density, respectively. The minimum Young's modulus *E*∗ <sup>0</sup> and the maximum value *E*<sup>∗</sup> <sup>1</sup> are related to the minimum shear modulus *μ*∗ <sup>0</sup> and the maximum value *μ*<sup>∗</sup> <sup>1</sup> according to *μ*<sup>∗</sup> *<sup>i</sup>* = *E*<sup>∗</sup> *<sup>i</sup>* /[2(1+*ν*)] (*i* = 0, 1), in which *ν* indicates Poisson's ratio.

For an open-cell metal foam, we have [75,76]:

$$\frac{E\_0^\*}{E\_1^\*} = \left(\frac{\rho\_0^\*}{\rho\_1^\*}\right)^2\tag{7}$$

Thus, the relation between *e*<sup>0</sup> and *em* is obtained as:

$$
\varepsilon\_m = 1 - \sqrt{1 - \varepsilon\_0} \tag{8}
$$

**Figure 1.** Schematic of a functionally graded (FG) nanoporous metal foam (NPMF) cylindrical nanoshell. (**a**) Coordinate system; (**b**) nanoporosity-1; (**c**) nanoporosity-2.

Figures 2 and 3 give the variations of mass density and Young's modulus, respectively, of the FG NPMF nanoshell along the thickness direction. Note that both kinds of nanoporosity distribution have the same minimum and maximum values of mass density and elasticity modulus. For the nanoporosity-1 nanoshell, it possesses the minimum values of mass density and Young's modulus on the middle surface (*z* = 0) of the nanoshell; while the maximum values are on the outer (*z* = *h*/2) and inner (*z* = −*h*/2) surfaces which are equal to the values of the nanoshell that consisted of solid metal. For the nanoporosity-2, mass density and elasticity modulus have the minimum values on the inner surface and gradually increase to the maximum values on the outer surface of the shell.

**Figure 2.** Variation of mass density of FG NPMF nanoshell: (**a**) nanoporosity-1; (**b**) nanoporosity-2.

**Figure 3.** Variation of Young's modulus of FG NPMF nanoshell: (**a**) nanoporosity-1; (**b**) nanoporosity-2.

#### **3. Theory and Formulation**

#### *3.1. General SGT*

ε

As we know, the strain energy density in the CT is described as the function of the strain tensor ε. The strain energy density in Mindlin's SGT, however, also incorporates the third-order strain gradient tensor ξ. Therefore, the strain energy density *W* has the most general form [40,77]:

$$\mathcal{W}(\mathfrak{e}, \mathfrak{k}) = \frac{1}{2}\lambda \varepsilon\_{\bar{i}\bar{i}}\varepsilon\_{\bar{j}\bar{j}} + \mu \varepsilon\_{\bar{i}\bar{j}}\varepsilon\_{\bar{i}\bar{j}} + a\_1 \mathfrak{z}\_{i\bar{k}k} \mathfrak{z}\_{i\bar{j}\bar{j}} + a\_2 \mathfrak{z}\_{k\bar{j}\bar{j}} \mathfrak{z}\_{i\bar{k}\bar{k}} + a\_3 \mathfrak{z}\_{i\bar{j}\bar{k}} \mathfrak{z}\_{i\bar{i}\bar{k}} + a\_4 \mathfrak{z}\_{i\bar{j}k} \mathfrak{z}\_{i\bar{j}\bar{k}} + a\_5 \mathfrak{z}\_{k\bar{j}i} \mathfrak{z}\_{i\bar{j}\bar{k}} \tag{9}$$

in which *ai* (*i* = 1, 2, ... , 5) are additional constants which can accommodate the small-scale effect of micro/nanostructures, and *λ* is Lame's first parameter defined as [78,79]:

$$
\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}\tag{10}
$$

ξ ε

In Equation (9), the third-order strain gradient tensor ξ and infinitesimal strain tensor ε are defined as [40]:

$$\boldsymbol{\varepsilon} = \frac{1}{2} \left( \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{\mathrm{T}} \right), \quad \boldsymbol{\varepsilon}\_{\mathrm{i}\bar{\jmath}} = \boldsymbol{\varepsilon}\_{\bar{\jmath}\bar{\imath}} = \frac{1}{2} \left( \boldsymbol{u}\_{\bar{\imath},\bar{\jmath}} + \boldsymbol{u}\_{\bar{\jmath},\bar{\imath}} \right) \tag{11}$$

$$\mathfrak{E} = \nabla \mathfrak{e}\_{\prime} \quad \mathfrak{J}\_{j|k} = \mathfrak{J}\_{i|\bar{k}\rangle} = \mathfrak{e}\_{\bar{\beta}k,i} = \frac{1}{2} \left( \mu\_{j,k} + \mu\_{k,\bar{j}} \right)\_{\ ,i} \tag{12}$$

in which **u** represents the displacement vector and ∇ is gradient operator.

The double stress tensor *τijk* and Cauchy stress tensor *σij* are written as [80]:

$$
\sigma\_{\bar{i}\bar{j}} = \sigma\_{\bar{j}\bar{i}} = \frac{\partial \mathcal{W}}{\partial \varepsilon\_{\bar{i}\bar{j}}} = \lambda \,\, \varepsilon\_{kk} \,\, \delta\_{\bar{i}\bar{j}} + 2\mu \,\, \varepsilon\_{\bar{i}\bar{j}} \tag{13}
$$

$$\begin{aligned} \mathfrak{T}\_{\mathrm{ijk}} = \mathfrak{T}\_{\mathrm{ikj}} = \frac{\partial W}{\partial \xi\_{\mathrm{ijk}}} = \frac{a\_{\mathrm{i}}}{2} \left( \mathfrak{T}\_{\mathrm{j}pp} \delta\_{\mathrm{ik}} + 2 \mathfrak{I}\_{ppp} \delta\_{\mathrm{jk}} + \mathfrak{I}\_{\mathrm{k}pp} \delta\_{\mathrm{ij}} \right) + 2a\_2 \mathfrak{I}\_{ipp} \delta\_{\mathrm{jk}} + a\_3 \left( \mathfrak{I}\_{ppp} \delta\_{\mathrm{ik}} + \mathfrak{I}\_{pp\mathrm{k}} \delta\_{\mathrm{ij}} \right) \\ &+ 2a\_4 \mathfrak{I}\_{ijk} + a\_5 \left( \mathfrak{I}\_{kji} + \mathfrak{I}\_{p\mathrm{ik}} \right) \end{aligned} \tag{14}$$

where *δij* represent the Kronecker delta.

#### *3.2. Constitutive Relations and Strain Energy*

The displacement field for the FG NPMF cylindrical nanoshell according to the FSD shell theory can be defined as [81–85]:

$$\left\{ \begin{array}{c} u\_x(x, y, z, t) \\ u\_y(x, y, z, t) \\ u\_z(x, y, z, t) \end{array} \right\} = \left\{ \begin{array}{c} u(x, y, t) \\ v(x, y, t) \\ w(x, y, t) \end{array} \right\} + z \left\{ \begin{array}{c} \psi\_x(x, y, t) \\ \psi\_y(x, y, t) \\ 0 \end{array} \right\} \tag{15}$$

In Equation (15), *ux*, *uy*, and *uz* stand for the displacements of any point in the nanoshell along the *x*, *y*, and *z* directions, respectively; *u*, *v*, and *w* denote displacement components of a point at the middle surface; *ψ<sup>x</sup>* and *ψ<sup>y</sup>* are the rotations of the transverse normals about the *y* and *x* axes, respectively; and *t* denotes time.

The nonzero constituents of strain tensor ε are given by: [86,87]

$$\begin{cases} \varepsilon\_{xx} = \phi\_0 + z\phi\_{1'}\\ \varepsilon\_{yy} = \phi\_0 + z\phi\_{1'}\\ \varepsilon\_{xy} = \varepsilon\_{yx} = (k\_0 + zk\_1)/2, \\ \varepsilon\_{yz} = \varepsilon\_{zy} = \gamma\_2/2, \\ \varepsilon\_{xz} = \varepsilon\_{zx} = \gamma\_1/2. \end{cases} \tag{16}$$

where *φ*0, *ϕ*0, and *k*<sup>0</sup> are the middle surface strains, *φ*1, *ϕ*1, and *k*<sup>1</sup> are changes in the curvature and torsion of the middle surface, and *γ*<sup>1</sup> and *γ*<sup>2</sup> are the transverse shear strains. They are given by:

$$\begin{cases} \begin{array}{l} \phi\_{0} = \frac{\partial \boldsymbol{\mu}}{\partial \mathbf{x}}, \phi\_{1} = \frac{\partial \boldsymbol{\psi}\_{x}}{\partial \mathbf{x}}, \varphi\_{0} = \frac{\boldsymbol{w}}{\mathbb{R}} + \frac{\partial \boldsymbol{\Gamma}}{\partial \boldsymbol{y}}, \varphi\_{1} = \frac{\partial \boldsymbol{\psi}\_{y}}{\partial \mathbf{y}}, k\_{0} = \frac{\partial \boldsymbol{\psi}}{\partial \mathbf{x}} + \frac{\partial \boldsymbol{\mu}}{\partial \mathbf{y}}, \\\ k\_{1} = \frac{\partial \boldsymbol{\psi}\_{y}}{\partial \mathbf{x}} + \frac{\partial \boldsymbol{\psi}\_{x}}{\partial \mathbf{y}}, \gamma\_{1} = \frac{\partial \boldsymbol{w}}{\partial \mathbf{x}} + \boldsymbol{\psi}\_{\mathbf{x}}, \gamma\_{2} = \boldsymbol{\psi}\_{y} - \frac{\boldsymbol{v}}{\mathbb{R}} + \frac{\partial \boldsymbol{w}}{\partial \mathbf{y}}. \end{array} \tag{17}$$

According to Equation (12), the following nonzero constituents of strain gradient tensor ξ are obtained:

⎧ ⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎩ *<sup>ξ</sup>xxx* = *∂φ*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>z</sup> ∂φ*<sup>1</sup> *<sup>∂</sup><sup>x</sup>* , *<sup>ξ</sup>yyy* <sup>=</sup> *∂ϕ*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>z</sup> ∂ϕ*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* , *<sup>ξ</sup>xyy* <sup>=</sup> *∂ϕ*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>z</sup> ∂ϕ*<sup>1</sup> *<sup>∂</sup><sup>x</sup>* , *<sup>ξ</sup>zyy* = *<sup>ϕ</sup>*1, *<sup>ξ</sup>yxx* = *∂φ*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>z</sup> ∂φ*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* , *<sup>ξ</sup>xxy* <sup>=</sup> *<sup>ξ</sup>xyx* <sup>=</sup> <sup>1</sup> 2 *∂k*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>z</sup> <sup>∂</sup>k*<sup>1</sup> *∂x* , *ξzxx* = *φ*1, *ξyxy* = *ξyyx* = <sup>1</sup> 2 *∂k*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>z</sup> <sup>∂</sup>k*<sup>1</sup> *∂y* , *ξzxy* = *ξzyx* = *<sup>k</sup>*<sup>1</sup> <sup>2</sup> , *<sup>ξ</sup>xxz* <sup>=</sup> *<sup>ξ</sup>xzx* <sup>=</sup> <sup>1</sup> 2 *∂γ*<sup>1</sup> *<sup>∂</sup><sup>x</sup>* , *ξyxz* = *ξyzx* = <sup>1</sup> 2 *∂γ*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* , *<sup>ξ</sup>xyz* <sup>=</sup> *<sup>ξ</sup>xzy* <sup>=</sup> <sup>1</sup> 2 *∂γ*<sup>2</sup> *<sup>∂</sup><sup>x</sup>* , *<sup>ξ</sup>yyz* <sup>=</sup> *<sup>ξ</sup>yzy* <sup>=</sup> <sup>1</sup> 2 *∂γ*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* . (18)

*Nanomaterials* **2019**, *9*, 271

By inserting Equations (16) and (18) into Equations (13) and (14), one can get the nonzero constituents of σ and τ as follows:

$$\begin{cases} \sigma\_{xx} = (\lambda + 2\mu)(\phi\_0 + \varphi g\_1) + \lambda(\phi\_0 + \varphi g\_1), \\ \sigma\_{yy} = (\lambda + 2\mu)(\rho\_0 + \varphi g\_1) + \lambda(\phi\_0 + \varphi g\_1), \\ \sigma\_{xx} = \sigma\_{xx} = \mu\tau\_{yy} \\ \sigma\_{xx} = \sigma\_{xx} = \mu\tau\_{zz} \\ \sigma\_{yz} = \sigma\_{zy} = \mu\tau\_{zz} \end{cases} \tag{19}$$

$$\begin{cases} \tau\_{xxx} = \beta\lambda\_1\xi\_{xxx} + \beta\lambda\_2\xi\_{xy} + \beta\lambda\_3\xi\_{xyy}, \\ \tau\_{yxx} = \beta\lambda\_2\xi\_{yy} + \beta\lambda\_3\xi\_{xxx} + \beta\lambda\_4\xi\_{xyy}, \\ \tau\_{xx} = \beta\phi\_2\xi\_{xx} + \beta\lambda\_4\xi\_{xx} + \alpha\lambda\_5\xi\_{yyz} + 2\alpha\xi\_{xyy}, \\ \tau\_{xy} = \beta\phi\_3\xi\_{xy} + \beta\lambda\_2\xi\_{xx} + \beta\lambda\_5\xi\_{yyy}, \\ \tau\_{xy} = \beta\lambda\_2\xi\_{xx} + \beta\lambda\_3\xi\_{xy} + \beta\lambda\_4\xi\_{yyy}, \\ \tau\_{xy} = \alpha\lambda\_1\xi\_{xxx} + 2\alpha\xi\_{xx} + \beta\lambda\_5\xi\_{yyy} + \beta\lambda\_4\xi\_{yyz}, \\ \tau\_{xy} = \tau\_{xy} = \tau\_{yy} = \beta\lambda\_2\xi\_{xy} + \beta\lambda\_4\xi\_{xyy} + \beta\lambda\_5\xi\_{xy}, \\ \tau\_{xy} = \tau\_{yx} = \tau\_{yy} = \frac{\beta\lambda\_2}{2}\xi\_{xyy} + \beta\lambda\_4\xi\_{xy} + \frac{\beta}{2}\lambda\_5\xi\_{xy}, \\ \tau\_{xy} = \tau\_{$$

in which

$$\begin{cases} \beta\_1 = 2(a\_1 + a\_2 + a\_3 + a\_4 + a\_5), & \beta\_2 = a\_1 + 2a\_2, \\ \beta\_3 = a\_1 + 2a\_3, & \beta\_4 = a\_1 + 2a\_5, \\ \beta\_5 = 2(a\_2 + a\_4), & \beta\_6 = a\_3 + 2a\_4 + a\_5. \end{cases} \tag{21}$$

Based on the general SGT, the stored strain energy, *ΠS*, in a linear elastic material occupying volume *V* can be given by [77]:

$$
\Pi\_S = \frac{1}{2} \int\_V \left( \sigma\_{i\bar{j}} \varepsilon\_{i\bar{j}} + \tau\_{i\bar{j}k} \xi\_{i\bar{j}k} \right) \mathrm{d}V \tag{22}
$$

If the strain energy is symbolized by classical part *Π<sup>C</sup>* and non-classical part *ΠNC*, the total strain energy is expressed as:

$$
\Pi\_{\mathbb{S}} = \Pi\_{\mathrm{NC}} + \Pi\_{\mathbb{C}} \tag{23}
$$

in which,

$$\Pi\_{\mathbb{C}} = \frac{1}{2} \int\_{A} \left( N\_{\text{xx}} \phi\_0 + M\_{\text{xx}} \phi\_1 + N\_{\text{xy}} k\_0 + M\_{\text{xy}} k\_1 + N\_{\text{yy}} \phi\_0 + M\_{\text{yy}} \phi\_1 + Q\_x \gamma\_1 + Q\_y \gamma\_2 \right) \mathrm{d}A \tag{24}$$

$$\begin{array}{rcl} II\_{NC} = & \frac{1}{2} f\_A \left( T\_{xxx} \frac{\partial \phi\_0}{\partial x} + M\_{xxx} \frac{\partial \phi\_1}{\partial x} + T\_{yxx} \frac{\partial \phi\_0}{\partial y} + M\_{yxx} \frac{\partial \phi\_1}{\partial y} + T\_{zxx} \phi\_1 \\ & + T\_{xyy} \frac{\partial \phi\_0}{\partial x} + M\_{xyy} \frac{\partial \phi\_1}{\partial x} + T\_{yyy} \frac{\partial \phi\_0}{\partial y} + M\_{yyy} \frac{\partial \phi\_1}{\partial y} \\ & + T\_{zyy} \phi\_1 + T\_{xxy} \frac{\partial \phi\_0}{\partial x} + M\_{xxy} \frac{\partial \phi\_1}{\partial x} + T\_{yxz} \frac{\partial \phi\_0}{\partial y} \\ & + M\_{yxx} \frac{\partial k\_1}{\partial y} + T\_{zxy} k\_1 + T\_{xxz} \frac{\partial \gamma\_1}{\partial x} + T\_{yxz} \frac{\partial \gamma\_1}{\partial y} \\ & + T\_{xyz} \frac{\partial \gamma\_2}{\partial x} + T\_{yzz} \frac{\partial \gamma\_2}{\partial y} \end{array} \tag{25}$$

*Nanomaterials* **2019**, *9*, 271

In Equations (24) and (25), the non-classical and classical resultant moments and forces are expressed as follows:

$$\begin{array}{ll} N\_{ij} = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \sigma\_{ij} \mathrm{d}z, & M\_{ij} = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \sigma\_{ij} z \mathrm{d}z, & Q\_i = K\_S \int\_{-\frac{h}{2}}^{\frac{h}{2}} \sigma\_{iz} \mathrm{d}z, \\\ T\_{ijk} = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \pi\_{ijk} \mathrm{d}z, & M\_{ijk} = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \pi\_{ijk} z \mathrm{d}z. \end{array} \tag{26}$$

where *KS* = 5/6 denotes the shear correction factor [88–91]; the non-classical and classical resultant moments and forces are given in Appendix A in detail.

#### *3.3. Kinetic Energy and External Work*

According to the FSD shell theory, the kinetic energy of the FG NPMF nanoshell, *ΠT*, is written as:

$$\begin{split} II\_{T} &= \frac{1}{2} \int\_{A} \int\_{-\frac{\Psi}{2}}^{\frac{\Psi}{2}} \rho(z) \left[ \left( \frac{\partial u}{\partial t} + z \frac{\partial \psi\_{x}}{\partial t} \right)^{2} + \left( \frac{\partial v}{\partial t} + z \frac{\partial \psi\_{y}}{\partial t} \right)^{2} + \left( \frac{\partial w}{\partial t} \right)^{2} \right] \mathrm{d}z \mathrm{d}A \\ &= \frac{1}{2} \int\_{A} \left[ I\_{0} \left( \frac{\partial u}{\partial t} \right)^{2} + 2I\_{1} \left( \frac{\partial u}{\partial t} \right) \left( \frac{\partial \psi\_{x}}{\partial t} \right) + I\_{0} \left( \frac{\partial v}{\partial t} \right)^{2} + I\_{0} \left( \frac{\partial w}{\partial t} \right)^{2} \\ &\qquad + 2I\_{1} \left( \frac{\partial v}{\partial t} \right) \left( \frac{\partial \psi\_{y}}{\partial t} \right) + I\_{2} \left( \frac{\partial \psi\_{x}}{\partial t} \right)^{2} + I\_{2} \left( \frac{\partial \psi\_{y}}{\partial t} \right)^{2} \right] \mathrm{d}A \end{split} \tag{27}$$

in which

$$I\_0 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \rho(z) \mathrm{d}z, \quad I\_1 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \rho(z) z \mathrm{d}z,\\ I\_2 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \rho(z) z^2 \mathrm{d}z. \tag{28}$$

Furthermore, the work *Π<sup>P</sup>* carried out by axial loads *N*<sup>0</sup> *xx* can be written as:

$$
\Pi\_p = \frac{1}{2} \int\_A \left[ N\_{xx}^0 \left( \frac{\partial w}{\partial x} \right)^2 \right] \text{d}A \tag{29}
$$

#### *3.4. Variational Formulation*

Using Hamilton's principle,

$$\delta \int\_{0}^{t} (\Pi\_T - \Pi\_S - \Pi\_P) \mathrm{d}t = 0 \tag{30}$$

Inserting Equations (23), (27) and (29) into Equation (30) yields the following governing equations:

$$
\delta u: \frac{\partial \overline{\mathcal{N}}\_{xx}}{\partial x} + \frac{\partial \overline{\mathcal{N}}\_{xy}}{\partial y} = I\_0 \frac{\partial^2 u}{\partial t^2} + I\_1 \frac{\partial^2 \psi\_x}{\partial t^2} \tag{31}
$$

$$\delta v: \frac{\partial \overline{N}\_{xy}}{\partial x} + \frac{\partial \overline{N}\_{yy}}{\partial y} + \frac{\overline{Q}\_y}{R} = I\_0 \frac{\partial^2 v}{\partial t^2} + I\_1 \frac{\partial^2 \psi\_y}{\partial t^2} \tag{32}$$

$$\delta w : \frac{\partial \overline{Q}\_x}{\partial x} + \frac{\partial \overline{Q}\_y}{\partial y} - \frac{\overline{N}\_{yy}}{R} + N\_{xx}^0 \frac{\partial^2 w}{\partial x^2} = I\_0 \frac{\partial^2 w}{\partial t^2} \tag{33}$$

$$\delta\psi\_x: \frac{\partial \overline{M}\_{xx}}{\partial x} + \frac{\partial \overline{M}\_{xy}}{\partial y} - \overline{Q}\_x = I\_2 \frac{\partial^2 \psi\_x}{\partial t^2} + I\_1 \frac{\partial^2 u}{\partial t^2} \tag{34}$$

$$\delta\psi\_{\mathcal{Y}} : \frac{\partial \overline{M}\_{xy}}{\partial \mathbf{x}} + \frac{\partial \overline{M}\_{yy}}{\partial y} - \overline{Q}\_{\mathcal{Y}} = l\_2 \frac{\partial^2 \psi\_{\mathcal{Y}}}{\partial t^2} + I\_1 \frac{\partial^2 v}{\partial t^2} \tag{35}$$

where,

$$\begin{cases} \overline{N}\_{xx} = N\_{xx} - \frac{\partial T\_{xx}}{\partial y} - \frac{\partial T\_{xx}}{\partial x}, \\ \overline{N}\_{xy} = N\_{xy} - \frac{\partial T\_{yy}}{\partial y} - \frac{\partial T\_{xy}}{\partial x}, \\ \overline{N}\_{yy} = N\_{yy} - \frac{\partial T\_{yy}}{\partial y} - \frac{\partial T\_{xy}}{\partial x} - \frac{T\_{yy}}{K}, \\ \overline{M}\_{xx} = M\_{xx} + T\_{zxx} - \frac{\partial M\_{yx}}{\partial x} - \frac{\partial M\_{xz}}{\partial x}, \\ \overline{M}\_{yy} = M\_{yy} + T\_{zy} - \frac{\partial M\_{yy}}{\partial y} - \frac{\partial M\_{xy}}{\partial x}, \\ \overline{M}\_{xy} = M\_{xy} + T\_{xy} - \frac{\partial M\_{yx}}{\partial y} - \frac{\partial M\_{xz}}{\partial x}, \\ \overline{Q}\_{x} = Q\_{x} - \frac{\partial T\_{yx}}{\partial y} - \frac{\partial T\_{xx}}{\partial x}, \\ \overline{Q}\_{y} = Q\_{y} - \frac{\partial T\_{yy}}{\partial y} - \frac{\partial T\_{xx}}{\partial x}. \end{cases} \tag{36}$$

Simultaneously, boundary conditions are derived as:

$$\begin{array}{lcl}\delta u = 0 & \text{or} & (\overline{\mathcal{N}}\_{xx})n\_x + (\overline{\mathcal{N}}\_{xy})n\_y = 0, \\\delta u\_{,x} = 0 & \text{or} & (T\_{xxx})n\_x + (T\_{yxx})n\_y = 0, \\\delta u\_{,y} = 0 & \text{or} & (T\_{xxy})n\_x + (T\_{yyx})n\_y = 0. \end{array} \tag{37}$$

$$\begin{array}{lcl}\delta v = 0 \text{ or } \left(\overline{\mathcal{N}}\_{xy} - \frac{T\_{xyx}}{\overline{\mathcal{R}}}\right) n\_x + \left(\overline{\mathcal{N}}\_{yy} - \frac{T\_{yx}}{\overline{\mathcal{R}}}\right) n\_y = 0, \\\delta v\_{,x} = 0 \text{ or } \left(T\_{xxy}\right) n\_x + \left(T\_{yyx}\right) n\_y = 0, \\\delta v\_{,y} = 0 \text{ or } \left(T\_{xyy}\right) n\_x + \left(T\_{yyy}\right) n\_y = 0. \end{array} \tag{38}$$

$$\begin{array}{ll}\delta w = 0 & \text{or} \quad \left(\overline{\mathbb{Q}}\_{x} + \frac{T\_{\text{xyz}}}{\overline{\mathbb{R}}}\right) n\_{x} + \left(\overline{\mathbb{Q}}\_{y} + \frac{T\_{\text{xyz}}}{\overline{\mathbb{R}}}\right) n\_{y} = 0, \\\delta w\_{,x} = 0 & \text{or} \quad \left(T\_{xxz}\right) n\_{x} + \left(T\_{yxz}\right) n\_{y} = 0, \\\delta w\_{,y} = 0 & \text{or} \quad \left(T\_{xyz}\right) n\_{x} + \left(T\_{yyz}\right) n\_{y} = 0. \end{array} \tag{39}$$

$$\begin{array}{lcl}\delta\psi\_{\mathbf{x}} = 0 & \text{or} & (\overline{M}\_{\mathbf{x}\mathbf{x}} + T\_{\mathbf{x}\mathbf{x}\mathbf{z}})n\_{\mathbf{x}} + (\overline{M}\_{\mathbf{x}\mathbf{y}} + T\_{\mathbf{y}\mathbf{z}\mathbf{z}})n\_{\mathbf{y}} = 0, \\\delta\psi\_{\mathbf{x},\mathbf{z}} = 0 & \text{or} & (M\_{\mathbf{x}\mathbf{x}\mathbf{z}})n\_{\mathbf{x}} + (M\_{\mathbf{y}\mathbf{x}\mathbf{z}})n\_{\mathbf{y}} = 0, \\\delta\psi\_{\mathbf{x},\mathbf{y}} = 0 & \text{or} & (M\_{\mathbf{x}\mathbf{y}})n\_{\mathbf{x}} + (M\_{\mathbf{y}\mathbf{x}\mathbf{z}})n\_{\mathbf{y}} = 0. \end{array} \tag{40}$$

$$\begin{cases} \delta \psi\_{\mathcal{Y}} = 0 \text{ or } \left( \overline{M}\_{xy} + T\_{xyz} \right) n\_x + \left( \overline{M}\_{yy} + T\_{yyz} \right) n\_{\mathcal{Y}} = 0, \\ \delta \psi\_{\mathcal{Y},x} = 0 \text{ or } \left( M\_{xxy} \right) n\_x + \left( M\_{yyx} \right) n\_{\mathcal{Y}} = 0, \\ \delta \psi\_{\mathcal{Y},\mathcal{Y}} = 0 \text{ or } \left( M\_{xyy} \right) n\_x + \left( M\_{yyy} \right) n\_{\mathcal{Y}} = 0. \end{cases} \tag{41}$$

where *nx* as well as *ny* indicate the direction cosines of the outward unit normal to the boundary of the mid-plane.

Substituting Equation (36) into Equations (31)–(35) and considering Equation (17) and Appendix A, it yields the governing equations in terms of *u*, *v*, *w*, *ψx*, and *ψy*:

*<sup>A</sup>*<sup>11</sup> *<sup>∂</sup>*2*<sup>u</sup> <sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *<sup>A</sup>*55*∂*2*<sup>u</sup> <sup>∂</sup>y*<sup>2</sup> <sup>+</sup> *<sup>∂</sup>*2*<sup>v</sup> <sup>∂</sup>x∂<sup>y</sup>* − *E*<sup>1</sup> *∂*4*u <sup>∂</sup>x*<sup>4</sup> <sup>−</sup> *<sup>E</sup>*<sup>6</sup> <sup>2</sup> *<sup>∂</sup>*4*<sup>u</sup> ∂y*<sup>4</sup> − *E*<sup>3</sup> + *E*<sup>4</sup> + *E*<sup>5</sup> + *<sup>E</sup>*<sup>6</sup> 2 *∂*4*u <sup>∂</sup>y*2*∂x*<sup>2</sup> <sup>−</sup> (2*E*2+*E*3+*E*4+*E*6) 2 *∂*4*v <sup>∂</sup>y∂x*<sup>3</sup> <sup>+</sup> *<sup>∂</sup>*4*<sup>v</sup> ∂y*3*∂x* +(*A*<sup>12</sup> <sup>+</sup> *<sup>A</sup>*55) *<sup>∂</sup>*2*<sup>v</sup> <sup>∂</sup>y∂<sup>x</sup>* <sup>+</sup> *<sup>A</sup>*<sup>12</sup> *<sup>R</sup> <sup>∂</sup><sup>w</sup> <sup>∂</sup><sup>x</sup>* <sup>−</sup> <sup>2</sup>*E*2+*E*3+*E*<sup>4</sup> <sup>2</sup>*<sup>R</sup> <sup>∂</sup>*3*<sup>w</sup> ∂y*2*∂x* <sup>−</sup>*E*<sup>2</sup> *<sup>R</sup> <sup>∂</sup>*3*<sup>w</sup> <sup>∂</sup>x*<sup>3</sup> <sup>+</sup> *<sup>B</sup>*<sup>11</sup> *<sup>∂</sup>*2*ψ<sup>x</sup> <sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *<sup>B</sup>*<sup>55</sup> *<sup>∂</sup>*2*ψ<sup>x</sup> <sup>∂</sup>y*<sup>2</sup> − *F*<sup>3</sup> + *F*<sup>4</sup> + *F*<sup>5</sup> + *<sup>F</sup>*<sup>6</sup> 2 *<sup>∂</sup>*4*ψ<sup>x</sup> <sup>∂</sup>y*2*∂x*<sup>2</sup> − *F*<sup>1</sup> *∂*4*ψ<sup>x</sup> ∂x*<sup>4</sup> <sup>−</sup>*F*<sup>6</sup> 2 *∂*4*ψ<sup>x</sup> <sup>∂</sup>y*<sup>4</sup> <sup>−</sup> <sup>2</sup>*F*2+*F*3+*F*4+*F*<sup>6</sup> 2 *<sup>∂</sup>*4*ψ<sup>y</sup> <sup>∂</sup>y∂x*<sup>3</sup> <sup>+</sup> *<sup>∂</sup>*4*ψ<sup>y</sup> ∂y*3*∂x* +(*B*<sup>12</sup> + *B*55) *∂*2*ψ<sup>y</sup> <sup>∂</sup>y∂<sup>x</sup>* = *I*<sup>0</sup> *∂*2*u <sup>∂</sup>t*<sup>2</sup> + *<sup>I</sup>*<sup>1</sup> *∂*2*ψ<sup>x</sup> ∂t*<sup>2</sup> (42)

+ 

+

−

−

+

<sup>−</sup>2*E*2+*E*3+*E*4+*E*<sup>6</sup> 2 *∂*4*u <sup>∂</sup>y∂x*<sup>3</sup> <sup>+</sup> *<sup>∂</sup>*4*<sup>u</sup> ∂y*3*∂x* <sup>+</sup> (*A*<sup>12</sup> <sup>+</sup> *<sup>A</sup>*55) *<sup>∂</sup>*2*<sup>u</sup> <sup>∂</sup>y∂<sup>x</sup>* − *E*<sup>3</sup> + *E*<sup>4</sup> + *E*<sup>5</sup> + *<sup>E</sup>*<sup>6</sup> 2 *∂*4*v <sup>∂</sup>y*2*∂x*<sup>2</sup> <sup>−</sup> *<sup>E</sup>*<sup>6</sup> <sup>2</sup> *<sup>∂</sup>*4*<sup>v</sup> <sup>∂</sup>x*<sup>4</sup> − *E*<sup>1</sup> *∂*4*v <sup>∂</sup>y*<sup>4</sup> <sup>+</sup> *<sup>A</sup>*<sup>55</sup> *<sup>∂</sup>*2*<sup>v</sup> <sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *<sup>A</sup>*<sup>11</sup> *<sup>∂</sup>*2*<sup>v</sup> <sup>∂</sup>y*<sup>2</sup> <sup>−</sup> *KSA*<sup>55</sup> *<sup>R</sup>*<sup>2</sup> *v* + <sup>1</sup> *<sup>R</sup>* (*A*<sup>11</sup> <sup>+</sup> *KSA*55) *<sup>∂</sup><sup>w</sup> <sup>∂</sup><sup>y</sup>* <sup>−</sup> <sup>1</sup> *R E*3+*E*4+2*E*<sup>5</sup> <sup>2</sup> <sup>+</sup> *<sup>A</sup>*<sup>3</sup> <sup>+</sup> *<sup>A</sup>*<sup>4</sup> <sup>+</sup> *<sup>A</sup>*<sup>5</sup> 2 *∂*3*w ∂y∂x*<sup>2</sup> +*E*1+*E*<sup>6</sup> *<sup>R</sup> <sup>∂</sup>*3*<sup>w</sup> <sup>∂</sup>y*<sup>3</sup> − 1 <sup>2</sup>*<sup>R</sup>* (2*A*<sup>1</sup> <sup>+</sup> <sup>2</sup>*A*<sup>3</sup> <sup>+</sup> <sup>2</sup>*A*5) <sup>+</sup> *<sup>B</sup>*<sup>12</sup> <sup>+</sup> *<sup>B</sup>*55 *<sup>∂</sup>*2*ψ<sup>x</sup> ∂y∂x* <sup>−</sup>2*F*2+*F*3+*F*4+*F*<sup>6</sup> 2 *<sup>∂</sup>*4*ψ<sup>x</sup> <sup>∂</sup>y*3*∂<sup>x</sup>* <sup>+</sup> *<sup>∂</sup>*4*ψ<sup>x</sup> ∂y∂x*<sup>3</sup> + *B*<sup>55</sup> + <sup>1</sup> <sup>2</sup>*<sup>R</sup>* (2*A*<sup>4</sup> + *A*5) *<sup>∂</sup>*2*ψ<sup>y</sup> ∂x*<sup>2</sup> + *<sup>B</sup>*<sup>11</sup> <sup>−</sup> *<sup>E</sup>*4+*E*<sup>6</sup> *R ∂*2*ψ<sup>y</sup> <sup>∂</sup>y*<sup>2</sup> − *F*<sup>1</sup> *∂*4*ψ<sup>y</sup> <sup>∂</sup>y*<sup>4</sup> <sup>−</sup> *<sup>F</sup>*<sup>6</sup> 2 *∂*4*ψ<sup>y</sup> ∂x*<sup>4</sup> − *F*<sup>3</sup> + *F*<sup>4</sup> + *F*<sup>5</sup> + *<sup>F</sup>*<sup>6</sup> 2 *<sup>∂</sup>*4*ψ<sup>y</sup> <sup>∂</sup>y*2*∂x*<sup>2</sup> <sup>+</sup> *KSA*<sup>55</sup> *<sup>R</sup> ψ<sup>y</sup>* = *I*<sup>0</sup> *∂*2*v <sup>∂</sup>t*<sup>2</sup> + *<sup>I</sup>*<sup>1</sup> *∂*2*ψ<sup>y</sup> ∂t*<sup>2</sup> (43) <sup>−</sup> *<sup>A</sup>*<sup>12</sup> *<sup>R</sup> <sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>x</sup>* <sup>−</sup> <sup>1</sup> <sup>2</sup>*<sup>R</sup>* (2*E*<sup>2</sup> <sup>+</sup> *<sup>E</sup>*<sup>3</sup> <sup>+</sup> *<sup>E</sup>*4) *<sup>∂</sup>*3*<sup>u</sup> <sup>∂</sup>y*2*∂<sup>x</sup>* <sup>+</sup> *<sup>E</sup>*<sup>2</sup> <sup>2</sup>*<sup>R</sup> <sup>∂</sup>*3*<sup>u</sup> ∂x*<sup>3</sup> − 1 *R A*<sup>11</sup> + *KSA*<sup>55</sup> + *<sup>E</sup>*<sup>6</sup> 2*R*<sup>2</sup> *∂v <sup>∂</sup><sup>y</sup>* <sup>+</sup> <sup>1</sup> *R E*3+*E*4+2*E*<sup>5</sup> <sup>2</sup> <sup>+</sup> *<sup>A</sup>*3+2*A*4+*A*<sup>5</sup> 2 *∂*3*v ∂y∂x*<sup>2</sup> + <sup>1</sup> *R* 2*E*1+*E*<sup>6</sup> 2 *∂*3*v <sup>∂</sup>y*<sup>3</sup> <sup>−</sup> (*A*<sup>3</sup> <sup>+</sup> <sup>2</sup>*A*<sup>4</sup> <sup>+</sup> *<sup>A</sup>*5) *<sup>∂</sup>*4*<sup>w</sup> ∂y*2*∂x*<sup>2</sup> <sup>−</sup>*E*<sup>6</sup> 2 *∂*4*w <sup>∂</sup>x*<sup>4</sup> <sup>+</sup> *<sup>∂</sup>*4*<sup>w</sup> ∂y*<sup>4</sup> + *<sup>A</sup>*55*KS* + <sup>2</sup>*E*5+*A*<sup>3</sup> 2*R*<sup>2</sup> *∂*2*w ∂x*<sup>2</sup> + *<sup>A</sup>*55*KS* + <sup>2</sup>*E*1+*E*<sup>6</sup> *R*2 *∂*2*w <sup>∂</sup>y*<sup>2</sup> <sup>−</sup> *<sup>A</sup>*<sup>11</sup> *<sup>R</sup>*<sup>2</sup> *w* − *<sup>A</sup>*1−*A*3−2*A*4−3*A*<sup>5</sup> 2 *<sup>∂</sup>*3*ψ<sup>x</sup> ∂y*2*∂x* − *E*4+*E*<sup>6</sup> <sup>2</sup> <sup>−</sup> *<sup>F</sup>*<sup>2</sup> *R ∂*3*ψ<sup>x</sup> <sup>∂</sup>x*<sup>3</sup> + *KSA*<sup>55</sup> <sup>−</sup> *<sup>B</sup>*<sup>12</sup> *<sup>R</sup>* <sup>+</sup> *<sup>A</sup>*1+*A*<sup>3</sup> *R*2 *∂ψ<sup>x</sup> ∂x* − *A*1+*A*3+2*A*4+3*A*<sup>5</sup> <sup>2</sup> <sup>+</sup> *<sup>F</sup>*3+*F*4+2*F*<sup>5</sup> *R <sup>∂</sup>*3*ψ<sup>y</sup> <sup>∂</sup>y∂x*<sup>2</sup> − *E*4+*E*<sup>6</sup> <sup>2</sup> <sup>+</sup> <sup>−</sup>2*F*<sup>1</sup> *R ∂*3*ψ<sup>y</sup> ∂y*<sup>3</sup> + *KSA*<sup>55</sup> <sup>−</sup> *<sup>B</sup>*<sup>11</sup> *<sup>R</sup>* <sup>+</sup> *<sup>E</sup>*4+*E*<sup>6</sup> 2*R*<sup>2</sup> *∂ψ<sup>y</sup> <sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>N</sup>*<sup>0</sup> *xx <sup>∂</sup>*2*<sup>w</sup> <sup>∂</sup>x*<sup>2</sup> = *<sup>I</sup>*<sup>0</sup> *∂*2*w ∂t*<sup>2</sup> (44) *<sup>B</sup>*<sup>11</sup> *<sup>∂</sup>*2*<sup>u</sup> <sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *<sup>B</sup>*<sup>55</sup> *<sup>∂</sup>*2*<sup>u</sup> <sup>∂</sup>y*<sup>2</sup> − *F*<sup>3</sup> + *F*<sup>4</sup> + *F*<sup>5</sup> + *<sup>F</sup>*<sup>6</sup> 2 *∂*4*u ∂y*2*∂x*<sup>2</sup> −*F*<sup>1</sup> *∂*4*u <sup>∂</sup>x*<sup>4</sup> <sup>−</sup> *<sup>F</sup>*<sup>6</sup> <sup>2</sup> *<sup>∂</sup>*4*<sup>u</sup> <sup>∂</sup>y*<sup>4</sup> + *<sup>B</sup>*<sup>12</sup> <sup>+</sup> *<sup>B</sup>*<sup>55</sup> <sup>−</sup> <sup>1</sup> <sup>2</sup>*<sup>R</sup>* (*A*<sup>1</sup> + *A*<sup>3</sup> + 2*A*5) *<sup>∂</sup>*2*<sup>v</sup> <sup>∂</sup>y∂<sup>x</sup>* <sup>−</sup>2*F*2+*F*3+*F*4+*F*<sup>6</sup> 2 *∂*4*v <sup>∂</sup>y*3*∂<sup>x</sup>* <sup>+</sup> *<sup>∂</sup>*4*<sup>v</sup> ∂y∂x*<sup>3</sup> − *KSA*<sup>55</sup> <sup>−</sup> *<sup>B</sup>*<sup>12</sup> *R ∂w ∂x* + *A*1+*A*3+2*A*4+3*A*<sup>5</sup> 2 *∂*3*w <sup>∂</sup>y*2*∂<sup>x</sup>* + *E*4+*E*<sup>6</sup> <sup>2</sup> <sup>−</sup> *<sup>F</sup>*<sup>2</sup> *R ∂*3*w ∂x*<sup>3</sup> −*G*<sup>1</sup> *∂*4*ψ<sup>x</sup> <sup>∂</sup>x*<sup>4</sup> <sup>−</sup> *<sup>G</sup>*<sup>6</sup> 2 *∂*4*ψ<sup>x</sup> <sup>∂</sup>y*<sup>4</sup> <sup>−</sup> (*G*<sup>3</sup> <sup>+</sup> *<sup>G</sup>*<sup>4</sup> <sup>+</sup> *<sup>G</sup>*<sup>5</sup> <sup>+</sup> *<sup>G</sup>*6) *<sup>∂</sup>*4*ψ<sup>x</sup> ∂y*2*∂x*<sup>2</sup> + *D*<sup>11</sup> + *E*<sup>4</sup> + *E*<sup>5</sup> + *<sup>E</sup>*<sup>6</sup> 2 *∂*2*ψ<sup>x</sup> <sup>∂</sup>x*<sup>2</sup> <sup>+</sup> (*D*<sup>55</sup> <sup>+</sup> <sup>2</sup>*A*<sup>4</sup> <sup>+</sup> *<sup>A</sup>*5) *<sup>∂</sup>*2*ψ<sup>x</sup> ∂y*<sup>2</sup> <sup>−</sup>*A*55*KSψ<sup>x</sup>* <sup>−</sup> *<sup>G</sup>*2+*G*3+*G*4+*G*<sup>6</sup> 2 *<sup>∂</sup>*4*ψ<sup>y</sup> <sup>∂</sup>y∂x*<sup>3</sup> <sup>+</sup> *<sup>∂</sup>*4*ψ<sup>y</sup> ∂y*3*∂x* + *D*<sup>12</sup> + *D*<sup>55</sup> + *A*<sup>1</sup> + 2*A*<sup>2</sup> + *<sup>A</sup>*<sup>3</sup> <sup>2</sup> <sup>+</sup> *<sup>A</sup>*<sup>4</sup> <sup>+</sup> *<sup>A</sup>*<sup>5</sup> 2 *∂*2*ψ<sup>y</sup> <sup>∂</sup>y∂<sup>x</sup>* = *I*<sup>2</sup> *∂*2*ψ<sup>x</sup> <sup>∂</sup>t*<sup>2</sup> + *<sup>I</sup>*<sup>1</sup> *∂*2*u ∂t*<sup>2</sup> (45) (*B*<sup>12</sup> <sup>+</sup> *<sup>B</sup>*55) *<sup>∂</sup>*2*<sup>u</sup> <sup>∂</sup>y∂<sup>x</sup>* <sup>−</sup> <sup>2</sup>*F*2+*F*3+*F*4+*F*<sup>6</sup> 2 *∂*4*u <sup>∂</sup>y∂x*<sup>3</sup> <sup>+</sup> *<sup>∂</sup>*4*<sup>u</sup> ∂y*3*∂x B*<sup>55</sup> + <sup>1</sup> <sup>2</sup>*<sup>R</sup>* (−2*A*<sup>4</sup> − *A*5) *∂*2*v <sup>∂</sup>x*<sup>2</sup> + *<sup>B</sup>*<sup>11</sup> <sup>−</sup> *<sup>E</sup>*4+*E*<sup>6</sup> 2*R ∂*2*v ∂y*<sup>2</sup> −*F*<sup>1</sup> *∂*4*v <sup>∂</sup>y*<sup>4</sup> <sup>−</sup> *<sup>F</sup>*<sup>6</sup> <sup>2</sup> *<sup>∂</sup>*4*<sup>v</sup> <sup>∂</sup>x*<sup>4</sup> − *F*<sup>3</sup> + *F*<sup>4</sup> + *F*<sup>5</sup> + *<sup>F</sup>*<sup>6</sup> 2 *∂*4*v <sup>∂</sup>y*2*∂x*<sup>2</sup> <sup>+</sup> *KSA*<sup>55</sup> *<sup>R</sup> v <sup>A</sup>*1+*A*3+2*A*4−3*A*5−*F*3−*F*4−2*F*<sup>5</sup> 2 *∂*3*w <sup>∂</sup>y∂x*<sup>2</sup> + *E*4+*E*<sup>6</sup> <sup>2</sup> <sup>−</sup> *<sup>F</sup>*<sup>1</sup> *R ∂*3*w ∂y*<sup>3</sup> *KSA*<sup>55</sup> <sup>−</sup> *<sup>B</sup>*<sup>11</sup> *R ∂w <sup>∂</sup><sup>y</sup>* <sup>−</sup> <sup>2</sup>*G*2+*G*3+*G*4+*G*<sup>6</sup> 2 *<sup>∂</sup>*4*ψ<sup>x</sup> <sup>∂</sup>y∂x*<sup>3</sup> <sup>+</sup> *<sup>∂</sup>*4*ψ<sup>x</sup> ∂y*3*∂x* <sup>−</sup> *<sup>G</sup>*<sup>6</sup> 2 *∂*4*ψ<sup>y</sup> <sup>∂</sup>x*<sup>4</sup> − *G*<sup>1</sup> *∂*4*ψ<sup>y</sup> ∂y*<sup>4</sup> *G*<sup>3</sup> + *G*<sup>4</sup> + *G*<sup>5</sup> + *<sup>G</sup>*<sup>6</sup> 2 *<sup>∂</sup>*4*ψ<sup>y</sup> <sup>∂</sup>y*2*∂x*<sup>2</sup> + (*D*<sup>55</sup> + <sup>2</sup>*A*<sup>4</sup> + *<sup>A</sup>*5) *∂*2*ψ<sup>y</sup> <sup>∂</sup>x*<sup>2</sup> + *D*<sup>11</sup> + *E*<sup>4</sup> + *E*<sup>5</sup> + *<sup>E</sup>*<sup>6</sup> 2 *∂*2*ψ<sup>y</sup> ∂y*<sup>2</sup> *D*<sup>12</sup> + *D*<sup>55</sup> + *A*<sup>1</sup> + 2*A*<sup>2</sup> + *<sup>A</sup>*<sup>3</sup> <sup>2</sup> <sup>+</sup> *<sup>A</sup>*<sup>4</sup> <sup>+</sup> <sup>3</sup>*A*<sup>5</sup> 2 *∂*2*ψ<sup>x</sup> <sup>∂</sup>y∂<sup>x</sup>* − *KSA*55*ψ<sup>y</sup>* = *I*<sup>2</sup> *∂*2*ψ<sup>y</sup> <sup>∂</sup>t*<sup>2</sup> + *<sup>I</sup>*<sup>1</sup> *∂*2*v ∂t*<sup>2</sup> (46)

in which *Aij*, *Bij*, *Dij*, *Ai*, *Bi*, *Ei*, *Fi*, and *Gi* (*i*, *j* = 1, 2, . . . , 6) are given in Appendix B.

It is worth mentioning that the present general strain gradient nanoshell model can reduce to those of MCST, MSGT, and CT. The MSGT model can be achieved if *ai* (*i* = 1, 2, ... , 5) are defined by three material length scale parameters as follows:

$$\begin{aligned} a\_1 &= \mu \left( l\_2^2 - \frac{4}{15} l\_1^2 \right), & a\_2 &= \mu \left( l\_0^2 - \frac{1}{15} l\_1^2 - \frac{1}{2} l\_2^2 \right), \\ a\_3 &= -\mu \left( \frac{4}{15} l\_1^2 + \frac{1}{2} l\_2^2 \right), & a\_4 &= \mu \left( \frac{1}{5} l\_1^2 + l\_2^2 \right), & a\_5 &= \mu \left( \frac{2}{5} l\_1^2 - l\_2^2 \right). \end{aligned} \tag{47}$$

where *l*0, *l*1, and *l*<sup>2</sup> are material length scale parameters corresponding to dilatation gradients, deviatoric stretch gradients and symmetric rotation gradients, respectively. In the following discussion, we assume that all the material length scale parameters are the same, namely, *l*<sup>0</sup> = *l*<sup>1</sup> = *l*<sup>2</sup> = *l*. In addition, by setting *a*<sup>1</sup> = *a*<sup>2</sup> = *a*<sup>3</sup> = *a*<sup>4</sup> = *a*<sup>5</sup> = 0, the present nanoshell model can be simplified to the CT-based model. Moreover, the MCST model [51] can be achieved if *ai* (*i* = 1, 2, . . . , 5) are set as:

$$a\_1 = a\_4 = -2a\_2 = -2a\_3 = -a\_5 = \mu l^2 \tag{48}$$

#### **4. Closed-Form Solution**

Herein, we employ the Navier solution technique to analyze the free vibration and axial buckling behaviors of an FG NPMF cylindrical nanoshell. Navier's method can obtain an analytical solution by introducing the double trigonometric series. Note that this method is only applicable to the simply supported boundary condition. For the other boundary conditions which are different from the simply supported boundary condition, other numerical methods such as the finite element method, differential quadrature method, finite difference method, meshless method, and wavelet method can be used. As an example, the boundary condition of the FG NPMF nanoshell considered in our study is simply supported at edges *x* = 0 as well as *x* = *L*, so one obtains:

$$\begin{cases} \boldsymbol{v} = \boldsymbol{w} = \boldsymbol{\upmu}\_{\boldsymbol{y}} = \overline{\boldsymbol{N}}\_{\boldsymbol{xx}} = \boldsymbol{0},\\ \frac{\partial \boldsymbol{\upmu}\_{\boldsymbol{y}}}{\partial \boldsymbol{y}} = \frac{\partial \boldsymbol{\upmu}\_{\boldsymbol{x}}}{\partial \boldsymbol{x}} = \frac{\partial \boldsymbol{w}}{\partial \boldsymbol{y}} = \frac{\partial \boldsymbol{v}}{\partial \boldsymbol{y}} = \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x}} = \boldsymbol{0},\\ \boldsymbol{T}\_{\boldsymbol{x}\boldsymbol{x}\boldsymbol{y}} = \boldsymbol{T}\_{\boldsymbol{x}\boldsymbol{x}\boldsymbol{z}} = \boldsymbol{M}\_{\boldsymbol{x}\boldsymbol{x}\boldsymbol{y}} = \overline{\boldsymbol{M}}\_{\boldsymbol{x}\boldsymbol{x}} + \boldsymbol{T}\_{\boldsymbol{x}\boldsymbol{x}\boldsymbol{z}} = \boldsymbol{0}. \end{cases} \tag{49}$$

The Navier procedure is used by assuming the displacements as follows:

$$\begin{Bmatrix} u(x,y,t) \\ v(x,y,t) \\ w(x,y,t) \\ \psi\_x(x,y,t) \\ \psi\_y(x,y,t) \end{Bmatrix} = \sum\_{n=1}^{\infty} \sum\_{m=1}^{\infty} \left\{ \begin{array}{l} u\_{mn}(t) \cos(a\_m x) \sin\left(\frac{my}{R}\right) \\ v\_{mn}(t) \sin(a\_m x) \cos\left(\frac{my}{R}\right) \\ w\_{mn}(t) \sin(a\_m x) \sin\left(\frac{my}{R}\right) \\ \psi\_{xmn}(t) \cos(a\_m x) \sin\left(\frac{my}{R}\right) \\ \psi\_{ymn}(t) \sin(a\_m x) \cos\left(\frac{my}{R}\right) \end{array} \right\} \tag{50}$$

in which *α<sup>m</sup>* = *mπ*/*L*, *n* is the circumferential wave number, and *m* is the axial half-wave number. Inserting Equation (50) into Equations (42)–(46) and then eliminating the trigonometric functions, the equations can be re-represented in the matrix form as:

$$\ddot{\mathbf{Md}} + \left(\mathbf{K} + \mathcal{N}\_{xx}^{0}\mathbf{K}\_{\mathbf{g}}\right)\mathbf{d} = 0\tag{51}$$

where the displacement vector **d**, mass matrix **M**, geometric stiffness matrix **Kg**, and stiffness matrix **K** are:

$$\mathbf{d} = \begin{bmatrix} \boldsymbol{u}\_{mm}, \boldsymbol{v}\_{mm}, \boldsymbol{w}\_{mm}, \boldsymbol{\psi}\_{xmn}, \boldsymbol{\psi}\_{ymn} \end{bmatrix}^{\mathrm{T}} \tag{52}$$

$$\mathbf{K} = \begin{bmatrix} K\_{11} & K\_{12} & K\_{13} & K\_{14} & K\_{15} \\ K\_{21} & K\_{22} & K\_{23} & K\_{24} & K\_{25} \\ K\_{31} & K\_{32} & K\_{33} & K\_{34} & K\_{35} \\ K\_{41} & K\_{42} & K\_{43} & K\_{44} & K\_{45} \\ \nu & \nu & \nu & \nu & \nu \end{bmatrix} \tag{53}$$

$$\mathbf{M} = \begin{bmatrix} & K\_{51} & K\_{52} & K\_{53} & K\_{54} & K\_{55} & \dots \end{bmatrix}$$

$$\mathbf{M} = \begin{bmatrix} & M\_{11} & 0 & 0 & M\_{14} & 0\\ 0 & & M\_{22} & 0 & 0 & M\_{25} \\ & 0 & 0 & M\_{33} & 0 & 0\\ & M\_{41} & 0 & 0 & M\_{44} & 0\\ & 0 & M\_{52} & 0 & 0 & M\_{55} \end{bmatrix} \tag{54}$$

$$\mathbf{K}\_{\mathfrak{F}} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & K\_{\mathfrak{F}} {}\_{3} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix} \tag{55}$$

in which the elements in these matrices are given in Appendix C.

If the dynamic displacement is considered, the form of the displacement vector **d** can be written as **d** = **d**\* ei*<sup>ω</sup><sup>t</sup>* . Once we ignore *N*<sup>0</sup> *xx*, the eigenvalue problem of free vibrating nanoshells can be obtained as:

$$\left(\mathbf{K} - \omega^2 \mathbf{M}\right)\mathbf{d}^\* = 0\tag{56}$$

where *ω* represents the natural frequency of the FG NPMF nanoshell. The non-trivial solution requires vanishing of the determinant of the coefficient matrix in Equation (56) [92–98].

Buckling loads of the FG NPMF nanoshell can be obtained by neglecting the inertia term in Equation (51). Letting *N*<sup>0</sup> *xx*= −*F*, one can get:

$$\left(\mathbf{K} - F\mathbf{K}\_{\mathcal{S}}\right)\mathbf{d} = 0\tag{57}$$

where *F* denotes the buckling load. For different combinations of *m* and *n*, there exists a minimum value which satisfies Equation (57). This minimum value is termed as the critical buckling load *F*cr.

#### **5. Validation**

Some comparative studies are first undertaken to prove the reliability of the present analysis.

#### *5.1. Example 1: Homogeneous Cylindrical Nanoshell Based on the MSGT*

In Table 1, the present results for a homogeneous simply supported cylindrical nanoshell are compared with those obtained by Zhang et al. [99]. The parameters used are: *E* = 1.06TPa, *ν* = 0.3, *ρ* = 2300 kg/m3, *R* = 2.32 nm, and *L*/*R* = 5. The frequency parameter *ω* = *ωR ρ*/*E* of the nanoshell is obtained based on the MSGT. One can see that the results from the current study coincide with those reported in Reference [99].

**Table 1.** Comparison of dimensionless natural frequency *ω* for a homogeneous nanoscale cylindrical shell.


*5.2. Example 2: Homogeneous Cylindrical Nanoshells Based on the MCST*

In Table 2, the comparison study is conducted for natural frequency Ω = *ωR ρ*/*E* of a homogeneous nanoscale cylindrical shell with a simply supported boundary condition by using the MCST. The adopted material properties are: *E* = 1.06 TPa, *ν* = 0.3, *ρ* = 2300 kg/m3. It is observed that the obtained results have a reasonable accordance with those reported [100].


**Table 2.** Comparison of dimensionless natural frequency Ω for a homogeneous cylindrical nanoshell (*R* = 2.32 nm and *L*/*R* = 5).

#### *5.3. Example 3: FG Cylindrical Shell*

Herein, a comparison study is conducted for a simply supported FG cylindrical shell without considering the size effect, as given in Table 3. The FG shell is made of the mixture of Stainless Steel (SS) and Nickel (Ni) with the following material parameters: *E*SS = 207.788 GPa, *ρ*SS = 8166 kg/m<sup>3</sup> and *ν*SS = 0.317756 for SS, and *E*Ni = 205.098 GPa, *ρ*Ni = 8900 kg/m<sup>3</sup> and *ν*NI = 0.31 for Ni. Our study yields an excellent agreement with Reference [101], bespeaking the correctness of the current research.

**Table 3.** Comparison of natural frequencies (Hz) for a simply supported FG cylindrical shell (*n* = 1, *R* = 1 m and *L*/*R* = 20).


#### **6. Results and Discussion**

In this section, size-dependent free vibration and axial buckling of an FG NPMF nanoshell simply supported at both ends are studied. The material properties of the nanoshell are *E*∗ <sup>1</sup> = 200 GPa, *ρ*∗ <sup>1</sup> <sup>=</sup> 7850 kg/m<sup>3</sup> , and *ν* = 1/3. The dimensionless natural frequency is defined as Ω = *ωR ρ*∗ 1/*E*<sup>∗</sup> 1 and the dimensionless buckling load is *F*= *F*/*A*110, where *A*<sup>110</sup> is the specific value of *A*<sup>11</sup> for the homogeneous nanoshell made of solid metal.

#### *6.1. Free Vibration Analysis*

Table 4 shows the variation of dimensionless natural frequency with the circumferential wave number for various length scale parameters. It is found that by increasing the dimensionless length scale parameter, the natural frequencies of the system decrease. Moreover, the fundamental natural frequency occurs at *n* = 2, independent of the length scale parameter. In the following studies, the mode (1, 2) is chosen as a representative mode.

**Table 4.** Effect of length scale parameter on dimensionless natural frequencies Ω based on the MSGT (nanoporosity-1, *m* = 1, *h* = 10 nm, *R* = 20*h*, *L*/*R* = 4, *e*<sup>0</sup> = 0.5).


The dimensionless natural frequency versus nanoporosity coefficient for different theories and nanoporosity distributions is illustrated in Figure 4. Results show that the natural frequency decreases by increasing the nanoporosity coefficient, indicating that the nano-pores decrease the effective stiffness of the nanoshell. Furthermore, the nanoporosity-2 nanoshell has a lower natural frequency than its nanoporosity-1 counterpart. It is observed that the natural frequencies predicted by the MCST and MSGT are greater than the natural frequency predicted by the CT. In other words, the additional length scale parameter makes the FG NPMF nanoshell stiffer. This is due to the extra stiffness introduced in the MCST and MSGT.

**Figure 4.** Dimensionless natural frequency versus nanoporosity coefficient with different theories and nanoporosity distributions (*m* = 1, *n* = 2, *h* = 10 nm, *R* = 20 *h*, *h* = 2*l*, *L* = 4*R*).

Depicted in Figure 5 is the variation of the dimensionless natural frequency against the dimensionless length scale parameter. It is seen that the size effect on natural frequency is more pronounced when the thickness of the nanoshell is comparable to the length scale parameter. The dimensionless natural frequencies from the MCST and MSGT converge to the results from the CT for a large value of the dimensionless length scale parameter, indicating that the larger dimensionless length scale parameter diminishes the size effect on the natural frequency of the FG NPMF nanoshell.

**Figure 5.** Dimensionless natural frequency versus dimensionless length scale parameter with different theories and nanoporosity distributions (*m* = 1, *n* = 2, *h* = 10 nm, *R* = 20 *h*, *L* = 4*R*, *e*<sup>0</sup> = 0.5).

Figure 6 plots the dimensionless natural frequency versus length-to-radius ratio with different theories and nanoporosity distributions. One can see that as the length-to-radius ratio increases, the dimensionless natural frequency decreases gradually. Compared to the MCST, the MSGT leads to more reasonable results due to the introduction of an additional deviatoric stretch gradient tensor and the dilatation gradient tensor in addition to the symmetric rotation gradient tensor.

**Figure 6.** Dimensionless natural frequency versus length-to-radius ratio with different theories and nanoporosity distributions (*m* = 1, *n* = 2, *h* = 10 nm, *R* = 20 *h*, *h* = 2*l*, *e*<sup>0</sup> = 0.5).

Figure 7 illustrates the effect of the thickness-to-radius ratio on the dimensionless natural frequency of the FG NPMF nanoshell. As expected, the natural frequency of the FG NPMF nanoshell increases with the rise of thickness-to-radius. This is because the larger thickness-to-radius ratio results in the enhancement of the nanoshell stiffness. Moreover, the difference among the results obtained from the MCST, MSGT, and CT becomes more and more notable as the ratio of thickness-to-radius increases, indicating that the size effect is more significant at the larger thickness-to-radius ratio.

**Figure 7.** Dimensionless natural frequency versus thickness-to-radius ratio with different theories and nanoporosity distributions (*m* = 1, *n* = 2, *h* = 10 nm, *L* = 4*R*, *h* = 2*l*, *e*<sup>0</sup> = 0.5).

#### *6.2. Buckling Analysis*

The effect of the length scale parameter on the dimensionless buckling load is shown in Table 5. It is revealed that by increasing the dimensionless length scale parameter, the buckling load of the system decreases. Additionally, with the increase of circumferential wave number, the buckling load first decreases and then increases. It is noted that the critical buckling load occurs at *n* = 2.

**Table 5.** Effect of the length scale parameter on dimensionless buckling load *F* based on the MSGT (nanoporosity-1, *m* = 1, *h* = 10 nm, *R* = 20*h*, *L*/*R* = 4, *e*<sup>0</sup> = 0.5).


Figure 8 plots the dimensionless critical buckling load versus nanoporosity coefficient for both nanoporosity distributions based on the MCST, MSGT, and CT. As can be observed, the larger nanoporosity coefficient results in a lower dimensionless critical buckling load. Moreover, the nanoporosity-1 nanoshell has a higher critical buckling load than its nanoporosity-2 counterpart. The difference between them tends to be significant with the increase of the nanoporosity coefficient. Furthermore, compared to the MCST, the MSGT leads to a more reasonable buckling load due to the introduction of an additional deviatoric stretch gradient tensor and dilatation gradient tensor in addition to the symmetric rotation gradient tensor.

**Figure 8.** Dimensionless critical buckling load versus nanoporosity coefficient (*m* = 1, *n* = 2, *h* = 10 nm, *R* = 20 *h*, *L* = 4*R*, *h* = 2*l*).

Figure 9 compares the variation of the dimensionless critical buckling load with the dimensionless length scale parameter based on classical and non-classical shell models. It is noted that the critical buckling load decreases with the increasing dimensionless length scale parameter. In addition, the difference among the results from the three models (MCST, MSGT, and CT) is diminishing when the dimensionless length scale parameter tends to large, indicating that the size effect is only significant when the thickness of the nanoshell is comparable to the length scale parameter. This phenomenon was also found in microplates and microbeams [41,47,102].

**Figure 9.** Dimensionless critical buckling load versus dimensionless length scale parameter (*m* = 1, *n* = 2, *h* = 10 nm, *R* = 20 *h*, *L* = 4*R*, *e*<sup>0</sup> = 0.5).

Depicted in Figure 10 is the variation of the dimensionless buckling load with the length-to-radius ratio for both kinds of nanoporosity distribution. It can be seen that with the increase of the length-to-radius ratio, the dimensionless buckling load first decreases and then increases. Moreover, the dimensionless buckling load obtained through the MSGT is greater than those predicted via the CT and MCST. The difference between the results obtained by the MCST, MSGT, and CT becomes more and more significant as the length-to-radius ratio rises.

**Figure 10.** Dimensionless buckling load versus length-to-radius ratio (*m* = 1, *n* = 2, *h* = 10 nm, *R* = 20 *h*, *h* = 2*l*, *e*<sup>0</sup> = 0.5).

Figure 11 plots the dimensionless buckling load with respect to the thickness-to-radius ratio for both kinds of nanoporosity distribution. As can be seen, the increase in the thickness-to-radius ratio contributes to the higher buckling load of the FG NPMF nanoshell. This is due to the fact that the larger thickness-to-radius ratio enhances the stiffness of the FG NPMF nanoshell.

**Figure 11.** Dimensionless buckling load versus thickness-to-radius ratio (*m* = 1, *n* = 2, *h* = 10 nm, *L* = 4*R*, *h* = 2*l*, *e*<sup>0</sup> = 0.5).

#### **7. Conclusions**

In this paper, size-dependent free vibration and buckling of FG NPMF cylindrical nanoshells are investigated based upon the FSD shell theory and general SGT. The symmetric and unsymmetric nanoporosity distributions are considered for the structural composition. Governing equations, as well as corresponding boundary conditions, are derived via Hamilton's principle. Moreover, the Navier solution technique is employed to derive the analytical solutions for FG NPMF nanoshells with a simply supported boundary condition. The conclusions can be summarized as follows:


**Author Contributions:** Conceptualization, Y.Z. and F.Z.; Methodology, F.Z.; Software, F.Z.; Validation, F.Z.; Formal Analysis, F.Z.; Investigation, F.Z.; Resources, Y.Z.; Data Curation, F.Z.; Writing-Original Draft Preparation, F.Z.; Writing-Review & Editing, Y.Z.; Visualization, F.Z.; Supervision, Y.Z.; Project Administration, Y.Z.; Funding Acquisition, Y.Z.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant no. 11672188). **Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

The non-classical and classical resultant moments and forces in Equation (26) are:

$$N\_{\rm xx} = B\_{11}\phi\_1 + A\_{11}\phi\_0 + B\_{12}\phi\_1 + A\_{12}\phi\_0 \tag{A1}$$

$$N\_{yy} = B\_{11}\varrho\_1 + A\_{11}\varrho\_0 + B\_{12}\phi\_1 + A\_{12}\phi\_0 \tag{A2}$$

$$N\_{xy} = B\_{55}k\_1 + A\_{55}k\_0\tag{A3}$$

$$Q\_{\mathbf{x}} = K\_{\mathbb{S}} A\_{\mathbb{S}\mathbf{5}} \gamma\_1, Q\_{\mathbf{y}} = K\_{\mathbb{S}} A\_{\mathbb{S}\mathbf{5}} \gamma\_2 \tag{A4}$$

$$M\_{\rm xx} = D\_{11}\phi\_1 + B\_{11}\phi\_0 + D\_{12}\phi\_1 + B\_{12}\phi\_0 \tag{A5}$$

$$M\_{yy} = D\_{11}\varrho\_1 + B\_{11}\varrho\_0 + D\_{12}\phi\_1 + B\_{12}\phi\_0 \tag{A6}$$

$$M\_{xy} = D\_{55}k\_1 + B\_{55}k\_0\tag{A7}$$

$$T\_{\rm xxx} = E\_1 \frac{\partial \phi\_0}{\partial \mathbf{x}} + E\_2 \frac{\partial \phi\_0}{\partial \mathbf{x}} + \frac{E\_3}{2} \frac{\partial k\_0}{\partial y} + F\_1 \frac{\partial \phi\_1}{\partial \mathbf{x}} + F\_2 \frac{\partial \varphi\_1}{\partial \mathbf{x}} + \frac{F\_3}{2} \frac{\partial k\_1}{\partial y} \tag{A8}$$

$$T\_{yxx} = E\_2 \frac{\partial \varphi\_0}{\partial y} + \frac{1}{2} \left( E\_4 \frac{\partial k\_0}{\partial x} + F\_4 \frac{\partial k\_1}{\partial x} \right) + E\_5 \frac{\partial \phi\_0}{\partial y} + F\_2 \frac{\partial \varphi\_1}{\partial y} + F\_5 \frac{\partial \phi\_1}{\partial y} \tag{A9}$$

$$T\_{zxx} = E\_5 \phi\_1 + \frac{E\_4}{2} \frac{\partial \gamma\_1}{\partial x} + \frac{A\_1}{2} \frac{\partial \gamma\_2}{\partial y} + 2A\_2 \phi\_1 \tag{A10}$$

$$T\_{xyy} = \frac{E\_4}{2} \frac{\partial k\_0}{\partial y} + E\_2 \frac{\partial \phi\_0}{\partial \mathbf{x}} + E\_5 \frac{\partial \varphi\_0}{\partial \mathbf{x}} + \frac{F\_4}{2} \frac{\partial k\_1}{\partial y} + F\_2 \frac{\partial \phi\_1}{\partial \mathbf{x}} + F\_5 \frac{\partial \varphi\_1}{\partial \mathbf{x}} \tag{A11}$$

$$T\_{yyy} = E\_2 \frac{\partial \phi\_0}{\partial y} + \frac{E\_3}{2} \frac{\partial k\_0}{\partial x} + E\_1 \frac{\partial \varphi\_0}{\partial y} + F\_2 \frac{\partial \phi\_1}{\partial y} + \frac{F\_3}{2} \frac{\partial k\_1}{\partial x} + F\_1 \frac{\partial \varphi\_1}{\partial y} \tag{A12}$$

$$T\_{zyy} = \frac{A\_1}{2} \frac{\partial \gamma\_1}{\partial x} + 2A\_2 \phi\_1 + E\_5 \phi\_1 + \frac{E\_4}{2} \frac{\partial \gamma\_2}{\partial y} \tag{A13}$$

$$T\_{xxy} = \frac{E\_3}{2} \frac{\partial \varphi\_0}{\partial y} + \frac{E\_6}{2} \frac{\partial k\_0}{\partial x} + \frac{E\_4}{2} \frac{\partial \phi\_0}{\partial y} + \frac{F\_3}{2} \frac{\partial \varphi\_1}{\partial y} + \frac{F\_6}{2} \frac{\partial k\_1}{\partial x} + \frac{F\_4}{2} \frac{\partial \phi\_1}{\partial y} \tag{A14}$$

$$T\_{yyx} = \frac{E\_4}{2} \frac{\partial \rho\_0}{\partial x} + \frac{E\_3}{2} \frac{\partial \phi\_0}{\partial x} + \frac{E\_6}{2} \frac{\partial k\_0}{\partial y} + \frac{F\_4}{2} \frac{\partial \rho\_1}{\partial x} + \frac{F\_3}{2} \frac{\partial \phi\_1}{\partial x} + \frac{F\_6}{2} \frac{\partial k\_1}{\partial y} \tag{A15}$$

$$T\_{zxy} = A\_4 k\_1 + \frac{A\_5}{2} \left(\frac{\partial \gamma\_2}{\partial x} + \frac{\partial \gamma\_1}{\partial y}\right) \tag{A16}$$

$$T\_{xxz} = \frac{A\_1}{2}\varphi\_1 + \frac{A\_3}{2}\frac{\partial \gamma\_2}{\partial y} + \frac{E\_6}{2}\frac{\partial \gamma\_1}{\partial x} + \frac{E\_4}{2}\phi\_1 \tag{A17}$$

$$T\_{yxz} = A\_4 \frac{\partial \gamma\_1}{\partial y} + \frac{A\_5}{2} \left(\frac{\partial \gamma\_2}{\partial x} + k\_1\right) \tag{A18}$$

$$T\_{xyz} = A\_4 \frac{\partial \gamma\_2}{\partial x} + \frac{A\_5}{2} \left(\frac{\partial \gamma\_1}{\partial y} + k\_1\right) \tag{A19}$$

$$T\_{yyz} = \frac{A\_1}{2}\phi\_1 + \frac{A\_3}{2}\frac{\partial \gamma\_1}{\partial x} + \frac{E\_6}{2}\frac{\partial \gamma\_2}{\partial y} + \frac{E\_4}{2}\varphi\_1 \tag{A20}$$

$$M\_{\rm xxx} = F\_1 \frac{\partial \phi\_0}{\partial x} + F\_2 \frac{\partial \phi\_0}{\partial x} + \frac{F\_3}{2} \frac{\partial k\_0}{\partial y} + G\_1 \frac{\partial \phi\_1}{\partial x} + G\_2 \frac{\partial \varphi\_1}{\partial x} + \frac{G\_3}{2} \frac{\partial k\_1}{\partial y} \tag{A21}$$

$$M\_{yxx} = F\_2 \frac{\partial \varphi\_0}{\partial y} + F\_5 \frac{\partial \phi\_0}{\partial y} + \frac{F\_4}{2} \frac{\partial k\_0}{\partial x} + G\_2 \frac{\partial \varphi\_1}{\partial y} + G\_5 \frac{\partial \phi\_1}{\partial y} + \frac{G\_4}{2} \frac{\partial k\_1}{\partial x} \tag{A22}$$

$$M\_{xyy} = \frac{F\_4}{2} \frac{\partial k\_0}{\partial y} + F\_2 \frac{\partial \phi\_0}{\partial x} + F\_5 \frac{\partial \varphi\_0}{\partial x} + \frac{G\_4}{2} \frac{\partial k\_1}{\partial y} + G\_2 \frac{\partial \phi\_1}{\partial x} + G\_5 \frac{\partial \varphi\_1}{\partial x} \tag{A23}$$

$$M\_{yyy} = F\_2 \frac{\partial \phi\_0}{\partial y} + \frac{F\_3}{2} \frac{\partial k\_0}{\partial x} + F\_1 \frac{\partial \phi\_0}{\partial y} + G\_2 \frac{\partial \phi\_1}{\partial y} + \frac{G\_3}{2} \frac{\partial k\_1}{\partial x} + G\_1 \frac{\partial \phi\_1}{\partial y} \tag{A24}$$

$$M\_{\mathbf{x}xy} = \frac{F\_3}{2} \frac{\partial \varphi\_0}{\partial y} + \frac{F\_6}{2} \frac{\partial k\_0}{\partial \mathbf{x}} + \frac{F\_4}{2} \frac{\partial \phi\_0}{\partial y} + \frac{G\_3}{2} \frac{\partial \varphi\_1}{\partial y} + \frac{G\_6}{2} \frac{\partial k\_1}{\partial \mathbf{x}} + \frac{G\_4}{2} \frac{\partial \phi\_1}{\partial y} \tag{A25}$$

*Nanomaterials* **2019**, *9*, 271

$$M\_{yyx} = \frac{F\_4}{2} \frac{\partial \varphi\_0}{\partial x} + \frac{F\_3}{2} \frac{\partial \phi\_0}{\partial x} + \frac{F\_6}{2} \frac{\partial k\_0}{\partial y} + \frac{G\_4}{2} \frac{\partial \varphi\_1}{\partial x} + \frac{G\_3}{2} \frac{\partial \phi\_1}{\partial x} + \frac{G\_6}{2} \frac{\partial k\_1}{\partial y} \tag{A26}$$

#### **Appendix B**

The parameters in Equations (42)–(46) are given by:

$$A\_{11} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} (\lambda + 2\mu) \mathrm{d}z , B\_{11} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} (\lambda + 2\mu) z \mathrm{d}z , \; D\_{11} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} (\lambda + 2\mu) z^2 \mathrm{d}z \tag{A27}$$

$$A\_{12} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \lambda \,\mathrm{d}z, B\_{12} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \lambda z \,\mathrm{d}z,\\ D\_{12} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \lambda z^2 \,\mathrm{d}z \tag{A28}$$

$$A\_{55} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \mu \mathrm{d}z \,\, \mathrm{d}z \,\, \mathrm{B}\_{55} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \mu z \mathrm{d}z \,\mathrm{d}z \,\mathrm{D}\_{55} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \mu z^2 \mathrm{d}z \,\tag{A29}$$

$$A\_1 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} a\_1 \,\mathrm{d}z \,\,\, A\_2 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} a\_2 \,\mathrm{d}z \,\,\,\, A\_3 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} a\_3 \,\mathrm{d}z \,\,\,\tag{A30}$$

$$A\_4 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} a\_4 \mathrm{d}z ,\\ A\_5 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} a\_5 \mathrm{d}z ,\\ B\_1 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} a\_1 z \mathrm{d}z \tag{A31}$$

$$B\_2 = \int\_{-\frac{h}{2}}^{\frac{h}{2}} a\_2 z \,\mathrm{d}z \,\, B\_3 = \int\_{-\frac{h}{2}}^{\frac{h}{2}} a\_3 z \,\mathrm{d}z \,\, B\_4 = \int\_{-\frac{h}{2}}^{\frac{h}{2}} a\_4 z \,\mathrm{d}z \tag{A32}$$

$$E\_1 = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \beta\_1 \mathrm{d}z,\\ E\_2 = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \beta\_2 \mathrm{d}z,\\ E\_3 = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \beta\_3 \mathrm{d}z \tag{A33}$$

$$E\_4 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_4 \mathrm{d}z, E\_5 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_5 \mathrm{d}z, E\_6 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_6 \mathrm{d}z \tag{A34}$$

$$F\_1 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_1 z \mathrm{d}z , F\_2 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_2 z \mathrm{d}z ,\ F\_3 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_3 z \mathrm{d}z \tag{A35}$$

$$F\_4 = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \beta\_4 z dz,\\ F\_5 = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \beta\_5 z dz,\\ F\_6 = \int\_{-\frac{h}{2}}^{\frac{h}{2}} \beta\_6 z dz \tag{A36}$$

$$\mathbf{G}\_1 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_1 z^2 \, \mathrm{d}z \, \mathrm{d}z = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_2 z^2 \, \mathrm{d}z \, \mathrm{d}z = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_3 z^2 \, \mathrm{d}z \tag{A37}$$

$$\mathbf{G}\_4 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_4 z^2 \, \mathrm{d}z \, \mathrm{d}z \, \mathrm{G}\_5 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_5 z^2 \, \mathrm{d}z \, \mathrm{G}\_6 = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \beta\_6 z^2 \, \mathrm{d}z \tag{A38}$$

#### **Appendix C**

The nonzero components in Equations (53)–(55) are:

$$K\_{11} = A\_{11}a\_m^2 + E\_1a\_m^4 + \left[A\_{55} + (E\_3 + E\_4 + E\_5 + E\_6)a\_m^2\right] \frac{n^2}{R^2} + \frac{E\_6n^4}{2R^4} \tag{A39}$$

$$K\_{12} = K\_{21} = \frac{a\_m n}{R} \left[ A\_{12} + A\_{55} + \left( \frac{2E\_2 + E\_3 + E\_4 + E\_6}{2} \right) \left( a\_m^2 + \frac{n^2}{R^2} \right) \right] \tag{A40}$$

$$K\_{13} = K\_{31} = -\frac{(E\_3 + E\_4)a\_m n^2}{2R^3} - \frac{E\_2 a\_m \left(a\_m^2 + \frac{n^2}{R^2}\right)}{R} - \frac{A\_{12} a\_m}{R} \tag{A41}$$

$$K\_{14} = K\_{41} = F\_1 a\_m^4 + \frac{F\_6}{2} \frac{n^4}{R^4} + \left(\frac{2F\_3 + 2F\_4 + 2F\_5 + F\_6}{2}\right) \frac{a\_m^2 n^2}{R^2} + B\_{11} a\_m^2 + \frac{B\_{55} n^2}{R^2} \tag{A42}$$

$$K\_{15} = K\_{51} = \left(\frac{2F\_2 + F\_3 + F\_4 + F\_6}{2}\right) \left(\frac{a\_m^3 n}{R} + \frac{a\_m n^3}{R^3}\right) + (B\_{12} + B\_{55})\frac{a\_m n}{R} \tag{A43}$$

$$\begin{array}{ll} \text{K2} = & \frac{E\_6 a\_m^4}{2} + \left( \frac{E\_6 + 2E\_3 + 2E\_4 + 2E\_5}{2} \right) \frac{a\_m^2 n^2}{R^2} + \frac{E\_1 n^4}{R^4} \\ & + \left( A\_{11} + \frac{E\_6}{R^2} \right) \frac{n^2}{R^2} + a\_m^2 \left( A\_{55} + \frac{A\_4}{R^2} \right) + \frac{K\_8 A\_{55}}{R^2} \end{array} \tag{A44}$$

$$\begin{array}{rcl} K\_{23} = K\_{32} = & -\frac{(2A\_3 + 2A\_4 + A\_5 + E\_3 + E\_4 + 2E\_5)a\_m^2 n}{2R^2} \\ & -\frac{(E\_1 + E\_6)n^3}{R^4} - \frac{\eta}{R^2}(A\_{11} + K\_S A\_{55}) \end{array} \tag{A45}$$

$$\begin{array}{l} K\_{24} = K\_{42} = \begin{array}{c} -\frac{a\_{\text{n}}u}{R^2} \left( A\_1 + A\_3 + A\_5 \right) \\ + \left[ \left( \frac{2F\_2 + F\_3 + F\_4 + F\_6}{2} \right) \left( \frac{a\_{\text{n}}^3 u}{R} + \frac{a\_{\text{n}} u^3}{R^3} \right) + \left( B\_{12} + B\_{55} \right) \frac{a\_{\text{n}} u}{R} \right] \end{array} \tag{A46}$$

$$\begin{array}{lcl} \text{K}\_{25} &= \text{K}\_{52} = -\frac{1}{2\overline{\mathcal{R}}} \Big[ (2A\_4 + A\_5)a\_m^2 + (E\_4 + E\_6)\frac{n^2}{\overline{\mathcal{R}}^2} \right] - \frac{\text{K}\_5 A\_{55}}{\overline{\mathcal{R}}}\\ &+ \left[ \frac{F\_6 a\_m^4}{2} + \frac{F\_1 n^4}{\overline{\mathcal{R}}^4} + \left( \frac{2F\_3 + 2F\_4 + 2F\_5 + F\_6}{2} \right) \frac{a\_m^2 n^2}{\overline{\mathcal{R}}^2} \right] + B\_{55} a\_m^2 + \frac{B\_{11} n^2}{R^2} \end{array} \tag{A47}$$

$$\begin{array}{l} K\_{33} = & (A\_3 + 2A\_4 + A\_5) \frac{a\_m^2 n^2}{R^2} + \frac{(2E\_1 + E\_6)n^2}{2R^4} + \frac{(2E\_5 + A\_3)a\_m^2}{2R^2} \\ & + K\_S A\_{55} \left(a\_m^2 + \frac{n^2}{R^2}\right) + \frac{E\_6}{2} \left(a\_m^4 + \frac{n^4}{R^4}\right) + \frac{A\_{11}}{R^2} \end{array} \tag{A48}$$

$$\begin{split} K\_{34} = K\_{43} &= \begin{array}{c} \frac{a\_m(A\_1 + A\_3)}{2R^2} + a\_m \left( K\_S A\_{55} - \frac{2F\_2 + F\_3 + F\_4}{2R} \frac{n^2}{R^2} + \frac{E\_6 + E\_4}{2} a\_m^2 \right) \\ + \frac{1}{2}(A\_1 + A\_3 + 2A\_4 + 3A\_5) \frac{a\_m n^2}{R^2} + \frac{a\_m^3 F\_2}{R} - \frac{B\_{12} a\_m}{R} \end{array} \tag{A49}$$

$$\begin{array}{ll} \mathbb{K}\_{55} = & \mathbb{K}\_{53} = \frac{1}{2} \Big[ (A\_1 + A\_3 + 2A\_4 + 3A\_5)a\_m^2 + (E\_4 + E\_6)\frac{\mu^2}{R^2} \Big] \frac{\mu}{R} \\ & -\frac{E\_3 - E\_4 - 2F\_5}{2R} \frac{a\_m^2 n}{R} - \frac{F\_1}{R} \frac{n^3}{R^3} - \left(\frac{B\_{11}}{R} - K\_S A\_{55}\right) \frac{\mu}{R} \end{array} \tag{A50}$$

$$\begin{array}{ll} \mathbb{K}\_{44} = & \mathbb{G}\_1 a\_m^4 + \left[ \left( \frac{2\mathbb{G}\_3 + 2\mathbb{G}\_4 + 2\mathbb{G}\_5 + \mathbb{G}\_6}{2} \right) \frac{n^2}{R^2} + D\_{11} + \frac{2E\_4 + 2E\_5 + E\_6}{2} \right] a\_m^2 \\ & + \frac{\mathbb{G}\_5 n^4}{2R^4} + \frac{n^2}{R^2} (D\_{55} + 2A\_4 + A\_5) + K\_S A\_{55} \end{array} \tag{A51}$$

$$\begin{array}{ll} K\_{45} = K\_{54} = & \frac{a\_{\overline{m}}}{\mathcal{R}} \left[ \frac{2G\_2 + G\_3 + G\_4 + G\_6}{2} \left( a\_m^2 + \frac{n^2}{\mathcal{R}^2} \right) \right] \\ & + \frac{a\_{\overline{m}}}{\mathcal{R}} \left( D\_{12} + D\_{55} + \frac{2A\_1 + 4A\_2 + A\_3 + 2A\_4 + 3A\_5}{2} \right) \end{array} \tag{A52}$$

$$\begin{array}{ll} K\_{55} = & \frac{G\_1 \mu^4}{R^4} + \frac{a\_m^4 G\_6}{2} + \left(\frac{2G\_3 + 2G\_4 + 2G\_5 + G\_6}{2}\right) \frac{a\_m^2 \mu^2}{R^2} \\ & + \frac{n^2}{R^2} \left(D\_{11} + \frac{2E\_4 + 2E\_5 + E\_6}{2}\right) \\ \end{array} \tag{A53}$$

$$K\_{\text{\\$33}} = -\frac{m^2 \pi^2}{L^2} \tag{A54}$$

$$M\_{11} = M\_{22} = M\_{33} = I\_0,\\ M\_{44} = M\_{55} = I\_2,\\ M\_{14} = M\_{25} = M\_{41} = M\_{52} = I\_1 \tag{A55}$$

#### **References**


#### *Nanomaterials* **2019**, *9*, 271


© 2019 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/).

## *Article* **Structural and Electrical Studies for Birnessite-Type Materials Synthesized by Solid-State Reactions**

**Nayda P. Arias 1,2, María E. Becerra 1,3,4,5,6 and Oscar Giraldo 1,3,5,\***


Received: 14 July 2019; Accepted: 8 August 2019; Published: 12 August 2019

**Abstract:** The focus of this paper is centered on the thermal reduction of KMnO4 at controlled temperatures of 400 and 800 ◦C. The materials under study were characterized by atomic absorption spectroscopy, thermogravimetric analysis, average oxidation state of manganese, nitrogen adsorption–desorption, and impedance spectroscopy. The structural formulas, found as a result of these analyses, were *K*<sup>+</sup> 0.29- *Mn*4<sup>+</sup> 0.84*Mn*3<sup>+</sup> 0.16 *<sup>O</sup>*2.07·0.61*H*2*<sup>O</sup>* and *<sup>K</sup>*<sup>+</sup> 0.48- *Mn*4<sup>+</sup> 0.64*Mn*3<sup>+</sup> 0.36 *O*2.06·0.50*H*2*O*. The N2 adsorption–desorption isotherms show the microporous and mesoporous nature of the structure. Structural analysis showed that synthesis temperature affects the crystal size and symmetry, varying their electrical properties. Impedance spectroscopy (IS) was used to measure the electrical properties of these materials. The measurements attained, as a result of IS, show that these materials have both electronic and ionic conductivity. The conductivity values obtained at 10 Hz were 4.1250 <sup>×</sup> <sup>10</sup>−<sup>6</sup> and 1.6870 <sup>×</sup> <sup>10</sup>−<sup>4</sup> <sup>Ω</sup>−1cm−<sup>1</sup> for Mn4 at 298 and 423 K respectively. For Mn8, the conductivity values at this frequency were 3.7074 <sup>×</sup> 10−<sup>7</sup> (298) and 3.9866 <sup>×</sup> 10−<sup>5</sup> <sup>Ω</sup>−1cm−<sup>1</sup> (423 K). The electrical behavior was associated with electron hopping at high frequencies, and protonic conduction and ionic movement of the K<sup>+</sup> species, in the interlayer region at low frequencies.

**Keywords:** Birnessite; nanoporous metal oxides; impedance spectroscopy

#### **1. Introduction**

From the 1960s, works of Wolkestein [1] state the importance of electron theory to elucidate the relation between the catalytic and electronic properties of a catalyst and also the semiconductor nature of this type of material. Swaminathan [2] produced research about the opportunities and ways to improve catalyst technology, specifically in their development, cost reduction, and field applications [2]. Among the catalysts, semiconductors such as TiO2, ZnO, and SnO2 play a pivotal role, as they can use the electromagnetic spectra to degrade contaminants [3,4]. It has been reported that these semiconductor oxides also have acidic and basic natures [2]. Therefore, they are useful as both solid acid and base catalysts. For this reason, efforts to develop new, cheaper, active, and selective catalysts are ongoing. Manganese oxides are among the most extensively studied raw, supported, or doped

catalysts [5–11]; however, their electrical properties have been less studied [12–14]. The family of the manganese oxide type of materials have variations in their structure depending on the way the Mn–O is linked together [15]. Therefore, layer, tunnel, spinel, and compact structures can be found. Manganese oxides have also been used as a cathode in secondary batteries [16,17]. Manganese, as a central atom in the octahedral coordination in these structures, mainly has the 4+ and 3+ oxidation states. However, average oxidations states of 3.5 and 3.8 are commonly found because of the presence of Mn4<sup>+</sup> and Mn3<sup>+</sup> in the same building blocks [18]. Birnessite is a special structure of this manganese oxide family. It is composed of Mn–O octahedra forming octahedral layers as clays, and it has monovalent or divalent ions surrounded by water molecules to compensate the electrical charge of its layers. Therefore, these monovalent or divalent cations are located in the region called the interlayer [9,12]. The cations in the interlayer region can be displaced by other species because of their mobile nature. Birnessite has been used in previous works as a catalyst in soot combustion processes and in methylene blue degradation, showing appreciable catalytic activity compared to traditional catalysts [9,11,19]. For this reason, the present work focuses on the charge transport mechanism for two birnessite material types, Mn4 and Mn8, synthesized at 400 and 800 ◦C, respectively, with the intention for a deeper understanding of the nature of this type of material for advanced applications.

#### **2. Results and Discussion**

#### *2.1. Chemical Composition, Thermogravimetric Analysis (TGA), and Average Oxidation State (AOS)*

The thermal reduction of KMnO4 has been reported by Kappestein [20] and Herbstein [21]. In the study conducted by Herbstein et al. [21] on the thermal decomposition of KMnO4 in a temperature range of 25 to 900 ◦C, in an atmosphere of air and nitrogen, it was found that the idealized equation for the decomposition of KMnO4 in air at 250 ◦C results in a soluble phase, the K2MnO4, and another phase that is insoluble, like this:

$$10\text{ }\text{KMnO}\_4 \text{ }\text{\textdegree 2.65 }\text{K}\_2\text{MnO}\_4 + [2.35\text{K}\_2\text{O}\text{ }\text{\textdegree 7.75}\text{MnO}\_2\text{\textdegree]} + 6\text{O}\_2\text{.}$$

Decomposition in air or nitrogen at higher temperatures (up to 540 ◦C) results in more considerable amounts of O2, corresponding to the change in the composition of the nonsoluble phase. At temperatures above 540 ◦C in both air and nitrogen, K2MnO4 decomposes into:

$$10\,\mathrm{K}\_{2}\mathrm{MnO}\_{4} \rightarrow 5.7\mathrm{K}\_{3}\mathrm{MnO}\_{4} + 0.5[2.9\,\mathrm{K}\_{2}\mathrm{O}\,\,8.6\,\mathrm{MnO}\_{2.1}] + 3.4\mathrm{O}\_{2}\mathrm{s}$$

These potassium manganates are stable within the temperature range of 25 to 900 ◦C, but they react quickly with water vapor. The thermal decomposition of KMnO4 is a redox reaction in which the oxoanion oxygen is oxidized to molecular oxygen, while the oxidation number of manganese in the oxoanion is reduced from Mn (VII) to Mn (VI), from *MnO*− <sup>4</sup> and *MnO*2<sup>−</sup> <sup>4</sup> . It has been reported that in the thermal decomposition of KMnO4 at 250 ◦C, an average of 1.6 Mn–O bonds are broken in 73% of the permanganate ions. The electrons are then transferred to the remaining 27% of permanganate ions, which are reduced to *MnO*2<sup>−</sup> <sup>4</sup> without presenting a significant change in its tetrahedral form. In the K2MnO4 at 600 ◦C, 1.6 Mn–O bonds are broken on average in 43% of *MnO*2<sup>−</sup> <sup>4</sup> ions, and the electrons are transferred to the remaining 57% of manganate ions, which are reduced to *MnO*3<sup>−</sup> 4 .

The manganates formed by this reaction are soluble salts. For this reason, the preparation of birnessite-type materials for this synthesis route involves a washing procedure after the calcination process to provide birnessite as the only crystallographic phase.

The K/Mn ratio, obtained by atomic absorption (AA) results (Table 1), shows a variation in the potassium and manganese content of the materials as the synthesis temperature increases [9].


**Table 1.** Structural formulas for birnessite-type materials.

Similarly, a decrease in the average oxidation state (AOS) (Table 1) indicates the presence of manganese in the 4+ and 3+ oxidation states. The percentage of Mn3<sup>+</sup> in these materials was 16%, and 36% for Mn4 and Mn8, respectively. The results of the AOS, for Mn4 and Mn8, are close to those reported in the literature [9,22] for birnessite produced by solid-state reactions at 400 and 800 ◦C. Thermal analysis in nitrogen and air atmospheres (Figure 1) showed thermal events typical of these layer materials [9,22,23]. Weight loss, up to 250 ◦C, is attributed to the evaporation of physiosorbed water and interlayer water. From 250 to about 800 ◦C, the materials lose oxygen because of the phase transformation from a layer to a tunnel structure [9]. When heated to temperatures above 250 ◦C, the Mn4 material loses 2.48% of its original mass in nitrogen atmosphere. Similarly, Mn8 experienced a weight loss of 1.69%. These results indicate that the materials synthesized at higher temperatures lose less oxygen, and show a greater thermal stability, compared to the ones synthesized by chemical routes [24]. The materials were more stable in nitrogen than in air, as it can be seen in Figure 1. The difference in weight losses at temperatures higher than 250 ◦C between both atmospheres indicates a release of oxygen from the framework, as it was reported by Suib et al. [25]. For Mn4, this difference was 0.53%, and it was 0.16% for Mn8. The combined results of atomic absorption, thermogravimetric analysis (TGA), and AOS were used to formulate the structural formula shown in Table 1.

**Figure 1.** Thermograms in air and N2 for synthesized birnessite-type materials. (**a**) birnessite type of material synthezised at 400◦C, (**b**) birnessite type of material synthezised at 800◦C.

#### *2.2. Structural Analysis and Rietveld Refinement*

The synthesis method produced a typical layered manganese oxide (birnessite) structure, as confirmed by X-ray diffraction (XRD) (Figure 2a). The d (001) basal spacing and estimated crystal size are shown in Table 1. The basal spacing of these layer materials is in accordance with the published works [18,22,26,27].

<sup>\*</sup> By Debye–Scherrer Equation, \*\* 21 ◦C.

**Figure 2.** Rietveld refinement results for (**a**) Mn4, and (**b**) Mn8. The insert shows the polyhedral structure of the triclinic and hexagonal birnessite for (**a**) Mn4, and (**b**) Mn8 respectively.

Analysis between 30◦ and 70◦ in 2θ showed the presence of structural disorders associated with the rotation and translational stacking faults of the layers, as demonstrated in phyllosilicate minerals [28]. The synthesis temperature and structural disorder ratio was also noted by Gaillot et al. [22,29]. The Mn4 material (Figure 2a), in the range mentioned above, had diffraction peaks located at 2.91, 2.47, 2.40, 2.29, 2.17, 1.81, 1.43, 1.39, and 1.35 Å, which are in agreement with those reported for triclinic birnessite [29,30]. For Mn8 it was observed that X-rays diffracted at angles of 12.42◦, 25.06◦, 36.48◦, 38.54◦, 41.36◦, 44.48◦, 44.80◦, and 53.46◦ 2θ, which correspond to interplanar spacing of 7.12, 3.55, 2.46, 2.33, 2.18, 2.03, 2.02, and 1.71 Å, respectively. These are characteristic of a hexagonal birnessite [29,30], an observation that also appears in Kim et al. [31].

The atomic position (fractional coordinates) can be shown in Appendix A, Table A1.

The structural symmetry of birnessite was confirmed by the Rietveld method. The structural parameters have been summarized in Table 2, and comparisons between the experimental and theoretical patterns are shown in Figure 2. The refined parameters confirm the triclinic structure for Mn4 and the hexagonal structure for Mn8 (Table 2). For Mn8 material, the crystallographic density value was higher than Mn4 (Table 2). The change in the crystallographic density for Mn8 was related to the variation in the crystal symmetry from triclinic to hexagonal. The c-cell parameter increased because of the rearrangement of the atoms in the interlayer region and the densification that resulted from temperature increase.

A typical Mn-O1 bond distance, in the layer of the analyzed materials, falls between 1.088–2.352 Å and 1.860 Å for Mn4 and Mn8 materials, respectively. The presence of Mn3<sup>+</sup>, in both the high- and low-spin states, as was evidenced by AOS (Table 1), was confirmed by the four short and two long Mn–O distances, which resulted from Jahn–Teller distortion. The Jahn–Teller distortion has important implications for the electrical conductivity of these materials. The Mn–O bond distances are in close agreement with those reported by Drits et al. [32], Gaillot et al. [22], Lopano et al. [33], and Post and Veblen et al. [34]. The crystal broadening, for Mn4, is isotropic, as can be seen in the aspect ratio index (Table 2), whereas it was anisotropic for Mn8. The crystal anisotropy affects the electrical properties, as will be discussed later.


**Table 2.** Summary of Rietveld refinement results for Mn4 and Mn8 materials.

Mn1: Mn layer, O1: O layer, and O2: O from interlayer water.

#### *2.3. Morphological Analysis*

Scanning electron microscopy analysis (Figure 3) shows that Mn4 (Figure 3a) was formed by particle aggregates small in size compared to Mn8 (Figure 3b). When subjected to calcination at 800 ◦C (Mn8), these agglomerates grew uniformly. These agglomerates tended to have a hexagonal plate morphology and had smoother surfaces, which is characteristic of these types of materials. These observations are consistent with the crystal size/temperature increment relation showed by XRD (Figure 2, Table 1) and Rietveld refinement (Table 2).

**Figure 3.** SEM images for Birnessite-type materials (**a**) Mn4 and (**b**) Mn8.

#### *2.4. N2 Adsorption–Desorption Analysis*

The N2 adsorption–desorption isotherms of the materials (Figure 4a) showed the existence of a type II isotherm, according to the International Union of Pure and Applied Chemistry (IUPAC) classification [35] and H3 hysteresis, characteristic of materials with slit-shaped and large pores [36]. The surface area, found by using the Brunauer-Emmet-Teller (BET) method (Table 1), decreased as the synthesis temperature increased. This trend is consistent with the crystal growth–smoothness relation, as evidenced by XRD (Figure 2) and SEM (Figure 3). The mesopore size distribution (Figure 4b) showed peaks centered at around 399.3 and 414.7 Å for Mn4 and Mn8, respectively, with a pore volume that decreased as the synthesis temperature increased. This observation is in accordance with the aggregate size increment of the birnessite particles. All materials showed microporosity, as observed in Figure 4c, and the data show that the microporosity changes with increasing synthesis temperature. The interaction energy of the porous solid surface with the N2, obtained by nonlocal density functional theory (NLDFT) [37], varied with the synthesis temperature for the birnessite-type material. The data presented in Figure 4d show values for the interaction energy around 48, 56, 60, 72, and 100 K, with a greater contribution at 72 K for the Mn4 material. The main contribution for Mn8 occurred at 28 K followed by interaction energies of 38, 70, and 100 K. Maddox et al. [38] reported that in mesoporous materials, like MCM41, values below 75 K correspond to the usual range of values reported for O–N2 interactions. Values between 125 and 75 K correspond to theoretical εsf/k values for the O–O interaction [38]. Therefore, the information presented asserts that the surfaces of the birnessite-type material are energetically heterogeneous, which probably contributes, in a positive manner, to its electrical processes.

**Figure 4.** N2 adsorption–desorption results: (**a**) N2 isotherms, (**b**) mesopore size distribution by Barret–Joyner–Halenda (BJH), (**c**) micropore size distribution by Hovarth–Kawazoe (HK), and (**d**) interaction energy by nonlocal density functional theory (NLDFT).

#### *2.5. Electrical Analysis of Nyquist Plots*

Nyquist plots from the impedance experimental data were conducted with the aim of studying the effect of temperature on the electrical response of the synthetized materials; the results are shown in Figure 5. The Nyquist plot exhibits a semicircle and a straight line at low frequency, but the magnitude of the diameter of the semicircle reduced as the temperature increased in both studied materials. This experimental evidence suggests an impedance reduction of the "bulk" and grain boundary processes. The presence of the semicircle in the Nyquist plot was associated with the electron transport processes on the "bulk", in which there are a contribution of inter-crystalline and intra-crystalline charge transport processes. Also, in this region, a possible additional contribution to the electrical response came from protonic conduction (H3O<sup>+</sup>) because of the interlayer water in the birnessite structure. The observations corresponding to the electron transport process are concurrent with the results of previous research [12,13]. However, the preparation methods reported here produced materials with less impedance compared to previous work [12]. Therefore, is possible to tune the electrical properties of manganese oxides. The straight line at low frequency was associated with the ionic diffusion of K<sup>+</sup> ions located in the interlayer region, which follow a Warburg diffusion process. The straight lines observed in Figure 5a,b indicated that the ionic diffusion was preserved up to 373 K, due to the satisfactory thermal stability of these materials (Figure 1), as opposed to soft synthesis routes [12,23].

**Figure 5.** Nyquist Plot as a function of temperature for (**a**) Mn4 and (**b**) Mn8. (**c**) equivalent circuit. R1 and R2: resistances; CPE-1 and CPE-2: constant phase element; and Wo1: Warburg open element. The insert in (**a**) corresponds to an increase in the scale between 0 to 8 kΩ, The insert in (**b**) corresponds to an increase in the scale between 0 to 50 kΩ.

The equivalent circuit, shown in Figure 5, that represents the experimental results, was constructed considering the morphology and structural issues for the studied materials, and it also takes into account the reports of the literature for the electrical behavior of porous materials [12,13,39–42]. In this circuit, each circuital element is composed of a constant phase element (CPE) connected in parallel with one resistor (R). CPE1 and R1 reflect the intra-/inter-crystalline electrical process, while CPE2 and R2 were assigned to fit the results concerning the charge transport in the grain boundary or aggregate boundary. These two circuital elements were put in series because it was the better circuit that better adjust the experimental results. Also, each of the circuital elements were connected in series because the Nyquist plot described consecutive processes [43], as reported in previous work [12]. The circuit was finally completed with a Warburg open element (W1). In these equivalent circuits, the constant phase elements (CPE1, CPE2) were selected based on SEM (Figure 3) and textural analysis (Figure 4). CPE describes a nonideal capacitive process, [12,44] which, for the purposes of this study, was used to model the porous structure of the material and its heterogeneous surface (Figure 4a–d) [45,46]. The numerical results of the equivalent circuits are summarized in Appendix A (Tables A2 and A3).

In the circuital model, R1 and R2 are the inter-/intra-crystalline resistance and grain boundary resistance, respectively; CPE1-T is the capacitance of inter-/intra-crystalline processes; CPE2-T corresponds to the capacitance of the grain boundary process; and CPE1-P and CPE2-P are the exponents of the impedance for CPE (Equation (1)) for the inter-/intra-crystalline and grain boundary processes.

$$Z\_{\rm CPE} = \frac{1}{\left(T(I \times \omega)^P\right)}\tag{1}$$

where *ZCPE* is the impedance of the CPE element, *T* corresponds to CPE-T, *P* corresponds to CPE-P, <sup>ω</sup> <sup>=</sup> angular frequency of alternating current signal, and *<sup>I</sup>* <sup>=</sup> <sup>√</sup> −1.

In the same model, Wo-R is the resistance to the ionic diffusion, and Wo-T and Wo-P are the coefficients according to the Warburg impedance equation (Equation (2)):

$$Z\_{Wo} = \frac{R\left(\text{ctg} h \left(I \times T \times \omega\right)^{P}\right)}{\left(I \times T \times \omega\right)^{P}},\tag{2}$$

where *Zwo* is the impedance for the Warburg open element, *R* describes the diffusion resistance (Ω), *P* is the exponential coefficient, and *T* is described through Equation (3):

$$\mathcal{W}\_0 - T = \frac{L^2}{D'} \tag{3}$$

where *W*<sup>0</sup> <sup>−</sup> *T* is the Warburg *T* coefficient, *L* is the effective diffusion thickness (m2), and *D* (m2s−1) is the effective diffusion coefficient of the charge carrier. The contribution of the charge transport processes, and their accumulation at the grain boundaries of the Birnessite particle's aggregates, overlap with the "bulk" contribution because of the proximity of the particles, which was produced by the sample preparation for the electrical experiments (Figure 3). Because of this reason, the CPE2 and R2 elements were introduced into the circuital model.

The increase of the resistance of the "bulk" process (intra-crystalline and inter-crystalline), in the order Mn4 < Mn8 (Tables A2 and A3, Appendix A), suggests an interplay between the crystal size, the average oxidation state of manganese, and the crystal symmetry (Tables 1 and 2). Therefore, the smaller the crystal size is, the higher the micropore and surface area. For instance, Mn4, compared to the other cases, allows the electrons to travel through the crystal more freely because of the shorter Mn–O bond distances, as discussed above.

The values for the relaxation time are presented in Figure 6. The increase in relaxation time with temperature is consistent with the decrease of the resistance at "bulk process", most likely due to higher K–O bond distances in Mn8, as was observed by the Rietveld refinement (Table 2). This relaxation time was almost linear with respect to the temperature for Mn8, but it decreased nearly exponentially for Mn4. This suggests that the relaxation time in the "bulk" is thermally activated and is lower for Mn4, indicating a fast electron mobility correlated to a high surface area and low crystal size (Table 1).

#### *2.6. Alternant Current (AC) conductivity*

AC conductivity measurements (Figure 7) unraveled a power law behavior at high frequencies [47], which suggests electron hopping [48,49] and short-range conductivity [50]. The power law was first reported by Jonscher [47] and referred to as the universal dynamic response. It is described as

$$
\sigma\_{\mathfrak{M}} = \sigma\_0 + A w^n,
$$

where σ*ac* is the AC conductivity, σ<sup>0</sup> is the DC conductivity, ω is the angular frequency of the applied electric field, ω = 2πf, and *A* and *n* are the fitting parameters. *n* generally falls between zero and unity.

**Figure 6.** Effect of temperature on bulk relaxation time (τb) for the synthesized birnessite-type material.

**Figure 7.** Real component of the complex conductivity as a function of temperature. (**a**) Mn4, birnessite type of material synthesized at 400◦C, (**b**) Mn8, birnessite type of material synthesized at 400◦C.

Such behavior has been observed in layered manganese oxides [12] and perovskite-type materials [51]. The sudden increment in the conductivity, above 10 Hz for Mn8 material (Figure 7), most likely is due to the high quantity of Mn3<sup>+</sup> (Table 1) compared to Mn4. However, Mn4 had the highest conductivity (Figure 7) because of its smaller crystal size and higher surface area (Table 1). The conductivity values obtained at 10 Hz were 4.1250 <sup>×</sup> <sup>10</sup>−<sup>6</sup> and 1.6870 <sup>×</sup> <sup>10</sup>−<sup>4</sup> <sup>Ω</sup>−1cm−<sup>1</sup> for Mn4 at 298 and 423 K, respectively. For Mn8, the conductivity values at this frequency were 3.7074 <sup>×</sup> 10−<sup>7</sup> (298 K) and 3.9866 <sup>×</sup> <sup>10</sup>−<sup>5</sup> <sup>Ω</sup>−1cm−<sup>1</sup> (423 K).

It is probable that the stacking faults in Mn4 provide crystallographic active sites for the movement of the charge carrier, as was evidenced in the surface energy plots (Figure 4c,d). Hence, it is probable that, at low frequencies, the ionic conduction process is the most favorable for Mn4, as represented by the appearance of the straight line in Figure 5a. Also, in this material, the K–O bond distances, the crystallographic position of K<sup>+</sup>, and the H2O interlayer allows the continuity in the channels for ionic conduction. Differences in the amount of potassium and water, and its position in the crystal cell of the studied materials, influence the electrical behavior in the low-frequency region. At higher contents of K<sup>+</sup> ions, the two positions in the crystal cell, for the water molecules in Mn8, and the proper alignment of the material´s layers limit the ionic conduction routes, such that this contributes to a diminished ionic conductivity at low frequency (Figure 7). These types of materials are ionic exchangers [52]; therefore, the ions in the interlayer space are mobile. Water is an important key for the mobility of the ions in the interlayer space [12,53]. A low amount of water in the interlayer region can diminish the ionic mobility of the ions [53]. Furthermore, the water´s ability to facilitate ion conduction is dependent on the charge density of these ion [12,53,54]. Another reason for the enhanced conductivity is that the stacking faults alter the electrical conduction pathways. It is probable that this heterogeneous surface creates interactions between the electric field and the charge carrier. The variations in surface energies (Figure 4c,d) support the heterogeneity of the material surfaces.

Figure 7 shows the AC conductivity for each material, by around two orders of magnitude, along with the induced temperature of the experiment. In Figure 7a, a low-frequency dispersion (LFD) can be seen for all temperatures. The region II (between 10 Hz up to 1 MHz) is dominated by DC conductivity, and region III (after 1 MHz) is characterized by the "universal Jonscher law" [47]. The results suggest that the charge transport mechanism is temperature-activated. Figure 8 shows the activation energy of each material. It is evident that, for Mn4 material (Figure 8a), the conductivity depends on the frequency and temperature. However, for Mn8 (Figure 8b), at temperatures higher than 348 K and frequencies below 10 MHz, the dependence is mainly associated with the temperature. The conductivity results suggest that the process is dominated by thermally activated "electron hopping" in the range of analyzed temperatures. When the electrons move forward between Mn sites, there is a change in the local structure caused by the change in the oxidation state of the MnO6 units that form the layer. These changes generate lattice vibrations. This phenomenon should be understood as a strong electron–phonon interaction or polaron mechanism [55], which is the other one present in the electrical conductivity process. The activation energies of each material (Figure 8) were lower than those reported for similar materials [12]. This is because of the higher average oxidation state of manganese (Table 1), which means a higher amount of Mn4+, as is reported for manganites obtained from a nickel permanganate precursor [56]. The lower activation energy of Mn4 compared to Mn6 and Mn8 could be associated with the difference in the quantity of Mn3<sup>+</sup> between the synthesized materials, the variation of the Mn–O–Mn bond lengths resulting from crystal symmetry (Table 2), the changes of the surface area that can lead to a redistribution of electron charge [57], "polaron hopping" [58], and the Jahn–Teller distortion of Mn3<sup>+</sup> ions [12]; thus, there are different electric behaviors with respect to the frequency of the applied electric field. Presumably, the increase of the conductivity involved in the transition zone, between Region I and Region II (Figure 8), was due to the release of physiosorbed water on the material's surface [12].

**Figure 8.** Arrhenius plot of the real component of the complex conductivity.(**a**) Mn4, (**b**) Mn8, (**c**) comparison of the activation energy between 298 to 348 K, (**d**) comparison of the activation energy between 373–423 K.

The comparison of conductivity values for manganese oxides and birnessites synthesized in this work is shown in the Table 3.


**Table 3.** Conductivity values for different manganese oxides.


**Table 3.** *Cont*.

#### **3. Materials and Methods**

#### *3.1. Synthesis*

Manganese oxides with birnessite-type structures (K-birnessite) were synthesized by the thermal reduction of a manganese salt, KMnO4 (Merck, Darmstand-Germany, 99%), as detailed in previous studies [9,22,31,67]. The procedure was as follows; the KMnO4 in powder form was homogeneously dispersed in a porcelain dish and exposed to calcination in a muffle furnace at temperatures of 400 and 800 ◦C for 6 h at 10 ◦C/min, to induce the thermal reduction, then they were cooled to room temperature. After calcination, each material was washed with deionized and distilled water to remove any soluble salts (potassium manganates) formed during the process and to reduce the pH from 12 to 9.5. The materials were dried at 60 ◦C for 24 h for further characterization. The materials were named as Mn4 and Mn8 (from thermal reduction at 400 and 800 ◦C, respectively).

#### *3.2. Elemental Analysis and Average Oxidation State of Manganese*

An elemental chemical analysis for the determination of K and Mn was performed by atomic absorption spectroscopy (AAS) in an atomic absorption spectrophotometer, PERKIN ELMER, model 3110, (Perkin Elmer Corporation, Waltham, MA, USA) and by using a standard procedure for sample digestion. The sample digestion process consisted of dissolving 100.0 mg of sample powder in 2.0 mL of 37% HCl and 1.0 mL of distilled and deionized water (DDW). The sample was then heated up to 100 ◦C, for about 30 min, until the solution became transparent and there were no undissolved solid particles. The Mn and K contents were determined at wavelengths of 279.5 and 766.5 nm, respectively. The average oxidation states (AOSs) of the manganese in the samples were determined by potentiometric and colorimetric titration [68]. In the first step, 40.0 mg of the sample was dissolved by heating in 18.5% v/v HCl solution, with the final volume being adjusted to 100 mL. After that, 10 mL of this solution was added to 100 mL of a saturated solution of Na2P2O7, the pH of the solution was adjusted to between 6.5–7.0, and then it was titrated with KMnO4 standard (around 1 <sup>×</sup> 10−<sup>3</sup> M). The final point for this titration was taken when the potential increased to more than 100–150 mV. In the second step, about 40 mg of sample was dissolved in 10.0 mL of ferrous ammonium sulfate (FAS) solution under constant nitrogen flow, and the volume was adjusted up to 100 mL with DDW in order to reduce Mn4<sup>+</sup> and Mn3<sup>+</sup> to Mn2<sup>+</sup>. For the final step, 10.0 mL of this solution was back titrated, with a permanganate standard, until a color change occurred. All measurements were repeated three times in order to obtain the standard deviation.

#### *3.3. X-ray Di*ff*raction (XRD) Data Collection and Structural Refinement*

All X-Ray Diffraction (XRD) patterns of the powder samples were performed at room temperature with a Brag–Brentano focusing geometry in a RIGAKU MINIFLEX II diffractometer, (Rigaku Company, Tokyo Japan), using CuKα radiation at 30 kV and 15 mA with a scan rate of 2◦ (2 θ min<sup>−</sup>1, sampling width of 0.01◦ (2θ), and between 3◦ and 70◦ (2θ). Rietveld refinement [69], for the obtained birnessite-type materials, was performed with a General Structure Analysis System (GSAS) [70] using the EXPGUI interface [71]. An experimental XRD pattern of a silicon powder was refined in order to obtain the instrumental parameter file. Because of the differences in the layer arrangements between each material, the initial structural parameters, triclinic [33] and hexagonal [22] symmetries, were used. The crystallographic information file (CIF) was obtained from Crystallography Open Database [72,73]. All polyhedral crystal structure drawings were made by using Vesta software.

#### *3.4. Thermal Analysis*

The thermal stability of the materials was studied by thermogravimetric analysis (TGA) on a TA instruments, model TGA Q500 thermogravimetric analyzer (TA Instrument, Delawere, DE, USA). Measurements were made on 10.0 mg of sample, using a high-resolution algorithm (sensitivity: 1, resolution: 5) with N2 flow (100 mL min<sup>−</sup>1) at a heating rate of 10 ◦C min<sup>−</sup>1. The measurement range was from 21 up to 800 ◦C.

#### *3.5. Surface Area and Porosity*

N2 adsorption–desorption isotherms were conducted to study the surface area, pore volume, and pore size distribution of the materials. The samples (100.0 mg) were degassed under vacuum at 120 ◦C for 24 h, and the adsorption–desorption isotherms were taken in a MICROMERITICS equipment model ASAP 2020, (Micromeritics, Norcross (Atlanta), Georgia, GA, USA), at 77 K with a pressure range from 1 <sup>×</sup> <sup>10</sup>−<sup>7</sup> <sup>P</sup>/Po to 0.99 P/Po. The specific area was calculated by the Brunauer–Emmett–Teller (BET) method. The Barret–Joyner–Halenda (BJH) method for mesopores and the Hovarth–Kawazoe (HK) method for micropores were used to determine the pore size distribution. To predict the surface energy of the synthesized materials, nonlocal density functional theory (NLDFT) was applied.

#### *3.6. Scanning Electron Microscopy (SEM)*

Micrographs were obtained in a JEOL JSM 5910LV microscope, (JEOL company Akishima, Tokyo, Japan), equipped with a secondary electron (SEI) detector operated at 15 kV in high-vacuum mode. For this analysis, around 3.0 mg of the synthesized material was placed onto a carbon ribbon, and then a thin layer of gold was deposited onto the surface of the sample by a sputter coating technique.

#### *3.7. Impedance Spectroscopy Analysis*

For impedance spectroscopy analysis, expressed as a function of temperature, the powder samples were uniaxially pressed at about 51 MPa to obtain pellets approximately 1.37 mm thick and 10 mm in diameter [12,14]. The pellets were kept between two platinum electrodes (the effective diameter for the working electrode was 10 mm) under spring-loaded pressure. The setup was put inside a furnace with temperature control from 298 to 423 K in steps of 25 K, and, in this case, measurements were taken using a dielectric interface SOLARTRON 1296 coupled to SOLARTRON 1260 analyzer (Solartron Analytical, Farnborough, UK). The cell was allowed to reach thermal stability for about 20 min before each measurement. Data were recorded in a frequency range of 10 MHz to 0.1 Hz with a voltage amplitude of 100 mVrms and 0 DC bias. The fitting circuit simulation of AC impedance data was performed with ZView® (Scribner Association, Southern Pines, North Carolina, NC, USA) software.

#### **4. Conclusions**

Changes in symmetry and microstructure of synthesized birnessite-type materials produced variations in their electrical properties and relaxation times for the electrical transport mechanism. Mn4 exhibited the highest conductivity, which was associated with turbostratic disorder, shorter Mn–O bond length, lower crystal size, higher Mn4<sup>+</sup> content, and higher surface area. For Mn8, the crystal size increment might have diminished the contact area between the particle agglomerates, making electron transport more difficult than in the Mn4 material. The enhanced conductivity with temperature suggests a thermally activated electron "hopping" and semiconductor behaviour for the synthesized materials. It is concluded that conductivity and other electronic-related phenomena can be tuned by varying the synthesis conditions.

**Author Contributions:** Synthesis conceptualization, O.G.; methodology, O.G., N.P.A., and M.E.B.; validation M.E.B.; formal analysis, N.P.A. and O.G.; investigation and data curation N.P.A. and M.E.B.; writing N.P.A. and O.G.; all authors reviewed and approved the manuscript.

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

**Acknowledgments:** The authors wish to thank the Dirección de Laboratorios de Sede, Laboratorio de Química (Atomic Absorption Analysis) and Laboratorio de Magnetismo y Materiales Avanzados (Thermal Analysis) at the Universidad Nacional de Colombia, Sede Manizales. The authors also thank the Laboratorio de Microscopia Avanzada (SEM) at the Universidad Nacional de Colombia, Sede Medellín. N.P.A acknowledge to Davián Martinez Buitrago for his advice in the Rietveld Refinement.

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

#### **Appendix A**


**Table A1.** Fractional coordinates for atom positions in Mn4 and Mn8.

**Table A2.** Circuital element values from equivalent circuit fitting of the experimental results for Mn4.


CPE1-T and CPE1-P are the coefficients of the constant phase element in Equation (1); Wo-R, Wo-T, and WoP are the coefficients for Warburg open element in Equation (2); and Rg and Rgb are the resistances in grain and grain boundary, respectively.


**Table A3.** Circuital element values from equivalent circuit fitting of the experimental results for Mn8.

CPE1-T and CPE1-P are the coefficients of the constant phase element in Equation (1); Wo-R, Wo-T, and WoP are the coefficients for Warburg open element in Equation (2); and Rg and Rgb are the resistances in grain and grain boundary, respectively.

#### **References**


**Sample Availability:** Samples of the compounds Mn4 and Mn8 are available from the authors.

© 2019 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/).

*Article*
