On Nonlocal Computation of Eigenfrequencies of Beams Using Finite Di®erence and Finite Element Methods

In this paper, we show that two numerical methods, speci¯cally the ¯nite di®erence method and the ¯nite element method applied to continuous beam dynamics problems, can be asymptotically investigated by some kind of enriched continuum approach (gradient elasticity or nonlocal elasticity). The analysis is restricted to the vibrations of elastic beams, and more speci¯cally the computation of the natural frequencies for each numerical method. The analogy between the ¯nite numerical approaches and the equivalent enriched continuum is demonstrated, using a continualization procedure, which converts the discrete numerical problem into a continuous one. It is shown that the ¯nite element problem can be transformed into a system of ¯nite di®erence equations. The convergence rate of the ¯nite numerical approaches is quanti¯ed by an

In this paper, we show that two numerical methods, speci¯cally the¯nite di®erence method and the¯nite element method applied to continuous beam dynamics problems, can be asymptotically investigated by some kind of enriched continuum approach (gradient elasticity or nonlocal elasticity). The analysis is restricted to the vibrations of elastic beams, and more speci¯cally the computation of the natural frequencies for each numerical method. The analogy between thē nite numerical approaches and the equivalent enriched continuum is demonstrated, using a continualization procedure, which converts the discrete numerical problem into a continuous one. It is shown that the¯nite element problem can be transformed into a system of¯nite di®erence equations. The convergence rate of the¯nite numerical approaches is quanti¯ed by an

