A high-order FEM formulation for free and forced vibration analysis of a nonlocal nonlinear graded Timoshenko nanobeam based on the weak form quadrature element method

The purpose of this paper is to provide a high-order ﬁnite element method (FEM) formulation of nonlocal nonlinear nonlocal graded Timoshenko based on the weak form quadrature element method (WQEM). This formulation offers the advantages and ﬂexibility of the FEM without its limiting low-order accuracy. The nanobeam theory accounts for the von Kármán geometric nonlinearity in addition to Eringen’s nonlocal constitutive models. For the sake of generality, a nonlinear foundation is included in the formulation. The proposed formulation generates high-order derivative terms that cannot be accounted for using regular ﬁrst- or second-order interpolation functions. Hamilton’s principle is used to derive the variational statement which is discretized using WQEM. The results of a WQEM free vibration study are assessed using data obtained from a similar problem solved by the differential quadrature method (DQM). The study shows that WQEM can offer the same accuracy as DQM with a reduced computational cost. Currently the literature describes a small number of high-order numerical forced vibration problems, the majority of which are limited to DQM. To obtain forced vibration solutions using WQEM, the authors propose two different methods to obtain frequency response curves. The obtained results indicate that the frequency response curves generated by either method closely match their DQM counterparts obtained from the literature, and this is despite the low mesh density used for the WQEM systems. and other small-scale structural elements constitute the building blocks of micro- and nanoelectromechanical systems (MEMS and NEMS), actuators, sensors and atomic force micro-M.

