Nonlinear Dynamic Behavior of Porous and Imperfect Bernoulli-Euler Functionally Graded Nanobeams Resting on Winkler Elastic Foundation

Abstract

Nonlinear free vibrations of functionally graded porous Bernoulli–Euler nano-beams resting on an elastic foundation through a stress-driven nonlocal elasticity model are studied taking into account von Kármán type nonlinearity and initial geometric imperfection. By using the Galerkin method, the governing equations are reduced to a nonlinear ordinary differential equation. The closed form analytical solution of the nonlinear natural flexural frequency is then established using the Hamiltonian approach to nonlinear oscillators. Several comparisons with existing models in the literature are performed to validate the accuracy and reliability of the proposed approach. Finally, a numerical investigation is developed in order to analyze the effects of the gradient index coefficient, porosity volume fraction, initial geometric imperfection, and the Winkler elastic foundation coefficient, on the nonlinear flexural vibrations of metal–ceramic FG porous Bernoulli–Euler nano-beams.

1. Introduction

Functionally graded materials (FGMs) are advanced composites designed and fabricated in a way that their physical and mechanical properties spatially vary in their structures. An overview of manufacturing methods for FGMs, their applications and future challenges, based on the available literature over 30 years, has been recently published by Saleh et al., in [1].

Due to the continuous variation in material properties, many problems related to the material discontinuities of conventional composite material (stress concentrations, residual stresses, delamination phenomena and damage growth) can be significantly reduced and high permeance requirements ensured [2,3,4].

In most cases, functionally graded materials are made of a combination of metal and ceramic to offer a wide range of applications for various equipment subject to extreme thermo-mechanical stresses [5,6,7,8,9,10,11,12,13].

Recent studies have also shown that, by managing some fabrication parameters during the manufacture of FGMs, different kinds of porosity distributions can be fabricated inside their structure to further improve the physical and mechanical characteristics of the material [14,15,16,17,18]. For example, the results of the investigations reported in [19,20,21,22] have highlighted the benefits of cellular metals and metal foams inside the structure of micro/nano-scale systems (MEMS/NEMS) in terms of electrical conductivity, thermal transport and energy absorption.

Various simple models have been widely adopted by many researchers to describe the interaction of the nano-beam with its foundation for a wide range of engineering applications. However, to the best of our knowledge, the analyses of the mechanical behavior of porous nano-beams resting on elastic foundations are rather rare and the available results are limited to FG nano-beams without porosities [23,24,25,26,27].

As is well-known, the mechanical response of nanostructures made of functionally graded materials is strongly size-dependent and therefore the small-scale effects on the mechanical behavior of FG nanostructures must be considered [28,29,30,31,32].

In order to overcome the enormous costs and computation times of molecular dynamics (MD) simulations, and the difficulties to conduct experimental investigations at the nanoscale, the size-dependent effects on the static and dynamic characteristics of small-scaled structural systems may be accurately captured by the so-called scale-dependent continuum mechanics-based theories, including the strain gradient theory of elasticity [33,34] and Eringen’s nonlocal elasticity theory [35,36,37]. In the first theory, the material response at a point of a continuum is assumed to be dependent on both the strain and the strain gradients of different orders, while, in the second one, the output field at a point of a continuum is assumed to be the integral convolution between the elastic source field and a suitable averaging kernel.

In the framework of nonlocal elasticity, the basic integral constitutive law is the strain-driven Eringen’s integral model (EIM), which was recommended, for the first time, by Peddieson [38] to analyze the mechanical behavior of nanoscale structures. For unbounded continua, when considering the bi-exponential nonlocal kernel function, due to tacit fulfillment of the vanishing boundary conditions at infinity, the integral strain-driven theory (EIM) can be replaced with the well-known differential type theory (EDM), by which it is much easier to model small-scale phenomena in nanostructures than with the corresponding integral formulation.

Recently, it was shown that the nonlocal differential-based and integral-based theories of elasticity models may be not equivalent to each other for boundary condition problems, since adequate higher-order homogeneous constitutive boundary conditions have to be prescribed.

Consequently, for bounded nanostructures, the strain-driven purely nonlocal elastic problems defined on bounded domains are mathematically ill-posed due to the incompatibility between the higher-order constitutive boundary conditions and equilibrium requirements (e.g., the paradox of a cantilever nano-beam subject to a transverse concentrated load at the free end) [39].

The ill-posed problem related to the pure nonlocal model may be circumvented by adopting the Eringen local-nonlocal mixture constitutive model [37] or by using coupled theories, such as the nonlocal strain gradient theory [40] and the combination of pure nonlocal theory with the surface theory of elasticity [41]. For example, some recent applications of the theories mentioned above are addressed in [42,43,44].

Furthermore, these difficulties can be overcome by adopting the stress-driven nonlocal integral model (SDM) recently proposed by Romano and Barretta [45], in which the roles of stress and elastic strain fields are swapped with respect to the strain-driven model. In addition, in this case, for a class of bi-exponential kernels, the integral form of the constitutive equations is shown to be mathematically equivalent to differential equations subjected to some higher order constitutive boundary conditions.

The stress-driven nonlocal theory of elasticity has been widely respected by the scientific community and has been successfully applied to investigate size-dependent behavior of elastic nano-beams [46,47,48,49,50,51,52,53,54] with bounded structural domains. In particular, closed form solutions for stress-driven integral models [55] have been provided for cantilever, simply supported, clamped-pinned and fully-clamped FG nano-beams subject to different load conditions.

Recently, the authors of this paper examined the nonlinear vibration behavior of geometrically imperfect functionally graded nano-beams on the basis of both the stress-driven nonlocal integral model (SDM) and the strain-driven model (EDM), especially in the presence of an axial pretension force [56].

In the present paper, the nonlinear free vibration of a metal-ceramic functionally graded Bernoulli–Euler nano-beam, resting on a Winkler elastic foundation, is examined taking into account von Kármán type nonlinearity and initial geometric imperfection. In particular, small scale effects are considered by using the stress-driven nonlocal elasticity. The equation of motion is determined by Hamilton’s principle and simplified by Galerkin’s method. The first order Hamiltonian approach [57] was then employed to solve the nonlinear vibrations problem.

