Nonlinear Vibration Analysis of Axially Functionally Graded Microbeams Based on Nonlinear Elastic Foundation Using Modified Couple Stress Theory

In this study, a non-classical approach was developed to analyze nonlinear free and forced vibration of an Axially Functionally Graded (AFG) microbeam by means of modified couple stress theory. The beam is considered as Euler-Bernoulli type supported on a threelayered elastic foundation with Von-Karman geometric nonlinearity. Small size effects included in the analysis by considering the length scale parameter. It is assumed that the mass density and elasticity modulus varies continuously in the axial direction according to the power law form. Hamilton's principle was implemented to derive the nonlinear governing partial differential equation concerning associated boundary conditions. The nonlinear partial differential equation was reduced to some nonlinear ordinary differential equations via Galerkin's discretization technique. He's Variational iteration methods were implemented to obtain approximate analytical expressions for the frequency response as well as the forced vibration response of the microbeam with doubly-clamped end conditions. In this study, some factors influencing the forced vibration response were investigated. Specifically, the influence of the length scale parameter, the length of the microbeam, the power index, the Winkler parameter, the Pasternak parameter, and the nonlinear parameter on the nonlinear natural frequency as well as the amplitude of forced response have been investigated.


Introduction
Experiments show that ignoring the internal length scale, which is the case in the classical continuum mechanics, in micro and nano structures can result in inaccurate structural predictions. To overcome this problem, a number of theories like couple stress theory and strain gradient theory [1][2][3][4][5] have been developed. In the couple stress theory, the effect of couple per unit area along with the effect of normal and shear forces is considered [6]. Yang et al. [7] developed a modified couple stress theory by use of a symmetrical couple stress tensor. Vibration analysis of micro-structures has attracted a lot of researchers in the last decade. Plenty of papers has been published on mircobeams [8][9][10][11][12][13], microbars [14,15] and micro plates [16,17] taking the length-scale parameters into account.
Functionally Graded Materials (FGMs), possessing spatial distributed properties, have been firstly developed for some special applications such as rocket engine components, aerospace structures, and turbine blades. The earliest FGMs were introduced by Japanese scientists in the 1980s, as high temperature resistant materials used in aerospace structures. Recently, these materials have got lots of applications in electrical devices, energy transformations, bio-engineering, and optics, among other fields [18]. Variation of mechanical properties through an FG beam makes its analysis difficult. Different techniques have been adopted to derive the equations and solving them to study FG beams. Chakraborty et al. [19] introduced a beam element to study the thermo-elastic behavior of FG beam structures. A unified approach was proposed by Li [20] for prediction of the static and dynamic behavior of FG beams with the rotary inertia and shear deformation included. Benatta et al. [21] presented a higher order theory for short FG symmetric beams under three-point bending. Aydogdu and Taskin [22] studied the free vibration of simply supported FG beams. They implemented Navier type solutions to obtain frequencies.
Şimşek [23] studied the nonlinear dynamics of an FG beam with a moving harmonic load. Kocaturk et al. [24] investigated on the non-linear static behavior of a cantilever Timoshenko FG beam under a non-follower transversal uniformly distributed load. Although FG structures are intensively studied by a lot of researchers, there exists limited number of publications on Axially Functionally Graded materials (AFG) in which the material constituents and the mechanical properties change axially. Huang and Li [25] studied the free vibration of non-uniform AFG beams. Shahba and Rajasekaran [26] studied the free vibration of Euler-Bernoulli AFG beams by means of differential transform method. Akgöz and Civalek [27] studied vibration behavior of a non-uniform AFG Euler-Bernoulli microbeam with various boundary conditions using the modified strain gradient and couple stress theories. Huang et al. [28] proposed a new method of analyzing the vibration of AFG Timoshenko beams with non-uniform cross-section. The size-dependent static and vibration behavior of micro-beams made of FGMs are analytically investigated on the basis of the modified couple stress theory in the elastic range by Asghari et al. [29]. The results showed that the static deflection and natural frequencies developed by the modified couple stress theory have a significant difference with those obtained by the classical beam theory when the ratio of the beam characteristic size to the internal material length scale parameter is small [29].
Recently, a number of researchers tried to study the nonlinear vibration analysis of FG structures. Moeenfard et al. [30] applied He's homotopic perturbation method to study the nonlinear free vibration of a Timoshenko microbeam. Ramezani [31] contributed to the field by examining the nonlinear free vibrations of a microscale Timoshenko beam employing the multiple scales method. Rajabi and Ramezani [32] introduced a microscale nonlinear beam model on the basis of strain gradient elasticity with surface energy. Bending response of non-homogenous microbeams embedded in an elastic medium was investigated by Akgöz and Civalek [33], based on modified strain gradient elasticity theory in conjunctions with various beam theories. It was observed that the effect of shear deformation may become more important for lower slenderness ratios. In addition, it can be said from the results that the small scale effect is more important when the thickness of the microbeam approaches to material length scale parameter [33].
In 2014, Bayat et al. [34] implemented an analytical variational approach for vibration analysis of an electrostatically-actuated microbeam. Şimşek [35] studied nonlinear free vibration of nano-beams based on Eringen's nonlocal elasticity theory using He's variational method. In another paper, he also studied the nonlinear static and free vibration of homogenous Euler-Bernoulli microbeams resting on a three-layered nonlinear elastic foundation within the framework of the modified couple stress theory [36]. A simple nonlocal finite element model was developed for buckling analysis of protein microtubules in [37]. Euler-Bernoulli beam model was used as continuum model of microtubules.
The results indicate that the effect of nonlocal parameter is to reduce the buckling load of microtubules for all type boundary conditions, in general [37].
Buckling analysis of linearly tapered micro-columns based on strain gradient elasticity studied by Akgoz and Civalek [38]. Bernoulli-Euler beam theory was used to model the non-uniform micro column. It was showed that the differences between critical buckling loads achieved by classical and those predicted by non-classical theories are considerable for smaller values of the ratio of the micro-column thickness (or diameter) at its bottom end to the additional material length scale parameters and the differences also increase due to increasing of the taper ratio [38].
Jia et al. [39] examined the effect of length-scale parameters on the nonlinear vibration behavior of an FG microbeam. Ebrahimi and Zia [40] studied the large-amplitude nonlinear vibration of FG Timoshenko beams made of porous material. Considering surface effect, Ansari et al. [41] conducted exact solution for the nonlinear forced vibration of FG nanobeams in the presence of surface effects. Additionally, Şimşek [42] studied the nonlinear free vibration of AFG microbeams with different boundary conditions based on the modified couple stress theory and Von-Karman's geometric nonlinearity. Farokhi et al. [43] studied the nonlinear dynamics of a geometrically imperfect microbeam numerically on the basis of the modified couple stress theory. They showed that the natural frequencies of the system increases with increasing either the initial imperfection or the dimensionless parameter η. The nonlinear analysis, on the other hand, revealed that the system exhibits both hardening and softening behaviour, depending on the amplitude of the initial curvature (imperfection) [43].
Asghari et al. [44] studied a size-dependent formulation for Timoshenko beams made of a FGM and showed that modeling beams on the basis of the couple stress theory causes more stiffness than modeling based on the classical continuum theory, such that for beams with small thickness, a significant difference between the results of these two theories is observed [44]. In [45] a microscale free vibration analysis of Composite Laminated Timoshenko Beam (CLTB) model was developed based on the newmodified couple stress theory. It was observed that the present beam model can capture the scale effects of micro-structures. The numerical results given by the CLTB model and the CLBB model showed that the natural frequency of the nonclassical beam model is always higher than that of the classical beam model [45].
Most recently, Shafiei et al. [46] investigated on nonlinear vibration of AFG tapered microbeams. For various boundary conditions, they derived the equations of motion and solved them via generalized differential quadrature method. Variation of first and second fundamental frequencies due to system parameters is studied. In [47] nonlinear vibration analysis of a simply supported AFG beam subjected to a moving harmonic load as an Euler-Bernoulli beam utilizing Green's strain tensor. The results indicate that these parameters have a considerable effect on both nonlinear natural frequency and response amplitude [47]. In another study, free and forced vibration analysis of FG doubly clamped micro-beams was investigated based on the third order shear deformation and modified couple stress theories [48]. The numerical results indicated that if the thickness of the beam is in the order of the material length scale parameter, size effects are more significant [48].
The present paper seems to be the first attempt to address the nonlinear free and forced vibration of an AFG microbeam on a nonlinear elastic foundation with doubly-clamped boundary condition, based on the modified couple stress theory and Von-Karman's geometric nonlinearity. A simple power law function is employed to predict the varying material properties of the AFG microbeam along its axis. The nonlinear governing equations and boundary conditions were derived using Hamilton's principle. Galerkin's method, He's variational method, and He's variational iteration method were employed to determine the nonlinear natural frequency and forced vibration response of the AFG microbeam.