scopes [1][2][3]. The choice of integrating small-scale components is related to exotic mechanical properties and size effects experimentally observed [4][5][6][7][8] at the nanoscale. While these size effects can be accurately captured and studied using molecular dynamics (MD) simulations, the computational cost of MD is generally prohibitive. Hence, higher-order continuum mechanics approaches have been widely adopted as an alternative in the modeling of small-scale structures. Several higher-order continuum theories have been developed, each of which was based on a different perspective of small-scale behavior. However, in general, most of these theories can be classified into three different categories, namely micro-continuum, strain gradient family and nonlocal elasticity theories.
The nonlocal elasticity theory postulates that the stress in a continuum at a given location depends not only on the strain at that location but also on the strains in a finite neighborhood of such point. This dependency on the nonlocal strain is captured by a size effect parameter called the nonlocal parameter. The nonlocal elasticity theory was first proposed by Kroner [9] and then later improved by Eringen and co-workers [10][11][12].
To simplify the implementation of the theory in practical problems, a differential form was developed [13] based on a specific kernel function. Lately, researchers have explored the possibility of combining nonlocal strain effects with strain gradient theory in a single higher-order theory [14][15][16] referred to as nonlocal strain gradient theory. These size-dependent theories were exploited to model nanorods [2], nanobeams [1,17,18] and nanoplates [19][20][21][22]. These simple structures are conventionally modeled based on Euler-Bernoulli beam theory (EBT) and classical plate theory (CPT), respectively. Other models like Timoshenko beam theory (TBT) and first-order shear deformation theory (FSDT) account for shear to model thick beams and plates, respectively, accurately [23,24].
To overcome the limitations of analytical solutions [24], methods like the finite element method (FEM), the differential quadrature method (DQM), the mesh-free method, the Ritz method, the Galerkin method, etc. were employed to solve small-scale problems and have become the most suitable methods for such problems. In general, numerical techniques are used to solve either the equation of motion or the variational statement. Although developing a solution to the former is generally simpler using collocation methods, solving the variational statement offers several advantages. For example, FEM has weaker regularity requirements (i.e., existence of high-order derivatives) and can easily handle complicated geometries and boundary conditions [25]. These advantages have made FEM the most commonly used method in the analysis of small-scale structures [24].
A few selected studies that use FEM in size-dependent beam problems are summarized in this paragraph, and a more comprehensive review can be found in [24]. Demir and Civalek used a linear nonlocal EBT element in two separate studies [1,18] . The EBT element is based on Hermite cubic interpolation with two nodes and two degrees of freedom per node. The effect of an elastic matrix was accounted for in both studies. Eltaher and his colleagues developed [26][27][28][29][30] nonlocal EBT elements for both functionally graded (FG) nanobeams [26][27][28] and homogeneous nanobeams [29,30]. The EBT element is a two-node element with three degrees of freedom per node: axial and transverse displacements in addition to rotation. The axial displacement is based on a Lagrange linear interpolation, while the transverse displacement is based on Hermite cubic interpolation. Eltaher et al. investigated free vibration problems of FG nanobeams on two separate occasions [26,28]. In the second paper, the authors reexamined the location of the nanobeam's neutral axis based on the physical neutral surface of FG beams [28]. In addition, Eltaher et al. [27] studied the buckling and bending response of graded nanobeams. Like their FG counterparts, homogeneous nanobeams have received considerable attention in the literature. The free vibration problem of homogeneous nanobeams was also examined by Eltaher et al. [29], while static bending of homogeneous nanobeams was considered by Alshorbagy et al. [30]. Nguyen et al. [31] proposed a mixed formulation consisting of developing a nonlocal mixed beam element to examine the static bending response of homogeneous nanobeams. This two-node element uses Lagrange interpolation for both deflection and bending. The literature also shows fewer nonlocal TBT elements. Reddy and El-Borgi [32] developed a finite element formulation for both nonlocal homogeneous EBT and TBT beams. Their models accounted for moderate rotations using the von Kármán strain nonlinearity. Hence, a nonlinear factor was added to the model. Similarly, the nonlinear EBT element relies on a mix of Lagrange and Hermite cubic interpolation for its axial and transverse displacements, while the nonlinear TBT element uses Lagrange interpolation for all its dependent variables. Later this work was extended to graded nanobeams [33]. Eltaher et al. [34] investigated the buckling and bending behavior of nonlocal graded Timoshenko nanobeams.
According to the literature, there have been several nonlocal beam element formulations. Each was tailored or designed to treat a specific problem. Though some elements were developed for TBT, EBT elements dominates the literature [24]. Technically, linear shape functions are sufficient to design an element model for nonlocal TBT. However, when nonlinear behaviors such as von Kármán strain nonlinearity and nonlinear elastic foundations are considered, even second-order elements may fall short of addressing all the high-order derivative continuity requirements in the variational statement. This problem was noted by Reddy and El-Borgi [32,33] where the authors chose to neglect high-order derivatives in the variational statement to be able to solve using FEM. Moreover, the literature shows that none of the cited FEM studies have been used to investigate force vibration response and generate frequency response curves (FRC). The lack of forced vibration response prediction in numerical studies can be traced to difficulties in obtaining steady-state responses for a system with a high number of degrees of freedom.
To address the shortcoming related to estimating higher-order derivatives, the problem was generally solved using high-order collocation methods like DQM [35,36] or the quadrature element method (QEM). This method is a high-order method used to solve FEM problems using a single or few high-order elements without the need to explicitly identify shape functions [37]. It relies on DQM matrices and a clever choice of the grid to simplify its implementation and hence may eliminate the need of an assembly subroutine [37]. QEM can be classified into two major families, namely the strong form quadrature element method (SQEM) [38,39] and the weak form quadrature element method (WQEM) [37,[40][41][42][43][44][45][46][47]. SQEM is also referred in the literature as the differential quadrature element method [37,48] or the strong formulation finite element method [49]. This approach is formulated similar to the regular DQM [37][38][39]48,49] with the additional freedom to subdivide the domain into few elements connected by their respective boundary conditions. This allows more flexibility and mitigates the weakness of DQM for discontinuous loading and geometries. On the other hand, WQEM can be formulated in a similar manner as FEM based on the minimum energy principle or the weak form of the integral or the variational statement. It has also been concluded that WQEM converges faster than FEM [50] and it is also more flexible than SQEM since it is essentially a higher-order FEM [37,40,43]. On another note, the applications to two-dimensional thin-plate problems by either DQM or SQEM have been mostly limited to simple domain shapes. Handling more complicated geometries, though theoretically possible [51], may come at the cost of accuracy and efficiency [52]. On the other hand, WQEM, and similar to FEM, can be employed to solve problems with any irregular shapes without any loss of accurately [40,43]. Furthermore, WQEM stiffness matrix is symmetrical, unlike that of SQEM which may have unstable complex eigenvalues.
In spite of the fact that WQEM is useful in estimating higher-order derivatives, few investigators have realized, however, its importance in solving size-dependent continuum mechanics problems [64] and most notably in the case of forced vibration. To fill this gap in the literature, the authors propose to develop a new FEM formulation based on WQEM to model the free and forced vibration response of a nonlocal TBT resting on a nonlinear elastic foundation accounting for moderate rotation through von Kármán strain. The foundation models the interaction between the beam and the medium in which the beam is embedded such as a protein microtubule embedded in a matrix [1] or a carbon nanotube (CNT) embedded in a foundation [3]. To model [65] the forced vibration response, the authors propose two new numerical methods to estimate the frequency response curves of the nanobeam which are validated based on results obtained by the main authors using the differential quadrature method [36]. The closest study to this work is the paper published recently by Jin and Wang [16] who investigated the free vibration response of a linear and classical Timoshenko graded beam using WQEM. As an extension of this paper, the authors added nonlocal and nonlinear effects in addition to forced vibration response.
Following this introduction, the size-dependent equations of motion and the corresponding variational statement for a nonlocal nonlinear graded TBT are established. The following section outlines how the variational statement is discretized using WQEM to obtain the free vibration response of the nanobeam. Section 4 summarizes the WQEM-based forced vibration solution using two different strategies to obtain the frequency response curves. Free and forced vibration results obtained by WQEM are presented and compared with DQM results in Sect. 5. Finally, a summary of this study and concluding remarks are given in Sect. 6.

