# Revisiting finite difference and finite element methods applied to structural mechanics within enriched continua

Noël Challamel, Vincent Picandet, Bernard Collet, Thomas Michelitsch, Isaac

Elishakoff, C. M. Wang

► To cite this version:

Noël Challamel, Vincent Picandet, Bernard Collet, Thomas Michelitsch, Isaac Elishakoff, et al. Revisiting finite difference and finite element methods applied to structural mechanics within enriched continua. European Journal of Mechanics - A/Solids, 2015, 53, pp.107-120. 10.1016/j.euromechsol.2015.03.003. hal-01459729

HAL Id: hal-01459729

https://hal.science/hal-01459729

Submitted on 20 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Revisiting finite difference and finite element methods applied
to structural mechanics within enriched continua

Noël Challamel<sup>a, *</sup>, Vincent Picandet<sup>a</sup>, Bernard Collet<sup>b</sup>, Thomas Michelitsch<sup>b</sup>,
Isaac Elishakoff<sup>c</sup>, C.M. Wang<sup>d</sup>

<sup>a</sup> Université de Bretagne Sud, EA 4250, LIMATB (Laboratoire d'Ingénierie des MATériaux de Bretagne), F-56100 Lorient, France

<sup>b</sup> Sorbonne Universités, Université Pierre et Marie Curie (Paris 6), Institut Jean le Rond d'Alembert, CNRS UMR 7190, Tour 55-65, 4 place Jussieu,
75252 Paris cedex 05, France

<sup>c</sup> Department of Ocean and Mechanical Engineering, Florida Atlantic University, Boca Raton, FL 33431-0991, USA

<sup>d</sup> Engineering Science Programme and Department of Civil and Environmental Engineering, National University of Singapore, Kent Ridge, Singapore 119260,
Singapore

In this paper, we revisit the capability of numerical approaches such as finite difference methods and finite element methods, in approximating exact one-dimensional continuous eigenvalue problems (such as lateral vibrations of a string, the axial or the torsional vibrations of a bar, and the buckling of elastic columns). The numerical methods analysed in this paper are converted into difference equations. Following a continualization procedure or the method of differential approximation, the difference op- erators are then expanded in differential operators via Taylor expansion or Padé approximants. Analogies between the finite numerical approaches and some equivalent enriched continuum are shown, using this continualization procedure. The finite difference methods (first-order or higher-order finite difference methods) are shown to behave as integral-based nonlocal media (or stress gradient media), while the finite element method is found to behave as gradient elasticity media (or strain gradient media). The length scale identification of each equivalent enriched continuum strongly depends on the order of the numerical method considered. For the finite difference methods, the length scale identification of the equivalent nonlocal medium depends on the static versus dynamic analysis, whereas this length scale appears to be independent of inertia effects for the finite element method. Some comparisons between the exact discrete eigenvalue problems and the approximated continuous ones show the efficiency of the continualization procedure. An equivalent enriched Rayleigh quotient can be defined for each numerical method: the integral-based nonlocal method gives a lower bound solution to the exact eigenvalue multiplier, whereas the gradient elasticity method furnishes an upper bound solution.

# 1. Introduction

The source of discreteness in mechanics or physics may come
from the inherent nature of matter which is composed of a discrete
(or a finite) number of local repetitive cells. One can say that the
matter composition is essentially discrete, even if some equivalent
continua may be used as efficient engineering models for many
practical problems. Discrete-based medium may be classified as a

microstructured medium, and the microstructure can be related to
the atomistic composition at a subscale, but also to the molecule,
crystal or grain composition at larger scales, depending on the
problems (Mindlin, 1964). The microstructure may depend on the
scale analysis, from nanoscale up to macroscale characterization.
One way to incorporate the microstructure into a continuum
modeling is the consideration of enriched continua. Enriched
continua, which can be classified as non-simple materials, have a
wide spectrum of applications for modeling the constitutive laws of
microstructured media, from atomistic scale to large repetitive
structures. Among these enriched continua, gradient elasticity
(where the constitutive law depends on the gradient of the strain
variables) and nonlocal integral-based elasticity (where the
constitutive law depends on an integral of the strain variables) can

* Corresponding author.
E-mail addresses: noel.challamel@univ-ubs.fr (N. Challamel), vincent.picandet@univ-ubs.fr (V. Picandet), bernard.collet@upmc.fr (B. Collet), michel@lmm.jussieu.fr (T. Michelitsch), elishako@fau.edu (I. Elishakoff), ceewcm@nus.edu.sg (C.M. Wang).
be utilized, and sometimes both enriched contributions are com-
bined. Nowadays, the interest of such constitutive laws stems from
the desire for accurately capturing some specific effects such as
scale factor, boundary layers or localization phenomena, although
the micromechanics foundation of these enriched laws is still
debated.

Another source of discreteness in mechanics is related to the
numerical discretization of a continuous engineering model. Based
on this viewpoint, the discrete problem is considered as an
approximation of the continuous one, which now plays the role of
the reference system. The aim of the characterization of a numerical
scheme is to understand how the discrete problem is able to
approximate the continuous one. This alternative perspective is the
numerical evaluation of approximate methods such as the Finite
Difference Method or the Finite Element Method for computing the
solutions of continuous problems (see for instance Greenwood,
1961; Zienkiewicz and Cheung, 1967; Cyrus and Fulton, 1968;
Walz et al., 1968; Strang and Fix, 1973; Gawain and Ball, 1978;
Subrahmanyam and Kaza, 1983).

The comparison of discrete methods with some equivalent
continuous ones is not new. In fact, it was investigated by many
researchers since the 1950's for beam vibration problems; see for
instance papers by Livesley (1955) and Leckie and Lindberg (1963).
Moreover, Lagrange (1788) and later Rayleigh (1894) determined
the exact vibration frequencies of a string with a finite number of
concentrated masses and compared the solution with the contin-
uous system asymptotically obtained for an infinite number of cells
(see also Livesley, 1955 or Gantmacher, 1970 on this topic). This
discrete string problem will be analysed herein by using a con-
tinualization procedure. Continualization procedures are based on
various approximations of the discrete operators by some contin-
uous ones via Taylor expansion or Padé approximants (see Collins,
1981; Rosenau, 1986, 1987; Wattis, 2000; Kevrekidis et al., 2002;
Andrianov et al., 2010). The so-called enriched continuum equiva-
lent to the discrete one is sometimes called a quasi-continuum (see
for instance Collins, 1981), and is generally dependent on the
truncated terms in the asymptotic expansion of the differential
operators. This method was pioneered by Kruskal and Zabusky
(1964) for a nonlinear discrete axial chain, and initially applied to
discrete wave equations. In fact, it can be shown that the discrete
chain problem is firmly linked to the investigation of the finite
difference scheme of the “local” continuous axial problem, as it can
be viewed as the physical support of the finite difference scheme.
Hence, the microstructured problem and the Finite Difference
method may be strongly related (and can be viewed as mathe-
matically similar). With respect to the numerical schemes, the
method of differential approximation which was developed in the
1950's by Zhukov (1957) (see also Shokin, 1983) is based on the
differential representation of a difference scheme. It was applied
successfully for discrete time-spatial schemes, to build new dif-
ference schemes with well defined properties, or to better analyse
the existing numerical schemes. In the following, we are concerned
by discrete spatiality with continuous time dependence. Therefore,
only spatial discrete operators will be asymptotically expanded for
building equivalent nonlocal continua.

It has been already shown by Kunin (1982) and more recently
by Charlotte and Truskinovsky (2012) that a discrete chain may
behave dynamically as an equivalent nonlocal continuum whose
kernel depends on the level of approximation of the reference
discrete medium. However, the link between the differential
format of Eringen's nonlocal model (Eringen, 1983) is more recent
(see Wang et al., 2013). It has been recently argued that discrete
elastic systems (including bending and shear interaction law)
behave as nonlocal structural elements, in the context of the so-
called continualization procedure (Wang et al., 2013; Challamel

et al., 2014a; 2014b, 2014c, 2015c; Zhang et al., 2014a). This
fundamental property has been shown on microstructured
bending structural systems, and the induced nonlocality may be
based on Eringen's classical model (Eringen, 1983) for the consti-
tutive law. The methodology has been also applied to discrete
shear systems which can be equivalently investigated within
nonlocal shear Timoshenko theories (Duan et al., 2013; Zhang
et al., 2013 or Challamel et al., 2014c). Eringen's nonlocality has
been shown to be an efficient engineering tool for capturing the
out-of-plane behavior of discrete systems (Challamel et al., 2015a).
A similar reasoning has been recently considered in Zhang et al.
(2014b) for justifying nonlocal elastic plate models from a
microstructured beam-grid model (see also Zhang et al., 2015). In
other words, Eringen's nonlocality model has a rigorous physical
support based on some repetitive discrete systems. In this way,
nonlocal approaches can be treated as a relevant enriched con-
tinuum for capturing scale effects in discrete systems. Within this
point of view, the enriched continuum is considered as an
approximate theory for approaching the main phenomena of the
reference discrete system.

In this paper, a similar continualization (or method of differ-
ential approximation) reasoning will be followed for approxi-
mating some finite numerical schemes (finite difference methods
and finite element methods) by some equivalent enriched con-
tinua. We have already shown for bending systems that finite
difference formulation (or its analogical lattice physical system)
has resemblance to nonlocal theory, whereas finite element
formulation may have a strain gradient continuum analogy (see
recently Challamel et al., 2015b for the analysis of beam vibra-
tions). We have used this analogy to calibrate the small length
terms of nonlocal theories with respect to the reference lattice
systems. Most of the papers devoted to this question in the past
were focused on bending systems, but some questions were
unanswered up to now:

- What about higher-order finite difference methods? Will they behave as higher-order nonlocal models (or higher-order stress gradient models)?

- What about finite element models? Will they behave as stress gradient or strain gradient models, and how will the enriched continuum depend on the choice of the shape function?

- Is the nonlocal length scale of the finite difference methods or the finite element methods inertia dependent (does it depend on the static versus dynamic analysis)?

- Is the equivalent enriched continuum dependent on the kine-
matics of the problem, for instance when considering the
bending, axial or torsional motion of a beam?

In the paper, we give some answers to these open questions,
from the analysis of exact one-dimensional discrete eigenvalue
problems (such as lateral vibrations of a string, the axial or the
torsional vibrations of a bar, and the buckling of elastic columns).
We show that the rate of convergence is strongly dependent on the
order of the finite discrete scheme: higher-order finite schemes
lead to higher-order enriched constitutive laws with a higher
convergence rate. For historical reasons, we first start from the
string problem and then investigate some column problems in
buckling by both the finite difference and finite element methods.

## 2. Discrete string problem

Starting from equilibrium considerations, the vibration fre-
quencies ω of a string of length L are obtained from the linear
second-order differential equation of the transverse deflection w
(see Lagrange 1788, Leissa and Qatu, 2011; Wang and Wang, 2013):
$$Tw'' + \mu\omega^2 w = 0 \qquad (1)$$

where $\mu = \rho S$ is the line density, $T$ the tensile force, $\rho$ the volumetric density, $S$ the cross sectional area and the prime denotes differentiation with respect to the axial coordinate $x$. For fixed end conditions, we have $w(0) = 0$ and $w(L) = 0$. By assuming a sinusoidal vibration mode $w = \bar{w} \sin(k\pi x/L)$, the vibration frequency is given by

$$\omega_{k,\infty} = \frac{k\pi}{L} \sqrt{\frac{T}{\mu}} \qquad (2)$$

where $k$ is an integer, the subscript “$\infty$” is related to the continuous string with an infinite number of degrees-of-freedom. If now, the string is assumed to compose a finite number $n + 1$ of equidistantly spaced masses ($L = n \times a$), the continuous problem is converted into a discrete problem governed by the following linear second-order finite difference equation

$$T \frac{w_{i+1} - 2w_i + w_{i-1}}{a^2} + \mu\omega^2 w_i = 0 \qquad (3)$$

where $a$ is the distance between the nodes. Equation (3) is the central finite difference approximation of the continuous problem defined by Eq. (1). On the other hand, Eq. (3) can be also considered as a true discrete problem, as introduced by Lagrange (1788) for the vibration frequencies of a system with a finite number of concentrated masses. It can be shown that the frequency equation of this discrete problem with the fixed–fixed boundary conditions $w_0 = w_n = 0$ can be analytically solved and the exact frequencies are given by (Lagrange, 1788):

$$\omega_{k,n} = \frac{2n}{L} \sin\left(\frac{k\pi}{2n}\right) \sqrt{\frac{T}{\mu}} \qquad (4)$$

In fact, the linear second-order difference equation can also be presented as

$$w_{i+1} + \left(\frac{\beta}{n^2} - 2\right)w_i + w_{i-1} = 0 \quad \text{with} \quad w_0 = 0 \quad \text{and} \quad w_n = 0 \qquad (5)$$

where $\beta = \mu\omega^2 L^2 / T$ is the dimensionless frequency parameter (or eigenvalue parameter). The solution of this finite difference boundary value problem is derived by Goldberg (1958). Goldberg (1958) also mentioned that this second-order finite difference boundary value problem arises in the mathematical theory of scale analysis (see also Guttman, 1950). We reproduce the main reasoning of the mathematical proof for the investigation of this discrete eigenvalue problem.

The characteristic equation is obtained by replacing $w_i = A\lambda^i$ in Eq. (5) which leads to

$$\lambda + \frac{1}{\lambda} = 2 - \frac{\beta}{n^2} \qquad (6)$$

Equation (6) is symmetrical with respect to interchanging $\lambda$ and $1/\lambda$. This equation admits the following two solutions for $|1 - \beta/2n^2| < 1$:

$$\lambda_{1,2} = \cos \phi \pm j \sin \phi \text{ with } \phi = \arccos\left(1 - \frac{\beta}{2n^2}\right) \quad \text{and} \quad j = \sqrt{-1} \qquad (7)$$

The solution of the finite difference equation can be expressed with the real basis as:

$$w_i = A \cos(\phi i) + B \sin(\phi i) \qquad (8)$$

The introduction of the two boundary conditions $w_0 = 0$ and $w_n = 0$ leads finally to the natural vibrations mode $w_i = B\sin(\phi i)$ for the fixed–fixed string with the natural frequency formulae:

$$\sin(\phi n) = 0 \Rightarrow \phi n = k\pi \Rightarrow \cos\frac{k\pi}{n} = 1 - \frac{\beta}{2n^2} \Rightarrow \beta = 4n^2 \sin^2\left(\frac{k\pi}{2n}\right) \qquad (9)$$

which is exactly the frequency value reported in Eq. (4). As detailed in Goldberg (1958), it can be shown that no nontrivial solution of the boundary value problem exists for $|1 - \beta/2n^2| \ge 1$.

Rayleigh (1894) and Livesley (1955) showed by using asymptotic expansion that the discrete equivalent system has slightly lower natural frequencies when compared to its continuum counterpart:

$$\frac{\omega_{k,n}}{\omega_{k,\infty}} = \frac{2n}{k\pi} \sin\left(\frac{k\pi}{2n}\right) = 1 - \frac{k^2\pi^2}{24n^2} + o\left(\frac{1}{n^4}\right) \le 1 \qquad (10)$$

This lower bound property of the discrete approach with respect to the continuous one could be argued a consequence of the nonlocal continuum property of the discrete system. In fact, the discrete equations can be 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 sufficiently smooth deflection function as:

$$w(x + a) = \sum_{k=0}^{\infty} \frac{a^k \partial_x^k}{k!} w(x) = e^{a \partial_x} w(x) \text{ with } \partial_x = \frac{\partial}{\partial x} \qquad (11)$$

The pseudo-differential operators can be introduced as (see Rosenau, 1986; Rosenau, 1987; Wattis, 2000; Kevrekidis et al., 2002; Andrianov et al., 2010):

$$\begin{aligned} \frac{w_{i-1} + w_{i+1} - 2w_i}{a^2} &= \frac{[e^{a\partial_x} + e^{-a\partial_x} - 2]}{a^2} w(x) \\ &= \frac{4}{a^2}\sinh^2\left(\frac{a}{2}\partial_x\right) w(x) \\ &= \frac{\partial_x^2}{1 - l_c^2\partial_x^2} w(x) + \dots \text{with } l_c^2 = \frac{a^2}{12} \end{aligned} \qquad (12)$$

As a consequence, the continualization approach based on Padé's approximant Eq. (12) and applied to the discrete system Eq. (3) gives a truncated nonlocal approach for the string dynamics (see Rosenau, 1986 or 1987 for the vibration model of the axial chain):

$$Tw'' + \mu\omega^2 w - \mu\omega^2 l_c^2 w'' = 0 \quad \text{with} \quad l_c^2 = \frac{a^2}{12} \qquad (13)$$

where $a = L/n$ is the size of the repetitive cell and the underlined term is related to the additional nonlocal effect. The vibration frequency of the fixed–fixed nonlocal string is then obtained by introducing the sinusoidal mode shape in Eq. (13), thereby leading to:

$$\frac{\omega_k^2}{\omega_{k,\infty}^2} = \frac{1}{1 + \frac{k^2 \pi^2 l_c^2}{L^2}} = 1 - \frac{k^2 \pi^2}{12n^2} + o\left(\frac{1}{n^4}\right) \qquad (14)$$

which is the same approximation as the one for the discrete problem given by Eq. (10). This does not contradict an assertion that
the discrete string behaves as a nonlocal string. Note that for this problem, no constitutive elasticity law is used for deriving the equations but only the equilibrium equations are employed. It is seen that the analysis of discrete systems within nonlocal continuous media can be handled through the string vibration model (which is formally analogous to the vibration of an axial rod) – an old problem exactly characterized by Lagrange for both the discrete version and its continuum limit. However, naturally, at that early time, the discussion on the properties of the discrete system in term of nonlocal continuum was not present.

The present analysis for the string problem can be easily extended for analysing the dynamics behavior of an axial bar or a torsional bar. In fact, the string problem is mathematically analogous to the axial chain problem, as observed from the second-order finite difference equation associated with the Born–von Kármán mono-atomic model (Born and von Kármán, 1912):

$$ES \frac{u_{i+1} - 2u_i + u_{i-1}}{a^2} + \omega^2 \rho S u_i = 0 \quad (15)$$

where $u_i$ is the axial displacement of the ith mass, E the Young modulus, S the equivalent cross sectional area and $\rho$ the density. Equation (15) can be seen as the finite difference form of the axial vibrations equations of a continuous bar, or the vibration equations of a lattice spring axial model. It is clear that Eq. (15) is analogous to Eq. (3). For similar reasons, the frequency equation of this discrete axial chain problem with the fixed–fixed boundary conditions $u_0 = 0$ and $u_n = 0$ can be analytically obtained as (see also Tong et al., 1971; Thomson, 1972 or Blevins, 2001 for this discrete axial problem):

$$\omega_{k,n} = \frac{2n}{L} \sin\left(\frac{k\pi}{2n}\right) \sqrt{\frac{E}{\rho}} \quad (16)$$

By using a nonlocal approach, the enriched vibration equation may be given in this axial case by Rosenau (1986, 1987):

$$ESu'' + \rho S \omega^2 u - \rho S \omega^2 \frac{l_c^2}{E} u'' = 0 \quad \text{with} \quad l_c^2 = \frac{a^2}{12} \quad (17)$$

This is the low frequency matching of the nonlocal approximated model with respect to the discrete one. Equations (10) and (14) are still valid for the axial case. The high frequency matching based on $k = n$ leads to the matching value $l_c^2/a^2 = 1/4 - 1/\pi^2$ given by Eringen (1983) – see Andrianov et al. (2010) or more recently Challamel et al. (2014b). Equation (17) is an enriched dynamics equation which may come from the nonlocal Eringen's differential formulation applied to the axial bar (Challamel et al., 2009; Aydogdu, 2009):

$$N - l_c^2 N'' = ESu' \quad \text{and} \quad N' = \rho S \bar{u} \quad (18)$$

The Eringen's type nonlocal model (Eringen, 1983) is recognized in the normal force N -strain $\epsilon = u'$ differential equation, and the nonlocal axial wave equation is written in this case as (Rosenau, 1986):

$$ESu'' + \frac{l_c^2}{E} \rho S \bar{u}'' = \rho S \bar{u} \quad (19)$$

This kind of “nonlocal” wave equation was already obtained by Rayleigh (1894) when including the kinetic energy due to lateral motion (the nonlocal length scale is related to the cross sectional properties for Rayleigh wave analogy). The same reasoning may be followed for the torsional vibration of a discrete elastic bar, governed by:

$$GJ \frac{\phi_{i+1} - 2\phi_i + \phi_{i-1}}{a^2} + \omega^2 \rho J \phi_i = 0 \quad (20)$$

where $\phi_i$ is the torsional angle, G the shear modulus, and J the torsional section constant. As previously developed, the frequency equation of this discrete torsional chain problem with the fixed–fixed boundary conditions $\phi_0 = 0$ and $\phi_n = 0$ can be given by:

$$\omega_{k,n} = \frac{2n}{L} \sin\left(\frac{k\pi}{2n}\right) \sqrt{\frac{G}{\rho}} \quad (21)$$

Similarly, the nonlocal torsional model which may be seen as an approximation of the discrete torsional problem may be written, using the same length scale calibration, from:

$$M_t - l_c^2 M_t'' = GJ\phi' \quad \text{and} \quad M_t' = \rho J \ddot{\phi} \quad (22)$$

One recognizes the nonlocal Eringen's differential formulation applied to the torsional problem, as investigated by Narendar (2011), Islam et al. (2014) or Arda and Aydogdu (2014). In summary, the wave equation governing the vibrations of a string is mathematically similar to the wave equation of an axial or a torsional elastic rod, whose discretized version using a finite difference scheme has been shown to be equivalent to a kind of nonlocal wave model. Both the discrete and the continuous problems of the string were already solved by Lagrange (1788), but the relationship between discrete elastic systems and enriched continuum mechanics is pointed out more recently, Rosenau (1986) already obtained the enriched wave equation for this one-dimensional problem, but the link between Eringen's nonlocality was not pointed out. Only recently, discrete microstructured systems which may be equivalent to a finite difference formulation of a continuous problem have been shown to behave as nonlocal elastic systems, in the Eringen sense (see for instance Challamel et al., 2014b).

In the sequel, the beam problem will be mainly investigated, and the specificities of the numerical schemes applied to an elastic buckling problem will be analysed. As already discussed by Silverman (1951), Hencky's chain (Hencky, 1920a, 1920b) – which is the lattice bending system (or microstructured model, also called discrete bending system) – is in fact strictly equivalent to the finite difference formulation of a continuous problem, i.e. the Euler–Bernoulli continuous beam problem. The performance of Finite Difference Method for solving buckling or vibrations eigenvalue problems has already been evaluated in the literature (see early studies by Salvadori, 1951; Wang, 1951, 1953; Seide, 1975, and more recent papers by Elishakoff, 1998; Santoro and Elishakoff, 2006), but without resorting to any nonlocal mechanics perspective (except recently in Challamel et al., 2014b). A consequence of the nonlocal equivalent principle for the modeling of discrete systems is that the finite difference system can be efficiently approached by nonlocal continuum mechanics tools. As it is known in the case of nonlocal mechanics behaviors, this result confirms the lower bound solution of such approximate Finite Difference 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.

### 3. Finite difference method – buckling

Consider the buckling problem of a local Euler–Bernoulli column axially loaded by an axial force P. The equilibrium equation is given by
$$M'' = -Pw'' \qquad (23)$$

where M is the bending moment, w the deflection of the beam, and the prime denotes differentiation with respect to x. The elastic bending constitutive law is given by the relationship:

$$M = EIw'' \qquad (24)$$

The finite difference approximation (central difference method) of the equilibrium equations is obtained from Eq. (23) as:

$$\frac{M_{i-1} - 2M_i + M_{i+1}}{a^2} + P \frac{w_{i-1} - 2w_i + w_{i+1}}{a^2} = 0 \qquad (25)$$

where a is the constant step of the finite difference approximation, and the elastic bending constitutive law is rewritten with this approximation as:

$$M_i = EI \frac{w_{i-1} - 2w_i + w_{i+1}}{a^2} \qquad (26)$$

By combining Eq. (25) and Eq. (26), we have the fourth-order finite difference equation:

$$EI \frac{w_{i+2} - 4w_{i+1} + 6w_i - 4w_{i-1} + w_{i-2}}{a^4} + P \frac{w_{i+1} - 2w_i + w_{i-1}}{a^2} = 0 \qquad (27)$$

For the buckling problem of a column with simply supported boundary conditions, it is possible to integrate this equation twice, thus leading to a second-order finite difference equation:

$$EI \frac{w_{i+1} - 2w_i + w_{i-1}}{a^2} + Pw_i = 0 \qquad (28)$$

This is the same finite difference equation as that obtained for axial vibrations. In fact, the linear second-order difference equation can also be presented in the form as in Eq. (5):

$$w_{i+2} + \left(\frac{\beta}{n^2} - 2\right)w_{i+1} + w_i = 0 \quad \text{with} \quad \beta = \frac{PL^2}{EI} \qquad (29)$$

By using the solution already developed in Eq. (9), the fundamental buckling load of this discrete system is obtained for $k = 1$ from the exact critical buckling load formulae (see also Wang, 1951):

$$\frac{P_{cr,n}}{P_{cr,\infty}} = \frac{\beta_n}{\beta_\infty} = \frac{4n^2}{\pi^2} \sin^2\left(\frac{\pi}{2n}\right) = 1 - \frac{\pi^2}{12n^2} + o\left(\frac{1}{n^4}\right) \qquad (30)$$

The continualization of the discrete equilibrium Eq. (25) leads to the unchanged “local” equilibrium Eq. (23), whereas the continualization of the discrete constitutive laws Eq. (26) leads to the Eringen's type nonlocal constitutive law:

$$M - l_c^2 M'' = EIw'' \quad \text{with} \quad l_c^2 = \frac{a^2}{12} \qquad (31)$$

By inserting Eq. (31) into Eq. (23), the equations of the uncoupled buckling problem reduce to

$$(EI - P l_c^2) w^{(4)} + P w'' = 0 \qquad (32)$$

Equation (32) is equivalent to considering an Eringen's based nonlocal model being applied to an Euler-Bernoulli beam kinematics (omitting the discussion on the boundary conditions - see Challamel et al., 2014d for the specific role of boundary conditions in nonlocal mechanics). As it is known, Eringen's nonlocal model leads to the lower bound solution of the buckling load with respect

to its local counterpart (Elishakoff et al., 2012; Challamel, 2013). This can be easily shown from the Rayleigh's quotient of Eringen's Euler-Bernoulli column model:

$$R = \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} w'^2 + l_c^2 w''^2 dx} \le \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} w'^2 dx} \qquad (33)$$