Backgrounds 2.1 Problem definition
Consider a functionally graded doubly-clamped microbeam with geometric dimensions of length L, width h and thickness h, as indicated in Fig. 1. The microbeam is composed of an Axially Functionally Graded material including metal and ceramic phases based on a nonlinear elastic foundation. The beam properties vary linearly in the axial direction according to a simple power law function. This problem aims in finding the characteristics of nonlinear free and forced vibration, including natural frequencies and nonlinear forced vibration response. The effect of the length scale parameter, the length of the microbeam, the power index, the Winkler parameter, the Pasternak parameter and the nonlinear parameter on results are also investigated. A uniform harmonic load of q(x, t) q 0 sin ωt is applied on the upper surface of the beam in forced vibration analysis.

Modified couple stress theory
In the modified couple stress theory proposed by Yang et al. [7], the strain energy in an isotropic linear body occupying a volume Ω can be written as: where σ is the stress tensor, ε is the strain tensor, m is the deviarotic part of of the couple stress tensor, and χ is the symmetric curvature tensor. These tensors are given by [7]: where u is the displacement vector, λ , l is the material length scale parameter and θ is the rotation vector that can be expressed as: The displacement component of an initially straight beam on the basis of Euler-Bernoulli beam theory can be written as: where u and w are the axial and transverse displacements of any point on the neutral axis, respectively. The Von-Karman's nonlinear strain-displacement relationship based on assumptions of large transverse displacements, moderate rotation, and small strain for straight beams are given by [47][48][49]: where ε xx is the longitudinal strain. The components of the rotation vector can be achieved by using Eq. (6) as: Inserting Eq. (12) into Eq. (5) yields the following expression for non-zero components of the symmetric curvature tensor: The strain energy induced by the nonlinear elastic medium can be written as: where k L , k p , k N are the spring constants of the Winkler elastic medium, Pasternak elastic medium, and nonlinear elastic medium, respectively. The kinetic energy is given by: where: The work done by external force q(x, t) is given by:

Modeling Axially Functionally Graded materials
It is assumed that material properties of the microbeam, such as Young's modulus (E) and mass density (ρ) vary linearly in the axial direction according to a power law function which can be described by [39]: where Λ R and Λ L are the corresponding material properties of the microbeam at the right and the left end of the microbeam and n is a non-negative number that defines the material variation profile along the length of the microbeam. It should be noted that the material length scale parameter l in this study.

Derivation of governing equations of dynamic equilibrium
Hamilton's principle can be expressed as: where T is the kinetic energy, U s is strain energy, U el is the strain energy induced by the nonlinear elastic medium, and V is the work done by external applied forces. Substituting Eq. (1), Eq. (14), Eq. (15) and Eq. (18) into Eq. (22), the governing equations and the corresponding boundary conditions can be achieved as: The corresponding boundary conditions are defined at x = 0 and x = L as follows: As can be seen from above equation, the governing equations are coupled with respect to the displacements u and w. To reduce the equations to a single equation in terms of w, the in-plane inertia can be neglected [15,16]. Further, Eq. (25) can be rewritten as: Integrating Eq. (28) with respect to x gives: where c is the integration constant to be calculated with respect to the boundary conditions. It is assumed that the microbeam has fixed supports in the axial direction. Thus, the boundary conditions related to the axial direction can be expressed as: Integrating both sides of Eq. (30) from 0 to L cosidering above boundary conditions gives: The governing equation is derived in terms of w by substituting Eq. (31) into Eq. (24) as follows:

Solution of governing equations
In this section, governing equations are solved by means of the semi-analytical method of Galerkin. The displacement function w(x, t) can be expanded into a finite series as follows: where α(t) is an unknown time-dependent coefficient to be determined, and w x ( ) is the basis function which must satisfy the kinematic boundary conditions. Transverse mode shapes of the doubly-clamped beam are as follows: where β is maximum amplitude. Substituting the approximate solution in Eq. (33) into Eq. (32), then multiplying both sides of the resulting equation by w x ( ) and integration over the domain (0, L) yields the following equation: where dot sign stands for derivation with respect to time. The coefficients a 1 , a 2 and Q in Eq. (35) can be expressed as: where w w w , and w 4 ( ) are the first, second, third, and fourth derivatives of w with respect to the axial coordinate x, respectively.

Approximate analytical solution for nonlinear free vibration
The microbeam is subjected to the following initial conditions: where β is the maximum vibration amplitude of the microbeam. For free vibration analysis, Eq. (35) can be rewritten in a compact form as: Based on He's method [50], using the semi-inverse method, Eq. (42) can be expressed as: where T is the period of the nonlinear oscillator. It is assumed that the approximate solution can be expressed as [42]: where β and ω n are the amplitude and natural frequency of the nonlinear oscillator, respectively. Substituting Eq. (44) in Eq. (43), and considering the transformation θ = ω n t, one can obtain: J a a n n n β ω ω β ω θ β θ θ π , sin c os cos According to the Ritz method, the stationary conditions After some mathematical amendment, Eq. (46) takes the following form: The nonlinear natural frequency ω NL can be found by performing the integral expression in Eq. (47), as follows: Using Eq. (48), the following approximate solution can be found for α(t):