Hamilton's principle
According to Eringen's nonlocal theory [10,13], the nonlocal stress is given bȳ where σ (x ) is the classical macroscopic Cauchy stress tensor at point x and K (|x − x|, τ ) is the kernel function of the nonlocal modulus, |x − x| being the distance and τ 0 being a material parameter that depends on internal and external characteristic lengths. Unlike classical mechanics, this relation stipulates that the stress at a given point in an elastic continuum depends on strains all over the body. An equivalent differential model, based on the exponential kernel, was proposed [13] as where e 0 is a material constant, ∇ 2 is the Laplacian operator, and a and are the internal and external characteristic lengths, respectively. It is usually assumed that the nonlocal size effect is only significant along the x-axis of the nanobeam which its along its longitudinal direction. In light of this assumption, Eq. (2) is reduced to the following: where ∇ 2 is reduced to ∂ 2 /∂ x 2 , E is the elastic modulus of the beam and G is its shear modulus. It is worth noting at this point that the transformation from the integral to the differential form of the nonlocal model comes with a paradox for beam bending problems with an exponential nonlocal kernel. In fact, Fernandez-Saez et al. [66] and Romano et al. [67] reported that this transformation yields a relationship that must be satisfied between the bending moment and the spatial derivative of the bending moment at the boundaries. The bending moment obtained from the solution of the differential equation should be checked to ensure the obtained solution is also a solution to the integral form of the model. This is easily done for problems with displacement-type boundary conditions, since the bending moment will be the solution of a second-order differential equation, and the constants of integration can be used to satisfy the bending moment boundary conditions. However, it should be noted that the integral form is incapable to model local effects at boundaries, which may result in some discrepancies between the actual and simulated bending moment at the boundary. Knowing that neither the integral form nor the differential form can solve all possible discrepancies at the boundaries, the differential form is selected in this study. These arguments were also used by the last two authors in a previous paper for choosing the nonlocal differential model [68].
A Timoshenko beam resting on a nonlinear foundation shown in Fig. 1 is considered in this study. Within the context of the small displacement and small deformation theory and only accounting for bending in the x-z plane, the components of the displacement field in a TBT model are assumed to be written as Accounting for von Kármán strain, the Green-Lagrange strain components can be expressed as where Hamilton's principle for the current nonlinear Timoshenko nanobeam can be written as where δ K is the variation of the kinetic energy, δU is the variation of strain energy and δW is the variation of the external work. These terms can be expressed as where q is the distributed transverse load and F v = μ fẇ represents the damping force assumed to be proportional to the velocityẇ wherein μ f is the damping coefficient. It is worth noting that there is no damping associated with rotation since the beam is considered elastic and not viscoelastic [69]. Finally, M