In order to show the efficiency and accuracy of the proposed approach, several comparisons with existing models in the literature are performed. Finally, the effects of initial geometric imperfection, the gradient index coefficient, the porosity volume fraction and the Winkler elastic foundation coefficient on nonlinear fundamental frequencies of porous FG simply supported nano-beams are presented and discussed.

2. Functionally Graded Porous Material

We consider a Bernoulli–Euler straight nano-beam, with length “L”, thickness “h” and width “b”, resting on a Winkler foundation, and denote by y’ and z’ the principal axes of geometric inertia of its rectangular cross-section (Σ), originating at the geometric center O (Figure 1).

The nano-beam is supposed to be made of a functionally graded (FG) material, composed of a mixture of metal and ceramic, whose distribution spatially varies from the bottom (z’ = −h/2) to the top (z’ = +h/2) surface. In this investigation, the top surface is ceramic-rich, whereas the bottom surface is metal-rich.

The nano-beam structure is assumed to have porosities inside the FG material generated during the manufacturing process of the two constituents. In particular, the porous FG nano-beam is assumed to have two kinds of porosity distributions across the thickness of the beam (Figure 1): in the first scenario, the porosity is evenly distributed among the metal and ceramic with widespread porosities; in the second scenario (uneven distribution), the porosities are distributed around the central area of the cross section of the FG nano-beam and tend to vanish both on the top surface (ceramic-rich) and on the bottom surface (metal-rich) of the nano-beam.

By denoting with k (k ≥ 0) the gradient index of the FG material, and with ζ the porosity volume fraction (ζ ≪ 1), in the case of an even distribution of porosity, the effective material properties, here described by the mass density, $\rho$(z’), and by Young’s modulus, E(z’), can be formulated by the following power-laws [58]:

$$\rho \left(z’\right)={\rho }_m+\left({\rho }_c-{\rho }_m\right){\left(\frac{1}{2}+\frac{z^{\prime }}{h}\right)}^k-\frac{\mathsf{\zeta }}{2} \left({\rho }_c+{\rho }_m\right),$$

(1)

$$E\left(z’\right)=E_m+\left(E_c-E_m\right){\left(\frac{1}{2}+\frac{z^{\prime }}{h}\right)}^k- \frac{\mathsf{\zeta }}{2} \left(E_c+E_m\right),$$

(2)

where $\rho$_c, $\rho$_m and E_c, E_m are the material densities and the Euler–Young moduli of ceramic and metal, respectively.

In the second scenario (uneven distribution of the porosity), the previous expressions are modified as follows

$$\rho \left(z’\right)={\rho }_m+\left({\rho }_c-{\rho }_m\right){\left(\frac{1}{2}+\frac{z^{\prime }}{h}\right)}^k-\frac{\mathsf{\zeta }}{2} \left({\rho }_c+{\rho }_m\right)\left(1-\frac{2\left|z^{\prime }\right|}{h}\right),$$

(3)

$$E\left(z’\right)=E_m+\left(E_c-E_m\right){\left(\frac{1}{2}+\frac{z^{\prime }}{h}\right)}^k- \frac{\mathsf{\zeta }}{2} \left(E_c+E_m\right)\left(1-\frac{2\left|z^{\prime }\right|}{h}\right).$$

(4)

As is well-known, in order to remove bending–stretching coupling caused by FG material variation, it is convenient to take as reference for the evaluation of the previous effective materials properties (Equations (1)–(4)) the elastic center C, whose position is shifted from the geometric center O of the following quantity

$$z_c^{’}=\frac{{\int }_{\mathsf{\Sigma }}E\left(z^{\prime }\right)z^{\prime }d\mathsf{\Sigma } }{{\int }_{\mathsf{\Sigma }}E\left(z^{\prime }\right)d\mathsf{\Sigma } }$$

(5)

Accordingly, in the new elastic Cartesian reference system at elastic center C, z = z’−$z_c^{’ }$ and y = y’ (Figure 1).

Moreover, determination of the bending stiffness, $EI$, and the axial stiffness, $EA$, of the FG porous nano-beam is useful in the subsequent relations. The first quantity is defined by the second moment of elastic area, weighted with the scalar field of the Euler–Young moduli, about the y axis ($I=I_y$), and is evaluated considering the bending abscissa z originating at C. Similarly, the axial stiffness is defined by the elastic cross-sectional area weighted with the scalar field of the Euler–Young moduli

$$\left\{EA,EI\right\}=b{\int }_{-\frac{h}{2}-z_c^{’}}^{\frac{h}{2}-z_c^{’}} E\left(z\right)\left\{1,z^2 \right\}d\mathrm{z}.$$

(6)

Moreover, the effective cross-sectional mass, $m_0$, and rotatory inertia, $m_2$, can be determined employing the previously introduced effective material density $\rho \left(z\right)$

$$m_0=b{{\displaystyle \int }}_{-\frac{h}{2}-z_c^{’}}^{\frac{h}{2}-z_c^{’}}\rho \left(z\right)d\mathrm{z},$$

(7)

$$m_2=b{{\displaystyle \int }}_{-\frac{h}{2}-z_c^{’}}^{\frac{h}{2}-z_c^{’}}\rho \left(z\right)z^{2}dz.$$

(8)

By substituting Equations (2) and (4) into Equation (5), we obtain the following coordinate of the elastic center C for the two distributions of porosity considered here

-

even distribution

$$z_{c_{even}}^{’}=\frac{\left(-E_c+E_m\right)hk}{\left(2+k\right)\left(-2\left(E_c+E_{m}k\right)+\left(E_c+E_m\right)\left(1+k\right)\zeta \right)},$$

(9)

-

uneven distribution

$$z_{c_{uneven}}^{’}=\frac{2\left(-E_c+E_m\right)hk}{\left(2+k\right)\left(-4\left(E_c+E_{m}k\right)+\left(E_c+E_m\right)\left(1+k\right)\zeta \right)}.$$

