Nonlocal and surface effects on nonlinear vibration response of a graded Timoshenko nanobeam

The free and forced vibration of a graded geometrically nonlinear Timoshenko nanobeam supported by on a nonlinear foundation is considered in this paper. The main contribution of this study is to propose a new formulation for the dynamic response of this beam by combining nonlocal and surface elasticity in addition to employing the physical neutral axis method which eliminates the quadratic nonlinearity from the equation of motion. Using the principle of virtual work, a fourth-order nonlinear partial differential equation is formulated and Galerkin technique is employed to yield a fourth-order ordinary differential equation with cubic nonlinearity in the temporal domain. The method of multiple scales is employed to obtain the analytical expression of the nonlinear frequency of the beam and its frequency response curve from a primary resonance analysis. To assess the accuracy of this analytical solution, it is compared with a numerical solution obtained using the differential quadrature method. The obtained analytical results are successfully validated for particular cases of the considered problem with results published by other authors. The effects of surface elasticity, nonlocality, the physical neutral axis, the beam aspect ratio, the power-law index and the elastic foundation coefficients on the free and forced vibration response of the graded Timoshenko nanobeam are thoroughly investigated for different types of boundary conditions .

does not have the capability to handle length scale effects. As a result, many models have been proposed during the last few decades that incorporated the size-dependency effects, namely nonlocal theory [5][6][7], modified coupled stress theory [8,9], strain gradient theory [10,11] or a combination of these theories such as the nonlocal strain gradient theory [12].
Functionally graded material (FGM) is a kind of composite material in which the volume fraction of the constituent materials changes from one surface to the other. Thus, FGM provides a smooth transition change in properties in a specified direction that leads to less stress concentration and consequently to fracture [13]. FGMs have been used in several engineering applications such as marine risers [14], biomedical applications [15], energy applications [16] and aerospace applications [17], and a detailed review can be found in [18]. In addition, FGMs have been recently employed in several micro/nanoapplications [19][20][21][22].
Using Eringen's nonlocal model, several investigators have focused on the linear dynamic response of functionally graded (FG) nanobeams. Eltaher et al. [23] performed a free vibration analysis of a nonlocal nanobeam using the finite element method (FEM). Navier's method was used by Uymaz [24] to investigate the dynamic response of a FG Euler-Bernoulli nanobeam. The same method was employed by Rahmani and Pedram [25] to examine the vibration behavior of a Timoshenko graded beam. A free vibration study was conducted by Nejad and Hadi [26] on a nanobeam graded in two directions in which they employed the generalized differential quadrature method (GDQM) to solve the equations of motion. Yao et al. [27] investigated the free vibration and wave propagation of an axially moving functionally graded Timoshenko microbeam. Generally, the nonlocal model induces a softening type of behavior; however, Li et al. [28] demonstrated that a cantilever Euler-Bernoulli beam can exhibit either a softening or a hardening type of behavior depending on the applied loading. Recently, Arefi et al. [29] conducted a free vibration analysis of graded polymer composite curved nanobeams reinforced with graphene nanoplatelets to predict their dynamic characteristics. Very recently, Luo et al. [30] conducted transverse vibration of an axisymmetric graded circular nanoplate with radial uniform loads. Esen et al. [31] investigated the vibrational characteristics of functionally graded (FG) cracked microbeam rested on elastic foundation and subjected to thermal and magnetic fields. On the other hand, the literature is abundant with studies using the strain gradient theory such as, for instance, the works among others of Arefi and Zenkour who performed dynamic and bending analyses for the sandwiched microbeams [32][33][34] and FG Timoshenko's sandwiched microbeams [35]. There are also several published papers which used the modified couple stress theory in studying the dynamic behavior of Timoshenko beams such as the recent work of Abdelrahman et al. [36] who investigated the dynamic behavior of perforated Timoshenko microbeams in a thermal environment and that of Abo-bakr et al. [37] who studied the static and dynamic behavior of axially graded microbeams with nonuniform cross section. The recent literature witnessed additional papers related to the nonlocal strain gradient theory. For instance, Daikh et al. [38] studied the static stability, free vibration and bending response of multilayer functionally graded carbon nanotubes reinforced using composite nanoplates.
Many authors studied the nonlinear dynamic response of graded nanobeams by considering the von Kármán geometric nonlinearity. He's variational method was employed by Simsek [39] to predict the nonlinear vibration response of a Euler-Bernoulli graded nanobeam. A similar method was used by Niknam and Aghdam [40] to determine the nonlinear response of a free vibration and buckling load of a Euler-Bernoulli graded beam supported by nonlinear elastic springs. Using the method of multiple scales (MMS), the nonlinear nonlocal free vibration response of a Euler-Bernoulli FG nanobeam was investigated by Nazemnezhad and Hosseini-Hashemi [41]. MMS was also used by El-Borgi et al. [42] who examined the nonlinear nonlocal free and forced vibration of a Euler-Bernoulli FG nanobeam supported by elastic foundation effects. A similar problem was solved by Trabelssi et al. [43,44] by employing both the locally adaptive differential quadrature method (LaDQM) and the regular differential quadrature method (DQM). These studies were extended by the same main authors to investigate the free and forced vibration response of a Timoshenko graded beam based on the weak quadrature element method [45].
In the last few years, several authors have been investigating the incorporation of surface effects on the behavior of nanostructures. These effects can become quite important due to abrupt changes in the surface area-to-volume ratio at the nanoscale dimensions. Overall, it was observed that the surface elasticity properties can influence the properties of the size-dependent materials at nanoscale [46][47][48][49][50]. The nonlinear response of the nonlinear free vibration of FG nanobeams while considering only surface effects was studied by [51]. The nonlinear vibration frequencies of a Timoshenko nanobeam were investigated by [52,53] in which the combined effects of nonlocal parameter, surface elasticity and residual surface stress were taken into account. The nonlinear vibration response of graded nanobeams was examined by Hosseini-Hashemi et al. [54] by emphasizing on the coupling between nonlocal and surface effects. In the above literature, it was assumed that the geometrical neutral axis coincides with the physical neutral axis for the undeformed plane of the beam. However, this assumption is valid only for uniform material properties along a given direction of a beam. This assumption was shown to result in an error up to 10% in the natural frequencies of a FG beam [55]. Recently, few researchers have recently adopted in their models the physical neutral axis in FG beams instead of using the geometrical central axis [55][56][57][58][59][60][61]. Lately, Arefi et al. [62,63] employed the concept of the physical neutral axis to study, respectively, the static behavior of graded piezoelectric plates and the vibration response of graded piezoelectric face-sheets. More recently, Shen et al. [64] investigated the vibration and stability response of graded nanoplates by considering physical neutral plane. The use of such technique can eliminate stretching-bending coupling effect due to nonuniform material distribution along the given direction in the beam.
From the above review and to the best of the authors' knowledge, it can be summarized that very few researchers have investigated the nonlinear free and forced vibration response of FG Timoshenko nanobeams by considering both nonlocal and surface effects [48] or surface effects only [49,51]. However, combining the aforementioned effects along with the physical neutral axis location for predicting the nonlinear dynamic behavior of FG beams has not been investigated thus far for deep beams. It was, however, investigated by the authors [65] for slender Euler-Bernoulli beams where shear effects were neglected. The effect of the neutral axis is expected to be more pronounced for a deep beam especially when shear effects are considered. For low aspect ratio, the effect of the neutral axis may cause further deviation from the standard model. To bridge these gaps, this study aims to introduce the physical neutral axis concept while considering both surface elasticity and nonlocality to study the nonlinear dynamic response of a deep Timoshenko beam by employing the method of multiple scales (MMS). Such a solution is based on the linear classical modeshape Galerkin technique. A DQM-based numerical solution proposed by Trabelssi et al. [44] is employed based on the actual nonclassical modeshape for the purpose of assessing the accuracy of the MMS solution. Therefore, the novelty of this study can be summarized as follows: (a) propose a new formulation for the dynamic response of a geometrically nonlinear TBT beam with a nonlinear foundation combining nonlocal and surface elasticity in addition to physical neutral axis, (b) develop an MMS-based analytical solution for the current problem and validate it using DQM, (c) explore low aspect ratios of the beam while remaining within the assumptions of the TBT beam theory, (d) study the effect of surface stress on the neutral axis position, (e) conduct a thorough and detailed parametric study covering all the nondimensional parameters governing the physics of the investigated problem. The proposed formulation can be potentially applied to several applications such as MEMS and NEMS devices [66] in addition to carbon nanotubes (CNTs) such as sensors, actuators, electromagnetics, electroacoustic and hydrogen storage [67] as well as protein microtubules surrounded by an elastic matrix [68].
This paper is structured as follows. Sects. 2 and 3 provide, respectively, the derivation of the motion equations for both the classical and the nonlocal Timoshenko beam theory by considering surface elasticity and the physical neutral axis location. Section 4 summarizes the solution formulation of vibration response of the nanobeam using MMS. Section 5 outlines the DQM formulation used to assess the accuracy of the MMS solution. In Sect. 6, the obtained analytical results are validated with published data for specific cases of the considered problem and then a detailed parametric study is presented to examine the effect of the various parameters on the free and forced vibration nonlinear response of the nanobeam. Finally, a summary of this investigation and main conclusions is provided in Sect. 7.