This Rayleigh's quotient R can be used to compute approximate upper bounds or exact solution of the variationally-based eigenvalue problem (see Rayleigh, 1894; Temple and Bickley, 1956 or more recently Bažant and Cedolin, 2003 or Ilanko and Monterrubio, 2015). For example, consider a pinned-pinned column, and introducing a sinusoidal shape function $\bar{w} = \bar{w} \sin(\pi x/L)$ as an admissible function into the Rayleigh's quotient, we get the exact buckling value of the pinned-pinned case as:

$$\tilde{w} = \bar{w} \sin\left(\frac{\pi x}{L}\right) \Rightarrow P_{cr} = \frac{\frac{\pi^2 EI}{L^2}}{1 + \frac{\pi^2 l_c^2}{L^2}} \le \frac{\pi^2 EI}{L^2} \quad \text{with} \quad l_c^2 = \frac{a^2}{12} \qquad (34)$$

It turns out that the buckling load of the finite difference system with respect to its continuous “local” counterpart can be expressed for the pinned-pinned column as:

$$\frac{P_{cr,n}}{P_{cr,\infty}} = \frac{1}{1 + \frac{\pi^2}{12n^2}} = 1 - \frac{\pi^2}{12n^2} + o\left(\frac{1}{n^4}\right) \qquad (35)$$