(10)

Utilizing the expressions of the effective Euler–Young modulus $E\left(z\right)$ introduced above for the two scenarios of porosity distribution taken into account (Equations (2) and (4)), bending stiffness, axial stiffness, cross-sectional mass and rotatory inertia can be determined as follows

-

even distribution

$$EI=\frac{\left(-48E_c^2-16E_c{E}_{m}k\left(7+k\left(4+k\right)\right)-4E_m^2{k}^2\left(7+k\left(4+k\right)\right)+4\left(E_c+E_m\right)\left(2+k\right)\left(2E_c\left(3+k\left(2+k\right)\right)+E_{m}k\left(7+k\left(4+k\right)\right)\right)\zeta -{\left(E_c+E_m\right)}^2\left(1+k\right){\left(2+k\right)}^2\left(3+k\right){\zeta }^2\right)}{2{\left(2+k\right)}^2\left(3+k\right)\left(-2\left(E_c+E_{m}k\right)+\left(E_c+E_m\right)\left(1+k\right)\zeta \right)} I{y}^{\prime },$$

(11)

$$EA= A\left(\frac{\left(E_c+k{E}_m\right)}{\left(1+k\right)}-\frac{1}{2}\left(E_c+E_m\right) \zeta \right),$$

(12)

$$m_0=A\left(\frac{\left({\rho }_c+k{\rho }_m\right)}{1+k}-\frac{1}{2}\left({\rho }_c+{\rho }_m\right)\zeta \right),$$

(13)

$$m_2=\frac{\left(-48{\rho }_c^2-16k\left(7+k\left(4+k\right)\right){\rho }_c{\rho }_m-4k^2\left(7+k\left(4+k\right)\right){\rho }_m^2-\left(1+k\right){\left(2+k\right)}^2\left(3+k\right){\zeta }^2{\left({\rho }_c+{\rho }_m\right)}^2+4\left(2+k\right)\zeta \left({\rho }_c+{\rho }_m\right)\left(2\left(3+k\left(2+k\right)\right){\rho }_c+k\left(7+k\left(4+k\right)\right){\rho }_m\right)\right)}{2{\left(2+k\right)}^2\left(3+k\right)\left(\left(1+k\right)\zeta \left({\rho }_c+{\rho }_m\right)-2\left({\rho }_c+k{\rho }_m\right)\right)} I{y}^{\prime };$$

(14)

-

uneven distribution

$$EI = \frac{\left(-32\left(12E_c^2+4E_c{E}_{m}k\left(7+k\left(4+k\right)\right)+E_m^2{k}^2\left(7+k\left(4+k\right)\right)\right)+4\left(E_c+E_m\right)\left(2+k\right)\left(E_{m}k\left(22+k\left(11+3k\right)\right)+E_c\left(18+k\left(11+7k\right)\right)\right)\zeta -{\left(E_c+E_m\right)}^2\left(1+k\right){\left(2+k\right)}^2\left(3+k\right){\zeta }^2\right)}{8 {\left(2+k\right)}^2\left(3+k\right)\left(-4\left(E_c+E_{m}k\right)+\left(E_c+E_m\right)\left(1+k\right)\zeta \right)} I{y}^{\prime },$$

(15)

$$EA=A \left(\frac{\left(E_c+E_{m}k\right)}{1+k}-\frac{1}{4}\left(E_c+E_m\right)\zeta \right),$$

(16)

$$m_0=A \left(\frac{\left({\rho }_c+{\rho }_{m}k\right)}{1+k}-\frac{1}{4}\left({\rho }_c+{\rho }_m\right)\zeta \right),$$

(17)

$$m_2=\frac{\left(-\left(1+k\right){\left(2+k\right)}^2\left(3+k\right){\zeta }^2{\left({\rho }_c+{\rho }_m\right)}^2+4\left(2+k\right)\zeta \left({\rho }_c+{\rho }_m\right)\left(\left(18+k\left(11+7k\right)\right){\rho }_c+k\left(22+k\left(11+3k\right)\right){\rho }_m\right)+32\left(-12{\rho }_c^2-4k\left(7+k\left(4+k\right)\right){\rho }_c{\rho }_m-k^2\left(7+k\left(4+k\right)\right){\rho }_m^2\right)\right)}{8{\left(2+k\right)}^2\left(3+k\right)\left(\left(1+k\right)\zeta \left({\rho }_c+{\rho }_m\right)-4\left({\rho }_c+k{\rho }_m\right)\right)} I{y}^{\prime },$$

(18)

where I_y’ = bh³/12 is the second geometric moment of area of the cross-section about the y’ axis about the geometric center and A = bh is the cross-sectional area.

If the porosity volume fraction ζ → 0, it is inferred from Equations (11)–(18) that, as k → 0 and k → ∞, the aforementioned material quantities approach the corresponding material quantities of pure ceramic and pure metal, respectively.

In this work, the nano-beam material is assumed to be made of silicon nitride-stainless steel (Si3N4-SuS3O4), whose properties are listed in the following

Table 1. Material properties of the FG nano-beam (silicon nitride-stainless steel).

Material	Young’s Modulus [GPa]	Poisson’s Ratio	Mass Density [kg/m3]
SuS3O4 (stainless steel)	207.8	0.3178	8166
Si3N4 (silicon nitride)	322.3	0.24	2370

[59,60].

Figure 2 shows the effect of the gradient index of the porous FG material ($k$) on the dimensionless position of the elastic center C (${\tilde{z}}_c=z_c^{\prime }/h)$of the rectangular cross section of the nano-beam.

Firstly, it can be noted that the dimensionless distance ${\tilde{z}}_c$increases as the gradient index $k$ increases, with a maximum attained at a value of $k$*, which depends on the value of $\zeta$; note that for $k>k$*, the distance ${\tilde{z}}_c$tends to vanish for $k$ → ∞. From Figure 2, one can also note that the curves corresponding to the second type of porosity distribution exhibit less of an increase than those of the first scenario of porosity.