Classical Timoshenko beam with surface effects
A graded nanobeam is considered in this study whose length, width and height are denoted L, b and h, respectively. The nanobeam is resting on a nonlinear elastic foundation as shown in Fig. 1. The considered boundary condition is either simply supported (S-S) or clamped-clamped (C-C). Due to the characteristics of the graded material, the beam's mid-plane does not match with the neutral axis of the beam, and hence, two separate frames of reference need to be considered in this study: • (x m ,y m ,z m ) is a reference coordinate system with the axis x m running through the geometrical mid-plane axis of the beam and y m and z m are two axes of symmetry defining the cross-sectional area of the beam, as depicted in Fig. 1. • (x,y,z) is a secondary reference coordinate system with the axis x = x m running through the neutral axis of the beam and y = y m and z = z m − h 0 are two axes of symmetry as shown in Fig. 1. Here, h 0 is the distance separating the neutral axis to the mid-plane of the nanobeam. For a Timoshenko beam with varying mechanical properties varying along z, it is assumed that the elastic and shear moduli, E and G, the surface elastic modulus E s and the density ρ follow a power-law function as given below where L and U denote, respectively, the lower and upper layers of the nanobeam. ν is the Poisson's coefficient assumed to be a constant [69]. The parameter n denotes the material gradation index. According to the Timoshenko beam theory (TBT), the displacement field can be expressed as follows: where u and w are, respectively, the displacements along x and z coordinates, φ is the shear angle, andê x andê z are, respectively, unit vectors along the x-and z-axes. Here, t designates the time variable. Within the assumptions of the small deformation and displacement theory, the strain tensor can be expressed as follows: where In a Hookean solid, the stress-strain law for a Timoshenko beam can be written as where K s is the shear correction factor with a value of 5/6 for a rectangular section. The beam surface elasticity along the longitudinal direction is represented by [48] It is worth emphasizing that the beam is assumed to be geometrically nonlinear (i.e., with van Karman strain nonlinearity appearing in Eq. (4)) but with a linear stress-strain behavior (i.e., Eq. (5)) and a nonlinear foundation.
Applying the principle of virtual work for the Timoshenko beam and integrating its expression by parts and separating the different components of the displacement field yield the following equations of motion: x x and M (0) xz are stress resultants and m 0 and m 2 are mass inertia which can be expressed as in which A = bh denotes the area of the beam's cross section. By employing the expressions from (2) to (6), the resultants of the stress components can be formulated as a function of the displacements The neutral axis location is derived by setting that the total axial force is equal to zero along the longitudinal axis. Therefore, Eq. (9a) can be rewritten as To determine the neutral axis location, the axial displacement is ignored and higher-order terms are neglected (i.e., first term of the above equation is neglected). Thus, Eq. (11) can be reduced to a closed-form expression for the position of the neutral axis Details about the calculation of h 0 can be found in Appendix A, and the expressions ofÃ,D andG appearing in (10) are written in a simpler form in Appendix B.
To combine the motion equations (8a) and (8b), a relation between the axial and transverse displacements is developed using the von Kármán nonlinearity. Based on a method proposed by [70], the displacement along the axial direction u is eliminated from the motion equations resulting in a constant von Kármán nonlinear term. To this end, the following assumptions are made in the study: (1) The inertia term (mü) is very small and can be ignored; (2) the beam is supported at x = 0 and x = L, and hence, the displacements along the longitudinal direction at the boundaries are zero, i.e., u(0, t) = 0 and u(L , t) = 0. The first assumption along with (8a) results in M x x being a constant. In light of that, the following equations of motion can be obtained after further mathematical manipulations: Here, dot ( . ) represents differentiation with respect to t and prime ( ) denotes differentiation with respect to x. Details of this simplification are provided in Appendix C. The force f z in Eq. (13a) comprises both the applied force and the elastic foundation reaction force, and its expression can be written as [71] where k L , k N L and k s represent, respectively, the linear, nonlinear and shear coefficients of the elastic foundation and F is the magnitude of the harmonic force. Finally, the following normalization is adopted: where r = I A is the radius of gyration of the cross section and I = bh 3 12 is the second moment of inertia. The equation of motion in dimensionless form in terms of the nondimensional parameters is given by