Equation (35) is numerically shown by Seide (1975) within the central finite difference method. Note that the Rayleigh's quotient can be derived from the Lagrangian equivalence obtained by continualization of the discrete Lagrangian, following a methodology introduced by Rosenau (2003) for axial systems.

The Rayleigh's quotient of the finite difference approximate system is given by

$$R = \frac{\sum_{i=0}^{n} EI \left( \frac{w_{i+1} - 2w_i + w_{i-1}}{a^2} \right)^2}{\sum_{i=0}^{n} \left( \frac{w_{i+1} - w_i}{a} \right)^2} \qquad (36)$$

Consider now the regular deflection variable w defined by

$$\partial_x^2 w = \frac{w_{i+1} - 2w_i + w_{i-1}}{a^2} \qquad (37)$$

This condition can be inverted, i.e.

$$w_i = \frac{a^2 \partial_x^2}{4 \sinh^2\left(\frac{a}{2} \partial_x\right)} w \qquad (38)$$

The potential energy can be expanded with this new variable as (see Appendix A):

$$\sum_{i=0}^{n} a \left( \frac{w_{i+1} - w_i}{a} \right)^2 = \int_{0}^{L} w'^2 + \frac{a^2}{12} w''^2 + o(a^4) dx - \frac{a^2}{12} [w'w'']_0^L \qquad (39)$$

which shows indeed, except for the boundary terms, the correspondence between the potential energy of the discrete system and the one of the nonlocal continuous system. This closes the proof from the Rayleigh's quotient of the discrete system to the one of the continuous enriched one. Note that the Rayleigh's quotient provides information on the rate of convergence of the discrete problem when compared to the exact solution. The length scale
calibration of the nonlocal model equivalent to the finite difference formulation in buckling is then equal to:

$$l_{c,\text{static}}^2 = \frac{a^2}{12} = \frac{l_{c,\text{dynamic}}^2}{2} \qquad (40)$$

As shown by Challamel et al. (2015b), a factor 2 affects the equivalent length scale of the finite difference formulation applied to the vibrations analysis, which shows that the nonlocal model associated with the finite difference formulation, depends on the static versus vibration analysis. Furthermore, using the generalized Rayleigh ratio, and for some other standard boundary conditions (including also the clamped-clamped, or clamped-free boundary conditions), the general formulae can be applied as well (see Challamel, 2013):

$$\frac{P_{cr,n}}{P_{cr,\infty}} = \frac{1}{1 + \frac{P_{cr,\infty}l_c^2}{EI}} \quad \text{with} \quad l_c^2 = \frac{a^2}{12} \qquad (41)$$

which is consistent with the asymptotic results of Seide (1975) for these other boundary conditions.

Some other discrete models can be investigated in a similar way, possibly introducing some numerically-based nonlocal equivalent laws. For instance, the model of Myklestad derived from the three-moment equation (see Myklestad, 1944; Leckie and Lindberg, 1963) leads to the finite difference equation:

$$\frac{M_{i+1} + 4M_i + M_{i-1}}{6} = EI \frac{w_{i+1} - 2w_i + w_{i-1}}{a^2} \qquad (42)$$

which leads to the continualization nonlocal law:

$$\frac{M}{1 - \frac{a^2 \partial_x^2}{6}} = EI \frac{\partial_x^2}{1 - \frac{a^2 \partial_x^2}{12}} w \Rightarrow M - \frac{a^2}{12} M'' = EI \left( w'' - \frac{a^2}{6} w^{(4)} \right) \qquad (43)$$

that can be reinterpreted as a mixed nonlocal law, as used for instance by Challamel and Wang (2008) at the beam scale.

Another approximation of Eq. (43) gives:

$$M = EI(w'' - l_c^2 w^{(4)}) \quad \text{with} \quad l_c^2 = \frac{a^2}{12} \qquad (44)$$

which is typically a gradient elasticity bending-curvature law (see Papargyri-Beskou et al., 2003) associated with the linear differential equation (the equilibrium equations are assumed to be unchanged at this stage):

$$-EI l_c^2 w^{(6)} + EIw^{(4)} + Pw'' = 0 \qquad (45)$$

As it is known, the gradient elasticity model leads to the upper bound solution of the buckling load with respect to its local counterpart (Elishakoff et al., 2012; Challamel, 2013). This can be easily shown from the Rayleigh's quotient of gradient elasticity model:

$$R = \frac{\int_{0}^{L} EI(w'' + l_c^2 w''^2) dx}{\int_{0}^{L} w'^2 dx} \geq \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} w'^2 dx} \qquad (46)$$

Consider again the example of a pinned-pinned column and introducing a sinusoidal shape function $\tilde{w} = \bar{w} \sin(\pi x / L)$ as an admissible function into the Rayleigh's quotient leads to the exact

buckling value of the pinned-pinned case (with the compatible higher-order boundary conditions):

$$\tilde{w} = \bar{w} \sin\left(\frac{\pi x}{L}\right) \Rightarrow P_{cr} = \frac{\pi^2 EI}{L^2} \left(1 + \pi^2 \frac{l_c^2}{L^2}\right) \geq \frac{\pi^2 EI}{L^2} \qquad (47)$$

Again, using the generalized Rayleigh ratio applied to the equivalent gradient elasticity system with clamped-clamped, or clamped-free boundary conditions, the general formulae can be applied as well (see Challamel, 2013):