Variations of the dimensionless axial stiffness ($\widetilde{EA}=\frac{EA}{E_{c}A})$and dimensionless bending stiffness ($\widetilde{EI}=\frac{EI}{E_{c}I}$) of the porous FG nano-beam, in terms of the gradient index, are reported in Figure 3, varying the porosity volume fraction of the material. As can be observed, the increase in the gradient index, as well as in the value of the porosity volume fraction, causes a decrease of the axial and bending stiffness of porous FG nano-beams, for both the two kinds of porosity distribution considered.

In particular, for each assigned value of $\zeta$ the dashed lines, corresponding to the uneven distributions of porosity, always have a greater value than the continuous ones, which correspond to the even distribution.

Figure 4 show the variations of the dimensionless cross-sectional mass ($m_o$) and rotatory inertia ($m_2)\mathrm{in}\mathrm{terms}$of the gradient index, respectively, varying the porosity volume fraction. From these graphs, one finds that$m_o$and$m_2$increase with increasing $\zeta$and with decreasing $k$.

The variations of the Euler–Young modulus across the thickness of a nano-beam are illustrated in Figure 5A, assuming $\zeta$ equal to zero and varying $k$ in the set {0, 0.3, 0.5, 1, 3, 5, → ∞}. As was to be expected, when $k$ = 0, the FG material reduces to pure ceramic (E = E_c = 322.3 GPa) while, on the contrary, for $k$ → ∞, the material properties tend to pure metal (E = E_m = 207.8 GPa).

The variations of the Euler–Young modulus through the thickness of the nano-beam cross section, assuming $\zeta$ = 0.2, are illustrated in Figure 5B,C for even and uneven porosity distributions, respectively. As can be noted, the curves of the variation of E relative to the even distribution (Figure 5B) have the same behavior of those illustrated in Figure 5A, corresponding to a non-porous material, but with a decrease in the values of Young’s moduli. Finally, Figure 5C shows that the maximum of Euler–Young’s modulus for the uneven distribution of porosity is reached at the top and bottom of the cross section and decreases in the direction of the middle zone.

3. Governing Equation

Using the Bernoulli–Euler theory, the displacement components $u_x,$ $u_y$ and $u_z$ along x, y and z directions, respectively, at any point of the nano-beam, can be written as

$$u_x\left(x,z,t\right)=u_0\left(x,t\right)+z {\phi }_y\left(x,t\right),$$

(19)

$$u_y\left(x,z,t\right)=0,$$

(20)

$$u_z\left(x,z,t\right)=w_0\left(x,t\right)+w^{*}\left(x\right),$$

(21)

where $u_0\left(x,t\right)$ and $w_0\left(x,t\right)$are the axial and transverse displacements of the elastic center C at time t, respectively; ${\phi }_y\left(x,t\right)=-\frac{\partial w_0}{\partial x}\left(x,t\right)$ is the rotation of the cross section about the y-axis, and $w^{*}\left(x\right)$is the initial assigned geometric imperfection in the transverse direction [61].

On the basis of the nonlinear Von-Kármán assumptions, with small strains and moderate rotation, the nonlinear strain–displacement relation may be written as [62]

$${\epsilon }_{xx}=\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2- z\frac{{\partial }^2{w}_0}{\partial x^2}+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}.$$

(22)

As is well-known, the nonlinear equations of motion can be obtained using Hamilton’s principle:

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

(23)

where the expression of $K$ (kinetic energy), $U$ (strain energy) and $W$ (work done by external forces) are given in the following.

The variation of kinetic energy $K$ is

$$\begin{gathered} \delta K= \underset{0}{\overset{L}{\int }}{\int }_{\mathsf{\Sigma }}\begin{array}{c}\rho \frac{\partial u_i}{\partial t} \delta \frac{\partial u_i}{\partial t} d\mathsf{\Sigma }dx=\\ \end{array} \underset{0}{\overset{L}{\int }}{\int }_{\mathsf{\Sigma }}\begin{array}{c}\rho \left[\left(\frac{\partial u_0}{\partial t}+z\frac{\partial {\phi }_y}{\partial t}\right)\left(\delta \frac{\partial u_0}{\partial t}+z \delta \frac{\partial {\phi }_y}{\partial t}\right)+\frac{\partial w_0}{\partial t}\delta \frac{\partial w_0}{\partial t}\right]d\mathsf{\Sigma }dx= \\ \end{array} \\ = {\int }_0^L (m_0 \frac{\partial u_0}{\partial t}\delta \frac{\partial u_0}{\partial t}+ m_2\frac{\partial {\phi }_y}{\partial t}\delta \frac{\partial {\phi }_y}{\partial t}+m_0\frac{\partial w_0}{\partial t}\delta \frac{\partial w_0}{\partial t}) dx. \end{gathered}$$

(24)

The expression for the variation of the strain energy $U$ is

