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

: 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.


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].
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 nanobeams 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 sizedependent 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 nanobeams [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 fullyclamped 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 nanobeams are presented and discussed.

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, (z'), and by Young's modulus, E(z'), can be formulated by the following power-laws [58]: where c, m and Ec, Em 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 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 Accordingly, in the new elastic Cartesian reference system at elastic center C, z = z'-and y = y' (Figure 1).
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 [59,60].  Figure 2 shows the effect of the gradient index of the porous FG material ( ) on the dimensionless position of the elastic center C ( ̃ = /ℎ) of the rectangular cross section of the nano-beam.
Firstly, it can be noted that the dimensionless distance ̃ increases as the gradient index increases, with a maximum attained at a value of *, which depends on the value of ; note that for > *, the distance ̃ tends to vanish for → ∞. 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 ( = ) and dimensionless bending stiffness ( = ) 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 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.   The variations of the Euler-Young modulus across the thickness of a nano-beam are illustrated in Figure 5A, assuming equal to zero and varying in the set {0, 0.3, 0.5, 1, 3, 5, → ∞}. As was to be expected, when = 0, the FG material reduces to pure ceramic (E = Ec = 322.3 GPa) while, on the contrary, for → ∞, the material properties tend to pure metal (E = Em = 207.8 GPa).
The variations of the Euler-Young modulus through the thickness of the nano-beam cross section, assuming = 0.2, are illustrated in Figure 5B and 5C 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.

Governing Equation
Using the Bernoulli-Euler theory, the displacement components , and along x, y and z directions, respectively, at any point of the nano-beam, can be written as ( , , ) = ( , ) + * ( ), where ( , ) and ( , ) are the axial and transverse displacements of the elastic center C at time t, respectively; ( , ) = − ( , ) is the rotation of the cross section about the y-axis, and * ( ) 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]: As is well-known, the nonlinear equations of motion can be obtained using Hamilton's principle: where the expression of (kinetic energy), (strain energy) and (work done by external forces) are given in the following. The variation of kinetic energy is The expression for the variation of the strain energy is where the stress resultants and introduced in Equation (25) are defined as: The expression of the virtual work of the external force can be expressed by where = ( ) and = ( ) are the axial and the transverse vertical distributed loads, respectively, and ( ) is the Winkler elastic foundation coefficient.
By substituting the expressions of , , and into Hamilton's principle (Equation (23)), performing integration-by-parts with respect to t as well as x to relieve the generalized displacements , and of any differentiations, and using the fundamental lemma of differential calculus, we obtain the following equations of motion: The boundary conditions involve specifying one element of each of the following three pairs: − at = 0, .

Stress-Driven Nonlocal Integral Model
In adapting the stress-driven integral formulation to the Bernoulli-Euler porous FG nano-beam, the nonlocal axial strain, , and the nonlocal bending curvature, χ, are obtained by the following convolutions: where , are position vectors of points of the domain of the Euclidean space occupied by the nanobeam at time t, and the stress resultants and fulfil the equilibrium conditions. By assuming the following bi-exponential function for the averaging kernel Φ with internal characteristic length, the stress-driven integrals (Equations (34) and (35)) admit the following set of solutions: with x ∈ [0, L], if and only if the following two pairs of constitutive boundary conditions are satisfied at the nano-beam ends

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: By substituting Equations (43,44) into Equations (29,30), the following stress-driven equations of motion can be derived: with the following natural boundary conditions at the nano-beam ends: where , and 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)(39)(40)(41)(42).
According to the Bernoulli-Euler kinematics, the geometric bending curvature is related to the transverse displacement by ( , ) = ( , ).
If we now neglect the axial inertia term, , from the first equation of motion, we obtain wherein is a constant.
Note that for a nano-beam with immovable ends ( | = | = 0 and | = | = 0), by integrating both sides of Equation (56) over the domain [0, L] yields the following expression: which is the so-called "mid-plane stretching effect" in SDM theory. Based on this assumption, from Equation (55), it follows: 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.

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) 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 in the set {0.00 + , 0.01, 0.03, 0.05, 0.1}.
In the second comparison study (Table 3), 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 = 50, and assuming that and are equal to zero. The third example (Table 4) concerns the evaluation of the ratios between the local nonlinear frequency (ωNL) and the local linear one (ωL), for a perfect simply-supported homogeneous nanobeams, varying the dimensionless amplitude of the nonlinear oscillator /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 Tables 5 and 6 [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.  Table 3. Linear dimensionless natural frequencies of homogeneous simply-supported (S-S) and clamped-clamped (C-C) nano-beams ( = = 0; = 50).

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 = 2 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 nanobeam, ωcL, are presented.

Gradient Index and Porosity Volume Fraction
Firstly, the combined effect of the gradient index, , and porosity volume fraction, , on the frequency ratio of a perfect FG nano-beam ( = 0) are plotted in Figures 6-8 assuming = 0.2 and = 0. From these figures it can be observed that the frequency ratios evaluated for < 1 are always greater than those obtained for > 1. Moreover, the effects of and , for the two kinds of porosity distributions considered, on the aforementioned frequency ratio are plotted in Figure 9, varying the dimensionless elastic foundation coefficient, , in the set {0, 20}, and the initial imperfection amplitude, , in the set {0.0, 0.2}. In particular, from Figure 9A it can be observed that, for a perfect FG nano-beam ( = 0) with an even distribution of porosity, the curves of the aforementioned ratios (continuous lines) tend to increase as the porosity volume fraction increases when < 1, while an opposite trend is obtained when > 1. Moreover, it is found that the nonlinear frequency ratios of imperfect nano-beams ( = 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 increases, no matter what values of gradient index and elastic foundation coefficient are adopted.

Nonlinear Oscillator Amplitude
The nonlinear frequency ratio versus the nonlinear oscillator amplitude, , is presented in Figure 10, varying the elastic foundation coefficient and the porosity volume fraction . As can be observed, the ratio ωNL/ωcL increases as the absolute value of 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 ( = 0.5). It is also found that the nonlinear frequency ratio increases as the porosity volume fraction increases and that, for a given value of , the ratio ωNL/ωcL also increases as the coefficient increases.

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, , 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 = 0.2 and =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.

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 < 1 ; an opposite trend was observed when > 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).
Author Contributions: Conceptualization, R.P. and L.F.; formal analysis, R.P. and L.F.; Methodology, R.P. and L.F.; Validation, R.P. and L.F.; Writing -original draft, R.P. and L.F.; Writing -review & editing, R.P. and L.F.; All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A. Solution Procedure
From Equation (24) where ( ) is the unknown i-th time-dependent coefficient and ( ) is the i-th test function, depending on the assigned boundary conditions.
In this study, we assume the test function form to be equal to the local modal shape (i = 1) where the time base function, ( ), is given by the following approximate cosine solution where is the first nonlinear vibration frequency, is the amplitude of the nonlinear oscillator, and ( ) is assumed to be equal to the linear spatial mode based on the stress-driven non-local integral model [64] ( ) = + + + + + . (A10) In Equation (A6), , , are the roots of the characteristic equation, and ( ,…, ) are six unknown constants to be determined by imposing suitable boundary conditions. Note that, from the numerical point of view, the evaluation of the linear fundamental natural frequencies of a porous FG nano-beam consists of solving the eigenvalue problem expressed in terms of a six dimensional array, = { , … , } of the aforementioned unknown constants.
Moreover, the initial imperfection, * ( ), is assumed to have the following expression * ( ) = ( ), where is the amplitude of the initial imperfection shape.
where the expressions of the four coefficients , , and are summarized in Table A1, in which ( ) ( ) = .
Finally, in agreement with the Hamiltonian approach to nonlinear oscillators, the following expression of the nonlinear fundamental vibration frequency of an imperfect porous FG nano-beam can be derived as = . (A13) Note that the linear vibration frequency of an imperfect porous FG nano-beam can be determined from the previous equation by setting = 0, while the nonlinear vibration frequency of a perfect nano-beam can be obtained by setting = 0. Moreover, by setting = = 0, we obtain the local linear frequency, L, of a perfect FG nano-beam.