Equations of motion
Substituting Eqs. (8a), (8b) and (8c) into Hamilton's principle (7), and then integrating by parts, yields the motion equations of the nanobeam which can be written as For a beam graded in the z direction, the elastic and shear moduli, appearing in Eqs. (3a) and (3b), are assumed to follow the power-law function below [27,46] where the subscripts U and L designate the upper and lower faces of the beam. Here, h is the thickness of the beam and n k is the material gradation index. In light of the above equations, the nanobeam is considered to be nonhomogeneous with an isotropic stress-strain law. Combining (3a) and (3b) with Eq. (9) yields: Next, using the technique developed by Nayfeh and Pai [70], the axial displacement u is eliminated from the equations of motion. To apply this technique, the following assumptions are adopted: (i) The beam is supported at its both boundary points such that u(0) = u(l) = 0 and (ii) the longitudinal acceleration m 0 ∂ 2 u/∂t 2 and the corresponding velocity are assumed to be very small and hence can be neglected. Applying these assumptions yields the following expression of M (0) x x which can be written as: Further details about this simplification can be found in [36]. With further manipulations, it can be shown that the nonlocal stress resultants can be written entirely in terms of displacements Finally, substituting the above equations into (10b) and (10c) yields the following reduced equations of motion: The nonlinear elastic foundation is assumed to be a transversely acting stiffness. Hence, in the case of a forced vibration load, q(x, t) is given by where k L , k N L and k s are, respectively, the linear, nonlinear and shear coefficients of the nonlinear medium in which the beam is embedded and therefore represent the effect of the surrounding medium. It is also worth noting that this model is a generalization of the linear models known as Winkler [71] and Pasternak [72] foundations, although a more complex model could be a viscoelastic foundation [73]. Previous papers published by the main authors [35,36] confirm that the effect of the surrounding material is crucial and the nonlinear stiffness parameter k N L appearing in Eq. (17) plays a dominant role in the response of the nanobeam. Therefore, it was decided to adopt the current nonlinear foundation in this study rather than the linear classical Winkler-type and Pasternak-type foundations. Finally, F(x) and ω, appearing in the above equation, designate, respectively, the forcing function amplitude and frequency. Finally, the amplitude is set to zero for free vibration case. For scaling purposes, the following normalization is utilized: This yields the following nondimensional equations of motion: For the hinged-hinged (HH) case, the following boundary conditions must be satisfied at the ends of the beam, i.e., at both ξ = 0 and ξ = 1:ŵ A clamped-clamped (CC) nanobeam must satisfy the following boundary conditions at ξ = 0 and at ξ = 1:

Variational statement
The aim of this study is to formulate a high-order variational method. To this end, Eq. (15) is substituted into Eqs. (8a), (8b) and (8c). The resulting expressions are then substituted into the expression of Hamilton's principle (7). Finally, integrating (8a) by parts, the variational formulation can be written as a function of displacements as follows: which can be rewritten in scalar product form as follows: Since only one element will be used to model the nanobeam, the external forces at the boundaries of the beam element are basically reaction forces. Hence, the work of external forces at the boundaries of the nanobeam is zero and does not need to be added to the variational statement [32].
Examining the above variational statement reveals that w is raised to the third derivative in several terms. One alternative used by Reddy et al. [33] is to neglect these terms and adopt a quadratic finite element model that does not account for all mechanical aspects of the system. A better alternative is to raise the order of the finite element model and one viable approach is the p-version of the finite element method. However, a simpler alternative adopted in this paper is to deploy WQEM to discretize the system, which in addition brings high-order accuracy.
Finally, utilizing the normalized variables in (18), the normalized variational statement can then be expressed as follows: The above nondimensional variational statement is the one subsequently solved in the free and forced vibration studies (Sects. 4 and 5). Furthermore, it is obvious that the highest derivativeŵ appearing in Eqs. (24) and (25) cannot be accounted for based on a regular FEM formulation even with second-order elements. The formulation of classical FEM requires the definition of shape functions whose order defines the element's order. Modifying the element order to accommodate higher derivatives (or to increase elements precision) requires the development of a brand new formulation. This locks the order FEM element at the formulation stage. WQEM, however, does not require an explicit computation of shape functions or their derivatives [37]. This allows the use of adaptive order of precision and hence can avoid any unnecessary approximations.