$$\begin{gathered} \delta U= \underset{0}{\overset{L}{\int }}{\int }_{\mathsf{\Sigma }}{\sigma }_{xx}\delta {\epsilon }_{xx} d\mathsf{\Sigma }dx= \underset{0}{\overset{L}{\int }}{\int }_{\mathsf{\Sigma }}{\sigma }_{xx}\delta \left(\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2+ z \frac{\partial {\phi }_y}{\partial x}+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)d\mathsf{\Sigma }dx= \\ ={\int }_0^L(N \left(\frac{\partial \delta u_0}{\partial x}+\frac{\partial w_0}{\partial x}\frac{\partial \delta w_0}{\partial x}+\frac{\partial \delta w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)+M\frac{\partial \delta {\phi }_y}{\partial x}) dx, \end{gathered}$$

(25)

where the stress resultants $N$ and $M$ introduced in Equation (25) are defined as:

$$N={\int }_{\mathsf{\Sigma }}{\sigma }_{xx} d\mathsf{\Sigma },$$

(26)

$$M={\int }_{\mathsf{\Sigma }}z {\sigma }_{xx} d\mathsf{\Sigma }.$$

(27)

The expression of the virtual work of the external force $W$can be expressed by

$$\delta W= -{\int }_0^L\left(q_x \delta u_0+q_z \delta w_0 - K_w\left(x\right)w_0\left(x,t\right)\delta w_0\right)dx,$$

(28)

where $q_x=q_x\left(x\right)$ and $q_z=q_z\left(x\right)$ are the axial and the transverse vertical distributed loads, respectively, and $K_w\left(x\right)$is the Winkler elastic foundation coefficient.

By substituting the expressions of $K$, $U$, and $W$ into Hamilton’s principle (Equation (23)), performing integration-by-parts with respect to t as well as x to relieve the generalized displacements $\delta u_0$, $\delta w_0$ and $\delta {\phi }_y$of any differentiations, and using the fundamental lemma of differential calculus, we obtain the following equations of motion:

$$\frac{\partial N}{\partial x}+q_x=m_0\frac{{\partial }^2{u}_0}{\partial t^2},$$

(29)

$$\frac{{\partial }^{2}M}{\partial x^2}+q_z- K_w\left(x\right) w_0+\frac{\partial }{\partial x}\left(N\frac{\partial w_0}{\partial x}\right)+\frac{\partial }{\partial x}\left(N\frac{d{w}^{*}}{dx}\right)=m_0\frac{{\partial }^2{w}_0}{\partial t^2}- m_2 \frac{{\partial }^4{w}_0}{\partial x^2 \partial t^2}.$$

(30)

The boundary conditions involve specifying one element of each of the following three pairs:

$$u_0 or N\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{at} x=0,L;$$

(31)

$$w_0 or \frac{\partial M}{\partial x}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{at}x=0,L;$$

(32)

$$-\frac{\partial w_0}{\partial x}\hspace{1em}or\hspace{1em}M\hspace{1em}\hspace{1em}\mathrm{at} x=0,L.$$

(33)

4. Stress-Driven Nonlocal Integral Model

In adapting the stress-driven integral formulation to the Bernoulli–Euler porous FG nano-beam, the nonlocal axial strain, ${\epsilon }_{xx}$, and the nonlocal bending curvature, $\chi$, are obtained by the following convolutions

$${\epsilon }_{xx}\left(x\right)={\int }_0^L{\mathsf{\Phi }}_{\lambda }\left(x-\xi \right) \frac{N\left(\xi \right)}{EA}d\xi ,$$

(34)

$$\chi \left(x\right)={\int }_0^{\mathrm{L}}{\mathsf{\Phi }}_{\lambda }\left(x-\xi \right)\frac{M\left(\xi \right)}{EI} d\xi ,$$

(35)

where $x$, $\xi$ are position vectors of points of the domain of the Euclidean space occupied by the nano-beam at time t, and the stress resultants $N$ and $M$ fulfil the equilibrium conditions.

By assuming the following bi-exponential function for the averaging kernel ${\mathsf{\Phi }}_{\lambda }$ with internal characteristic length, $L_c$

$${\mathsf{\Phi }}_{\lambda }\left(x,L_c\right)=\frac{1}{2L_c}\exp (-\frac{\left|x\right|}{L_c}),$$

(36)

the stress-driven integrals (Equations (34) and (35)) admit the following set of solutions

$${\epsilon }_{xx}\left(x\right)-L_c{}^2 \frac{{\partial }^2{\epsilon }_{xx}\left(x\right)}{\partial x^2}=\frac{N\left(x\right)}{EA} ,$$

(37)

$$\chi \left(x\right)-L_c{}^2 \frac{{\partial }^2\chi \left(x\right)}{\partial x^2}=\frac{M\left(x\right)}{EI}{,}_{}$$

(38)

with x ∈ [0, L], if and only if the following two pairs of constitutive boundary conditions are satisfied at the nano-beam ends

$$\frac{\partial {\epsilon }_{xx} }{\partial x}\left(0\right)=\frac{1}{L_c}{\epsilon }_{xx}\left(0\right),$$

(39)

$$\frac{\partial {\epsilon }_{xx} }{\partial x}\left(\mathrm{L}\right)=-\frac{1}{L_c}{\epsilon }_{xx}\left(\mathrm{L}\right),$$

(40)

$$\frac{\partial \chi }{\partial x}\left(0\right)=\frac{1}{L_c}\chi \left(0\right){,}_{}$$

(41)

$$\frac{\partial \chi }{\partial x}\left(\mathrm{L}\right)=-\frac{1}{L_c}\chi \left(\mathrm{L}\right){.}_{}$$

(42)

4.1. Equation of Motion (SDM)

According to the stress-driven differential model (SDM), the nonlocal axial force and moment resultants can be described explicitly in terms of generalized strain components as follows

$$\begin{gathered} N\left(x,t\right)={\int }_{\mathsf{\Sigma }}{\sigma }_{xx} d\mathsf{\Sigma }= {\int }_{\mathsf{\Sigma }}E \left({\epsilon }_{xx}\left(x,t\right)-L_c{}^2 \frac{{\partial }^2{\epsilon }_{xx}\left(x,t\right)}{\partial x^2}\right)d\mathsf{\Sigma }= \\ = {\int }_{\mathsf{\Sigma }}E \left(\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2- z \frac{{\partial }^2{w}_0}{\partial x^2}+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right) \\ -E L_c{}^2 \frac{{\partial }^2}{\partial x^2}\left(\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2- z \frac{{\partial }^2{w}_0}{\partial x^2}+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)d\mathsf{\Sigma }= \\ = EA \left(\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right) \\ -EA{L}_c{}^2 \frac{{\partial }^2}{\partial x^2}\left(\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right); \end{gathered}$$

(43)

$$\begin{gathered} M\left(x,t\right)={\int }_{\mathsf{\Sigma }}{\sigma }_{xx} z d\mathsf{\Sigma }= {\int }_{\mathsf{\Sigma }}E \left({\epsilon }_{xx}\left(x,t\right)-L_c{}^2 \frac{{\partial }^2{\epsilon }_{xx}\left(x,t\right)}{\partial x^2}\right)z d\mathsf{\Sigma }= \\ = {\int }_{\mathsf{\Sigma }}E \left(\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2- z \frac{{\partial }^2{w}_0}{\partial x^2}+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)z \\ -E{L}_c{}^2 \frac{{\partial }^2}{\partial x^2}\left(\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2- z \frac{{\partial }^2{w}_0}{\partial x^2}+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)z d\mathsf{\Sigma }= \\ =- EI \left(\frac{{\partial }^2{w}_0}{\partial x^2}-L_c{}^2\frac{{\partial }^4{w}_0}{\partial x^4} \right). \end{gathered}$$

(44)

By substituting Equations (43) and (44) into Equations (29) and (30), the following stress-driven equations of motion can be derived

$$EA \frac{\partial }{\partial x}\left(\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)-EA L_c{}^2 \frac{{\partial }^3}{\partial x^3}\left(\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)+q_x=m_0\frac{{\partial }^2{u}_0}{\partial t^2},$$

(45)

$$EI \left(\frac{{\partial }^4{w}_0}{\partial x^4}-L_c{}^2\frac{{\partial }^6{w}_0}{\partial x^6} \right)- q_z+K_w\left(x\right) w_0-\frac{\partial }{\partial x}\left(N\frac{\partial w_0}{\partial x}\right)-\frac{\partial }{\partial x}\left(N\frac{d{w}^{*}}{dx}\right)=m_2 \frac{{\partial }^4{w}_0}{\partial x^2 \partial t^2} - m_0\frac{{\partial }^2{w}_0}{\partial t^2},$$

(46)

with the following natural boundary conditions at the nano-beam ends

$$N= \overline{N}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{at}\hspace{1em}x=0,L;$$

(47)

$$m_2 \frac{{\partial }^3{w}_0}{\partial x \partial t^2}+ \frac{\partial M}{\partial x}=\overline{V}\hspace{1em} \mathrm{at}\hspace{1em}x=0,L;$$

(48)

$$M=\overline{M}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{at}\hspace{1em}x=0,L;$$

(49)

where $\overline{N}$, $\overline{M}$ and$\overline{V}$ the assigned generalized forces acting at the nano-beam ends together and with the constitutive boundary conditions at the nano-beam ends given by Equations (38)–(42).

According to the Bernoulli–Euler kinematics, the geometric bending curvature is related to the transverse displacement $w$ by

$$\chi \left(x,t\right)=\frac{{\partial }^{2}w}{\partial x^2}\left(x,t\right).$$

(50)

Therefore, the constitutive boundary conditions (Equations (38)–(42)) can be rewritten as a function of the displacement components

$$\frac{{\partial }^2{u}_0}{\partial x^2}\left(0,t\right)=\frac{1}{L_c}\frac{\partial u_0}{\partial x}\left(0,t\right),$$

(51)

$$\frac{{\partial }^2{u}_0}{\partial x^2}\left(L,t\right)=-\frac{1}{L_c}\frac{\partial u_0}{\partial x}\left(L,t\right);$$

(52)

$$\frac{{\partial }^3{w}_0}{\partial x^3}\left(0,t\right)=\frac{1}{L_c}\frac{{\partial }^2{w}_0}{\partial x^2}\left(0,t\right),$$

(53)

$$\frac{{\partial }^3{w}_0}{\partial x^3}\left(L,t\right)=-\frac{1}{L_c}\frac{{\partial }^2{w}_0}{\partial x^2}\left(L,t\right).$$

(54)

Nonlinear Transverse Free Vibrations (SDM)

The nonlinear transverse free vibrations ($q_x=$ $q_z=0)$of an imperfect porous FG nano-beam, based on the stress-driven model, are derived assuming that the rotary inertia term, $m_2 \frac{{\partial }^4{w}_0}{\partial x^2 \partial t^2}$, is very small compared to the other terms. Hence, from Equation (46), it follows

$$EI\frac{{\partial }^4{w}_0}{\partial x^4 }- EI L_c{}^2\frac{{\partial }^6{w}_0}{\partial x^6} -\frac{\partial }{\partial x}\left(N\frac{\partial w_0}{\partial x}\right)-\frac{\partial }{\partial x}\left(N\frac{d{w}^{*}}{dx}\right)+ m_0\frac{{\partial }^2{w}_0}{\partial t^2}=- K_w\left(x\right) w_0.$$

(55)

If we now neglect the axial inertia term, $m_0\frac{{\partial }^2{u}_0}{\partial t^2}$, from the first equation of motion, we obtain

$$N\left(x\right) = EA \left(\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)-EA L_c{}^2 \frac{{\partial }^2}{\partial x^2}\left(\frac{\partial u_0}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)=\widehat{N},$$

(56)

wherein $\widehat{N}$ is a constant.

Note that for a nano-beam with immovable ends ($u_0{|}_{x=0}$= $u_0{|}_{x=L}$= 0 and $w_0{|}_{x=0}$= $w_0{|}_{x=L}$= 0), by integrating both sides of Equation (56) over the domain [0, L] yields the following expression

$$\widehat{N}=\frac{EA}{L} {\int }_0^L\left(\left(\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)- L_c{}^2 \frac{{\partial }^2}{\partial x^2} \left(\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)\right)dx,$$

(57)

which is the so-called “mid-plane stretching effect” in SDM theory.

Based on this assumption, from Equation (55), it follows:

$$\begin{gathered} EI{L}_c{}^2\frac{{\partial }^6{w}_0}{\partial x^6}-EI\frac{{\partial }^4{w}_0}{\partial x^4 } - K_w\left(x\right) w_0+ \\ +(\frac{EA}{L} {\int }_0^L\left(\left(\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)- L_c{}^2 \frac{{\partial }^2}{\partial x^2} \left(\frac{1}{2}{\left(\frac{\partial w_0}{\partial x}\right)}^2+\frac{\partial w_0}{\partial x}\frac{d{w}^{*}}{dx}\right)\right)dx)\left(\frac{{\partial }^2{w}_0}{\partial x^2}+\frac{{\mathrm{d}}^2{w}^{*}}{\mathrm{d}{x}^2}\right)={\mathrm{m}}_0\frac{{\partial }^2{w}_0}{\partial t^2}. \end{gathered}$$

(58)

Equation (58) describes the nonlinear transverse free oscillations of a functionally graded (FG) porous and imperfect nano-beam resting on an elastic foundation. The solution procedure is reported in Appendix A.

5. Convergence and Comparison Study

Some numerical examples are presented in this paragraph to validate the accuracy and reliability of the proposed approach. In the first comparison example (

Table 2. Linear dimensionless natural frequencies of homogeneous simply-supported (S-S) and clamped-clamped (C-C) nano-beams (${\mathcal{A}}_w={\mathcal{A}}_0=0;{\tilde{K}}_w=0).$

λ	S-S		C-C
	Present Approach	Ref. [55]	Present Approach	Ref. [55]
0.00+	9.8696	9.8696	22.3733	22.3733
0.01	9.8744	9.8744	22.8518	22.8518
0.03	9.9107	9.9107	23.9976	23.9976
0.05	9.9786	9.9786	25.3918	25.3918
0.1	10.2534	10.2534	29.8303	29.8303

) geometrical nonlinearity, initial imperfection and the elastic foundation coefficient are neglected in order to compare the results obtained with the present approach. The comparison, in terms of linear dimensionless natural frequencies of homogeneous simply-supported (S-S) and clamped-clamped (C-C) FG nano-beams, is made with the results of Apuzzo et al., who in [55] obtained results while varying the dimensionless nonlocal parameter $\lambda$ in the set {0.00⁺, 0.01, 0.03, 0.05, 0.1}.

In the second comparison study (

Table 3. Linear dimensionless natural frequencies of homogeneous simply-supported (S-S) and clamped-clamped (C-C) nano-beams (${\mathcal{A}}_w={\mathcal{A}}_0=0;{\tilde{K}}_w=50).$

λ	S-S		C-C
	Present Approach	Ref. [59]	Present Approach	Ref. [59]
0.00+	12.1412	12.1412	23.4641	23.4641
0.01	12.1451	12.1451	23.9208	23.9208
0.03	12.1747	12.1747	25.0177	25.0177
0.05	12.2300	12.2300	26.3579	26.3579
0.1	12.4552	12.4552	30.6569	30.6569

), the present approach is compared with the model of Mahmoudpou et al., [59] for homogenous S-S and C-C nano-beams resting on an elastic foundation with ${\tilde{K}}_w=50$, and assuming that ${\mathcal{A}}_w$ and ${\mathcal{A}}_0$ are equal to zero.

The third example (

Table 4. Local nonlinear frequency ratios of perfect simply-supported (S-S) nano-beams ($\lambda =0;{\mathcal{A}}_0=0;{\tilde{K}}_w=0).$

${\mathcal{A}}_{\mathit{w} }/\mathit{L}$	Present Approach			RGM Ref. [63]	HAM Ref. [59]
	ωNL	ωL	ωNL/ωL	ωNL/ωL	ωNL/ωL
1.0	10.7552	9.8696	1.0897	1.0897	1.0892
2.0	13.0563	9.8696	1.3229	1.3229	1.3178
3.0	16.1798	9.8696	1.6394	1.6394	1.6259
4.0	19.7392	9.8696	2.0000	-	1.9766

) concerns the evaluation of the ratios between the local nonlinear frequency (ω_NL) and the local linear one (ω_L), for a perfect simply-supported homogeneous nano-beams, varying the dimensionless amplitude of the nonlinear oscillator ${\mathcal{A}}_w/\mathrm{L}$, and neglecting the elastic foundation coefficient. The results obtained with the present approach are compared with the results of Singh et al., [63], based on the Ritz–Galerkin method (RGM), and with those of Mahmoudpou et al., [59], based on the first-order homotopy analysis method (HAM).

Finally, in

Table 5. Nonlocal dimensionless natural frequencies of a homogeneous simply-supported (S-S) nano-beam resting on an elastic foundation for various dimensionless amplitude ratios ($\lambda =0.1;{\mathcal{A}}_0=0;{\tilde{K}}_w=50,100).$

${\mathcal{A}}_{\mathit{w} }/\mathit{L}$	${\tilde{\mathit{K}}}_{\mathit{w}}=50$			${\tilde{\mathit{K}}}_{\mathit{w}}=100$
	Present Approach	Ref. [59]	Δ1 [%]	Present Approach	Ref. [59]	Δ2 [%]
0.0	12.4552	12.4552	0.00	14.3225	14.3225	0.00
0.1	12.4694	-	-	14.3348	-	-
0.2	12.5120	-	-	14.3719	-	-
0.3	12.5827	12.5639	0.15	14.4334	14.4171	0.11
0.4	12.6810	-	-	14.5192	-	-
0.5	12.8062	-	-	14.6287	-	-
0.6	12.9576	12.8835	0.57	14.7614	14.6967	0.44

and

Table 6. Nonlocal dimensionless natural frequencies of a homogeneous clamped-clamped (C-C) nano-beam resting on an elastic foundation for various dimensionless amplitude ratios ($\lambda =0.1;{\mathcal{A}}_0=0;{\tilde{K}}_w=50,100).$

${\mathcal{A}}_{\mathit{w} }/\mathit{L}$	${\tilde{\mathit{K}}}_{\mathit{w}}=50$			${\tilde{\mathit{K}}}_{\mathit{w}}=100$
	Present Approach	Ref. [59]	Δ1 [%]	Present Approach	Ref. [59]	Δ2 [%]
0.0	30.6591	30.6575	0.01	31.4640	31.4624	0.01
0.1	30.6606	-	-	31.4654	-	-
0.2	30.6651	-	-	31.4698	-	-
0.3	30.6725	30.7440	0.23	31.4770	31.5467	0.22
0.4	30.6830	-	-	31.4872	-	-
0.5	30.6964	-	-	31.5003	-	-
0.6	30.7128	31.0017	0.93	31.5163	31.7980	0.89

the values of the nonlocal nonlinear dimensionless frequencies, evaluated with the present approach varying the dimensionless amplitude of the nonlinear oscillator ${\mathcal{A}}_w/L$ in the set {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6}, are summarized and compared with the results by Mahmoudpou et al., in [59] for two values of the elastic foundation coefficient: 50 and 100.

From these comparison studies, the accuracy of the Hamiltonian approach to nonlinear oscillators here employed is validated.

6. Results and Discussion

The nonlinear dynamic behavior of a simply-supported Bernoulli–Euler imperfect and porous FG nano-beam, resting on a Winkler elastic foundation, whose material properties are listed in Table 1, is considered as a case study in this section. The nano-beam has a length L = 20 nm, internal characteristic length $L_c=2\mathrm{nm}$, and a squared cross-section (b = h = 0.05 L).

In particular, the effects of the gradient index, porosity distribution, nonlinear oscillator amplitude and the dimensionless elastic foundation coefficient, on the frequency ratio between the non-linear frequency of porous FG nano-beam, ω_NL, and the linear frequency of a pure ceramic nano-beam, ω_cL, are presented.

6.1. Gradient Index and Porosity Volume Fraction

Firstly, the combined effect of the gradient index, $k$, and porosity volume fraction, $\zeta ,$ on the frequency ratio of a perfect FG nano-beam (${\mathcal{A}}_0=$ 0) are plotted in Figure 6, Figure 7 and Figure 8 assuming${\mathcal{A}}_w$= 0.2 and ${\tilde{K}}_w=$ 0. From these figures it can be observed that the frequency ratios evaluated for $k$ < 1 are always greater than those obtained for$k$ > 1. Moreover, the effects of $k$ and $\zeta ,$ for the two kinds of porosity distributions considered, on the aforementioned frequency ratio are plotted in Figure 9, varying the dimensionless elastic foundation coefficient, ${\tilde{K}}_w$, in the set {0, 20}, and the initial imperfection amplitude,${\mathcal{A}}_0$, in the set {0.0, 0.2}.

In particular, from Figure 9A it can be observed that, for a perfect FG nano-beam (${\mathcal{A}}_0$= 0) with an even distribution of porosity, the curves of the aforementioned ratios (continuous lines) tend to increase as the porosity volume fraction $\zeta$ increases when $k<1$, while an opposite trend is obtained when $k>1$. Moreover, it is found that the nonlinear frequency ratios of imperfect nano-beams (${\mathcal{A}}_0$ = 0.2) have values which are always greater than those of perfect nano-beams. On the contrary, for the second scenario of porosity distribution (Figure 9B), both in the case of perfect and imperfect FG nano-beams, the frequency ratios always increase as $\zeta$ increases, no matter what values of gradient index $k$ and elastic foundation coefficient ${\tilde{K}}_w$are adopted.

6.2. Nonlinear Oscillator Amplitude

The nonlinear frequency ratio versus the nonlinear oscillator amplitude,${\mathcal{A}}_w$, is presented in Figure 10, varying the elastic foundation coefficient ${\tilde{K}}_w$ and the porosity volume fraction $\zeta$.

As can be observed, the ratio ω_NL/ω_cL increases as the absolute value of ${\mathcal{A}}_w$increases. Moreover, the curves exhibit a symmetric behavior when a perfect FG nano-beam is considered, while a different trend can be observed for an imperfect FG nano-beam (${\mathcal{A}}_0$ = 0.5). It is also found that the nonlinear frequency ratio increases as the porosity volume fraction $\zeta$ increases and that, for a given value of $\zeta$, the ratio ω_NL/ω_cL also increases as the coefficient ${\tilde{K}}_w$increases.

6.3. Winkler Elastic Foundation Coefficient

The effect of the Winkler elastic foundation coefficient on the non-linear frequency ratio was investigated by varying the initial imperfection amplitude,${\mathcal{A}}_0$, in the set {0, 0.1, 0.3} for the first (Figure 11A) and the second (Figure 11B) kind of porosity distribution here considered.

The curves illustrated in these Figures, obtained for${\mathcal{A}}_w=0.2$and $k$ = 1, show that the non-linear frequency ratio tends to increase with increasing elastic foundation coefficient and increasing initial imperfection amplitude. Moreover, from the comparison of graphs shown in Figure 11, one finds that the nonlinear frequency ratios associated with the even distribution of porosity are always the largest for all considered values of the initial imperfection amplitude.

7. Conclusions

In this paper, the nonlinear dynamic behavior of a porous FG Bernoulli–Euler nano-beam, resting on a Winkler elastic foundation, with von Kármán type nonlinearity and initial geometric imperfection, has been studied, employing the stress-driven nonlocal integral model.

The governing equations have been reduced to a nonlinear ordinary differential equation by using the Galerkin method. Then, the Hamiltonian approach to nonlinear oscillators was employed to obtain a closed form analytical solution of the nonlinear natural frequency.

In view of the numerical results obtained in the present study, the following conclusions may be formulated:

(1)

the increase in the gradient index, as well as in the porosity volume fraction, cause a decrease in the values of the axial and bending stiffness of porous FG nano-beams;

(2)

an increase in the porosity volume fraction of perfect FG nano-beams causes an increase in the nonlinear frequency ratio when $k<1;$an opposite trend was observed when $k>1$;

(3)

the nonlinear frequencies of imperfect porous FG nano-beams are always greater than those obtained for perfect nano-beams;

(4)

the nonlinear frequency ratio always increases as the porosity volume fraction increases in the case of an uneven distribution of porosity;

(5)

an increase in the elastic foundation coefficient and in the initial imperfection amplitude causes an increase in the nonlinear frequency ratio.

In conclusion, the proposed approach is able to capture the nonlinear dynamic behavior of porous and imperfect Bernoulli–Euler functionally graded nano-beams resting on a Winkler elastic foundation and the approach provides a cost-effective method for the design and optimization of a wide range of nano-scaled beam-like components of nano-electro-mechanical-systems (NEMS).