A novel formulation for the weak quadrature element method for solving vibration of strain gradient graded nonlinear nanobeams

A novel formulation of the weak form quadrature element method, referred to as the locally adaptive weak quadrature element method, is proposed to develop elements for nonlinear graded strain gradient Timoshenko and Euler–Bernoulli nanobeams. The equations of motion are obtained based on Hamilton principle while accounting for the position of the physical neutral axis. The proposed elements use Gauss quadrature points to ensure full integration of the variational statement. The proposed formulation develops matrices based on the differential quadrature method which employs Lagrange-based polynomials. These matrices can be modified to accommodate any number of extra derivative degrees of freedom including third-order beams and higher-order strain gradient beams without requiring an entirely new formulation. The performance of the proposed method is evaluated based on the free vibration response of the linear and nonlinear strain gradient Timoshenko and Euler–Bernoulli nanobeams. Both linear and nonlinear frequencies are evaluated for a large number of configurations and boundary conditions. It is shown that the proposed formulation results in good accuracy and an improved convergence speed as compared to the locally adaptive quadrature element method and other weak quadrature element methods available in the literature.


Introduction
Nano and microstructure research has seen an increasing interest recently [1]. At such a small scale, mechanical devices exhibit exotic mechanical properties related to atomic-scale interactions, called length scale mechanical properties [2][3][4][5][6] . These effects can be computationally captured using molecular dynamics (MD) simulations. This approach though accurate is very computationally intensive; thus, higher-order continuum mechanics is usually used as a more practical alternative. The strain gradient theory (SGT) is a higher-order continuum theory that is extensively used in small-scale mechanics. This theory introduces different forms of strain energy gradients that can be included in the expression of strain energy. The choice of the strain gradient form yields a different version of SGT [1,[7][8][9][10][11][12] commonly named as SGT family, a detailed list of which can be found in [1,13]. Most SGTs generate additional high-order derivative degrees of freedom (DOFs) at the boundaries.
To model key nanostructures like carbon nanotubes and graphene sheets [14][15][16][17][18][19][20][21], researchers couple SGT constitutive equations with common beam and plate theories such as Euler-Bernoulli beam theories (EBT), Timoshenko beam (TBT), Reddy beam (RBT), and other beam and plate theories. With the addition of highorder DOFs, deriving shape functions for a finite element method (FEM) has generated a lot of interest in the literature [1,[21][22][23][24][25][26][27][28][29][30][31][32]. In fact, either utilizing a higher-order SGT or raising the order of the element significantly complicates the process of computing the shape functions [22,25,31]. Hence, developing high-order elements for nano-mechanics applications is challenging and difficult. However, this type of element provides several advantages including faster convergence and reducing potential numerical problems. To derive high-order solutions, differential quadrature method (DQM) and locally adaptive differential quadrature method (LaDQM) [33,34] constitute attractive alternatives to FEM. However, these methods lack the geometric flexibility of FEM and suffer from few limitations related to concentrated loads, implementation of boundary conditions (BCs), in addition to spurious eigenvalues [35]. The introduction of the strong form quadrature-element method (SQEM) [36] helps alleviate a few of the conventional DQM deficiencies. However, unlike FEM, both DQM and SQEM do not solve the weak form formulation. Lately, higher-order FEMs [37] and weak quadrature element method (WQEM) [38] were proposed as weak form high-order alternatives to FEM. These methods still offer the same advantages as the regular FEM while providing a high-order convergence rate. WQEM is a weak-form method that uses DQM differentiation matrices and quadrature integral to compute the mass and stiffness matrices without the need to compute explicitly the shape functions [39].
Unlike DQM and LaDQM, SQEM and WQEM explicitly require the presence of derivative degrees of freedom at the boundaries. This limits the number of methods used to derive the DQM differentiation matrices for the latter methods. In fact, to the best of the authors' knowledge, there is no available method to implement WQEM formulation for more than two extra derivative DOFs at the boundaries [36,[39][40][41]. Besides, computing these matrices is both challenging and complicated. In fact, Hermite-based DQM matrices used for the second strain gradient (SSG) EBT are totally different from the ones used for SSG TBT which complicates its implementation. Moreover, few SSG beam and plate models may require more than two extra derivative DOFs at the boundaries which is still not supported by the current DQM matrix implementations. This makes implementing WQEM difficult for SGT EBT, TBT, and other third-order beams [42]. To cover this gap in the literature, a new WQEM formulation is proposed in this paper referred to as the locally adaptive weak quadrature element method (LaWQEM). The proposed formulation uses a transfer matrix-based technique to build DQM matrices with derivative degrees of freedom at the boundaries. Hermite polynomial shape functions are also developed using the same technique. This formulation is employed to develop high-order elements for nonlinear SSG TBT and SSG EBT. The proposed elements are designed to ensure full integration capability and account for the nonlinear von Kármán strain.
This manuscript is organized as follows. Following the Introduction, variational statements for a nonlinear SSG TBT and a nonlinear SSG EBT are developed in Sect. 2. Then, the implementation of standard WQEM is presented for both SSG beam models in Sect. 3. The proposed LaWQEM, as well as the proposed LaWQEM SSG TBT and SSG EBT element formulations, are detailed in Sect. 4. Section 5 reports the performance of the proposed formulation based on data from the literature and results from well-established methods. Finally, concluding remarks are provided in Sect. 6.

Weak integral for SSG EBT and TBT beam theory
A graded nanobeam of rectangular cross section with a width b, a height h and a length l, shown in Fig. 1, is considered in this study. For a beam graded in the depth direction, the elastic and shear moduli, E and G, and the density ρ are assumed to vary following a power-law function of the transverse coordinate as follows: where U and L represent, respectively, the upper and lower layers of the functionally graded material, and z G is measured from the mid-plane of the beam. Here, n is the material gradation index, and h is the height of the beam. For a homogeneous beam, the neutral axis is located at the mid-plane of the beam. However, due to the uneven distribution of stiffness in the FGM, the position of the neutral axis is located based on the following equation [43,44]: where (x G , z G ) and (x, z) are the coordinates at the mid-plane and the neutral axis, respectively. Such a choice should remove stretching and bending couplings caused by the FG material variation and reduce the number of nonlinear terms in the equation of motion. Accounting for the von Kármán strain, the Green-Lagrange strain tensor for the nanobeam can be written as where u and w are, respectively, the axial and the transverse displacement at the beam's neutral axis and φ is the angle of rotation of the normal to the neutral axis of the beam. SGT is selected to account for the length scale effect in the nanobeam. According to Jafari and Azzati [31], the different SGTs can generally fall into two different categories, namely the first strain gradient (FSG) and the second strain gradient (SSG). The deformation energy for SSG is given by [31] U ssc= where in which l s is the first non-classical constant and ε i j is the strain tensor. The notation of the derivative convention used in this paper is given by The weak integral form for the nonlinear nanobeam can be derived using Hamilton's principle, where δ K and δU denote, respectively, the variation of the kinetic and strain energy and which can be written as For an EBT where m 0 = A ρd A, m 2 = A z 2 ρd A, and the stress resultants appearing in the above equations are defined as The axial stress can be removed following the same process used in [44,45] based on the following assumptions [45,46]: 1. m 0ü is very small and can be neglected. 2. The beam is supported at x = 0 and x = l such that Then, the following nonlinear term appearing in the expression of δU can be simplified: The process of determining the expression of M x x as a function of the displacement is derived in [47], and the stress resultants (10) can be expressed in terms of displacements as follows: The following non-dimensional parameters are introduced: WQEM requires the expression of the variational statement as a function of the displacement. Thus, utilizing the normalized variables in (16), applying Hamilton's principle (8) whose components are given by (9.1), then writing the stress resultants as a function of displacements (14.2) yields the following normalized variational statement: For an EBT : Integrating by parts, the variational statement yields the boundary conditions of the problem [47]. For a hingedhinged beam (H H) and a clamped-clamped beam (CC), the classical boundary conditions are, respectively, given by These conditions should be coupled with the following non-classical boundary conditions:

Proposed LaWQEM formulation
This Section aims to describe the newly proposed WQEM formulation for SSG beams, which can be employed to develop WQEM elements for any other theory requiring several degrees of freedom at the boundary such as non-local strain gradient theory [13,47,48] or other versions of SGT [13,31]. First, generalities about the WQEM are described including the required degrees of freedom for each beam theory, the matrix form of the integration quadrature, and the derivation of the matrix form of the variational statement. Then LaWQEM transfer matrix-based formulation is presented to facilitate the mapping between the LaWQEM degrees of freedom and LaDQM degrees of freedom. Such mapping ensures the compatibility of the required degrees of freedom with the order of the polynomial shape functions. Furthermore, this mapping allows the use of regular WQEM tools with Hermite-type SSG beams' degrees of freedom and enables the development of various types of Hermite polynomials using linear combinations of Lagrange polynomials.

Standard weak form Quadrature Element Method
The goal of this Section is to discretize the variational integral statements (17) and (18). To this end, the following mesh distribution with n grid points is generally used with DQM and WQEM: Let Y 1 and Y 2 be the discretized displacement vectors relative to w(x, t) and φ(x, t) and Y be the general displacement vector such as Y = Y 1 Y 2 for the TBT problem. For the EBT case, Y = Y 1 is the nodal displacement vector for w which can be written as Derivative degrees of freedom for the TBT problem, Derivative degrees of freedom Internal degrees of freedom w n w n Derivative degrees of freedom for the EBT problem. (23.2) and the velocity and acceleration vectors are, respectively, designated byẎ andŸ . It is worth noting that, similar to FEM, SSG EBT WQEM requires continuity to the second derivatives across elements boundaries. This is due to the fact that the weak form of the SSG EBT has third-order derivatives. Such requirement cannot be satisfied with Lagrange-based WQEM elements. The proposed LaWQEM formulation can generate derivative DOFs at the boundaries to any prescribed derivative order. It is, therefore, capable of ensuring continuity for the SSG EBT system. LaDQM, however, is not designed to handle such a task despite possessing the same number of DOFs. To simplify the implementation of WQEM, the discretized system can be written in a matrix form. The discretization of the normalized variational statement (17) yields the WQEM matrix form for an SSG TBT which can be written as On the other hand, the WQEM matrix form of the normalized variational statement (18) of an SSG EBT can be expressed as follows: where ω ξ is the integral quadrature weight vector, [M k ] is the DQM k th derivative matrix whose expression can be found in [39,42,49] It is important to note that the following relation is employed to obtain the WQEM matrix form of both normalized variational statements: where in this particular context ξ i , (i = 1, ..., n), are the quadrature points coordinates. Note that ω ξ , the integral quadrature weight coefficient vector, should be compatible with the selected grid. For example, a Gauss-Lobatto-Legendre (GLL) grid can be used in conjunction with compatible integral quadrature weight coefficients. However, in the current cases, the presence of derivative degrees of freedom raises the order of the interpolation polynomial. Thus, the accurate evaluation of the discretized variational statement requires a far higher accuracy than offered by GLL. In fact, using GLL, the accuracy of the integration is limited to an order of 2n − 3, which is enough to achieve reduced integration accuracy for a problem without derivative degrees of freedom. However, unlike non-local TBT [39,46], neither the mass nor the stiffness matrices are fully integrated. To improve the accuracy of the solution, another transfer matrix is required [39].
To obtain a fully integrated system, a higher number of integration points and a high-order integration scheme must be used such as Gauss quadrature integration. Let, for example, f (x) be an arbitrary function discretized in the SSG EBT Hermite grid. Using Hermite interpolation, the k th derivative of f (x) can be interpolated as follows: where H i (x) are the Hermite interpolation functions based on grid coordinates ξ i and f ξ is f (x) DOFs vector. Applying the Gauss quadrature rule [38,39] gives  (28) yields Thus, all DQM derivative matrices [M k ] and quadrature weights ω ξ in the discretized variational statement should be now replaced by [ς M H,k ] and ω ς , respectively. Further details about this formulation can be found in [46].

Derivation of the proposed LaWQEM formulation
The proposed LaWQEM transfer matrix-based formulation is presented in this Section to facilitate the mapping between the LaWQEM degrees of freedom and LaDQM degrees of freedom, thereby ensuring the compatibility of the required degrees of freedom with the order of the polynomial shape functions. As a result, it is possible to employ regular WQEM tools with Hermite-type SSG beams' degrees of freedom and build various types of Hermite polynomials using linear combinations of Lagrange polynomials.

Development of LaWQEM DOFs
The proposed LaWQEM formulation relies on an LaDQM compatible mesh. The coordinates of the LaDQM mesh adopted in this paper are given by [39] where δ is a scaling factor, n is the number of internal nodes, and n R and n L denote, respectively, the number of external nodes at the right and the left sides. Based on this mesh, a vector of n + n R + n L degrees of freedom is required for each dependent variable. Unlike the Hermite-based DQM matrices, the basic LaWQEM matrices do not provide a direct interpolation polynomial basis. The only known interpolation polynomials are the Lagrange polynomial associated with the LaDQM grid. Thus, applying a similar approach as in Sect. (3.1) is not straightforward.
Assuming that the normalized displacement w is discretized on an LaDQM grid, then Y 1 (i.e., the vector of nodal displacements for w) can be written as For the EBT case, Y = Y 1 is the displacement vector such as n R = n L = 2, The DOFs of the LaDQM displacement vector in (32) are related to the DOFs of the standard WQEM displacement vector in (23.2) using the following relationship: where M i is the DQM i th derivative matrix relative to the LaDQM mesh and [T ] is the LaWQEM mapping matrix. Further details about the derivation of [T ] can be found in Appendix A.

LaWQEM transfer matrix-based formulation
Equation (33) establishes a relation between the LaWQEM and the LaDQM DOFs. However, it does not allow to generate a direct interpolation polynomial basis associated with LaWQEM matrices. The closest useful interpolation polynomials are the Lagrange polynomial basis relative to the original LaDQM grid. This means that using the relation established in Eqs. (28) and (29) is difficult since building [i T H ] is not possible without the explicit expression of the Hermite basis. To get around this problem, consider the integration of a function f (x), discretized in the LaDQM grid, using the Gauss quadrature rule [38,39], where ς M L,k is the LaWQEM differentiation matrix relative to ς i and ξ i . Comparing (28) and (35)

LaWQEM polynomial-based new formulation
It was stated previously that the proposed formulation does not come with an explicit polynomial basis. However, it is possible to obtain a suitable basis using the same transfer matrix method. Consider a function expressed in the current LaDQM mesh using the following scalar product: Transitioning to LaWQEM coordinates can be obtained using Eq. (37), for the current SSG EBT formulation. Technically, the product of the LaDQM Lagrange base with [T ] −1 appearing in (37) produces a Hermite basis with two derivative DOFs at the boundaries. Hence, [HL(x)] is considered to be the analytical basis associated with LaWQEM coordinates, and it is, therefore, possible to create any Hermite basis using a combination of Lagrange polynomial and a DQ rule. Figure 2 shows a plot of these polynomials for a mesh of 7 internal nodes and 11 DOFs. It is clear from the plots that the [HL(x)] verifies a number of required properties at the mesh points. The mesh along with the nodal location is plotted using short vertical lines along the horizontal axis of all the plots illustrated in Fig. 2. It is also worth mentioning that it is possible to perform the same procedure for Hermite basis with single derivative DOFs at the boundaries. Using the new [HL(x)] basis, it is possible to calculate b a f (k) (x)dx similar to the Hermite interpolation case in Eq. (28), . . . Details about applying the boundary conditions for an LaWQEM system can be found in [46].

Validation of the proposed LaWQEM formulation
The aim of this Section is to validate the proposed LaWQEM approach using data from the literature in addition to analytical results and well-established numerical methods. As detailed in Section 3.2, LaWQEM relies on modified LaDQM matrices obtained using the transfer matrix approach. The solutions of free vibration of SSG TBT and SSG EBT were developed using LaWQEM. SSG usually results in a large number of possible boundary conditions. In this study, only a symmetric set of boundary conditions is selected, a list of which is presented in Table 1. Each boundary condition will be referred to by a reference number as stated in Table 1. However, to assess the validity of the proposed formulation, an analytical solution to a simple vibration problem is treated beforehand. High-order methods rely generally on raising the order of the elements to achieve numerical convergence. The number of element is, thus, kept to a minimum [39,41,46]. Elements are introduced at discontinuities such as an abrupt loading, geometric or material property variation. Therefore, one single element discretization is used though this study unless indicated otherwise.

Analytical free vibration of classical EBT validation
EBT beam vibration is one of the most common classical mechanics problems that requires imposing two BCs at each boundary. This problem is certainly one major contributor that led to the development of several methods based on DQM. Since an analytical solution to this problem is available, it can be used to assess the validity of the proposed LaWQEM formulation. This problem was used before in the literature for a similar purpose [35]. The normalized version of the equation of motion of the free vibration of an EBT beam is given by [35] w (4) where w(x) is the transverse displacement. The weak integral form is given by The convergence error with respect to the analytical solution of the first five nonzero linear natural frequencies of a classical EBT is reported in Fig. 3. The data are plotted for LaWQEM EBT using six different mesh densities. Similar data are also reported using LaDQM to help assess and compare the convergence speed. The results indicate that LaWQEM rapidly converges to the analytical solution faster than LaDQM which proves the validity of the proposed approach to deriving DQM matrices. It is well known that variational methods like LaWQEM generally converge faster than collocation methods like LaDQM. A similar observation was also made by Trabelssi et al. [46] when comparing DQM and QEM convergence speed. The modeshapes generated using LaWQEM are also plotted in Fig. 4. It can be seen that the generated modeshapes are similar to the conventional EBT modeshapes

Hermite-based WQEM SSG EBT validation
The free vibration of a simply supported linear SSG EBT study by Ishaquddin and Gopalakrishnan [41] is used as the second validation benchmark for the proposed LaWQEM formulation. Both analytical and numerical solutions of this problem are available in the literature. The numerical solution was obtained using Hermitebased WQEM with the GLL grid for numerical integration. Assuming the rotary inertia to be negligible, the following parameters are selected in consistent units for this validation: L = 1, Young's modulus E = 3×10 6 , Poisson's ratio ν = 0.3, density ρ = 1. The simply supported boundary conditions for the linear SSG EBT used in this study are given by at both boundaries of the beam. The results of this study are presented in Table 2. It can be seen that LaWQEM gives an accuracy similar to the Hermite-based WQEM and to the analytical results. However, the use of full integration in the LaWQEM formulation provided, in general, a faster convergence rate compared to Hermite-based WQEM.

Results and discussion
The goal of this Section is to study the accuracy of the proposed LaWQEM formulation for linear and nonlinear vibration problems based on SSG nanobeams. First, the proposed method is applied to a linear vibration problem. Then, the general procedure for obtaining the nonlinear fundamental vibration mode is detailed. Finally, a study of the accuracy of the proposed method in comparison with the regular LaDQM approach is performed.

Linear natural frequency results and mesh convergence study
In this Section, the beam is set to the following configuration l s = 0.2, L h = 10, n = 3 while the remaining material properties are provided in Table 3. The first set of data is shown in Figure 5. It summarizes a convergence study of the first three linear natural frequencies using different mesh densities. To limit the amount of data, only four boundary conditions among the ones listed in Table 1 are considered for each beam theory. It is generally known that variational methods converge faster than their collocation counterparts [39,41,46]. Examining the LaWQEM natural frequencies in Fig. 5 seems to confirm this statement for both SSG TBT and SSG EBT systems. In fact, for the TBT case, the first mode can be obtained accurately within a mesh of only 9 nodes. A mesh of 9 nodes is also required to obtain the same accuracy for the EBT problem. Given that the EBT case requires two more derivative DOFs at the boundaries, it can be concluded that for both systems LaWQEM converges at a comparable mesh density. To better assess the convergence error, the percentage difference between the current mesh and the highest mesh density is presented in Fig. 5. It is clear that the rest of the eigenvalues data converges within 9 to 11 nodes for the TBT problem, while 11 nodes are required to deliver sufficient accuracy for the EBT problem. In general, the accuracy of LaWQEM with 11 nodes is comparable to LaDQM with 15 nodes. Therefore, it can be concluded that the proposed method was successfully employed to calculate the first three linear frequencies of both SSG TBT and SSG EBT problems.

Nonlinear natural frequency and mesh convergence study
To solve for the eigenvalues, Y is set as Y = Y e iωt . The discretized variational statement, given by (24) and (25) for the TBT and EBT beams, respectively, can be written in a unified manner as follows: The linear natural frequencies can be obtained by replacing K Sys (Y ) with K Sys ({0}) in (43). The nonlinear frequency is, however, generally obtained by an iterative process [45,46,[51][52][53] or using the harmonic balance method for highly nonlinear problems [54].
In this Section, a study of the mesh convergence of the fundamental nonlinear mode is performed. To this end, the following parameters are selected: l s = 0.3, A = 1, L h = 10, n = 3. The same boundary conditions employed in Section 5.1 are applied for each beam theory. The results of this study are provided in Fig. 6. Examination of this Figure reveals few interesting observations. Despite possible rounding errors due to the presence of the transfer matrix given by (35) and (37), the LaWQEM data show rapid convergence compared to LaDQM. In fact, even for a mesh of 9 nodes the results are below 0.01% error for most cases for both SSG TBT and SSG EBT systems. As far as mesh convergence is concerned, the data in Fig. 6 indicate that LaDQM converges at about 13 to 15 nodes. In addition, a mesh of 9 nodes is sufficient to have an error less than 0.001% for both SSG TBT and SSG EBT beams. Therefore, a mesh of 11 nodes is selected for LaWQEM subsequently.

Nonlinear natural frequencies: results and discussion
To further demonstrate the accuracy of the proposed method, an extensive parametric study is performed based on the selected mesh configuration for each beam theory. In this Section, a combination of different material properties such as n andl s is used while the same set BCs kept. The vibration amplitude A ranges from 0 to 1, which basically covers the linear and most nonlinear cases. The aspect ratio L h of the beam ranges between 10 and 100 for the SSG EBT and between 5 to 100 for the SSG TBT. The choice of adding a lower aspect ratio for the SSG TBT is to ensure that the proposed formulation still provides accurate results for low aspect ratios and does not experience shear locking [46]. The data are presented in Tables 4, 5 , 6, 7, 8, 9, 10, 11, each of which covers one single boundary condition configuration. Similar to the previous Section, LaDQM results are included for assessing the accuracy of each method.
The results show similar trends as those mentioned in the mesh convergence Section (Sect. 5.2). The LaWQEM results seem to be unaffected by the presence of a transfer matrix. In fact, using only a mesh of 11 nodes, the results are as accurate as the LaDQM results obtained with a 15-node mesh. This should offset the complexity of developing a discretized variational statement including the integration procedure. Besides, despite the low aspect ratio used in the TBT results, comparing both LaDQM and LaWQEM shows that the  proposed method provides the same results as LaDQM even when the vibration amplitude is set to 1. Since LaDQM does not experience shear locking, these results show that the proposed LaWQEM SSG TBT element does not experience shear locking either [46]. It can be concluded that the proposed LaWQEM formulation offers similar performance for both nonlinear SSG TBT and EBT systems, despite the differences in the DQM matrix formulation.

Conclusions
A novel formulation of the weak form quadrature element method referred to as the locally adaptive weak quadrature element method was proposed to develop elements for nonlinear graded strain gradient Timoshenko and Euler-Bernoulli nanobeams. The proposed formulation develops a Hermite-based polynomial basis using the Lagrange-based differential quadrature method. These matrices can be customized to accommodate any number of extra derivative degrees of freedom encountered in higher-order strain gradient beam problems without requiring an entirely new formulation. The proposed DQM matrix formulation combines Similarly, the proposed strain gradient Euler-Bernoulli beam element formulation accounts for the von Kármán strain and covers the classical and non-classical EBT degrees of freedom. In the latter case, two extra derivative degrees of freedom at each boundary are required. Unlike the Hermite-based DQM, these elements were developed using similar transfer matrices with minimal coding effort while providing support for both full quadrature integration and GLL reduced quadrature integration. The variational statements were derived using Hamilton's principle while accounting for the von Kármán strain and the position of the beam neutral axis. The validation first task was performed based on the analytical results of the free vibration problem of a classical beam. A second validation was performed using numerical data from an SSG EBT Hermite-based WQEM element formulation whose natural frequencies are available in the literature. The results showed a good agreement with the analytical solution of the classical beam and the numerical solution of the SSG EBT Hermite-based WQEM element. The use of Gauss quadrature in the proposed formulation yields faster convergence than the SSG EBT Hermite-based WQEM element which utilizes reduced GLL quadrature.
The mesh convergence study was carried out for both linear and nonlinear free vibration versions of the strain gradient Timoshenko and Euler-Bernoulli beams. LaDQM results were also computed to compare the accuracy and the convergence rate with the proposed method. The proposed LaWQEM results converged consistently faster than the LaDQM results for the same mesh density. These results were consistent for both linear and nonlinear tests.
Finally, several nonlinear natural frequencies were computed for a large number of material properties combinations, vibration amplitudes, and boundary conditions. A low-aspect-ratio case was also included in SSG TBT results. Based on the obtained results, it was concluded that the proposed LaWQEM method can produce accurate solutions of strain gradient Timoshenko and Euler-Bernoulli beams for various boundary conditions and material properties. The formulation may eventually be extended to third-order nanobeams and higher-order strain gradient beams that may require more than two derivative degrees of freedom at each boundary.
Appendix B Discretizing the 1 0 w 2 dξ integral and the expression of iW (Y 1 ) Conventionally 1 0 w 2 dξ is integrated using the Newton-Cotes formula [55], which is given by the following:  which can be equally approximated using the Newton-Cotes formula [55]. Newton-Cotes formula is generally not the most accurate method to evaluate this integral. In fact, applying WQEM/LaWQEM full integration method developed earlier should give more accurate results.