Introduction
In this paper, the vibration behavior of elastic beams will be investigated from a numerical point of view. It is known that numerical approaches such as the¯nite di®erence or the¯nite element methods convert a continuous vibrations problem into a discrete one. In fact, the discretization methods lead to the resolution of an algebraic problem for an initial continuous eigenvalue problem. The possibility to solve automatically the algebraic problem using a computer makes the discretization approach advantageous as compared to the initial continuous one. However, the reduction process inherent to the discretization may incur the loss of some fundamental mathematical properties of the initial continuous system. There appears to be a need to better explore the e±ciency of the numerical schemes with respect to the continuous problem.
A possible way to handle these numerical discrete problems is to de¯ne a kind of an equivalent continuum that is a representative of the discrete problem. The de¯nition of an equivalent continuum from a discrete one may be labeled as a continualization procedure. Continualization procedures are based on various approximations of the discrete operators by some continuous ones via Taylor expansion or Pad e approximants. [1][2][3][4][5][6][7] The so-called enriched continuum equivalent to the discrete one is sometimes called a quasi-continuum. 2 It is generally dependent on the truncated terms in the asymptotic expansion of the di®erence operators. This method was pioneered by Kruskal and Zabusky 1 and initially applied to discrete wave equations, with applications to the Fermi-Pasta-Ulam model, 8 an axial lattice with nonlinear interaction (see also the analysis of Zabusky and Kruskal within the dynamics of solitons 9 ). The reader can refer to Rosenau, 3 Palais 10 and Maugin 11 for an historical perspective on the link between the Fermi-Pasta-Ulam lattice model and the continualized wave propagation equation. Kruskal and Zabusky 1 used a Taylor expansion of the second-order¯nite di®erence operator arising in the discrete lattice up to the fourth-order spatial derivative. Collins 2 proposed to use the inverse of the second-order¯nite di®erence operator, thereby avoiding the use of fourth-order spatial operators. Pad e approximants of the¯nite di®erence operators were introduced by Rosenau 3 and are shown to be e±cient for capturing the wave propagation in the dynamics of axial lattice. [3][4][5][6][7] It is worth noting that some continualization processes have been already applied to¯nite di®erence structural problems by Cyrus and Fulton,12,13 or to¯nite element problems by Walz et al. 14 In this paper, a similar continualization reasoning will be followed for approximating some¯nite numerical schemes (¯rst-order and higher-order¯nite di®erence methods and¯nite element methods) in beam problems by some equivalent enriched continua. We show that the rate of convergence is strongly dependent on the order of the¯nite discrete scheme: higher-order¯nite schemes lead to higher-order enriched constitutive laws with a higher convergence rate.
The comparison of discrete methods with the equivalent continuum is not new. In fact, it was conducted by Livesley, 15 Greenwood, 16 Leckie and Lindberg 17 for beam vibration problems using¯nite di®erence methodology for instance. Livesley 15 or Leckie and Lindberg 17 gave exact solutions of the vibration frequencies for the¯nite di®erence beam vibration problems. The results have been recently generalized by Zhang et al. 18 for general boundary conditions. In fact, already Lagrange 19 and Rayleigh 20 had determined the exact vibration frequencies of a string with a¯nite number of concentrated masses and compared the solution with the continuous system that was asymptotically obtained as a limit (see also Refs. 15 or 21 on this topic). It can be shown that the discrete string problem is equivalent to the¯nite di®erence formulation of the continuous string problem. The same analogy is also valid for elastic beams. As already discussed by Silverman, 22 Hencky's bar chain 23 which is the bending discrete system (or microstructured beam model) -is in fact equivalent to the central¯nite di®erence formulation of a continuous problem, i.e., the Euler-Bernoulli continuous beam problem. Therefore, the¯rst-order central nite di®erence formulation of an Euler-Bernoulli continuous beam problem is strictly equivalent to the Hencky's microstructured chain. The performance of Finite Di®erence Method for solving buckling or vibration eigenvalue problems has already been evaluated in the literature (see early studies in Refs. 12,13,[15][16][17][24][25][26][27][28][29][30][31], but generally without resorting to any nonlocal mechanics perspective (except in recent paperssee Ref. 32). Furthermore, the nonlocal equivalence will be extended here for higher-order¯nite di®erence schemes. A consequence of the nonlocal equivalency principle for the modeling of discrete systems is that the¯nite di®erence system can be e±ciently approached by nonlocal continuum mechanics tools. As it is known in the case of nonlocal mechanics behaviors, this result con¯rms the lower bound solution of such approximate Finite Di®erence Methods, at least for homogeneous structures (with respect to both convergence and rate of convergence arguments). We extend such a result for approximate Finite Element Methods using gradient elasticity constitutive law, which shows the upper bound solution of Finite Element results based on the work-energy formulation. Excellent monographs are available for the computation of eigenfrequencies of structural members using Finite Element Methods (see for instance Refs. [33][34][35][36][37]. Exact solutions of the Finite Element formulation of the dynamics of Euler-Bernoulli beams have been given by Tong et al., 38 Belytschko and Mindle 39 or Xie and Steven 40 using a cubic-based interpolation function for the displacement. These last authors also discussed the possibility to have a lumped mass matrix or a consistent mass matrix with respect to the displacement interpolation¯eld. Belytschko and Mindle 39 and Xie and Steven 40 derived the exact frequency solutions using the consistent mass matrix with cubic-based interpolation function for the displacement. Tong et al. 38 also obtained the asymptotic solution of the frequency parameter with respect to the size of the¯nite element. These numerical results are revisited in this paper using an equivalent gradient elasticity model.

Continuous Problem
The vibration problem of a continuous elastic Euler-Bernoulli beam is investigated herein. The continuous reference problem is¯rst brie°y presented. The elastic bending momentcurvature constitutive law is expressed for the Euler-Bernoulli kinematics by where M is the bending moment, E the Young modulus, I the area moment of inertia, w the de°ection of the beam, and the prime denotes di®erentiation with respect to x. The equation of motion including the inertia forces is given by where the superdot denotes the di®erentiation with respect to time, and is the mass per unit length of the beam material. By combining Eqs. (1) and (2), one obtains the fourth-order partial di®erential equation of motion: EIw ð4Þ þ w :: ¼ 0; By assuming harmonic motion with ! as the angular frequency of vibration, the free vibration problem of the Euler-Bernoulli beam is governed by the linear di®erential equation: which can be equivalently handled in a weak format using the principle of virtual work, i.e., The Rayleigh's quotient R½w can be introduced for the computation of ! 2 : For simply supported boundary conditions, the exact natural kth vibration mode is the trigonometric solution given by: The substitution of Eq. (7) in the Rayleigh's quotient Eq. (6) leads to the wellknown solution of the continuous Euler-Bernoulli beam problem: where the subscript k refers to the kth mode considered, and the in¯nite character refers to the continuous reference problem which possesses an in¯nite number of degrees-of-freedom.

First-Order Finite Di®erence Method
The¯rst-order¯nite di®erence formulation of this vibration problem will now be presented, and exactly solved for the simply supported boundary condition. The constitutive law Eq. (1) is now written in this¯nite di®erence formulation as: a is the uniform grid spacing also de¯ned by a = L/n where n is the number of grid spacing (and n + 1 is the number of uniformly spaced grid points). The equilibrium Eq. (2) is expressed with the second-order¯nite di®erence equation: :: By using again both the equilibrium equation and the constitutive law expressed in the¯nite di®erence formulation, one obtains the discretized equation of motion : By assuming harmonic motion, the free vibration equation of motion is given by An exact solution to this problem can now be investigated. The same methodology can be followed for the vibration equation, as detailed for instance in Ref. 31, from the fourth-order linear¯nite di®erence Eq. (12) restricted to the vibration terms: The characteristic equation is obtained by replacing w i ¼ A i in Eq. (13) which leads to: where A is a constant. Equation (14) is symmetrical with respect to interchanging and 1=,a n d admits four solutions written as (see also Ref. 31): For the simply supported discrete beam system, the natural vibration modes are obtained from the trigonometric shape function w i ¼ W sinðiÞ; thus leading to the natural vibration frequency n ¼ k as obtained in Refs. 17 or 31: We then obtain for the square of the natural frequencies: where ! 2 E ¼ EI ð k L Þ 4 is the natural frequency parameter of the continuous beam. The discrete equations are extended to an equivalent continuum via a continualization method. The following relation between the discrete and the equivalent continuous system w i ¼ wðx ¼ iaÞ holds for a su±ciently smooth de°ection function as: The pseudo-di®erential operators can be introduced as: The pseudo-di®erential operator can be e±ciently approximated by the Pad e's approximant (see for instance Refs. 3-7): By using such a pseudo-di®erential operator, the constitutive law Eq. (9) using Eqs. (19) and (20) may be continualized as: One recognizes an Eringen's type di®erential equation 41 applied at the beam scale, thereby leading to a nonlocal bending momentcurvature constitutive law. By using the same methodology, the equilibrium equations [Eq. (10)] may be also continualized as: In view of Eqs. (21) and (22) and neglecting the terms in l 4 c , one obtains the modi¯ed nonlocal bending wave equation as: Equation (23) can be also directly obtained from the expansion of the discrete operator in Eq. (12) as: The continualized problem governed by Eq. (23) is also equivalent to the following nonlocal model, where the nonlocal operator is applied only to the constitutive law and the equilibrium equations remain local: M À 2l 2 c M 00 ¼ EIw 00 and M 00 ¼Àw Within this point of view, it is worth mentioning that a factor 2 a®ects the length scale calibration in the nonlocal law. The Rayleigh's quotient for the computation of ! 2 of the nonlocal problem can be presented as: For the simply supported beam problem, the substitution of Eq. (7) in the Rayleigh's quotient Eq. (26) leads to the nonlocal solution of the continuous Euler-Bernoulli beam problem: Equation (27) is consistent with the asymptotic expansion of the exact solution given by Eq. (17), and it also coincides with the asymptotic expansion given by Leckie and Lindberg. 17 We note that there is a factor 2 between the equivalent length scale for the buckling problem and the one of the dynamics problem 32 :

Higher-Order Finite Di®erence Method
The same continuous problem will now be handled using a higher-order¯nite difference scheme. An improved¯nite di®erence analysis may be based on the introduction of the second-order central di®erence for the expressions of the¯rst and the second derivatives of the displacement (see Refs. 16,[28][29][30]. The higher-order¯nite di®erence formulation of the constitutive equation [Eq. (1)] now reads as The higher-order¯nite di®erence equilibrium equation [Eq. (2)] is presented as: The consideration of harmonic motion and in view of Eqs. (29) and (30) lead to the eight-order¯nite di®erence equation: It is also possible to get an exact analytical solution to this numerical scheme, by solving the linear¯nite di®erence equation. The characteristic equation is obtained by replacing w i ¼ A i in Eq. (31) which leads, with 2 ¼ ! 2 L 4 EI to: which can also be presented as: It is possible to factorize in the following way: This equation admits eight solutions for . By solving the¯rst term, one gets The¯rst group of solutions are obtained as: The second group of solutions comes from the second-order polynomial equation: which also possesses four solutions: It is also possible to rewrite the set of solutions 3;4 using the trigonometric functions: 3;4 ¼ cos AE j sin with ¼ arccos 4 À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 33 þ n 2 s "# and j 2 ¼À1: ð39Þ For the simply supported discrete system, the natural vibration modes are obtained from the trigonometric shape function w i ¼ W sinðiÞ; thus leading to the natural vibration frequency n ¼ k, which can also be expressed by: which has hitherto not been reported to the best of the authors' knowledge. It is possible to show from Eq. (40) by using an asymptotic expansion that: The higher-order pseudo-di®erential operator in Eq. (29) or in Eq. (30) can also be expanded as: The constitutive law Eq. (29) can then be continualized such as: The equilibrium equations [Eq. (30)] can be continualized as well: M 00 ¼Àw :: À l 4 c w :: ð4Þ with l 4 c ¼ By combining Eq. (43) with Eq. (44), one obtains the nonlocal bending wave equation as Equation (45) could have been obtained directly from the asymptotic expansion: Equation (46) is also equivalent to the nonlocal model: M þ 2l 4 c M ð4Þ ¼ EIw 00 and M 00 ¼Àw As observed for the¯rst-order¯nite di®erence method, a factor 2 also a®ects the length scale calibration in the nonlocal law for the higher-order¯nite di®erence method. The Rayleigh's quotient for the computation of ! 2 of the new nonlocal problem can be presented as: For the simply supported beam problem, the substitution of Eq. (7) in the Rayleigh's quotient Eq. (48) leads to the nonlocal solution of the continuous Euler-Bernoulli beam problem: We also note that there is a factor 2 between the equivalent length scale for the (static) buckling problem and the one of the dynamics problem:

Finite Element Method
The Finite Element Method is now applied to the beam vibrations problem, using Hermitian cubic-based interpolation functions. The Hermitian cubic functions can be used for the interpolation function of the displacement¯eld: where ¼ x=a. By substituting this Hermitian cubic function into the Rayleigh's quotient, one obtains which is now calculated for the cubic-based Hermitian interpolation function: ðw iÀ1 a iÀ1 À w i a i Þþ 13 6 ðw i a iÀ1 À w iÀ1 a i Þ : ð53Þ By taking the stationarity conditions of the Rayleigh's quotient R½w;¼0 for the two-variable¯eld ðw i ; i Þ, we obtain the coupled system of¯nite di®erence equations: which would have been equivalently obtained by using the weak formulation Eq. (5) of the problem with the Hermitian cubic functions. It is possible to de¯ne the following¯nite di®erence operators The¯nite di®erence system associated with the Finite Element Method can then be presented using the¯nite di®erence operators: The¯nite di®erence equation is then obtained in a single format as: which can also be written as 2 ! 4 a 6 63ðEIÞ 2 w i þ We then obtain for the square of the natural frequencies: where ! 2 E ¼ EI ð k L Þ 4 is the natural frequency parameter of the continuous beam. This result Eq. (64) with the factor 1/720 has been already obtained by Tong et al. 38 Now, by using a continualization procedure, an asymptotic expansion of each di®erence operator in Eq. (57) gives: which can be e±ciently approximated by the following sixth-order di®erential equation, when collecting the terms up to the fourth-order in a 4 : The di®erential equation [Eq. (66)] can also be factorized as which means that the cubic-based¯nite element model can be equivalently reduced to the eight-order di®erential equation: EI a 4 720 w ð8Þ þ EIw ð4Þ À ! 2 w ¼ 0: Walz et al. 14 also obtained an eight-order di®erential equation for the continualized bending problem which was investigated by the Hermitian-based Finite Element model, with the correct coe±cient 1/720 but with a di®erent sign. For the¯nite element model considered herein, the associated Rayleigh's quotient can then be expressed by: leading to the gradient elasticity solution Eq. (68), associated with the gradient elasticity constitutive law: Considering again a simply supported beam, and introducing a sinusoidal shape function wðxÞ¼W sinðx=LÞ as a test function into the Rayleigh's quotient leads to the natural frequencies from this approximated gradient elasticity solution. Equation (69) shows that the Finite Element model gives an upper bound of the \local" problem asymptotically found for n tending towards in¯nite.  Hermitian interpolation function appears to be the most e±cient numerical approach for this problem in comparison to both the¯rst-order and the higher-order¯nite di®erence approaches. Of course, the e±ciency of the¯nite element method strongly depends on the adopted displacement interpolation¯eld.

Concluding Remarks
The vibration behavior of elastic beams is studied using standard numerical methods such as¯nite di®erence or¯nite element methods. Some exact solutions of the¯nitedimensional discretized problems are presented for archetypal boundary conditions. The¯nite numerical approaches are then approximated by some enriched nonlocal or gradient elastic problems. The analogy between the¯nite numerical approaches and the equivalent enriched continuum is demonstrated, using a continualization procedure, which converts the discrete numerical problem into a continuous nonlocal or gradient one. For di®erent orders of¯nite di®erence (FD) schemes, di®erent nonlocal length scales are obtained, with respect to the step size a of the numerical scheme, i.e. l c ¼ a= ffiffi ffi 6 p % 0:408a for¯rst-order FD scheme as shown in Eq. (27), and l c ¼ a= ffiffiffiffiffi 45 4 p % a= ffiffiffiffiffiffiffiffiffiffiffi 6:708 p % 0:386a for higher-order FD scheme as shown in Eq. (49). In both cases, the nonlocal length scale of the dynamics problem is larger than the one of the statics problem. The¯nite di®erence schemes leads to a lower bound of the \local" continuous problem.
For the cubic Hermitian¯nite element, the discretized system is e±ciently approximated by a gradient elasticity law. In this case, the gradient elasticity length scale is the same than the one of the statics problem. The¯nite element formulation leads to an upper bound of the \local" continuous problem. The comparison between the exact numerical solution and the approximated nonlocal approaches shows the e±ciency of the continualization methodology. These analogies between enriched continuum and¯nite numerical schemes provide a new attractive framework for potential applications of enriched continua in computational mechanics.