$$\frac{P_{cr,n}}{P_{cr,\infty}} = 1 + \frac{P_{cr,\infty} l_c^2}{EI} \quad \text{with} \quad l_c^2 = \frac{a^2}{12} \qquad (48)$$

An improved finite difference analysis may be based on the introduction of the second-order central difference for the expressions of the first and the second derivatives of the displacement (see Greenwood, 1961; Subrahmanyam and Kaza, 1983; Subrahmanyam and Leissa, 1985 or more recently Elishakoff and Santoro, 2005). The higher-order finite difference formulation of the equilibrium equations now reads as

$$\frac{-M_{i-2} + 16M_{i-1} - 30M_i + 16M_{i+1} - M_{i+2}}{12a^2} + P \frac{-w_{i-2} + 16w_{i-1} - 30w_i + 16w_{i+1} - w_{i+2}}{12a^2} = 0 \qquad (49)$$

where a is the constant step of the finite difference approximation, and the constitutive law is rewritten with this approximation as:

$$M_i = EI \frac{-w_{i-2} + 16w_{i-1} - 30w_i + 16w_{i+1} - w_{i+2}}{12a^2} \qquad (50)$$

For the pinned-pinned column, the finite difference equation to be solved with this improved finite difference analysis is given by

$$EI \frac{-w_{i-2} + 16w_{i-1} - 30w_i + 16w_{i+1} - w_{i+2}}{12a^2} + Pw_i = 0 \qquad (51)$$

The characteristic equation is obtained by replacing $w_i = A\lambda^i$ in Eq. (51) which leads to

$$-(\lambda + \frac{1}{\lambda})^2 + 16(\lambda + \frac{1}{\lambda}) - 28 + \frac{12\beta}{n^2} = 0 \quad \text{with} \quad \beta = \frac{PL^2}{EI} \qquad (52)$$

that admits the following four solutions

$$\begin{align*}
\lambda_{1,2} &= \cos \phi \pm j \sin \phi \quad \text{and} \\
\lambda_{3,4} &= 8 - \cos \phi \pm \sqrt{(8 - \cos \phi)^2 - 1} \quad \text{with} \\
\phi &= \arccos \left(4 - \sqrt{9 + \frac{3\beta}{n^2}}\right) \qquad (53)
\end{align*}$$

For the pinned-pinned column, the buckling mode is obtained from the trigonometric shape function $w_i = B \sin(\phi i)$ thus furnishing the exact buckling load formulae:

$$\begin{align*}
\cos \frac{\pi}{n} &= 4 - \sqrt{9 + \frac{3\beta}{n^2}} \Rightarrow \frac{P_{cr,n}}{P_{cr,\infty}} = \frac{\beta_n}{\beta_\infty} = \frac{n^2}{3\pi^2} \left[ 7 - 8 \cos \frac{\pi}{n} + \cos^2 \frac{\pi}{n} \right] \\
&= 1 - \frac{\pi^4}{90n^4} + o\left(\frac{1}{n^6}\right) \qquad (54)
\end{align*}$$
The higher-order pseudo-differential operator can also be expanded as:

$$ 
\begin{aligned}
& \frac{-w_{i-2} + 16w_{i-1} - 30w_i + 16w_{i+1} - w_{i+2}}{12a^2} \\
&= \frac{\left[-e^{-2a\partial_x} + 16e^{-a\partial_x} - 30 + 16e^{a\partial_x} - e^{2a\partial_x}\right]}{12a^2} w(x) \\
&= \partial_x^2 \left(1 - \frac{a^4 \partial_x^4}{90}\right) w(x) + \dots
\end{aligned}
\quad (55) $$

In other words, the bending moment – curvature enriched law associated with the higher-order finite difference scheme is written as:

$$ M + l_c^4 M^{(4)} = EIw'' \quad \text{with} \quad l_c^4 = \frac{a^4}{90} \quad (56) $$

To the authors' knowledge, this nonlocal law has not been reported before. By inserting Eq. (56) into Eq. (23), the equations of the uncoupled buckling problem reduce to

$$ P l_c^4 w^{(6)} + E I w^{(4)} + P w'' = 0 \quad (57) $$

Equation (57) is equivalent to considering a nonlocal model being applied to an Euler-Bernoulli beam kinematics; the nonlocal model leads to the lower bound status of the buckling load with respect to its local counterpart. This can be easily shown again from the Rayleigh's quotient of Eringen model:

$$ R = \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} w'^2 + l_c^4 w'''^2 dx} \le \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} w'^2 dx} \quad (58) $$

Consider again the pinned-pinned column, and introducing a sinusoidal shape function $\tilde{w} = \bar{w} \sin(\pi x / L)$ as an admissible function into the Rayleigh's quotient leads to the exact buckling value of the pinned-pinned case:

$$ \tilde{w} = w_0 \sin\left(\frac{\pi x}{L}\right) \Rightarrow P_{cr} = \frac{\frac{\pi^2 EI}{L^2}}{1 + \pi^4 \frac{l_c^4}{L^4}} \le \frac{\pi^2 EI}{L^2} \quad \text{with} \quad l_c^4 = \frac{a^4}{90} \quad (59) $$

In this case, the buckling load of the higher-order finite difference system with respect to its continuous “local” counterpart can be expressed for the pinned-pinned column as:

$$ \frac{P_{cr,n}}{P_{cr,\infty}} = \frac{1}{1 + \frac{\pi^4}{90n^4}} = 1 - \frac{\pi^4}{90n^4} + o\left(\frac{1}{n^6}\right) \quad (60) $$

which, of course, gives a higher rate of convergence than for the central difference method. The length scale calibration of the nonlocal model equivalent to the higher-order finite difference formulation in buckling is then equal to:

$$ l_{c,\text{static}}^4 = \frac{a^4}{90} = \frac{l_{c,\text{dynamic}}^4}{2} \quad (61) $$

As shown by Challamel et al. (2015b), a factor 2 affects the equivalent length scale of the higher-order finite difference formulation applied to the vibrations analysis, which shows that the nonlocal model associated with the finite difference formulation, depends on the static versus vibration analysis. Furthermore, using the generalized Rayleigh ratio, a new general formulae can be

applied as well for some other boundary conditions (including the clamped-free and the clamped-clamped boundary conditions):

$$ \frac{P_{cr,n}}{P_{cr,\infty}} = \frac{1}{1 + \left(\frac{P_{cr,\infty}}{EI}\right)^2} \quad \text{with} \quad l_c^4 = \frac{a^4}{90} \quad (62) $$

It is clearly shown that the buckling load is more sensitive to the cell size, which is affected by a fourth-order power in the fraction. This result is not surprising as the efficiency of higher-order finite difference schemes is already known in the literature (see Greenwood, 1961; Subrahmanyam and Kaza, 1983; Subrahmanyam and Leissa, 1985 or more recently Elishakoff and Santoro, 2005). The analysis of these higher-order schemes in terms of continualized enriched continuum contributes here to a new understanding of this efficiency. This is clearly related to the nonlocal nature of the equivalent continuous medium, involving the fourth-order derivative of the generalized stress variable.

A similar result is obtained from the following discrete equations issued of a mixed variational principle, as detailed by Seide (1975):