Nonlocal Timoshenko beam theory with surface effects
Eringen [5,6] introduced the nonlocal theory which states that the stress σ at a given point x in an elastic domain is not only a function of the local strain at that point ε but is a function of the strains at all points in the domain. Based on this theory, the nonlocal stress can be expressed as follows: where σ (x ) is the local stress at point x and K (|x − x|, τ ) is the nonlocal modulus, |x − x| denotes the distance between x and x', and τ is a material constant. Eringen [6] also introduced an equivalent form of the differential model which can be written as where e 0 is a material constant, a is an internal characteristic length, is an external characteristic length and μ 0 is the nonlocal parameter, and ∇ 2 is the Laplacian operator. In addition, for a beam-type structure, it is reasonable to assume that the size-dependent effect is dominant along the axial direction. Thus, Eq. (19) can be simplified to the following:σ where the operator ∇ 2 is substituted with ∂ 2 /∂ x 2 andÃ,D andG are given by (B.1) and the expressions of A andD incorporate surface effects. Integration of the above relations over the section of the beam yields the nonlocal stress resultants as a function of the strain components which can be written as Substituting (2) and (4) into the above equations yields the nonlocal stress resultants as a function of displace-mentsM By employing the principle of virtual displacement, the following nonlocal equations of motion are obtained [44]: Taking the derivative of the above equations with respect to x and substituting the expressions The assumptions listed in Sect. 2 are used to remove the axial displacement u from (24a) and yield the stress resultants as a function of the vertical deflection w and the shear angle φ. The above equations reduce to the following:M where the expression of c 1 (t) is provided in (C.3). Substituting the above equations into (23b) and (23c) yields the nonlocal equations of motion of the nanobeam The dimensionless parameters introduced in equation (15) are substituted in the above equations to yield the following nondimensional equations of motion: wheref z =k Lŵ +k N Lŵ 3 −k s w +F cos(θt) andμ 0 = μ 0 L 2 is the dimensionless nonlocal parameter and the remaining quantities have been defined in (17). Whenμ 0 = 0, then the nonlocal equations of motion of the nanobeam reduce to the equations of motion given by (16).