Forced vibration analysis
As it is stated in problem description, the beam is only under harmonic side loading in forced vibration analysis. where ω is the forced vibration frequency. According to the variational iteration method [50,51], the correction functional can be written as follows: where λ is a Lagrange multiplier [52], L is a linear operator, N is a nonlinear operator, and g is a known function.
The multiplier can be obtained as λ τ The following form is assumed as an approximate solution [37]: where γ( a 2 ) is a non-zero unknown function of a 2 with γ(0) = 1 [51]. The substitution of Eq. (52) into Eq. (35) results in the following residual: By the variational iteration formula (51), we have: Generally speaking, the residual Eq. (53) is not equal to zero. The right-hand side of this equation would vanish, if α 0 (t) is a solution of Eq. (35). We may ensure the vanishing of the factor of cos a t 1 γ ( ) by setting [51]: By implementing formula Eq. (47), we obtain:

Numerical results
In this section, free and forced vibration of a slender microbeam of height h = 10 −6 , width b = 3 × h, length L = 10 × h and l = 1 × h is presented. Mechanical properties of the desired FG material are shown in Table 1.
In order to verify the extracted formulation the FGM variability, foundation effect as well as size effect was set to zero in the developed code. The first natural frequency had a difference about 1 % compared to a classical doubly clamped aluminum beam.

Free vibration
The nonlinear natural frequency ω NL versus l for different values of power index is shown in Fig. 2. The default values of k NL , k P , k L are assumed 3, 2 and 1 N/m 2 , respectively. By increasing the value of n, the overall nonlinear natural frequency ω NL is decreased. Despite of the value of power index, an increasing trend is observed for the ω NL for higher length scale parameters. The variation of nonlinear natural frequency ω NL versus l for different values of L is shown in Fig. 3. As the value of L increases, the nonlinear natural frequency of ω NL is decreased. Fig. 4 depicts the influence of k L on the nonlinear natural frequency for the selected values of the length scale parameter. As seen from this figure, an increase in the value of k L leads to an increase in the nonlinear natural frequency. Fig. 5 shows the effect of pasternak medium stiffness on the nonlinear natural  frequency according to length scale parameter l. As the value of k P increases, the nonlinear natural frequency of ω NL is also increased. The variation of nonlinear natural frequency ω NL versus l is shown for different values of k NL in Fig. 6. Again, increasing the value of k NL makes the nonlinear natural frequency to increase.

Forced vibration
Time history of the maximum deflection W Max is shown for different values of l in Figs. 7 and 8. As the value of l increases, period times of system is reduced. Fig. 9  . As the value of n , is shown in Fig. 11.
As the value of l increases, the maximum deflection of the system is decreased.

Results
As shown in Figs. 2, 3, and 5, the length of the microbeam l, the power index n, and the spring constant k P had significant effect on the natural frequency of the nonlinear system. With small changes in these parameters, the nonlinear natural frequency of a nonlinear system underwent major changes. On the other hand, Figs. 4 and 6 showed that the effect of parameters k L (linear spring constant) and k NL (nonlinear spring constant) on the nonlinear natural frequency was not too much. Figs. 7 and 8 showed that increasing the material length scale parameter l reduced the periodicity of the system. The shapes in Figs. 9, 10 and 11 revealed that increasing the parameters n and L increased the maximum deflection, while increasing l reduced the maximum deflection.

Conclusion
In this study, a non-classical approach was developed to analyze nonlinear free and forced vibration of an Axially Functionally Graded (AFG) microbeam by means of modified couple stress theory. The beam is considered as Euler-Bernoulli type supported on a three-layered elastic foundation with Von-Karman geometric nonlinearity. Small size effects included in the analysis by considering the length scale parameter. Hamilton's principle was implemented to derive the nonlinear governing partial differential equation concerning associated boundary conditions. The nonlinear partial differential equation was reduced to some nonlinear ordinary differential equations via Galerkin's discretization technique. He's Variational iteration methods were implemented to obtain approximate analytical expressions for the frequency response as well as the forced vibration response of the microbeam with doubly-clamped end conditions. In this study, some factors influencing the forced vibration response were investigated. Specifically, the influence of the length scale parameter, the length of the microbeam, the power index, the Winkler parameter, the Pasternak parameter, and the nonlinear parameter on the nonlinear natural frequency as well as the amplitude of forced response have been investigated. The results show that: • The length of the microbeam L, the power index n, and the spring constant k P had significant effect on the natural frequency of the nonlinear system.
• The effect of parameters k L (linear spring constant) and k NL (nonlinear spring constant) on the nonlinear natural frequency was not too much. • Increasing the parameters n and L increased the maximum deflection, while increasing l reduced the maximum deflection.