$$ \left\{ \begin{aligned} EI \frac{\theta_{i+1} - \theta_{i-1}}{2a} - \frac{M_{i+1} + 4M_i + M_{i-1}}{6} &= 0 \\ \frac{M_{i+1} - M_{i-1}}{2a} + P \frac{\theta_{i+1} + 4\theta_i + \theta_{i-1}}{6} &= 0 \end{aligned} \right. \quad (63) $$

where $\theta_i = w'_i$ is the discrete rotation field. This two-dimensional discrete system can also be expressed with the difference operators $\delta_0$ and $\delta_1$:

$$ \left\{ \begin{aligned} EI\delta_1\theta - \delta_0M &= 0 \\ \delta_1M + P\delta_0\theta &= 0 \end{aligned} \right. \quad (64) $$

where the difference operators are defined by

$$ \delta_0 = \frac{1}{6} \left[ e^{-a\partial_x} + 4 + e^{a\partial_x} \right] \quad \text{and} \quad \delta_1 = \frac{1}{2a} \left[ -e^{-a\partial_x} + e^{a\partial_x} \right] \quad (65) $$

It is possible to use present the finite difference system in terms of the kinematics variable:

$$ \begin{aligned} EI\delta_1^2\theta + P\delta_0^2\theta &= 0 \Rightarrow EI \frac{\theta_{i+2} - 2\theta_i + \theta_{i-2}}{4a^2} \\ & \quad + P \frac{\theta_{i+2} + 8\theta_{i+1} + 18\theta_i + 8\theta_{i-1} + \theta_{i-2}}{36} \\ &= 0 \end{aligned} \quad (66) $$

An asymptotic expansion of the two difference finite operators gives:

$$ \begin{aligned} \delta_0^2 &= 1 + \frac{(a\partial_x)^2}{3} + \frac{(a\partial_x)^4}{18} + o(a^6) \quad \text{and} \quad \delta_1^2 \\ &= \left[ 1 + \frac{(a\partial_x)^2}{3} + \frac{2(a\partial_x)^4}{45} + o(a^6) \right] \partial_x^2 \end{aligned} \quad (67) $$

Finally, the finite difference system can be approximated by the following differential equation:

$$ EI\theta'' + P\left[1 + \frac{(a\partial_x)^4}{90}\right]\theta = 0 \quad (68) $$

which is equivalent to the differential equation Eq. (57) by setting $\theta = w'$. In this case, the mixed variational principle and the higher-order finite difference method lead to equivalent first-order
numerical results, associated with an equivalent linear differential equation. This continualization analysis of the mixed variational principle is consistent with the numerical results presented by Seide (1975) for the buckling of pinned-pinned columns, or other standard boundary conditions already mentioned.

### 4. Finite element method – buckling

In this section, the status of the finite element method will be investigated by using higher-order constitutive laws. For the buckling problem, and following the reasoning of Seide (1975) (called finite element method la in Seide, 1975), we will show that the linear interpolation shape form of the rotation function leads to a gradient elasticity constitutive law. The linear interpolation field is assumed as:

$$ egin{cases} \frac{dw}{d\xi} = \theta_{i-1}(1-\xi) + \theta_i\xi \\ \frac{d^2w}{d\xi^2} = -\theta_{i-1} + \theta_i \end{cases} \quad (69) $$

where $\xi = x/a$ is the dimensionless abscissa and the notation $\theta_i = w'_i$ is used for the discrete rotation field. The buckling load can be computed from the Rayleigh's quotient (see also Seide, 1975):

$$ R = \frac{\sum_{i=1}^{n} \int_{0}^{1} \frac{EI}{a^4} \left( \frac{d^2w}{d\xi^2} \right)^2 d\xi}{\sum_{i=1}^{n} \int_{0}^{1} \frac{1}{a^2} \left( \frac{dw}{d\xi} \right)^2 d\xi} \quad (70) $$

For the linear interpolation shape function considered, this Rayleigh's quotient is calculated as:

$$ R = \frac{EI \sum_{i=1}^{n} (\theta_i - \theta_{i-1})^2}{a^2 \sum_{i=1}^{n} \frac{\theta_i^2 + \theta_i \theta_{i-1} + \theta_{i-1}^2}{3}} \quad (71) $$

Consider the regular deflection variable w defined by

$$ (\partial_x w)^2 = \frac{\theta_i^2 + \theta_i \theta_{i-1} + \theta_{i-1}^2}{3} \quad (72) $$

This condition can be inverted (see Appendix B), i.e.

$$ \theta_i = \left[ 1 - \frac{(a \partial_x)^2}{12} \right] \partial_x w \quad (73) $$

The strain energy can be expanded with this new variable as (see Appendix B):

$$ \sum_{i=1}^{n} a \left( \frac{\theta_i - \theta_{i-1}}{a} \right)^2 = \int_{0}^{L} w''^2 + \frac{a^2}{12} w'''^2 + o(a^4) dx - \frac{a^2}{12} \left[ w'' w''' \right]_0^L \quad (74) $$

For this finite element model, the Rayleigh's quotient can then be continualized as:

$$ R = \frac{\int_{0}^{L} EI(w'' + l_c^2 w'''^2) dx}{\int_{0}^{L} w'^2 dx} \ge \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} w'^2 dx} \quad \text{with} \quad l_c^2 = \frac{a^2}{12} \quad (75) $$

leading to the gradient elasticity solution:

$$ -EI l_c^2 w^{(6)} + EIw^{(4)} + Pw'' = 0 \quad \text{with} \quad l_c^2 = \frac{a^2}{12} \quad (76) $$

associated with the gradient elasticity constitutive law:

$$ M = EI(w'' - l_c^2 w^{(4)}) \quad \text{with} \quad l_c^2 = \frac{a^2}{12} \quad (77) $$

In this case, the buckling load of the finite element system (with linear interpolation of the rotation field) with respect to its continuous “local” counterpart can be expressed as:

$$ \frac{P_{cr,n}}{P_{cr,\infty}} = 1 + \frac{\pi^2}{12n^2} + o\left(\frac{1}{n^4}\right) \quad (78) $$

as numerically shown by Seide (1975) for the finite element method (with linear interpolation of the rotation field).

The discrete equations of this finite element method are obtained from taking the stationarity condition of Rayleigh's quotient $\delta R = 0$ as defined in Eq. (71), thereby leading to (see also the discussion in Bažant and Cedolin, 2003 on Rayleigh's quotient properties):

$$ EI \frac{\theta_{i+1} - 2\theta_i + \theta_{i-1}}{a^2} + P \frac{\theta_{i+1} + 4\theta_i + \theta_{i-1}}{6} = 0 \quad (79) $$

The characteristic equation is obtained by replacing $\theta_i = A\lambda^i$ in Eq. (79) which leads, with $\beta = PL^2/EI$ to:

$$ \left(1 + \frac{\beta}{6n^2}\right)\lambda^2 + \left(-2 + \frac{4\beta}{6n^2}\right)\lambda + 1 + \frac{\beta}{6n^2} = 0 \quad (80) $$

For the pinned-pinned column, the buckling load of this finite element problem is obtained by solving Eq. (80):

$$ \lambda_{1,2} = \cos \phi \pm j \sin \phi \quad \text{with} \quad \cos \phi = \frac{1 - \frac{2\beta}{3n^2}}{1 + \frac{\beta}{6n^2}} \quad (81) $$

The buckling load is then given by $\phi = \pi/n$ and is expressed as (Seide, 1975):

$$ \begin{aligned} \beta &= 6n^2 \frac{1 - \cos \frac{\pi}{n}}{2 + \cos \frac{\pi}{n}} \Rightarrow \frac{P_{cr,n}}{P_{cr,\infty}} = \frac{3}{2 + \cos \frac{\pi}{n}} \left( \frac{\sin \frac{\pi}{2n}}{\frac{\pi}{2n}} \right)^2 \\ &= 1 + \frac{\pi^2}{12n^2} + o\left(\frac{1}{n^4}\right) \end{aligned} \quad (82) $$

The discrete equations can be directly continualized by:

$$ EI \frac{\partial_x^2 \theta}{1 - \frac{(a \partial_x)^2}{12}} + P \frac{\theta}{1 - \frac{(a \partial_x)^2}{6}} = 0 \Rightarrow EI \left[ 1 - \frac{(a \partial_x)^2}{12} \right] \partial_x^2 \theta + P \theta = 0 \quad (83) $$

and then Eq. (76) is found again by differentiation.

Returning to the axial, torsional and string vibrations problem, the finite element formulation would also lead to a gradient elasticity-type solution. For instance, when studying the axial bar problem, the linear displacement interpolation field can be used:

$$ \left\{ \begin{array}{l} u = u_{i-1}(1-\xi) + u_i\xi \\ \frac{du}{d\xi} = -u_{i-1} + u_i \end{array} \right. \quad \text{where} \quad \xi = \frac{x}{a} \quad (84) $$
The Rayleigh's quotient associated with the free vibrations problem of the axial bar is given by:

$$R = \frac{\sum_{i=1}^{n} \int_{0}^{1} \frac{ES}{a^2} \left(\frac{du}{d\xi}\right)^2 d\xi}{\sum_{i=1}^{n} \int_{0}^{1} \mu u^2 d\xi} \quad (85)$$

The discrete equations of this finite element method are obtained from taking the stationarity condition of the Rayleigh's quotient $\delta R = 0$ as defined in Eq. (85),

$$w = w_{i-1}(1 - 3\xi^2 + 2\xi^3) + w_i\xi^2(3 - 2\xi) + \theta_{i-1}a\xi(1 - \xi)^2 - \theta_i a\xi^2(1 - \xi) \quad \text{where} \quad \xi = \frac{x}{a} \quad (90)$$

Using again the definition of the discrete Rayleigh's quotient, Eq. (70) gives the Rayleigh's quotient of the cubic-based Hermitian interpolation function:

$$R = \frac{\sum_{i=1}^{n} \int_{0}^{1} \frac{EI}{a^4} \left(\frac{d^2 w}{d\xi^2}\right)^2 d\xi}{\sum_{i=1}^{n} \int_{0}^{1} \frac{1}{a^2} \left(\frac{dw}{d\xi}\right)^2 d\xi} = \frac{10EI}{a^2} \frac{\sum_{i=1}^{n} \left[ (6w_i - 6w_{i-1})^2 - (6w_i - 6w_{i-1})(6a\theta_{i-1} + 6a\theta_i) + 12a^2(\theta_{i-1}^2 + \theta_{i-1}\theta_i + \theta_i^2) \right]}{\sum_{i=1}^{n} \left[ (6w_i - 6w_{i-1})^2 - (6w_i - 6w_{i-1})(a\theta_{i-1} + a\theta_i) + a^2(4\theta_{i-1}^2 - 2\theta_{i-1}\theta_i + 4\theta_i^2) \right]} \quad (91)$$

$$ES \frac{u_{i+1} - 2u_i + u_{i-1}}{a^2} + \mu\omega^2 \frac{u_{i+1} + 4u_i + u_{i-1}}{6} = 0 \quad (86)$$

By taking the stationarity conditions of the Rayleigh's quotient $\delta R = 0$ for the two-variable field ($w_i, \theta_i$), we obtain the coupled system of finite difference equations:

$$\begin{cases} (w_{i+1} - 2w_i + w_{i-1}) \left(4 - \frac{2Pa^2}{5EI}\right) - a(\theta_{i+1} - \theta_{i-1}) \left(2 - \frac{Pa^2}{30EI}\right) = 0 \\ (w_{i+1} - w_{i-1}) \left(6 - \frac{Pa^2}{10EI}\right) - a(\theta_{i+1} + \theta_{i-1}) \left(2 + \frac{Pa^2}{30EI}\right) - a\theta_i \left(8 - \frac{4Pa^2}{15EI}\right) = 0 \end{cases} \quad (92)$$

This difference equation is also obtained by Tong et al. (1971) for the so-called consistent mass matrix. These difference equations can be continually, leading to the following enriched equation:

$$ESu'' + \mu\omega^2 u - ES I_c^2 u^{(4)} = 0 \quad \text{with} \quad I_c^2 = \frac{a^2}{12} \quad (87)$$

For fixed-fixed boundary conditions, one finds from the fourth-order differential equation Eq. (87):

$$\frac{\omega_k^2}{\omega_{k,\infty}^2} = 1 + \frac{k^2 \pi^2 l_c^2}{L^2} \quad \text{with} \quad l_c^2 = \frac{a^2}{12} \quad (88)$$

Equation (87) is an enriched dynamics equation which may come from the gradient elasticity formalism (see Tsepoura et al., 2002 or Challamel, 2013) as:

$$N = ES(u' - l_c^2 u'') \quad \text{and} \quad N' = \rho S \bar{u} \quad (89)$$

Now, we will show that a higher-order displacement field in the Finite Element Method is associated with a higher-order gradient elasticity law. Going back to the buckling problem, the Hermitian cubic functions can also be used for the interpolation function of the displacement field:

which were also obtained by Seide (1975). This system looks like the coupled system of finite difference equations considered for the discrete Timoshenko beam with some additional discrete shear effect (see Zhang et al., 2013), even if the structure of the governed equations differ between both problems. This can be easily seen from the expression of the Rayleigh's quotient of this two field problem, which is clearly different from the one of the discrete Timoshenko equations.

The finite difference system can be presented using the finite difference operators:

$$\begin{cases} \left(4 - \frac{2Pa^2}{5EI}\right) \delta_2 w - \left(4 - \frac{Pa^2}{15EI}\right) \delta_1 \theta = 0 \\ \left(12 - \frac{Pa^2}{5EI}\right) \delta_1 w - \left[12\delta_0 + \frac{Pa^4}{30EI}\delta_2 - \frac{Pa^2}{5EI}\right] \theta = 0 \end{cases} \quad (93)$$

with the finite difference operator defined in Eq. (65) and

$$\delta_2 = \frac{1}{a^2} \left[ e^{-a\bar{\alpha}_x} - 2 + e^{a\bar{\alpha}_x} \right] \quad (94)$$

The finite difference equation is then obtained as:

$$\left[ - \left( 12\delta_0 + \frac{Pa^4}{30EI}\delta_2 - \frac{Pa^2}{5EI}\right) \left( 4 - \frac{2Pa^2}{5EI} \right) \delta_2 + \left( 4 - \frac{Pa^2}{15EI} \right) \left( 12 - \frac{Pa^2}{5EI} \right) \delta_1^2 \right] w = 0 \quad (95)$$
which can be also presented as:

$$ -12\left(4 - \frac{2Pa^2}{5EI}\right) \frac{w_{i-2} + 2w_{i-1} - 6w_i + 2w_{i+1} + w_{i+2}}{6a^2} \\ - \frac{Pa^4}{30EI}\left(4 - \frac{2Pa^2}{5EI}\right) \frac{w_{i-2} - 4w_{i-1} + 6w_i - 4w_{i+1} + w_{i+2}}{a^4} \\ + \frac{Pa^2}{5EI}\left(4 - \frac{2Pa^2}{5EI}\right) \frac{w_{i-1} - 2w_i + w_{i+1}}{a^2} \\ + \left(4 - \frac{Pa^2}{15EI}\right)\left(12 - \frac{Pa^2}{5EI}\right) \frac{w_{i-2} - 2w_i + w_{i+2}}{4a^2} = 0 \qquad (96) $$

The characteristic equation is obtained by replacing $w_i = A\lambda^i$ in Eq. (96) which leads, with $\beta = PL^2/EI$ to:

$$ 4 + \frac{4\beta}{15n^2} + \frac{\beta^2}{60n^4} \left(\lambda + \frac{1}{\lambda}\right)^2 + -16 + \frac{44\beta}{15n^2} - \frac{2\beta^2}{15n^4} \left(\lambda + \frac{1}{\lambda}\right) \\ + 16 - \frac{104\beta}{15n^2} + \frac{\beta^2}{5n^4} = 0 \qquad (97) $$

that admits the following four solutions:

$$ -12\left(4 - \frac{2Pa^2}{5EI}\right) \left[1 + \frac{(a\partial_x)^2}{4} + \frac{11}{360}(a\partial_x)^4 + \frac{43}{20160}(a\partial_x)^6\right. \\ + o(a^8)\left]w'' - \frac{Pa^4}{30EI}\left(4 - \frac{2Pa^2}{5EI}\right) \left[1 + \frac{(a\partial_x)^2}{6} + \frac{1}{80}(a\partial_x)^4\right. \\ + o(a^6)\left]w^{(4)} + \frac{Pa^2}{5EI}\left(4 - \frac{2Pa^2}{5EI}\right) \left[1 + \frac{(a\partial_x)^2}{12} + \frac{1}{360}(a\partial_x)^4\right. \\ + \frac{1}{20160}(a\partial_x)^6 + o(a^8)\left]w'' + \left(4 - \frac{Pa^2}{15EI}\right)\left(12 - \frac{Pa^2}{5EI}\right) \left[1 \\ + \frac{(a\partial_x)^2}{3} + \frac{2}{45}(a\partial_x)^4 + \frac{1}{315}(a\partial_x)^6 + o(a^8)\right]w'' \\ = 0 \qquad (101) $$

which can be efficiently approximated by the following sixth-order differential equation, when collecting the terms up to the fourth-order in $a^4$:

$$ EI \frac{a^2}{6} w^{(6)} + \left(EI + \frac{3Pa^2}{20}\right) w^{(4)} + P \left(1 - \frac{Pa^2}{60EI}\right) w'' = 0 \qquad (102) $$

$$ \lambda_{1,2} = \cos \phi \pm j \sin \phi \quad \text{and} \quad \lambda_{3,4} = \frac{16 - \frac{44\beta}{15n^2} + \frac{2\beta^2}{15n^4}}{8 + \frac{8\beta}{15n^2} + \frac{\beta^2}{30n^4}} - \cos \phi \pm \sqrt{\left(\frac{16 - \frac{44\beta}{15n^2} + \frac{2\beta^2}{15n^4}}{8 + \frac{8\beta}{15n^2} + \frac{\beta^2}{30n^4}} - \cos \phi\right)^2 - 1} \quad \text{with} \\ \phi = \arccos \left[ \frac{16 - \frac{44\beta}{15n^2} + \frac{2\beta^2}{15n^4}}{4 + \frac{4\beta}{15n^2} + \frac{\beta^2}{60n^4}} - \sqrt{\left(16 - \frac{44\beta}{15n^2} + \frac{2\beta^2}{15n^4}\right)^2 - 4\left(4 + \frac{4\beta}{15n^2} + \frac{\beta^2}{60n^4}\right)\left(16 - \frac{104\beta}{15n^2} + \frac{\beta^2}{5n^4}\right)} \right] \qquad (98) $$

The buckling load is then given by $\phi = \pi/n$ and is expressed as:

$$ \frac{\cos \frac{\pi}{n} - 3}{60n^4} \beta^2 + \frac{4 \cos \frac{\pi}{n} + 26}{15n^2} \beta + 4 \left( \cos \frac{\pi}{n} - 1 \right) = 0 \qquad (99) $$

For the pinned-pinned column, the buckling load of this finite element problem is then obtained from:

$$ \frac{P_{cr,n}}{P_{cr,\infty}} = \frac{13 + 2 \cos \frac{\pi}{n} - \sqrt{124 + 112 \cos \frac{\pi}{n} - 11 \cos^2 \frac{\pi}{n}}}{3 - \cos \frac{\pi}{n}} \left(\frac{2n}{\pi}\right)^2 \\ = 1 + \frac{\pi^4}{720n^4} + o\left(\frac{1}{n^6}\right) \qquad (100) $$

which is also the value reported by Seide (1975). Now, by using a continualization procedure, an asymptotic expansion of each difference operator in Eq. (96) gives:

This differential equation can be factorized by:

$$ \left[ \frac{(a\partial_x)^2}{6} + 1 - \frac{Pa^2}{60EI} \right] \left[ EI\partial_x^4 + P\partial_x^2 \right] w = 0 \qquad (103) $$

Equation (103) shows the product of differential operators, with the so-called local differential operator. As a consequence, the local solution $P_{cr,n} = P_{cr,\infty}$ is a solution of this system, which means that a higher-order asymptotic expansion is needed for developing the finite difference operators in Eq. (101) with higher-order terms. Now collecting the terms up to the sixth-order in $a^6$ leads to the enriched equation:

$$ EI \frac{a^4}{80} w^{(8)} + \left(EI \frac{a^2}{6} + \frac{Pa^4}{72}\right) w^{(6)} + \left(EI + \frac{3Pa^2}{20} + \frac{P^2 a^4}{360EI}\right) w^{(4)} \\ + P \left(1 - \frac{Pa^2}{60EI}\right) w'' = 0 \qquad (104) $$

By introducing the sinusoidal shape function $\tilde{w} = \bar{w} \sin(\pi x/L)$ in this eight-order differential equation, we obtain the following buckling load
$$ \frac{P_{cr,n}}{P_{cr,\infty}} = \frac{1 - \frac{3\pi^2}{20n^2} + \frac{\pi^4}{72n^4} - \sqrt{\left(1 - \frac{3\pi^2}{20n^2} + \frac{\pi^4}{72n^4}\right)^2 - 4\left(1 - \frac{\pi^2}{6n^2} + \frac{\pi^4}{80n^4}\right)\left(\frac{\pi^2}{60n^2} + \frac{\pi^4}{360n^4}\right)}}{\frac{\pi^2}{30n^2} + \frac{\pi^4}{180n^4}} \quad (105) $$

An asymptotic expansion shows that:

$$ \frac{P_{cr,n}}{P_{cr,\infty}} = 1 + \frac{\pi^4}{720n^4} + o\left(\frac{1}{n^6}\right) \quad (106) $$

which is the value reported by Seide (1975) based on the cubic-based finite element method. This coefficient 1/720 in the asymptotic expansion of the eigenvalue problem is also found by Tong et al. (1971) for the FEM solution of the eigenfrequencies of a pinned-pinned column. The differential equation Eq. (104) can also be factorized as

$$ EI \frac{a^4}{80} \partial_x^8 + \left( EI \frac{a^2}{6} + \frac{Pa^4}{72} \right) \partial_x^6 + \left( EI + \frac{3Pa^2}{20} + \frac{P^2a^4}{360EI} \right) \partial_x^4 + P \left( 1 - \frac{Pa^2}{60EI} \right) \partial_x^2 = \left[ 1 + \frac{(a\partial_x)^2}{6} - \frac{Pa^2}{60EI} + \frac{(a\partial_x)^4}{90} + \frac{Pa^2}{360EI} (a\partial_x)^2 \right] \times \left[ \frac{EIa^4}{720} \partial_x^8 + EI\partial_x^4 + P\partial_x^2 \right] + o(a^6) \quad (107) $$

which means that the cubic-based finite element model can be equivalently reduced to the eight-order differential equation:

$$ EI \frac{a^4}{720} w^{(8)} + EIw^{(4)} + Pw'' = 0 \quad (108) $$

Walz et al. (1968) also obtained a corrected eight-order differential equation for the continualized bending problem which was investigated by the Hermitian-based Finite Element model, with the correct coefficient 1/720 but with a different sign. For the finite element model considered herein, the associated Rayleigh's quotient can then be expressed by:

$$ R = \frac{\int_{0}^{L} EI (w''^2 + l_c^4 w^{(4)2}) dx}{\int_{0}^{L} w'^2 dx} \geq \frac{\int_{0}^{L} EIw''^2 dx}{\int_{0}^{L} w''^2 dx} \quad \text{with} \quad l_c^4 = \frac{a^4}{720} \quad (109) $$

leading to the gradient elasticity solution Eq. (109), associated with the gradient elasticity constitutive law:

$$ M = EI(w'' + l_c^4 w^{(6)}) \quad \text{with} \quad l_c^4 = \frac{a^4}{720} \quad (110) $$

Consider again a pinned-pinned column and introducing a sinusoidal shape function $\tilde{w} = \bar{w} \sin(\pi x/L)$ as an admissible function into the Rayleigh's quotient leads to the exact buckling value of the pinned-pinned case, also given in Eq. (106). Equation (106) shows that the Finite Element column model gives an upper bound of the "local" problem asymptotically found for n tending towards infinite.

The length scale calibration of the gradient elasticity model equivalent to the cubic hermitian finite element formulation in buckling is then equal to:

$$ l_{c,static}^4 = \frac{a^4}{720} = l_{c,dynamic}^4 \quad (111) $$

As shown by Challamel et al. (2015b), one also recognizes the equivalent length scale of the finite element formulation applied to the vibrations analysis, which shows that the gradient elasticity model associated with the finite element formulation, is independent of the static versus vibration analysis (a property which is different from the one associated with the finite difference schemes). Furthermore, using the generalized Rayleigh ratio, a new general formulae can be applied as well for some other boundary conditions (including the clamped-free and the clamped-clamped boundary conditions):

$$ \frac{P_{cr,n}}{P_{cr,\infty}} = 1 + \left(\frac{P_{cr,\infty}}{EI}\right)^2 l_c^4 \quad \text{with} \quad l_c^4 = \frac{a^4}{720} \quad (112) $$

which is again consistent with the asymptotic results of Seide (1975).

In other words, the central finite difference method can be efficiently approximated by a nonlocal Eringen's model based on the second derivative of the generalized stress variable (the bending moment for the considered beam model) as a corrected term, whereas the higher-order finite difference method can be captured by a higher-order nonlocal model based on the fourth-order derivative of the generalized stress variable. Both finite difference based models can be classified as stress gradient models. In both cases, a lower bound of the buckling load is obtained as shown by the associated nonlocal Rayleigh's quotient (see Fig. 1 or Fig. 2). On the other hand, the finite element method with a linear interpolation field of the rotation

<figure><img src="image_2.png" /><figcaption>Fig. 1. Finite difference method (or some equivalent nonlocal elastic systems) leads to a lower bound of the buckling load, whereas finite element method (or some equivalent gradient elasticity systems) leads to an upper bound.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Comparison of finite difference and finite element approaches for buckling of simply supported column; Efficiency of enriched continua with respect to exact numerical results.</figcaption></figure>

can be efficiently approximated by a gradient elasticity model based on the second derivative of the generalized strain variable (the curvature variable for the considered beam model) as a corrected term, whereas the cubic-based Finite Element method can be captured by a higher-order gradient elasticity model based on the fourth-order derivative of the generalized strain variable (see Figs. 1 and 2). Both finite element based models can be classified as strain gradient models. As it is also shown from the equivalent formulation of the Rayleigh's quotient, finite element methods give an upper bound of the buckling load of the local model obtained for an infinite number of elements.

## 5. Conclusions

Numerical approaches often convert the mathematical equations of continuous systems into a finite number of variables associated with some discrete equivalent systems in order to reduce the complexity of the mathematical problem for solution. It is shown herein, from the vibrations and the buckling analyses of one-dimensional systems, that the obtained discrete equations behave as a kind of enriched nonlocal continua with respect to the local reference continua. The analogy between the finite numerical approaches and some equivalent enriched continuum is shown by using a continualization procedure (or equivalently by using the method of differential approximation). The finite difference method is investigated using an integral-based nonlocal method, whereas the finite element method appears to be equivalent to a gradient elasticity method. As a consequence, the integral-based nonlocal method gives a lower bound solution to the exact eigenvalue multiplier whereas the gradient elasticity method furnishes an upper bound solution.

The convergence rate of the finite numerical approaches may be quantified through the use of Rayleigh's quotient for the continualization approach. The rate of convergence is strongly dependent on the order of the finite discrete scheme: higher-order finite schemes lead to higher-order enriched constitutive law with a higher convergence rate. These analogies between enriched continuum and finite numerical schemes give a new attractive framework for potential applications of enriched continua in computational mechanics. Surprisingly, these higher-order numerical schemes give the opportunity to build some new higher-order nonlocal or gradient constitutive laws.

# Appendix A. Rayleigh's quotient for finite difference system

The potential energy depends on the difference:

$$ \frac{w_{i+1} - w_i}{a} = \frac{[\exp(a\partial_x) - 1]}{a} w_i \quad \text{with} \quad w_i = \frac{a^2 \partial_x^2}{4 \sinh^2 \frac{a\partial_x}{2}} w \quad (A1) $$

A new differential operator can be introduced as:

$$ \begin{aligned} \frac{w_{i+1} - w_i}{a} &= Q_1(a\partial_x)\partial_x w \quad \text{with} \quad Q_1(a\partial_x) = \frac{[\exp(a\partial_x) - 1]a\partial_x}{4 \sinh^2 \frac{a\partial_x}{2}} \\ &= \frac{\exp\left(\frac{a\partial_x}{2}\right) a\partial_x}{2 \sinh\left(\frac{a\partial_x}{2}\right)} \end{aligned} \quad (A2) $$

The last differential operator can be expanded as:

$$ \begin{aligned} Q_1(a\partial_x) &= \frac{\exp\left(\frac{a\partial_x}{2}\right) a\partial_x}{2 \sinh\left(\frac{a\partial_x}{2}\right)} = \frac{-a\partial_x}{\exp(-a\partial_x) - 1} = \sum_{k=0}^{\infty} B_k (-1)^k (a\partial_x)^k \\ &= 1 + \frac{a\partial_x}{2} + \frac{1}{12}(a\partial_x)^2 + \dots \end{aligned} \quad (A3) $$

where $B_k$ being the Bernoulli numbers defined by the expansion of $\xi/\exp(\xi) - 1 = \sum_{k=0}^{\infty} B_k/k!\xi^k$. The potential energy can then be obtained from:

$$ \begin{aligned} \int_0^L \left( \frac{w_{i+1} - w_i}{a} \right)^2 dx &= \int_0^L Q_1(a\partial_x)(\partial_x w)Q_1(a\partial_x)(\partial_x w)dx \\ &= \int_0^L (\partial_x w)Q_1^*(a\partial_x)Q_1(a\partial_x)(\partial_x w)dx \end{aligned} \quad (A4) $$

where $Q_1^*(a\partial_x) = Q_1(-a\partial_x)$. For infinite medium, $Q_1^*(a\partial_x) = Q_1(-a\partial_x)$ is the adjoint operator whereas for finite systems as the ones considered in this paper, the adjointness of $Q_1^*(a\partial_x) = Q_1(-a\partial_x)$ depends on the boundary conditions of the system. It is possible to expand this differential operator as:

$$ Q_1^*(a\partial_x)Q_1(a\partial_x) = \frac{a^2 \partial_x^2}{4 \sinh^2\left(\frac{a\partial_x}{2}\right)} = \left[ 1 - \frac{(a\partial_x)^2}{12} + o(a^4) \right] \quad (A5) $$

The potential energy term can then be rewritten as:

$$ \int_0^L \left( \frac{w_{i+1} - w_i}{a} \right)^2 dx = \int_0^L w' \left( w' - \frac{a^2}{12} w'' + o(a^4) \right) dx \quad (A6) $$

An integration by part finally gives Eq. (A7).

$$ \int_0^L \left( \frac{w_{i+1} - w_i}{a} \right)^2 dx = \int_0^L w'^2 + \frac{a^2}{12} w''^2 + o(a^4) dx - \frac{a^2}{12} \left[ w'w'' \right]_0^L \quad (A7) $$
## Appendix B. Rayleigh's quotient for finite element system

The discrete potential energy is written in the following form, with the new continuous variable:

$$ \frac{\theta_i^2 + \theta_i \theta_{i-1} + \theta_{i-1}^2}{3} = (\partial_x w)^2 \qquad (B1) $$

The discrete left hand-side term can be also expressed by:

$$ \begin{aligned} \frac{\theta_i^2 + \theta_i \theta_{i-1} + \theta_{i-1}^2}{3} &= \frac{1}{3} (\theta_i \quad \theta_{i-1}) \begin{pmatrix} 1 & \frac{1}{2} \\ \frac{1}{2} & 1 \end{pmatrix} \begin{pmatrix} \theta_i \\ \theta_{i-1} \end{pmatrix} \\ &= \frac{1}{3} (1 \quad e^{-a\partial_x}) \theta_i \begin{pmatrix} 1 & \frac{1}{2} \\ \frac{1}{2} & 1 \end{pmatrix} \begin{pmatrix} 1 \\ e^{-a\partial_x} \end{pmatrix} \theta_i \end{aligned} \qquad (B2) $$

and then by integration, the following expansion can be performed:

$$ \begin{aligned} \int_0^L \frac{\theta_i^2 + \theta_i \theta_{i-1} + \theta_{i-1}^2}{3} dx &= \int_0^L \frac{1}{3} \theta_i [2 + \cosh(a\partial_x)] \theta_i dx \\ &= \int_0^L \theta_i \left[ 1 + \frac{(a\partial_x)^2}{6} + o(a^4) \right] \theta_i dx \end{aligned} \qquad (B3) $$

Using the definition Eq. (B1), the continuous rotation can then be expressed by:

$$ \partial_x w = \left[ 1 + \frac{(a \partial_x)^2}{12} + \dots \right] \theta_i \qquad (B4) $$

This condition can be inverted, i.e.

$$ \theta_i = \left[ 1 - \frac{(a \partial_x)^2}{12} + \dots \right] \partial_x w \qquad (B5) $$

A new differential operator can be introduced as:

$$ \begin{aligned} \theta_i - \theta_{i-1} &= Q_2(a\partial_x)\partial_x w \quad \text{with} \quad Q_2(a\partial_x) \\ &= \left[ 1 - e^{-a\partial_x} \right] \left[ 1 - \frac{(a\partial_x)^2}{12} + \dots \right] \end{aligned} \qquad (B6) $$

The differential operator can be expanded as:

$$ Q_2(a\partial_x) = 1 - \frac{a\partial_x}{2} + \frac{(a\partial_x)^2}{12} + \dots a\partial_x \qquad (B7) $$

The strain energy can then be expressed with this differential operator $Q_2$ or with a new differential operator $Q_3$:

$$ \begin{aligned} \int_0^L (\theta_i - \theta_{i-1})^2 dx &= \int_0^L Q_2(a\partial_x)(\partial_x w)Q_2(a\partial_x)(\partial_x w)dx \\ &= a^2 \int_0^L Q_3(a\partial_x)(\partial_x^2 w)Q_3(a\partial_x)(\partial_x^2 w)dx \quad \text{with} \\ & \quad \times Q_3(a\partial_x) = 1 - \frac{a\partial_x}{2} + \frac{(a\partial_x)^2}{12} + \dots \end{aligned} \qquad (B8) $$

It is then possible to expand the strain energy in the following form:

$$ \begin{aligned} \int_0^L (\theta_i - \theta_{i-1})^2 dx &= a^2 \int_0^L Q_3(a\partial_x)(\partial_x^2 w)Q_3(a\partial_x)(\partial_x^2 w)dx \\ &= a^2 \int_0^L (\partial_x^2 w)Q_3^*(a\partial_x)Q_3(a\partial_x)(\partial_x^2 w)dx \end{aligned} \qquad (B9) $$

where $Q_3^*(a\partial_x) = Q_3(-a\partial_x)$. It is easy to notice that:

$$ Q_3^*(a\partial_x)Q_3(a\partial_x) = Q_1^*(a\partial_x)Q_1(a\partial_x) = 1 - \frac{(a\partial_x)^2}{12} + o(a^4) \qquad (B10) $$

The strain energy can be expanded with this new variable as:

$$ \int_0^L \left( \frac{\theta_i - \theta_{i-1}}{a} \right)^2 dx = \int_0^L \partial_x^2 w \left[ 1 - \frac{(a\partial_x)^2}{12} \right] \partial_x^2 w dx \qquad (B11) $$

Finally, the strain energy can be presented in its gradient elasticity format, by integration by parts:

$$ \int_0^L \left( \frac{\theta_i - \theta_{i-1}}{a} \right)^2 dx = \int_0^L w''^2 + \frac{a^2}{12} w'''^2 + o(a^4) dx - \frac{a^2}{12} [w'' w']_0^L \qquad (B12) $$

## References

Andrianov, I.V., Awrejcewicz, J., Weichert, D., 2010. Improved continuous models for discrete media. Math. Probl. Eng. 1–35, 986242.

Arda, M., Aydogdu, M., 2014. Torsional statics and dynamics of nanotubes embedded in an elastic medium. Compos. Struct. 114, 80–91.

Aydogdu, M., 2009. Axial vibration of the nanorods with the nonlocal continuum rod model. Phys. E 41, 861–864.

Bažant, Z.P., Cedolin, L., 2003. Stability of Structures – Elastic, Inelastic, Fracture, and Damage Theories. Dover Publications, Inc., New-York.

Blevins, R.D., 2001. Formulas for Natural Frequency and Mode Shape. Krieger Publishing Company.

Born, M., von Kármán, T., 1912. On fluctuations in spatial grids. Phys. Z. 13, 297–309.

Challamel, N., Wang, C.M., 2008. The small length scale effect for a non-local cantilever beam: a paradox solved. Nanotechnology 19, 345703.

Challamel, N., Rakotomanana, L., Le Marrec, L., 2009. A dispersive wave equation using non-local elasticity. C. R. Méc. 337, 591–595.

Challamel, N., 2013. Variational formulation of gradient or/and nonlocal higher-order shear elasticity beams. Compos. Struct. 105, 351–368.

Challamel, N., Lerbet, J., Wang, C.M., Zhang, Z., 2014a. Analytical length scale calibration of nonlocal continuum from a microstructured buckling model. Z. Angew. Math. Mech. 94 (5), 402–413.

Challamel, N., Wang, C.M., Elishakoff, I., 2014b. Discrete systems behave as nonlocal structural elements: bending, buckling and vibration analysis. Eur. J. Mech. A/ Solids 44, 125–135.

Challamel, N., Lerbet, J., Wang, C.M., 2014c. On buckling of granular columns with shear interaction: discrete versus nonlocal approaches. J. Appl. Phys. 115, 234902.

Challamel, N., Zhang, Z., Wang, C.M., Reddy, J.N., Wang, Q., Michelitsch, T., Collet, B., 2014d. On non-conservativeness of Eringen's nonlocal elasticity in beam mechanics: correction from a discrete-based approach. Arch. Appl. Mech. 84 (9), 1275–1292.

Challamel, N., Camotim, D., Wang, C.M., Zhang, Z., 2015a. On lateral-torsional buckling of discrete elastic systems: a non-local approach. Eur. J. Mech. A/ Solids 49, 106–113.

Challamel, N., Picandet, V., Elishakoff, I., Wang, C.M., Collet, B., Michelitsch, T., 2015b. On nonlocal computation of eigenfrequencies of beams using finite difference and finite element methods. Int. J. Struct. Stab. Dyn. 15 (7) (in press), Micromech. ASCE 5, A4014004.

Challamel, N., Zhang, Z., Wang, C.M., 2015c. Nonlocal equivalent continua for buckling and vibration analyses of microstructured beams. J. Nanomech. Micromech. ASCE 5, A4014004.

Charlotte, M., Truskinovsky, L., 2012. Lattice dynamics from a continuum viewpoint. J. Mech. Phys. Solids 60, 1508–1544.

Collins, M.A., 1981. A quasicontinuum approximation for solitons in an atomic chain. Chem. Phys. Lett. 77 (2), 342–347.
Cyrus, N.J., Fulton, R.E., 1968. Accuracy Study of Finite Difference Methods. NASA Technical Note, NASA-TN D-4372. NASA Langley Research Center, pp. 1-29.

Duan, W., Challamel, N., Wang, C.M., Ding, Z., 2013. Development of analytical vibration solutions for microstructured beam model to calibrate length scale coefficient in nonlocal Timoshenko beams. J. Appl. Phys. 114, 1-11, 104312.

Elishakoff, I., 1998. Lecture Notes in Stability of Structures. Florida Atlantic University, Boca Raton, USA.

Elishakoff, I., Santoro, R., 2005. Error in the finite difference based probabilistic dynamic analysis: analytical evaluation. J. Sound Vib. 281, 1195-1206.

Elishakoff, I., Pentaras, D., Dujat, K., Versaci, C., Muscolino, G., Storch, J., Bucas, S., Challamel, N., Natsuki, T., Zhang, Y.Y., Wang, C.M., Ghyselinck, G., 2012. Carbon Nanotubes and Nanosensors: Vibrations, Buckling and Ballistic Impact. Wiley - ISTE, London.

Eringen, A.C., 1983. On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves. J. Appl. Phys. 54, 4703-4710.

Gantmacher, F., 1970. Lectures in Analytical Mechanics. Mir Publishers, Moscow.

Gawain, T.H., Ball, R.E., 1978. Improved finite difference formulas for boundary value problems. Int. J. Num. Meth. Eng. 12 (7), 1151-1160.

Goldberg, S., 1958. Introduction to Difference Equations with Illustrative Examples from Economics, Psychology and Sociology. Dover Publications, New-York, p. 184.

Greenwood, D.T., 1961. The Use of Higher-order Finite Differences in Beam Vibrations. NASA TND-964.

Guttman, L., 1950. The principal components of scale analysis (Chapter 9). In: Stouffer, S.A. (Ed.), Measurement and Prediction, Studies in Social Psychology in World War II, vol. 4. Princeton University Press, Princeton.

Hencky, H., 1920a. Habilitationschrift, D. Sc. Thesis. Technical University of Darmstadt, Germany (in German).

Hencky, H., 1920b. Über die angenaherte Lösung von Stabilitätsproblemen im Raummittels der elastischen Gelenkkette. In: Der Eisenbau, vol. 11, pp. 437-452 (in German).

Ilanko, S., Monterrubio, L.E., 2015. The Rayleigh-Ritz Method for Structural Analysis. ISTE-Wiley.

Islam, Z.M., Jia, P., Lim, C.W., 2014. Torsional wave propagation and vibration of circular nanostructures based on nonlocal elasticity theory. Int. J. Appl. Mech. 6 (2), 1450011.

Kevrekidis, P.G., Kevrekidis, I.G., Bishop, A.R., Titi, E.S., 2002. Continuum approach to discreteness. Phys. Rev. E 65, 1-13, 046613.

Kruskal, M.D., Zabusky, N.J., 1964. Stroboscopic perturbation for treating a class of nonlinear wave equations. J. Math. Phys. 5, 231-244.

Kunin, I.A., 1982. Elastic Media with Microstructure. I - One-dimensional Models. Springer, Berlin.

Lagrange, J.L., 1788. Mécanique Analytique, Paris, third ed. In: Mallet-Bachelier, Gendre et successeur de Bachelier, Imprimeur-Libraire du bureau des longitudes, de l'école Polytechnique de l'école centrale des arts et manufactures, Paris, 1853, p. 367.

Leckie, F.A., Lindberg, G.M., 1963. The effect of lumped parameters on beam frequencies. Aeronaut. Q. 14, 224-240.

Leissa, A.W., Qatu, M.S., 2011. Vibrations of Continuous Systems. McGraw-Hill, New York.

Livesley, R.K., 1955. The equivalence of continuous and discrete mass distributions in certain vibration problems. Q. J. Mech. Appl. Math. 8 (3), 353-360.

Mindlin, R.D., 1964. Micro-structure in linear elasticity. Arch. Rat. Mech. Anal. 16, 51-78.

Myklestad, N.O., 1944. A new method of calculating natural modes of uncoupled bending vibration of airplane wings and other types of beams. J. Aeronaut. Sci. 11.

Narendar, S., 2011. Nonlocal torsional vibration of nanorods. J. Nanosci. Nanoeng. Appl. 1 (2), 36-51.

Papargyri-Beskou, S., Tsepoura, K.G., Polyzos, D., Beskos, D.E., 2003. Bending and stability analysis of gradient elastic beams. Int. J. Solids Struct. 40, 385-400.

Rayleigh, Lord, 1894. The Theory of Sound, second ed., p. 172 London.

Rosenau, P., 1986. Dynamics of nonlinear mass-spring chains near the continuum limit. Phys. Lett. A 118 (5), 222-227.

Rosenau, P., 1987. Dynamics of dense lattices. Phys. Rev. B 36 (11), 5868-5876.

Rosenau, P., 2003. Hamiltonian dynamics of dense chains and lattices: or how to correct the continuum. Phys. Lett. A 311, 39-52.

Salvadori, M.G., 1951. Numerical computation of buckling loads by finite differences. Trans. ASCE 116, 590-624 (590-636 with the discussion).

Santoro, R., Elishakoff, I., 2006. Accuracy of the finite difference method in stochastic setting. J. Sound Vib. 291, 275-284.

Seide, P., 1975. Accuracy of some numerical methods for column buckling. J. Eng. Mech. 101 (5), 549-560.

Shokin, Yu.I., 1983. The Method of Differential Approximation. Springer.

Silverman, I.K., 1951. Discussion on the paper of "Salvadori M.G., numerical computation of buckling loads by finite differences. Trans. ASCE 116, 590-636. Trans. ASCE, 116, 625-626, 1951.

Strang, G., Fix, G.J., 1973. An Analysis of the Finite Element Method. Prentice-Hall, Englewood Cliffs, NJ.

Subrahmanyam, K.B., Kaza, K.R.V., 1983. An Improved Finite-difference Analysis of Uncoupled Vibrations of Tapered Cantilever Beams. In: NASA Technical Memorandum 83495, September, pp. 1-46.

Subrahmanyam, K.B., Leissa, A.W., 1985. An improved finite difference analysis of uncoupled vibrations of cantilevered beams. J. Sound Vib. 98 (1), 1-11.

Temple, G., Bickley, W.G., 1956. Rayleigh's Principle and its Applications to Engineering. Dover Publications.

Thomson, W.T., 1972. Theory of Vibration with Applications. Prentice-Hall, Englewood Cliffs, NJ.

Tong, P., Pian, T.H.H., Bucciarelli, L.L., 1971. Mode shapes and frequencies by finite element method using consistent and limped masses. Comput. Struct. 1, 623-638.

Tsepoura, K.G., Papargyri-Beskou, S., Polyzos, D., Beskos, D.E., 2002. Static and dynamic analysis of a gradient elastic bar in tension. Arch. Appl. Mech. 72, 483-497.

Walz, J.E., Fulton, R.E., Cyrus, N.J., 1968. Accuracy and Convergence of Finite Element Approximations. Report AFFDL-TR-68-150. Air Force Flight Dynamics Laboratory, pp. 996-1027. see also AD0703685, Proceedings of the Conference on Matrix Methods in Structural Mechanics (2nd) Held at Wright-Patterson Air Force Base, Ohio, on 15-17 October 1968.

Wang, C.M., Zhang, Z., Challamel, N., Duan, W.H., 2013. Calibration of Eringen's small length scale coefficient for initially stressed vibrating nonlocal Euler beams based on microstructured beam model. J. Phys. D. Appl. Phys. 46, 345501.

Wang, C.T., 1951. Discussion on the paper of "Salvadori M.G., numerical computation of buckling loads by finite differences". Trans. ASCE 116, 629-631.

Wang, C.T., 1953. Applied Elasticity. McGraw-Hill, New-York.

Wang, C.Y., Wang, C.M., 2013. Structural Vibration – Exact Solutions for Strings, Membranes, Beams and Plates. CRC Press, Boca Raton,FL.

Wattis, J.A.D., 2000. Quasi-continuum approximations to lattice equations arising from the discrete non-linear telegraph equation. J. Phys. A: Math. General 33, 5925-5944.

Zhang, Z., Challamel, N., Wang, C.M., 2013. Eringen's small length scale coefficient for buckling of nonlocal Timoshenko beam based on a microstructured beam model. J. Appl. Phys. 114, 1-6, 114902.

Zhang, Z., Wang, C.M., Challamel, N., Elishakoff, I., 2014a. Obtaining Eringen's length scale coefficient for vibrating nonlocal beams via continualization method. J. Sound. Vib. 333, 4977-4990.

Zhang, Z., Wang, C.M., Challamel, N., 2014b. Eringen's length scale coefficient for buckling of nonlocal rectangular plates from microstructured beam-grid model. Int. J. Solids Struct. 51, 25-26, 4307-4315.

Zhang, Z., Wang, C.M., Challamel, N., 2015. Eringen's length scale coefficient for vibration and buckling of nonlocal rectangular plates with simply supported edges. J. Eng. Mech. 141 (2), 04014117.

Zhukov, A.I.O., 1957. Convergence of the solution of a difference equation to the solution of the differential equation. Dokl. Akad. Nauk. SSSR Ser. Mat. Fiz. 117 (2), 174-176 (in Russian).

Zienkiewicz, O.C., Cheung, Y.K., 1967. The Finite Element Method in Structural and Continuum Mechanics: Numerical Solution of Problems in Structural and Continuum Mechanics. Mc Graw-Hill, London.