Galerkin discretization
It is possible to combine the above equations of motion system into a single partial differential equation in terms of onlyŵ. To this end, the expression ofφ is extracted from (27a) and then injected into thex derivative of Eq. (27b). The resulting equation can be written as By employing Galerkin technique, the transverse displacement is assumed to be written asŵ(x,t) = ψ(x)q(t) and the resulting equation (28) can be transformed into an ordinary differential equation with time being the independent variable. The term ψ(x) is expressed as a classical linear beam modeshape. Table 1 lists these linear modeshapes corresponding to the either simply supported (S-S) or clamped-clamped (C-C). The aforementioned assumption is substituted into (28), and the resulting equation is integrated over the beam domain Table 1 Linear classical modeshape functions of a uniform beam for different boundary conditions where the coefficients β i (i = 1, 5) and μ i (i = 1, 4) are defined in Appendix D and μ 1 and μ 2 contain the damping coefficientĉ. In the above equation, is the amplitude ofF(x), and ψ F (x) is a function used to describe the applied force distribution along the x-axis. However, in general, this distribution is set as ψ F (x) = 1.

Free vibration solution
To study the free vibration response of the nanobeam, the damping term and the external force are removed from (29). Adding the initial conditions gives where A represents the amplitude of vibration. An approximate analytical solution for equation (30a) is achieved using the method of multiple scales [70]. The solution is obtained through an expansion of both the time variablet and the dependent variable q(t). This expansion requires that the system response is expressed as a function of several independent variables called time scales. The solution is expressed in terms of this expansion as follows: where ε is a small scaling term. The fast time scale T 0 =t is defined as the time scale of the primary oscillation, and the slow time scales T n = ε nt , n ≥ 1, are the times used to account for amplitude and phase modulation. Substituting (31) into (30a) and then gathering the powers of ε result in a system of three differential equations where D i (i = 0, 1, 2) indicates a differentiation with respect to time and the subscript i denotes the scale of time. The solution of Eq. (32a) may be written as where cc designates the complex conjugate of the first two terms and ω 1 and ω 2 are linear frequencies which are given by Substituting the solution given by (33) into (32b) and then setting the secular terms e iω 1 T 0 , e iω 2 T 0 to zero lead to the following conditions: which simplify to Then, substituting the solutions of (32a) and (32b) and applying the above conditions given by (35) into the differential equation (32c) then dropping the secular terms e iω 1 T 0 , e iω 2 T 0 , the following equations are obtained: Here, A 1 and A 2 are given in polar forms as where a j and θ j ( j = 1, 2) are reals. Substituting A 1 and A 2 into (37a) and (37b) and then separating into real and imaginary parts, it can be concluded that a 1 and a 2 are constants and θ 1 and θ 2 are given by In the above equations, θ 10 and θ 20 are constants and δ 1 and δ 2 are given by Substituting (39) into (38) yields the following expressions of A 1 and A 2 : Here, T 2 = ε 2t . Substituting A 1 and A 2 into (33), the overall solution given by Eq. (31) becomes where the nonlinear frequencies j ( j = 1, 2) are given by To obtain the unknown constants a 1 , a 2 , θ 10 and θ 20 , solving Eq. (42) by applying the initial conditions provided in (30b) and gathering like powers of ε ε (a 1 cos (θ 1 ) + a 2 cos ( The solution of the above system of equations can be written as Finally, the nonlinear frequencies j ( j = 1, 2) can be expressed as Both eigenvalues 1 and 2 exist to be relative to one single modeshape, and this was also observed by Majkut [72].

Forced vibration solution
To study the forced vibration response of the nanobeam, the equation of motion (28) is considered with both external harmonic forcing function and damping coefficient and its solution is sought using the MMS. For the primary resonance, i.e.,θ ≈ ω 1 , assume that the excitation frequencyθ and the system linear frequency ω l are approximately the same. The following equation can be expressed aŝ where σ represents the detuning parameter. An approximate solution is obtained for the ordinary differential equation (ODE) given by (29) using Eq. (31) to derive the n th -order differential equations. The only difference with free vibration analysis is that the damping terms μ 1 and μ 2 and the dimensionless amplitude of the forcē F 1 are scaled to the third order. Substituting (31) into (29) and gathering like powers of ε yield the a system of three differential equations Evidently, Eqs. (48a) and (48b) are, respectively, analogous to the case of free vibration equations (32a) and (32b) and hence have identical solutions. Substituting the solutions of equations (48a) and (48b) into the third-order equation (48c) results in the following secular terms: Using the expression of A 1 and A 2 provided by Eq. (38) and splitting into imaginary and real components yield a set of equations in which the unknowns are a 1 and θ 1 A similar system of differential equations can be obtained for a 2 and θ 2 which can be written as , a steady-state condition is sought. Hence, the frequency response equations can be expressed by setting all derivative terms, in Eqs. (50a), (50b), (51a) and (51b), to zero. The resulting equations can be written as From Eq. (53c), it can be concluded that a 2 = 0. In light of this, Eqs. (53a) and (53b) can be further simplified By taking the square and summation of Eqs. (54a) and (54b), the following frequency response equation is obtained:

Numerical solution using DQM
The MMS solution was developed based on linear classical modeshapes defined in Table 1. To assess the accuracy of such assumption, a numerical differential quadrature method (DQM) solution proposed by Trabelssi et al. [44] based on the actual nonclassical modeshapes is employed. The DQM formulation is briefly reviewed in this section. The dynamic response of the nanobeam is governed by two nonlinear partial differential equations (PDEs) (27a) and (27b), each subjected to two boundary conditions. DQM is a high-order method which leverages all the mesh points in the system to interpolate and calculate the derivatives involved in the PDE system. DQM and its related methods have been used several times to analyze the nonlinear response of structures at the nanoscale [43,44,[73][74][75][76]. Because of its efficiency, DQM was used to discretize the equation of motion in a weight optimization of axially functionally graded nonlinear microbeams [76]. Differentialintegral quadrature method (DIQM), another variant of DQM, was exploited to investigate the buckling and post-buckling behaviors of nonlinear carbon nanotubes [75]. Lately, Trabelsi at al. developed a high-order finite element for a nanobeam using DQM matrices to study the dynamic response of a nonlinear nanobeam [45]. In this section, these PDEs are discretized using DQM based on the shifted Chebyshev-Gauss-Lobatto grid points [77] where n x is the number of nodes. Let Y 1 and Y 2 be the nodal displacement vectors where the element of each vector is The velocity and acceleration vectors are denotedŸ andẎ , respectively. The discretized equations of motion (27a) and (27b) is given by following system of ODEs: where the m th DQM derivative matrices M m are given by [77]. Note that M 0 is the identity matrix. iW is a vector integral operator, and its expression can be given in [44]. Here, the matrix product of M 2 and Y is denoted M 2 .Y and the scalar product of two vectors is iW.Y . M 2 M 1 or M 2 1 is the Hadamard product. Only S-S and C-C boundary conditions are considered in this study. Their discretized forms are, respectively, given by  6 Numerical results and discussion

Validation studies
The first validation study is performed to check the correctness of the nonlocal nonlinear frequency ratio, which is defined as follows: Nonlocal nonlinear frequency ratio = Nonlocal nonlinear natural frequency Classical nonlinear natural frequency (60) For the S-S case, the present study results are given in Table 2 along with those reported in references [78] and [79] for different values of amplitude parameters A and gradation index n. The geometrical configuration of the beam such as length L, breadth b and height h is 10nm, 1nm and 1nm, respectively. The beam material properties used in this study are the same as those of [78]. This validation study does not consider the elastic foundation effects. The obtained results are in agreement with the published results as reported in Table 2.
In the previous validation, it is assumed that the neutral axis and the geometrical axis coincide with each other. The second validation is performed to check the correctness of the linear frequency based on the above assumption. The natural frequency was estimated for various values of the ratio E U /E L and for two values of the material gradation index n = 0 and 4. The properties of the beam are E L = 210 GPa and ρ = 7800kg/m 3 and are the same used by Eltaher et al. [55]. The beam length is L = 100nm, and the width and height of its cross section are chosen such that b = 0.1L and h = 0.01L. The present study results are reported in Table 3 and found to be in good agreement as compared with those of reference [55].
The previous validations were conducted on EBT nanobeams with a length-to-thickness ratio equal to 10. The purpose of the third validation is to assess the accuracy of the nondimensional nonlinear nonlocal frequency for a deep S-S TBT nanobeam with a length-to-thickness ratio varying between 5 and 8. The results obtained using MMS are tabulated in Table 4 and are compared with those published by Refs. [44,45] using DQM for different values of the inhomogeneity index n, the vibration amplitude A and the nonlocal parameter μ 0 . Excellent agreement is achieved between the obtained and published results.

Free vibration response
In this section, the effect of varying the aspect ratio L/ h, material inhomogeneity n, nonlocal parameter μ 0 , elastic modulus ratio E U /E L , amplitude A, surface elastic modulus ratio E s U /E s L and elastic foundation stiffness coefficientsk L ,k s andk N L on the nonlinear frequency of the beam is studied in this section. The results are obtained for S-S and C-C boundary conditions using MMS and accounting for the physical neutral axis and using both MMS and DQM. The parametric study is conducted by considering the graded nanobeam with a beam length of 10 nm. Based on the power-law model given by (1), the bottom surface is enriched with aluminum material properties with an elastic modulus E = 70 GPa, density ρ = 2700 kg/m 3 and surface elastic modulus E s =5.1882 N/m. The modeshape functions used to obtain the results are provided in Table 1. Table 3 Comparison of the nondimensional linear natural frequency obtained for the simply supported beam with shifted neutral axis (NA) with results obtained by Eltaher et al. [55] for different elastic modulus ratios (μ 0 = 0, L/ h = 100, ρ U /ρ L = 1 and  Figure 2a, b depicts the influence of varying the inhomogeneity index n from 1 to 100 and the amplitude A from 0 to 1 on the nonlocal nonlinear frequency versus the beam aspect ratio for the cases where the normalized nonlocal parameterμ 0 is equal to 0.04. The elastic modulus ratio E U /E L was chosen to be 2 and 10 which correspond to the critical cases to investigate the influence of the aspect ratio of the beam on the nondimensional frequency. The upper and lower two plots correspond, respectively, to the S-S and C-C boundary conditions. A substantial variation of the nondimensional frequency in the low aspect ratio range can be observed especially below the value of 10 and regardless of the value of the ratio E U /E L . As expected, a large value of E U /E L amplifies the effect of the inhomogeneity index n. Furthermore, it can be seen that the higher the aspect ratio, the lower the impact of the inhomogeneity index. Such observation is expected since low aspect ratio beams provide higher moment contribution for the unevenly distributed material. In addition, it can be noted that the higher the value of the amplitude, the higher the effect of the inhomogeneity. Moreover, the effect of the inhomogeneity index n seems to be strongly related to E U /E L ratio and this is true for both S-S and C-C cases. Indeed, given the nature of the material distribution, it is expected that raising n does not always move the nonlinear frequency toward the same direction. This is clearly illustrated when observing the case of E U /E L = 2 and A = 1 where the frequency curves order does not follow their inhomogeneity index. Furthermore, when E U /E L is raised to a value of 10, the curves presetting the same amplitude but different inhomogeneity index n are affected differently. This can probably be attributed to the fact that for each value of E U /E L , there is a different value of n that maximizes or minimizes the position of h 0 to yield a larger effect of E U /E L . In the linear case, where A = 0, the frequency curves split at a value of 10. However, increasing the amplitude drives the splitting value further away. In general, the effect of the amplitude appears to be highly nonlinear. It is also interesting to explore the influence of considering the physical neutral axis. For this study, the data used to generate Fig. 2a, b, are recreated without considering the physical neutral axis (i.e., the neutral axis is positioned in the center line of the beam). The resulting frequencies are presented in Table 5 as percentage of their corresponding frequencies had the physical neutral axis been adopted. Furthermore, the same data are regenerated while considering a high surface stress ratio. Given the results in Table 5, it can be observed that for a low value of E U /E L , the effect of the position of the neutral axis is small. In fact, this variation reaches a maximum at 99.36% for the S-S and 97.39% for the C-C case for the lowest aspect ratio. However, when E U /E L ratio is high, the difference reaches 9% for the C-C case. Interestingly, it can be seen that raising the amplitude or lowering the aspect ratio of the beam widens the variation of the natural frequency which emphasizes the importance of considering the neutral axis especially for the nonlinear results. When exploring the data where the surface stress is accounted, up to a 12% deviation is recorded between frequencies computed with and without accounting for the physical neutral axis . This further highlights the importance of accounting for surface stress when calculating the position of the neutral axis. Figure 3a,b depicts the effect of varying the amplitude A from 0 to 1 and the nonlocal parameterμ 0 from 0 to 2 on the nonlocal nonlinear frequency versus the inhomogeneity index n for the cases where the elastic moduli ratio E U /E L is equal to 2 and 10, respectively. The beam aspect ratio L/ h was chosen to be 6 to represent a deep beam. The upper and lower two plots correspond, respectively, to the S-S and C-C boundary conditions. Unlike the amplitude A,μ 0 does not seem to alter significantly the shape of the frequency curves. It is also observed that the variation of the frequencies registers a minimum especially for high values of the amplitude A. A higher value of E U /E L not only boosts the effect of the inhomogeneity index but also moves the position of the minimum. It can also be noted that raising the value ofμ 0 not only results in a softening but also increases the effect of the amplitude A as the plot of the same line type is more spread for higher values ofμ 0 . The frequency curves of C-C case seem to be less affected by the variation of the inhomogeneity index.
Here, the correlation between the effect of n and E U /E L can clearly be seen when observing the location of the observed minima which is slightly affected byμ 0 and heavily affected by the elastic moduli ratio E U /E L . These results are in agreement with the observations presented earlier. Figure 4a, b describes the effect of varying the inhomogeneity index n from 5 to 20 and the amplitude A from 0 to 1 on the nonlocal nonlinear frequency versus the elastic moduli ratio for the cases where the beam aspect ratio L/ h is equal to 4 and 8, respectively. The nonlocal parameterμ 0 was chosen to be equal to 2. The upper and lower two plots correspond, respectively, to the S-S and C-C boundary conditions. An aspect ratio L/ h ranging between 4 and 8 is chosen. This should lead to a more pronounced effect of E U /E L . The reduction of the value of L/ h does not seem to affect the general evolution of the frequency plots. However, it does highlight that thicker beams are more sensitive to the variation of E U /E L , especially for the C-C cases. Both the ratio E U /E L and the inhomogeneity index n seem to alter the values of the natural frequency in different ways resulting in the occurrence of minima in few curves. Raising the value of the amplitude A mainly exaggerates these effects. This is clearly illustrated in Fig. 5 which shows a rapid increase in the frequencies as A increases especially for the S-S case. Figure 6a, b represents the effect of varying the inhomogeneity index n from 5 to 20 and the beam aspect ratio L/ h from 4 to 8 on the nonlocal nonlinear frequency versus the elastic surface moduli ratio E s U /E s L for the cases where the amplitude A is equal to 0 and 1 and the elastic moduli ratio E U /E L is equal to 2 and 10. Modifying the value of E s U /E s L changes the stiffness distribution across the beam in a different way than E U /E L . As indicated in Appendix B, while modifying E U /E L affects the shear stiffnessD, modifying the surface stress ratio E s U /E s L does not. On the other hand, it can be seen that the expression of κ 0 given by (17) indicates that the geometric nonlinearity is directly affected by the surface stress. Figure 6a,b shows significant variation in the frequency as the ratio E s U /E s L increases. This effect tends to be more pronounced for low aspect ratio regardless of the boundary conditions. This effect is slightly influenced by the vibration amplitude especially for the S-S case. In Fig. 7a,b, the natural frequencies for different surface stiffness ratios, amplitudes, inhomogeneity indices are plotted as a function of the aspect ratio. It can be observed that a low aspect ratio greatly enhances the effect of the surface stress ratio for both C-C and S-S cases. It is important to note that the surfaces stresses have higher lever arm to the neutral axis in deep beams, thereby a higher contribution to the moment. Furthermore, the expression of h 0 is dependent on the surface stresses values. These two facts may explain the importance of surface stress in deep beams. A crossing in the frequency is also observed as the aspect ratio is raised and its location seems to be affected by the vibration amplitude A. As stated earlier, this can be traced to the expression of κ 0 . This study shows the importance of modeling surface stresses in an FGM beam using Timoshenko beam theory. In fact, since Timoshenko beam theory accounts for shear effects, it is now clearer that the effect of surface stress gradient and gradient elasticity is different than for the Euler-Bernoulli beam theory. This effect is more pronounced for deep beams which further highlights the importance of using Timoshenko beam theory.
The effect of the elastic foundation under different boundary conditions on the nonlinear frequency is illustrated in Fig. 8a, b with and without the nonlinear stiffness coefficientk N L , respectively. The plotted results are based on the current MMS formulation in addition to the DQM formulation proposed by Trablelssi et al. [44]. It can be observed that there is a slight increase in the frequency with the increase in amplitude. The net effect of the shear stiffness coefficients dominates the nonlinear frequency. The set of values chosen fork L andk S in Fig. 8 is chosen based on the expression of β 2 and β 3 in Appendix D.k L andk S only appear in the expression of β 2 and β 3 . It can be seen that the effects ofk L andk S are comparable except thatk S is multiplied by the second space derivative of the weighted sum of the shape function and few of its space derivative. A rough calculation using the linear modeshapes for both S-S and C-C case shows that the effect ofk S should be about 10 times the effect ofk L for S-S case and around 20 times for the C-C case. This observation can be easily detected in the MMS results since linear modeshape are used for the evaluation of the results. However, while the DQM and MMS produce comparable results for the S-S case, the results for the C-C case are vastly  different for both methods. This is most likely due to the fact that for the C-C case the modeshape can no longer be approximated with a linear modeshape. A comparable observation was made by Trabelssi et al. [43] on a similar problem. It can be concluded that these errors occur mainly because MMS assumes the first modeshape function for both S-S and C-C cases which are classical modeshape functions as given in Table 1, but DQM employs a nonclassical modeshape function. The data presented in Fig. 8 can be used to examine the effect of the nonlinearities on the beam response which is either caused by the van Karman strain nonlinearity or the nonlinear stiffness coefficientk N L or a (a) (b) Fig. 7 Nonlocal nonlinear frequency versus aspect ratio for different values of the surface elastic moduli ratio, amplitude, elastic moduli ratio and inhomogeneity index a n = 5 and b n = 20 for both S-S and C-C cases combination of both. As stated earlier, this can be achieved by comparing the frequencies obtained for an amplitude equal to and different from zero corresponding to the linear and nonlinear behavior, respectively. Whenk N L is set to zero, only the geometric nonlinearity is responsible for the dependence of the frequency on the amplitude. Compared to the linear case, the geometric nonlinearity results in a slight increase in the natural frequency when the vibration amplitude is raised. However, whenk N L is set to 100, the effect of the nonlinear foundation increases the dependence of the natural frequencies on the amplitude making it further depart from the linear case.

Forced vibration response
The proposed forced vibration study investigates the effects of the aspect ratio, the material inhomogeneity, the nonlocal parameter, the foundation coefficients while considering the effect of surface elasticity on the frequency response curve (FRC) of the nanobeam. The damping coefficient is taken asĉ = 0.1, and the force amplitude is selected asF 1 = 1 for the S-S case. The effect of the aspect ratio L/ h on the frequency response of the nanobeam under different elastic modulus ratios without and with surface elasticity effects is shown in Fig. 9. It is observed that the effect of surface stress, though small, is noticeable on the FRCs. It is noted that setting E S U /E S L = 10 enhances the effect of the aspect ratio. These observations agree with the free vibration results shown in Fig. 2 where low aspect ratio amplifies the surface stress contribution due to higher moment contribution for the unevenly distributed material.
The influence of the inhomogeneity index n and the nonlocal parameter of the beam on the frequency response of the beam is explored in Fig. 10. Different elastic modulus ratios are considered, with or without surface elasticity effects. A relatively low aspect ratio of L/ h = 6 is selected, while elasticity ratio is set to E U /E L = 10. As established in the free vibration response, such a choice should help highlight the effect of the inhomogeneity of the beam material. It can be seen from Fig. 10a that raising the inhomogeneity index n does not always produce the same effect. This is expected since raising n does not always make the beam less homogeneous. Looking back at Figs. 2 and 3 in the free vibration study, a similar behavior was observed where raising the value does not produce a monotonous effect. The effect of the nonlocal parameter on the frequency response under different elastic modulus ratios is depicted in Fig. 10b. As expected, regardless of surface elasticity, raising the nonlocal parameter induces a hardening of the FRC. A similar observation was made by Trabelssi et al. [44].
The effect of the nonlinear foundation is considered in Fig. 11. It can be seen that increasing either the linear stiffness or the shear coefficient softens the FRCs. However, as discussed in the free vibration section, given the expression of β 2 and β 3 , the effect of the shear coefficient is more pronounced than the linear stiffness coefficient and is shape function dependent.
The nonlinear stiffness coefficient appears to be the major contributor to the nonlinear response of the FRCs of the nanobeam, as shown in Fig. 11. In fact, for low values ofk N L , the FRCs show a very limited hardening. This hardening can be rapidly increased by raising the valuek N L .

Conclusions
This paper investigates the effects of combining surface elasticity, nonlocality and physical neutral axis location on the dynamic performance of a graded nonlinear Timoshenko deep beam. The principle of virtual displacements was utilized to obtain the equations of motion by considering the surface effects. By employing the method of multiple scales (MMS), a fourth-order nonlinear ordinary differential equation was obtained to yield the nondimensional nonlinear frequencies and the frequency response curves of the beam. Since MMS was derived using linear classical modeshapes, a DQM numerical solution which uses the nonclassical modeshapes was employed to assess the accuracy of the MMS solution. The present study results were validated based on published results. A parametric study was performed to examine the effects of the aspect ratio, surface modulus ratio, the elastic modulus ratio, the nonlocal parameter, the inhomogeneity index, the elastic foundation coefficients on the dynamic behavior of the graded nanobeam. The major findings in this work are outlined as given below • The nonlinear natural frequency reduces rapidly for low aspect ratios for both S-S and C-C cases.
• The nonlocal parameter reduces the natural frequency, while the amplitude is dependent on n and the ratio E U /E L . • For both C-C and S-S cases, the inhomogeneity index slightly decreases and then increases the natural frequency for low values of E U /E L and remains insensitive otherwise. The inhomogeneity index effect is also slightly sensitive to the amplitude and μ 0 . • Deeper beams are found to be more sensitive to the variation of E U /E L especially for the C-C cases. A low aspect ratio also enhances the effect of the surface stiffness ratio E s U /E s L . • Unlike the S-S case, a study, based on the expressions of β 2 and β 3 , revealed that for low aspect ratio the C-C nanobeam modeshape deviates from the classical one. • A significant variation in the natural frequency is observed as the ratio E s U /E s L increases especially for low aspect ratio. • The forced vibration study showed that raising to E s U /E s L = 10 boosts the effect of the aspect ratio. • The forced vibration study showed that raising the inhomogeneity index n does not always produce the same effect on FRC. This behavior was also observed in the free vibration study. • As expected, in the forced vibration study, raising the nonlocal parameter results in a hardening of the FRC. • Given the expression of β 2 and β 3 , the effect of the shear coefficient on the FRCs is more pronounced than the linear stiffness coefficient and is shape function dependent • In the forced vibration study, it was observed that nonlinear stiffness coefficient is responsible for most of the nonlinear response of the FRCs of the nanobeam. • For the Timoshenko beam, the ratio E s U /E s L affects its bending stiffness and does not contribute to its shear stiffness, while the ratio E U /E L alters both the bending and shear stiffness. For the Euler-Bernoulli beam, there are no shear effects accounted in the model. • Surface stress gradient affects the position of the beam neutral axis and possesses a higher lever arm resulting in a larger moment contribution especially in deep beams.
• The expressions of the nonlinear frequencies j ( j = 1, 2) given by Eqs. (46a) and (46b) and that of the frequency response curve given by Eq. (55) are validated using the differential quadrature method and can be potentially employed in a design problem of a nanobeam.

A Expression of the position of the neutral axis appearing in (12):
To determine the neutral axis location, the axial displacement is ignored and higher-order terms are neglected. Thus, Eq. (11) reduces toB In what follows,B is equated to zero. Simplifying (A.1) results in the following expression: Making the following change of variables z = z m − h 0 , Eq. (A.2) can be written as follows: From Eq. (A.3), the position of the neutral axis can be determined n + 1 (B.1)