Free vibration WQEM formulation
To simplify the computation of the discretized system and later write the variational statement in matrix form, the following integral is first approximated Here, f (ξ ) and g(ξ ) are arbitrary functions that can be interpolated using a Lagrange polynomial basis of order n and f (ξ ) (m) and g(ξ ) (k) are, respectively, the mth-order and kth-order derivatives of f (ξ ) and g(ξ ) with respect to ξ . To this end, an n-node mesh has to be selected and the integral can be evaluated as where where [M m ] is the mth-order differentiation DQM matrix whose expression is given in "Appendix" 7. M 0 designates the identity matrix which has the same order as the DQM matrices.
In addition, δY 1 and δY 2 are vectors of virtual displacements such as [ Finally, the general displacement vector Y is defined as Y = Y 1 Y 2 , while the velocity and acceleration vectors are denoted byẎ andŸ , respectively. A WQEM discretization is utilized to obtain the free vibration solution of the nanobeam. The mesh coordinates ξ i (i = 1, . . . , n) for n nodes are chosen based on the Gauss-Lobatto-Legendre (GLL) quadrature grid which yields an integration accuracy up to a polynomial of degree (2n -3) [37,45]. For a general linear TBT, this should result in a fully integrated stiffness matrix and reduced integrated mass matrix [37,45]. Applying the spatial discretization in (27) to the variational statement (25) yields the following: and ω ξ is an integral quadrature weight coefficient vector compatible with the GLL grid. M 2 M 1 or M 2 1 is an element by element product (i.e., the Hadamard product) of either two vectors or two matrices. The matrix form in Eqs. (28a) and (28b) produces a system of 2n coupled differential equations, n for each degree of freedom. Hinged-hinged (HH) and clamped-clamped (CC) beams are the only boundary conditions considered herein. Since all boundary conditions present in this study are homogeneous, only essential boundary conditions need to be explicitly stated for variational methods. These boundary conditions are given by where k b is either 1 or n. The system of Eqs. (28a) and (28b) is then reduced to its basic degrees of freedom using the procedure outlined in [35,64,65,74] Here, the superscript {R} denotes the reduced formulation of the system, M Sys is the mass matrix and K Sys (Y ) is the nonlinear stiffness matrix. To obtain the eigenvalues, Y is assumed to have the following form Y = Y e iωt . Then, Eq. (30) is rewritten in the following form: To obtain the linear natural eigen-system, (31). The ith nonlinear natural frequency is obtained through an iterative process which starts with the ith linear eigenvector to evaluate K Since the nonlinear frequency is amplitude dependent, the ith eigenvector must always be scaled relative to the mode shape of w (x, t) in order to keep the mode shape's amplitude constant.

Forced vibration WQEM formulation
The free vibration studies related to the problems similar to the one under investigation largely outnumbered their force vibration counterparts in the literature. In fact, there are limited number studies of forced vibration involving DQM [36,65,74,75], especially in size-dependent mechanics. A number of nonclassical mechanics WQEM studies are extremely rare and are focused on the free vibration response [64]. To the best of the authors knowledge, there have been no forced vibration WQEM studies in the literature similar to the one present with DQM. To fill this gap, two force vibration methods are proposed in this section. Each of the proposed methods aims at finding the periodic steady-state solution of the system for different excitation frequencies.
For this aim, the time is discretized using a periodic grid and periodic derivation matrices. However, adding a time discretization is equivalent to adding another dimension to the problem with consequential important computational cost. Hence, knowing that the forcing term should only excite a few modes, it is necessary to reduce the spatial degrees of freedom. In this section, each proposed method introduces a different approach to perform this task.
Adding the forcing and the damping terms to (28a) and (28b) yields the following WQEM formulation of the variational statement for the forced vibration case: in whichμ f is the normalized damping coefficient andF denotes the normalized discretized force distribution.

WQEM formulation using a mode shape interpolation basis
The mode shape-based forced vibration approach follows three main steps [35,74]: 1. switching the interpolation basis from a Lagrange basis to a modal basis to reduce the number of degrees of freedom; 2. discretizing time using a periodic method, such as the spectral method (SM) or the harmonic quadrature method (HQM); 3. solving the discretized system for a different forcing frequency at the vicinity of its first eigen-frequency and plotting the frequency response curve.

Switching the interpolation basis
As explained earlier, WQEM is a high-order FEM that relies on DQM to express the derivatives of the shape functions at the integration points [64]. This is technically, equivalent to making the following assumptions: where L(ξ ) is a vector of the Lagrange basis relative to {ξ } and the dimensions of each term are specified as a subscript in parentheses. Note that this discretization applies for both the displacementsŵ (ξ, τ ) andφ (ξ, τ ), and virtual displacements δŵ (ξ, τ ) and δφ (ξ, τ ). To reduce the size of the problem,ŵ (ξ, τ ) andφ (ξ, τ ), as well as δŵ (ξ, τ ) and δφ (ξ, τ ), are expressed using a reduced number of mode shapes m. Terms related to the reduced coordinates will use a double script font or will be underlined as indicated beloŵ The mode shape approximation basis forŵ(ξ,τ ) (34c) The mode shape approximation basis forφ(ξ,τ ) in which Y 1 (τ ) and Y 2 (τ ) denote, respectively, the reduced generalized coordinate vectors forŵ (ξ, τ ) and φ (ξ, τ ). The virtual displacements δŵ (ξ, τ ) and δφ (ξ, τ ) are also interpolated in a similar manner. and are a collection of m columns representing linear eigenvectors relative toŵ andφ, respectively, such that { Y w,i } and { Y φ,i } are the ith linear eigenvectors with respect toŵ andφ, respectively. The basic concept here is to interpolate the dependent variables and the virtual displacements using the limited number of dominant linear mode shapes. Hence, the reduced generalized coordinates in Y 1 (τ ) and Y 2 (τ ) are simply the amplitudes of each dominant linear mode. Substituting into (25) gives a system of equations similar to Eqs. (32a) and (32b) as Switching to a mode shape interpolation basis reduces the problem size from 2n equations (n is number of nodes in the mesh) to 2m. This step is driven by the fact that a limited number of first mode shapes generally dominate the dynamic behavior of the nanobeam. In addition, the forcing term is focused to mostly excite the first w modes. Thus,F orF is selected such that whereF 1 denotes the amplitude of the forcing term and μ F represents a scaling factor.Finally, let us define Y = Y 1 Y 2 , so that Eqs. (35a) and (35b) can be expressed in the following form:

Time discretization
Introducing a new timescaleτ such thatτ = τ T = τ 2πω , Eq. (37) becomes The choice of this timescale eliminates the need to update the forcing term for various values ofω. A periodic steady-state solution must be reached to compute the frequency response of the system. This can be expressed with the following periodic initial conditions: A numerical solution for the forced vibration problem requires an adequate discretization of the time dimension. Both spectral method (SM) and harmonic differential quadrature method (HQM) can be utilized for this aim since both methods implicitly implement periodic initial conditions with additional accuracy of high-order methods. Hence, the corresponding compatible mesh is adopted where n τ designates the time increment number. Here, the SM and HQM differentiation matrices are, respectively, provided in Appendices B and C and the discretized time space coordinate matrix is defined as Here, the columns and lines of [Q] correspond to the discretized space and time, respectively. This approach yields the following discretized equation of motion: is the kth-order time derivative matrix and [A] is a (1 × n τ ) line matrix such as [A] i = cos(2πτ i ). Finally, solving (43) for different values ofω in the neighborhood of the first linear frequency is required to obtain the frequency response curve.

Frequency response curve
A frequency response curve can be generated in the neighborhood of each resonance. This applies to both system's dependent variables, namelyŵ andφ. The amplitudes ofŵ's mode shapes as a function of time are obtained through the following transformation: ⎡ where [Q w ] is defined in (42). Finally, [L(τ )], the time discretization basis is given by [76][77][78] for HQM (45) The frequency response curve was established such that the forcing term excited only the fundamental natural frequency. This corresponds to finding max(w 1 (τ )) of the QEM system (32) at several values ofω in the neighborhood of the first natural frequency.

WQEM formulation using Galerkin technique
It is also possible to write (32a) and (32b) as Based on this form, the forced vibration problem can be treated similarly to the DQM case in [74,75,79,80]. Despite the similarity, there are few differences to be noted 1. In this case, a Galerkin projection is applied to discretize the variational statement. 2. Two nested integrals have to be calculated: the WQEM integral and the Galerkin integral.
3. Thanks to the careful mesh choice, it is possible to use the same high-order integration scheme for the WQEM and Galerkin integrals. 4. This method is computationally more intensive than the method presented in the previous section. 5. The main difference between this approach and the one presented in the previous section is the procedure of reducing the number of the degrees of freedom.
Going back to (46), all matrices present in this equation along with [F Sys ] are established using WQEM, i.e., the high-order variational statement. Now that the system looks like a time-dependent differential equation, the Galerkin technique is adopted to limit the size of the WQEM system to 2m equations. As stated in the previous section, only a limited number of mode shapes dominate the nanobeam's vibrational response. Hence, the 2m first linear nonlocal mode shapes are selected as the basis in applying the Galerkin technique [74,75,79,80]. Consequently, the following change of variables has to be made: . This system looks similar to (37), although the resulting matrices are not the same. Nevertheless, the solution procedure is exactly the same from this point onward.

Numerical results and discussion
Various aspects of the graded nanobeam vibration are presented herein including nonlocal linear and nonlinear frequencies in addition to force vibration frequency response curves. A schematic of the beam is presented in Fig. 1. It is assumed that the beam has a square cross section such that b = h = 1 10 L. It is further assumed that its material distribution is graded in the z-direction according to Eq. (11) where − h 2 ≤ z ≤ h 2 and P L and P U designate aluminum and silicon properties, respectively. In this study, the utilized material properties are reported in Table 1 and the considered boundary conditions include hinged-hinged (HH) and clamped-clamped (CC).

Performance of WQEM
A mesh convergence study is the first step. To give a better assessment of the convergence performance of WQEM, it was decided to conduct a mesh convergence study on all cases treated later in Table 5. Generally, the CC case is the slowest to converge. A selection of these difficult CC cases are presented in Table 2 along with their HH analogs. Table 2 shows that the errors for the CC cases fall below 0.5% for as few as 7 nodes. A choice of 7 nodes is totally acceptable, although an 11-node grid is selected for better accuracy. In fact, the error for the HH cases drops below 0.05% using just 7 nodes. This rapid convergence is one of the major advantages for using a high-order variational statement. A similar system would have required 15 nodes to converge using DQM [36]. It is important to note that these results are obtained despite the fact that the nanobeam is highly nonlinear. In this study, a choice of 11 nodes is adopted. The DQM results in this paper were provided from a study by Trabelssi et al. [36] which used a 15-node grid.

Free vibration response
The aim of this section is to replicate the results obtained by Trabelssi et al. [36] for a similar problem using WQEM. The results obtained herein can be divided into two tables. First, the effect of different parameters including the material inhomogeneity index n k , the amplitude of the free vibration A and the stiffness parameters of the elastic foundation, on the free vibration of the nanobeam is investigated in Table 3 for HH and CC nanobeams. The results reported in Table 3 were generated based on the WQEM discretization. For the sake of comparison, Table 3 contains DQM data obtained by Tarbelssi et al. [36], which helps to assess the performance of the proposed approach. The present data utilize a range of values of the amplitude and the nonlocal parameter A andμ 0 , while the inhomogeneity index n k varies from 1 to 4. Both HH and CC configurations are included in this study, while L/ h varies between 10 to 100. The foundation stiffness configuration is described with the following parametersk s = 5,k L = 50 andk N L = 50. Table 3 shows that the frequencies obtained using WQEM match their DQM counterparts up to the fourth digit regardless of the configuration of the   Table 2 WQEM convergence performance ( L h = 10,μ 2 0 = 5,n k = 0.5,A = 1): slowest converging cases selected from the study performed in Table 5 WQEM frequencies WQEM Error to the finest grid (n x = 15) k LkN LkS n x 5  nanobeam, the vibration amplitude or the boundary conditions. Knowing that WQEM results were obtained with a significantly lower mesh density, this truly assesses the accuracy of the WQEM data.
To assess the sensitivity of the proposed formulation to shear locking, the same configuration used to generate the data in Table 3 is used to recompute the nonlocal nonlinear frequencies for thick nanobeams where the aspect ratio L h is kept below 10. To the best of the authors' knowledge, there has been no study that confirmed the presence of shear locking in DQM. In light of this, DQM data were also generated to assess the accuracy of the WQEM results. According to Table 4, the data generated using WQEM are in good agreement with DQM data, indicating the absence of shear locking. This is in agreement with Jin and Wang [82] who also found no shear locking in their linear classical WQEM TBT model.
The effect of the nonlinear foundation is investigated in Table 5. Using the same set of boundary conditions, the nonlinear nonlocal frequencies were computed using WQEM along with DQM data [36]. The values of the different parameters are chosen to underline the effect of the nonlinear elastic foundation which is the highest nonlinear element in the system. It is also the only controllable nonlinearity in the system. The foundation's linear and nonlinear coefficientsk L andk N L are set to vary between 0 and 100 for the former and between 10 and 100 for the latter. The shear coefficientk s varies between 0 and 50, while the inhomogeneity index n k and the amplitude A are set to a fixed value of 0.5 and 1, respectively. Bothk S andk N L are known to affect considerably the behavior of the nanobeam [36]. Despite the wide range of the selected variables, Table 5 shows that the WQEM and DQM results still show great consistency regardless of the selected configuration. In fact, the WQEM low density mesh still performs as good as DQM's higher-density mesh. Technically, variational methods converge faster than their collocation counterparts indicating that these results were expected. The use of the GLL grid helped improve the integration accuracy although the mass matrix is still not fully integrated. It is still possible to improve the accuracy of WQEM with a higher-order integration technique, although this may increase the complexity of the implementation [37].

Forced vibration response
The forced vibration of the nanobeam is examined in this section. This investigation covers both WQEM force vibration approaches detailed in Sects. 4.1 and 4.2 in addition to the previously validated DQM [36]. For each method, the FRC is plotted individually for various values ofk s and for n k = 2,μ 2 0 = 2,k N L = 50,k L = 10,F 1 = 0.75 and L h = 10. This configuration is chosen due to its highly nonlinear behavior and its high dependency onk s . All WQEM FRCs were generated using only 9 nodes, while the DQM FRCs were generated based on a 15-node grid. Figure 2 shows all of the FRC plots for each method as well as a plot of all methods together. These plots are generated for the HH nanobeam. These results show that, despite the lower mesh density, both WQEM methods were able to achieve the same level of convergence as the 15-node DQM results. In fact, Fig. 2 shows that the FRCs totally overlap. The same FRCs are plotted for the CC case in Fig. 3. For this configuration,F 1 had to be raised toF 1 = 1.5 in order to obtain a comparable deformation to the HH case. The mesh density is left unchanged and the FRCs are plotted in a similar manner. Figure 3 shows similar convergence of all methods since again all of the FRCs overlap. Based on the reported results, it is clear that WQEM offers a significant computational advantage over DQM. The method requires fewer nodes than DQM, and WQEM solves this FEM problem using a single high-order element without the need to explicitly identify the shape functions. This simplifies its implementation and eliminates the need of element  [36] assembly. It is important to note that the first approach requires less computational effort than the second one since only one integration is required, while the second approach requires two integrations.

Conclusion
A general formulation of the WQEM with arbitrary element order is presented for an FG nonlocal nonlinear nanobeam based on Timoshenko beam theory. Eringen's nonlocal elasticity was employed to capture size effects of the nanobeam, and a power-law function was utilized to model the material property distribution in the nanobeam. The use of DQ rule simplifies the implementation of the proposed WQEM elements and easily allows an increase in the element order, while the GLL mesh guarantees a significant integration accuracy. For the sake of generality, the formulation accounts for the nonlinear von Kármán strain as well as the contribution of the nonlinear elastic foundation.
The suitability and computational efficiency of the proposed quadrature elements for the vibration analysis of FG beams are demonstrated. A free vibration study was carried out for several mesh densities as well as a set of nonlinear configurations. The study shows an improved convergence rate compared to DQM. The free vibration results indicate that the proposed quadrature FG Timoshenko nonlinear nanobeam element is  highly accurate and efficient. A study of the nonlinear free vibration of thick nanobeams revealed that the proposed element is shear locking free and can yield accurate solutions with a small number of nodes for both thin and moderately thick nanobeams. This approach offers the precision of high-order methods such as DQM without sacrificing the flexibility of variational methods such as FEM. In addition, due to the presence of von Kármán strain and the nonlinear foundation, high-order derivatives are required to evaluate accurately the system response. Such behavior is hard to capture using conventional low-order methods and high-order collocation methods can be used to overcome this problem. WQEM solves this by offering high-order accuracy in a variational method. The study also aims to establish a standard procedure to solve forced vibration WQEM problems. In general, plotting FRCs using numerical methods is challenging due to limitations related to time integration and transient response. Various types of periodic time discretization are generally used which allows time to be treated similarly to a space dimension. However, it consequently adds an extra dimension to the computational problem. Generally DQM force vibration systems resort to Galerkin techniques to reduce the computational cost. To achieve a similar goal for WQEM problems, the authors proposed two different methods to reduce the number degrees of freedom of the WQEM system. The first approach replaces the Lagrange interpolation shape functions generally used in WQEM with mode shape-based functions. The second method utilizes a similar approach to the one used in DQM. Both methods were validated based on a previous DQM study performed by the authors. Despite the lower WQEM mesh density, the FRCs generated using either of these methods overlapped with the DQM FRCs. Obtaining FRCs generally requires that time is discretized like a space dimension which relatively adds a considerable computational cost. The fact that regardless of the approach, WQEM offers comparable results to DQM using a much smaller mesh and this highlights the accuracy and efficiency of WQEM.