# Geometrically Nonlinear Response of a Fractional-Order Nonlocal Model of
Elasticity

Sai Sidhardh  School of Mechanical Engineering, Ray W. Herrick Laboratories,
Purdue University, West Lafayette, IN 47907  Sansit Patnaik  School of
Mechanical Engineering, Ray W. Herrick Laboratories, Purdue University, West
Lafayette, IN 47907  Fabio Semperlotti† School of Mechanical Engineering, Ray
W. Herrick Laboratories, Purdue University, West Lafayette, IN 47907

###### Abstract

This study presents the analytical and finite element formulation of a
geometrically nonlinear and fractional-order nonlocal model of an Euler-
Bernoulli beam. The finite nonlocal strains in the Euler-Bernoulli beam are
obtained from a frame-invariant and dimensionally consistent fractional-order
(nonlocal) continuum formulation. The finite fractional strain theory provides
a positive definite formulation that results in a mathematically well-posed
formulation which is consistent across loading and boundary conditions. The
governing equations and the corresponding boundary conditions of the
geometrically nonlinear and nonlocal Euler-Bernoulli beam are obtained using
variational principles. Further, a nonlinear finite element model for the
fractional-order system is developed in order to achieve the numerical
solution of the integro-differential nonlinear governing equations. Following
a thorough validation with benchmark problems, the fractional finite element
model (f-FEM) is used to study the geometrically nonlinear response of a
nonlocal beam subject to various loading and boundary conditions. Although
presented in the context of a 1D beam, this nonlinear f-FEM formulation can be
extended to higher dimensional fractional-order boundary value problems.

† All correspondence should be addressed to<EMAIL_ADDRESS>

##  1 Introduction

Several theoretical and experimental studies have shown that nonlinear and
scale-dependent effects are prominent in the response of several structures
ranging from sandwich composites, to layered and porous media, to random and
fractal media, to media with damage and cracks, to biological materials like
tissues and bones [1, 2, 3, 4]. In all these classes of problems, the
coexistence of different material spatial scales renders the response fully
nonlocal [5, 6]. In many applications involving slender structures, large and
rapidly varying loads typically induce nonlinear response. These nonlinear
size-dependent effects are particularly prominent in microstructures and
nanostructures that have far-reaching applications in medical devices, atomic
devices, micro/nano-electromechanical devices and sensors [7, 8, 9]. These
devices are primarily made from a combination of slender structures such as
beams, plates and shells. Clearly, accurate modeling of the nonlocal and
nonlinear effects in the response of these structures is paramount in many
engineering as well as biomedical applications. However, it appears that there
is a limited amount of studies focusing on the modeling of geometrically
nonlinear and nonlocal slender structures. In the following, we briefly review
the main characteristics of these studies as well as discuss key limitations.

Nonlocal continuum theories use different strategies to enrich the classical
governing equations, that describe the response of the medium at a point, with
information on the behavior of distant points contained in a prescribed area
of influence. Depending on the field in which this nonlocal description of the
medium is applied, the same concept can be referred to with different
terminologies such as action at a distance, long range interactions, or scale
effects. The key concept at the basis of nonlocal continuum theories is the
horizon of influence or horizon of nonlocality which describes areas of the
medium where distant points can still influence one another via long range
cohesive forces [5, 6]. Several theories based on integral methods, gradient
methods, and very recently, the peridynamic approach have been developed to
capture these long range interactions and analyze their effect on the response
of either media or structures. Gradient elasticity theories [10, 11, 12, 13,
14] account for the nonlocal behavior by introducing strain gradient dependent
terms in the stress-strain constitutive law. Integral methods [15, 16, 17]
define constitutive laws in the form of convolution integrals between the
strain and the spatially dependent elastic properties over the horizon of
nonlocality.

As emphasized earlier, several applications involving nonlocal slender
structures also experience geometrical nonlinearities. Nonlocal effects on the
nonlinear transverse response of shear deformable beams were modeled in [18,
19] using von Kármán nonlinearity and Eringen’s nonlocal elasticity. Reddy et
al. [20] modeled a nonlocal plate starting from the differential constitutive
relations of Eringen and the von Kármán nonlinear strains. The differential
model of Eringen’s nonlocal elasticity was also used in [21, 22] to analyse
nonlocal and nonlinear beams. Similar models for beams and plates have also
been developed in [23, 24]. Several studies were also carried out to simulate
the multiphysics response of beams such as, for example, the nonlinear static
and dynamic response of a capacitive electrostatic nanoactuator [25] or the
effect of nonlocal elasticity over the electroelastic response of
flexoelectric nanobeams [17].

Although the above nonlinear studies based on classical approaches have been
able to capture key aspects of nonlocal elasticity, they are subject to
important drawbacks. As shown in [26, 27], the integral nonlocal methods do
not allow defining a self-adjoint quadratic potential energy which results in
paradoxical predictions such as hardening and absence of nonlocal effects for
certain combinations of boundary conditions [27, 28]. Note that hardening
effects per se are not inaccurate or paradoxical from a physical perspective
and have been observed in several micro- and nanostructures [29]. However,
hardening effects are not expected in classical integral approaches. This
specific inaccuracy is the result of the ill-posed nature of the integral
governing equations [29]. Also, given the complex nature of the nonlocal and
nonlinear differential equations, research on the corresponding numerical
methods is limited so far. It appears that only quadrature based methods [25,
30] have been developed for the solution of nonlinear and nonlocal boundary
value problems (BVPs).

In recent years, fractional calculus has emerged as a powerful mathematical
tool to model a variety of nonlocal and multiscale phenomena. Fractional
derivatives, which are a differ-integral class of operators, are intrinsically
multiscale and provide a natural way to account for nonlocal effects. In fact,
the order of the fractional derivative dictates the shape of the influence
function (kernel of the fractional derivative) while its domain of integration
defines the horizon of influence, that is the distance beyond which
information is no longer accounted for in the derivative. It follows that time
fractional operators enable memory effects (i.e. the response of a system is a
function of its past history) while space fractional operators account for
nonlocal and scale effects. The unique set of properties of fractional
operators have determined a surge of interest in the exploration of possible
applications. Among the areas that have rapidly developed, we find the
modeling of viscoelastic materials [31, 32, 33], transport processes in
complex media [34, 35, 36], and model-order reduction of lumped parameter
systems [37, 38]. Given the multiscale nature of fractional operators,
fractional calculus has also found several applications in nonlocal
elasticity. Space-fractional derivatives have been used to formulate nonlocal
constitutive laws [39, 40, 41, 36] as well as to account for microscopic
interaction forces [42, 43, 44]. Very recently fractional-order theories have
been extended to model and analyze the static response, buckling
characteristics, as well as the dynamic response of nonlocal beams [7, 45, 9,
46, 47, 48]. However, all these studies focused on the linear geometric
analysis of nonlocal structures.

Previous works conducted on the development of nonlocal continuum theories
based on fractional calculus have highlighted that the differ-integral nature
of the fractional operators allows them to combine the strengths of both
gradient and integral based methods while at the same time addressing a few
important shortcomings of these integer-order formulations which have been
highlighted earlier [39, 40, 41, 47, 49]. We highlight here that, while the
fractional continuum formulation developed in this study builds upon [40, 41],
in our study we have extended and generalized the formulation in order to
account for asymmetric horizon conditions that typically occur near
geometric/material boundaries or discontinuities. Further, note that while
[41] focused primarily on anisotropic nonlocality, our model is for isotropic
materials with possible asymmetric horizon conditions. In this study, we
address the shortcomings of the classical nonlocal formulations by developing
a consistent nonlinear and fractional-order nonlocal beam model. Further, we
develop a nonlinear fractional-order finite element method (f-FEM) to obtain
the numerical solution of nonlinear and nonlocal BVPs.

The overall goal of this study is two-fold. First, we derive the governing
equations for the geometrically nonlinear and nonlocal beam in strong form
using variational principles. We build upon the fractional-order nonlocal
continuum model proposed in [47] to develop a fractional-order constitutive
relation for a geometrically nonlinear Euler-Bernoulli beam. The model
proposed in [47] was shown to be self-adjoint, positive definite and well-
posed, thus overcoming one of the important limitations mentioned above.
Further, note that the fractional-order nonlocal continuum model used in [47]
is frame-invariant, capable of dealing with asymmetric horizons, and it admits
integer-order boundary conditions. The thermodynamic consistency of the same
formulation has also been established in [50]. Second, we formulate a fully
consistent and highly accurate nonlinear f-FEM to numerically investigate the
response of the geometrically nonlinear fractional-order nonlocal beams.
Although several FE formulations for fractional-order equations have been
proposed in the literature, they are based on Galerkin or Petrov-Galerkin
methods that are capable of solving only linear hyperbolic and parabolic
differential equations [51, 52]. We develop a Ritz FEM that is capable of
obtaining the numerical solution of the nonlinear fractional-order elliptic
boundary value problem that describes the static response of the fractional-
order nonlocal and nonlinear beam. Further, by using the f-FEM we show that,
independently from the boundary conditions, the fractional-order theory
predicts a consistent softening behavior for the geometrically nonlinear
fractional-order beam as the degree of nonlocality increases. This is a direct
result of the well-posed nature of the fractional-order formulation which
accepts a unique solution [47]. With these results, we explain the paradoxical
predictions of hardening and the absence of nonlocal effects for certain
combinations of boundary conditions, as predicted by classical integral
approaches to nonlocal elasticity [27, 28].

The remainder of the paper is structured as follows: we begin with the
development of the nonlinear fractional-order beam theory based on the
fractional-order von Kármán strain-displacement relations. Next, we derive the
governing equations of the nonlinear fractional-order beam in strong form
using variational principles. Then, we derive a strategy for obtaining the
numerical solution to the nonlinear beam governing equation using fractional-
order FEM. Finally, we validate the nonlinear f-FEM, establish its
convergence, and then use it to analyze the effect of the fractional-order
nonlocality on the geometrically nonlinear static response of the beam under
different loading conditions.

##  2 Nonlinear and nonlocal fractional constitutive model

In this section, the geometrically nonlinear response of a nonlocal Euler-
Bernoulli beam is modeled using the fractional-order nonlocal continuum
formulation from [47]. Starting from this continuum theory, the nonlinear
constitutive relations for the nonlocal Euler-Bernoulli beam are derived.
Finally, the governing differential equations and the corresponding boundary
conditions for the fractional-order beam theory are derived using variational
principles.

###  2.1 Fundamentals of the fractional-order nonlocal continuum formulation

In analogy with the traditional approach to continuum modeling, we perform the
deformation analysis of a nonlocal solid by introducing two stationary
configurations, namely, the reference (undeformed) and the current (deformed)
configurations. The motion of the body from the reference configuration
(denoted as X) to the current configuration (denoted as x) is assumed as:

$\textbf{x}=\bm{\Psi}(\textbf{X})$ (1)

such that $\bm{\Psi}(\textbf{X})$ is a bijective mapping operation. The
relative position vector of two point particles located at $P_{1}$ and $P_{2}$
in the reference configuration of the nonlocal medium is denoted by
${\mathrm{d}\tilde{\textbf{X}}}$ (see Fig. (1)). After deformation due to the
motion $\bm{\Psi}(\textbf{X})$, the particles move to new positions $p_{1}$
and $p_{2}$, such that the relative position vector between them is
${\mathrm{d}\tilde{\textbf{x}}}$. It appears that
$\mathrm{d}\tilde{\textbf{X}}$ and $\mathrm{d}\tilde{\textbf{x}}$ are the
material and spatial differential line elements in the nonlocal medium,
conceptually analogous to the classical differential line elements
$\mathrm{d}{\textbf{X}}$ and $\mathrm{d}{\textbf{x}}$.

Figure 1:  (a) Schematic indicating the infinitesimal material
$\mathrm{d}\tilde{\textbf{X}}$ and spatial $\mathrm{d}\tilde{\textbf{x}}$ line
elements in the nonlocal medium under the displacement field u. (b) Horizon of
nonlocality and length scales at different material points. The current
nonlocal model can account for an asymmetric horizon condition that occurs at
points X close to a boundary or interface. The schematic is adapted from [47].

The mapping operation described in Eq. (1) is used to model the differential
line elements in the nonlocal solid using fractional calculus in order to
introduce nonlocality into the system [47]. More specifically, the
differential line elements of the nonlocal medium are modeled by imposing a
fractional-order transformation on the classical differential line elements as
follows:

$\mathrm{d}\tilde{\textbf{x}}=\big{[}D^{\alpha}_{\textbf{X}}\bm{\Psi}(\textbf{X})\big{]}\mathrm{d}\textbf{X}=\big{[}\tilde{\textbf{F}}_{X}(\textbf{X})\big{]}\mathrm{d}\textbf{X}$
(2a)
$\mathrm{d}\tilde{\textbf{X}}=\big{[}D^{\alpha}_{\textbf{x}}\bm{\Psi}^{-1}(\textbf{x})\big{]}\mathrm{d}\textbf{x}=\big{[}\tilde{\textbf{F}}_{x}(\textbf{x})\big{]}\mathrm{d}\textbf{x}$
(2b)

where $D^{\alpha}_{\square}\square$ is a space-fractional derivative whose
details will be presented below. Given the differ-integral nature of the
space-fractional derivative, the differential line elements
$\mathrm{d}\tilde{\textbf{X}}$ and $\mathrm{d}\tilde{\textbf{x}}$ have a
nonlocal character. Using the definitions for $\mathrm{d}\tilde{\textbf{X}}$
and $\mathrm{d}\tilde{\textbf{x}}$, the fractional deformation gradient tensor
$\overset{\alpha}{\textbf{F}}$ with respect to the nonlocal coordinates is
obtained in [47] as:

$\frac{\mathrm{d}\tilde{\textbf{x}}}{\mathrm{d}\tilde{\textbf{X}}}=\mathop{\textbf{F}}^{\alpha}=\tilde{\textbf{F}}_{X}\textbf{F}^{-1}\tilde{\textbf{F}}_{x}^{-1}$
(3)

where F is the classical deformation gradient tensor given as
$\textbf{F}=\mathrm{d}\textbf{x}/\mathrm{d}\textbf{X}$ given in local and
integer order form.

The space-fractional derivative
$D^{\alpha}_{\textbf{X}}\bm{\Psi}(\textbf{X},t)$ is taken here according to a
Riesz-Caputo (RC) definition with order $\alpha\in(0,1)$ defined on the
interval
$\textbf{X}\in(\textbf{X}_{A},\textbf{X}_{B})\subseteq\mathbb{R}^{3}$, and
given by:

$D^{\alpha}_{\textbf{X}}\bm{\Psi}(\textbf{X},t)=\frac{1}{2}\Gamma(2-\alpha)\big{[}\textbf{L}_{A}^{\alpha-1}~{}{}^{C}_{\textbf{X}_{A}}D^{\alpha}_{\textbf{X}}\bm{\Psi}(\textbf{X},t)-\textbf{L}_{B}^{\alpha-1}~{}{}^{C}_{\textbf{X}}D^{\alpha}_{\textbf{X}_{B}}\bm{\Psi}(\textbf{X},t)\big{]}$
(4)

where $\Gamma(\cdot)$ is the Gamma function, and
${}^{C}_{\textbf{X}_{A}}D^{\alpha}_{\textbf{X}}\bm{\Psi}$ and
${}^{C}_{\textbf{X}}D^{\alpha}_{\textbf{X}_{B}}\bm{\Psi}$ are the left- and
right-handed Caputo derivatives of $\bm{\Psi}$, respectively. The length
scales introduced in the above definition are given as
$\textbf{L}_{A}=\textbf{X}-\textbf{X}_{A}$ and
$\textbf{L}_{B}=\textbf{X}_{B}-\textbf{X}$ (see Fig. (1)). Before proceeding,
we highlight certain implications of the above definition of the fractional-
order derivative. The interval of the fractional derivative
$(\textbf{X}_{A},\textbf{X}_{B})$ defines the horizon of nonlocality (also
called attenuation range in classical nonlocal elasticity) which is
schematically shown in Fig. (1) for a generic point
$\textbf{X}\in\mathbb{R}^{2}$. The length scale parameters
$L_{A_{j}}^{\alpha-1}$ and $L_{B_{j}}^{\alpha-1}$ ensure the dimensional
consistency of the deformation gradient tensor, and along with the term
$\frac{1}{2}\Gamma(2-\alpha)$ ensure the frame invariance of the constitutive
relations. As shown in Appendix 1, the deformation gradient tensor introduced
via Eq. (3) enables an efficient and accurate treatment of the frame
invariance and completeness of kernel in the presence of asymmetric horizons,
material boundaries, and interfaces (Fig. (1)). More specifically, the use of
the different length scales $L_{A_{j}}^{\alpha-1}$ and $L_{B_{j}}^{\alpha-1}$
in Eq. (4) allows dealing with possible asymmetries in the horizon of
nonlocality (e.g. resulting from a truncation of the horizon, when approaching
a boundary point or an interface). We emphasize that this should not be
confused with the general concept of anisotropy wherein the length-scales
would potentially be different at each point in the nonlocal continuum (see,
for example, [41]).

The differential line elements are used to define strain measures in the
nonlocal medium. Here, we provide only the definition of strain measures in
Lagrangian coordinates. In analogy with the classical strain measures, the
fractional strain can be defined using the differential line elements as
$d\tilde{\textbf{x}}d\tilde{\textbf{x}}-d\tilde{\textbf{X}}d\tilde{\textbf{X}}$.
Using Eq. (3), the Lagrangian strain tensor in the nonlocal medium is obtained
as:

$\mathop{\textbf{E}}^{\alpha}=\frac{1}{2}(\mathop{\textbf{F}}^{\alpha}{}^{T}\mathop{\textbf{F}}^{\alpha}-\textbf{I})$
(5)

where I is the second-order identity tensor. Using kinematic position-
displacement relations, the expressions of the strains can be obtained in
terms of the displacement gradients. The fractional displacement gradient is
obtained using the definition of the fractional deformation gradient tensor
from Eq. (1) and the displacement field
$\textbf{U}(\textbf{X})=\textbf{x}(\textbf{X})-\textbf{X}$ as:

$\nabla^{\alpha}{\textbf{U}}_{X}=\tilde{\textbf{F}}_{X}-\textbf{I}$ (6)

The fractional gradient denoted by $\nabla^{\alpha}\textbf{U}_{X}$ is given as
$\nabla^{\alpha}\textbf{U}_{X_{ij}}=D^{\alpha}_{X_{j}}U_{i}$. Using the
nonlocal strain defined in Eq. (5) and the fractional deformation gradient
tensor $\overset{\alpha}{\textbf{F}}$ given in Eq. (3) together with Eq. (6),
the relationship between the strain and displacement gradient tensors is
obtained as:

$\mathop{\textbf{E}}^{\alpha}=\frac{1}{2}\bigr{(}\nabla^{\alpha}{\textbf{U}}_{X}+\nabla^{\alpha}{\textbf{U}}_{X}^{T}+\nabla^{\alpha}{\textbf{U}}_{X}^{T}\nabla^{\alpha}{\textbf{U}}_{X}\bigl{)}$
(7)

From the strain-displacement relations, assuming moderate rotations but small
strains, the fractional-order von-Kármán strain-displacement relations are
defined analogously to classical elasticity theories [53]. The expressions for
the full 3D nonlinear von-Kármán strain-displacement relations are:

$\tilde{\epsilon}_{ij}={\frac{1}{2}\left(D^{\alpha}_{X_{j}}u_{i}+D^{\alpha}_{X_{i}}u_{j}\right)}+{\frac{1}{2}\left(D^{\alpha}_{X_{i}}u_{3}~{}D^{\alpha}_{X_{j}}u_{3}\right)},~{}~{}~{}i,j=1,2$
(8)

where $u_{3}(\textbf{x})$ is the transverse displacement field and
$u_{k}(\textbf{x})$ $(k=1,2)$ are the in-plane variables. The transverse
strains (normal and shear) are simply the linearized forms of the respective
expressions derived from Eq. (7). Complete expressions for the components of
the fractional-order von-Kármán strains can be written from the above equation
to be:

$\tilde{\epsilon}_{11}=D_{x_{1}}^{\alpha}u_{1}+\frac{1}{2}(D_{x_{1}}^{\alpha}u_{3})^{2}$
(9a)
$\tilde{\epsilon}_{22}=D_{x_{2}}^{\alpha}u_{2}+\frac{1}{2}(D_{x_{2}}^{\alpha}u_{3})^{2}$
(9b)
$\tilde{\epsilon}_{12}=\tilde{\epsilon}_{21}=\frac{1}{2}\left(D_{x_{1}}^{\alpha}u_{2}+D_{x_{2}}^{\alpha}u_{1}\right)+\frac{1}{2}(D_{x_{1}}^{\alpha}u_{3}~{}D_{x_{2}}^{\alpha}u_{3})$
(9c) $\tilde{\epsilon}_{33}=D_{x_{3}}^{\alpha}u_{3}$ (9d)
$\tilde{\epsilon}_{13}=\tilde{\epsilon}_{31}=\frac{1}{2}\left(D_{x_{1}}^{\alpha}u_{3}+D_{x_{3}}^{\alpha}u_{1}\right)$
(9e)
$\tilde{\epsilon}_{23}=\tilde{\epsilon}_{32}=\frac{1}{2}\left(D_{x_{2}}^{\alpha}u_{3}+D_{x_{3}}^{\alpha}u_{2}\right)$
(9f)

which agree with the integer-order counterparts given in [53] for
$\alpha=1.0$.

In this study, we have analyzed the effect of the above nonlinear strain-
displacement relations on the response of an isotropic solid. The stress in
this isotropic solid is given as:

$\tilde{\sigma}_{ij}=C_{ijkl}\tilde{\epsilon}_{kl}$ (10)

where $C_{ijkl}$ is the constitutive matrix for the isotropic solid and
$\tilde{\epsilon}_{kl}$ is the fractional-order strain. We emphasize that the
stress defined through the above equation is nonlocal in nature. This follows
from the fractional-order definition of the deformation gradient tensor which
is then reflected in the nonlocal strain as evident from Eq. (5). The
thermodynamic consistency of the above constitutive relations has been
addressed in detail in [50]. As expected, classical continuum mechanics
relations are recovered when the order of the fractional derivative is set as
$\alpha=1$.

###  2.2 Nonlinear constitutive relations

We use the fractional-order nonlinear strain-displacement formulation
presented above to derive the constitutive model for geometrically nonlinear
nonlocal response of the beam. A schematic of a beam subject to transverse
loads is shown in Fig. (2). The Cartesian coordinates are chosen such that
$x_{3}=\pm h/2$ coincide with the top and bottom surfaces of the beam, and
$x_{1}=0$ and $x_{1}=L$ correspond to the two ends of the beam in the
longitudinal direction. The origin of the reference frame is chosen at the
left-end of the beam and $x_{3}=0$ is coincident with the mid-plane.

Figure 2: Schematic of an elastic beam subject to a distributed transverse
load $F_{t}(x_{1})$.

In the above defined coordinate system, the axial and transverse components of
the displacement field $\textbf{u}(x_{1},x_{3})$, denoted by
$u_{1}(x_{1},x_{3})$ and $u_{3}(x_{1},x_{3})$ at any spatial location
$\textbf{x}(x_{1},x_{3})$, are related to the mid-plane displacements
according to the Euler-Bernoulli assumptions:

$u_{1}(x_{1},x_{3})=u_{0}(x_{1})-x_{3}\left[\frac{\mathrm{d}w_{0}(x_{1})}{\mathrm{d}x_{1}}\right]$
(11a) $u_{3}(x_{1},x_{3})=w_{0}(x_{1})$ (11b)

where $u_{0}(x_{1})$ and $w_{0}(x_{1})$ are the axial and transverse
displacements of a point on the mid-plane, respectively. The fractional-order
nonlinear axial strain is obtained from the fractional-order von-Kármán
strain-displacement relations defined in Eq. (8). For the above displacement
field, the non-zero strain $\tilde{\epsilon}_{11}(x_{1},x_{3})$ is obtained
as:

$\tilde{\epsilon}_{11}(x_{1},x_{3})=D_{x_{1}}^{\alpha}u_{1}(x_{1},x_{3})+\frac{1}{2}(D_{x_{1}}^{\alpha}u_{3}(x_{1}))^{2}$
(12)

where $u_{1}$ and $u_{3}$ denote the displacement fields in the
$\textbf{x}_{1}$ and $\textbf{x}_{3}$ directions. For the Euler-Bernoulli
displacement field assumptions in Eq. (11), we obtain:

$\tilde{\epsilon}_{11}(x_{1},x_{3})=\tilde{\epsilon}_{0}(x_{1})+x_{3}\tilde{\kappa}(x_{1})$
(13)

where the fractional-order axial strain $\tilde{\epsilon}_{0}(x_{1})$ and
bending strain $\tilde{\kappa}(x_{1})$ are given as:

$\tilde{\epsilon}_{0}(x_{1})=D_{x_{1}}^{\alpha}u_{0}(x_{1})+\frac{1}{2}\left[D_{x_{1}}^{\alpha}w_{0}(x_{1})\right]^{2}$
(14a)
$\tilde{\kappa}(x_{1})=-D_{x_{1}}^{\alpha}\left[\frac{\mathrm{d}w_{0}(x_{1})}{\mathrm{d}x_{1}}\right]$
(14b)

In the above equation, $D^{\alpha}_{x_{1}}(\cdot)$ denotes the fractional-
order RC derivative given in Eq. (4) such that $\alpha\in(0,1)$.
$\textbf{x}_{A}(x_{A_{1}},0)$ and $\textbf{x}_{B}(x_{B_{1}},0)$ are the
terminals of the left- and right-handed Caputo derivatives within the RC
derivative. Consequently, the domain $(\textbf{x}_{A},\textbf{x}_{B})$ is the
horizon of nonlocal interactions for $\textbf{x}(x_{1},0)$, along the mid-
plane. Further, $l_{A}=x_{1}-x_{A_{1}}$ and $l_{B}=x_{B_{1}}-x_{1}$ are the
length scales to the left and right hand sides of the point
$\textbf{x}(x_{1},0)$ along the direction $\hat{\textbf{x}}_{1}$. As
emphasized in §2.1, these length scales are truncated for points near the
boundaries to account for the truncation of the nonlocal horizon. We have
illustrated this effect on the length scales $l_{A}$ and $l_{B}$ in Fig. (3),
where we have considered three points (P,Q, and R) at different distances from
the beam boundaries. At points P and R close to the boundaries, the asymmetry
in the nonlocal horizons is addressed by an appropriate truncation of the
relevant length scales ($l_{A}^{P}$ and $l_{B}^{R}$). At point Q, which is
sufficiently within the domain, a full horizon is captured by the complete
length scales on either side of the point. Note that equal values for
$l_{A}^{Q}$ and $l_{B}^{Q}$ represent an isotropic horizon.

Figure 3: A schematic presenting the variable nature of the length scale
corresponding to the horizon of nonlocality for different points along the
length of the beam. The length scale on the left side of the point P,
$l_{A}^{P}$ is truncated such that $l^{P}_{A}<l^{Q}_{A}$. Similarly,
$l^{R}_{B}<l^{Q}_{B}$.

The axial stress $\tilde{\sigma}_{11}$ resulting from the fractional-order
strain is evaluated using Eq. (10) as:

$\tilde{\sigma}_{11}(x_{1},x_{3})=E\tilde{\epsilon}_{11}(x_{1},x_{3})$ (15)

where $E$ denotes the Young’s Modulus of the isotropic elastic beam. The
deformation energy $\mathcal{U}$ accumulated in the nonlocal elastic beam due
to the above stress-strain distribution is obtained as:

$\mathcal{U}=\frac{1}{2}\int_{\Omega}\tilde{\sigma}_{11}(x_{1},x_{3})\tilde{\epsilon}_{11}(x_{1},x_{3})\mathrm{d}V$
(16)

where $\Omega$ denotes the volume occupied by the beam. Note that, by using
the definition for the nonlocal strain in Eq. (9e) in the Euler-Bernoulli
displacement field given in Eq. (11), a non-zero transverse shear strain would
be obtained. The remaining strains are trivially zero for the chosen
displacement field, analogous to classical formulations. The non-zero
transverse shear strains are not considered in the deformation energy in Eq.
(16) because for under slender beam assumptions, the rigidity to transverse
shear forces is much higher when compared to the bending rigidity. Hence, the
contribution of the transverse shear deformation towards the deformation
energy of the solid can be neglected.

The total potential energy functional, assuming that the beam is subjected to
distributed forces $F_{a}(x_{1})$ and $F_{t}(x_{1})$ acting on the mid-plane
along the axial and transverse directions, is given by:

$\Pi[\textbf{u}(\textbf{x})]=\mathcal{U}-\int_{0}^{L}F_{a}(x_{1})u_{0}(x_{1})\mathrm{d}x_{1}-\int_{0}^{L}F_{t}(x_{1})w_{0}(x_{1})\mathrm{d}x_{1}$
(17)

To obtain the total potential energy of the system, we have assumed no body
forces to be applied. Using the above constitutive model for the fractional-
order nonlocal and nonlinear elastic beam, we can derive the governing
differential equations in their strong form and the associated boundary
conditions by imposing optimality conditions on the total potential energy
functional $\Pi[\textbf{u}(\textbf{x})]$.

###  2.3 Governing equations

Using a variational approach, the governing differential equations and the
associated boundary conditions for the fractional-order beam are:

${}^{R-RL}_{x_{1}-l_{B}}D^{\alpha}_{x_{1}+l_{A}}\left[N(x_{1})\right]+F_{a}(x_{1})=0~{}~{}\forall~{}~{}x_{1}\in(0,L)$
(18a)
$\frac{\mathrm{d}}{\mathrm{d}x_{1}}\big{[}{}^{R-RL}_{x_{1}-l_{B}}D^{\alpha}_{x_{1}+l_{A}}\left[M(x_{1})\right]\big{]}+{}^{R-RL}_{x_{1}-l_{B}}D^{\alpha}_{x_{1}+l_{A}}\left[N(x_{1})D_{x_{1}}^{\alpha}\left[w_{0}(x_{1})\right]\right]+F_{t}(x_{1})=0~{}~{}\forall~{}~{}x_{1}\in(0,L)$
(18b) and subject to the boundary conditions:
$N(x_{1})=0~{}~{}\text{or}~{}~{}\delta
u_{0}(x_{1})=0~{}~{}\text{at}~{}~{}x_{1}\in\\{0,L\\}$ (18c)
$M(x_{1})=0~{}~{}\text{or}~{}~{}\delta\Bigg{[}\frac{\mathrm{d}w_{0}(x_{1})}{\mathrm{d}x_{1}}\Bigg{]}=0~{}~{}\text{at}~{}~{}x_{1}\in\\{0,L\\}$
(18d)
$\frac{\mathrm{d}M(x_{1})}{\mathrm{d}x_{1}}+N(x_{1})\frac{\mathrm{d}w_{0}(x_{1})}{\mathrm{d}x_{1}}=0~{}~{}\text{or}~{}~{}\delta
w_{0}(x_{1})=0~{}~{}\text{at}~{}~{}x_{1}\in\\{0,L\\}$ (18e)

In the above, $N(x_{1})$ and $M(x_{1})$ are the axial and bending stress
resultants, and are given by:

$N(x_{1})=\int_{-b/2}^{b/2}\int_{-h/2}^{h/2}\tilde{\sigma}_{11}~{}\mathrm{d}x_{3}~{}\mathrm{d}x_{2}$
(19a)
$M(x_{1})=\int_{-b/2}^{b/2}\int_{-h/2}^{h/2}x_{3}~{}\tilde{\sigma}_{11}~{}\mathrm{d}x_{3}~{}\mathrm{d}x_{2}$
(19b)

where $b$ denotes the width of the beam. In the fractional-order governing
equations provided in Eq. (18),
${}^{R-RL}_{x_{1}-l_{B}}D^{\alpha}_{x_{1}+l_{A}}[\cdot]$ is a Riesz-type
Riemann-Liouville (R-RL) derivative of order $\alpha\in(0,1)$. The R-RL
fractional derivative of an arbitrary function $f(x_{1})$ is defined analogous
to the RC derivative given in Eq. (4) as:

${}^{R-RL}_{x_{1}-l_{B}}D^{\alpha}_{x_{1}+l_{A}}f(x_{1})=\frac{1}{2}\Gamma(2-\alpha)\left[l_{B}^{\alpha-1}\left({}^{RL}_{x_{1}-l_{B}}D^{\alpha}_{x_{1}}f(x_{1})\right)-l_{A}^{\alpha-1}\left({}^{RL}_{x_{1}}D^{\alpha}_{x_{1}+l_{A}}f(x_{1})\right)\right]$
(20)

where ${}^{RL}_{x_{1}-l_{B}}D^{\alpha}_{x_{1}}f(x_{1})$ and
${}^{RL}_{x_{1}}D^{\alpha}_{x_{1}+l_{A}}f(x_{1})$ are the left- and right-
handed Riemann Liouville (RL) derivatives of an arbitrary function $f(x_{1})$
to the order $\alpha$, respectively. The detailed steps to obtain the strong
form of the governing equations by minimization of the total potential energy
are reported in Appendix 2.

Note that the governing equations for the axial and transverse displacements
are coupled, similar to what is seen in the classical nonlinear Euler-
Bernoulli beam formulation. This is unlike the linear elastic fractional-order
response of a nonlocal beam, where the governing equations for the axial and
transverse displacements are uncoupled [47]. Further, as expected, the
classical nonlinear Euler-Bernoulli beam governing equations and boundary
conditions are recovered for $\alpha=1$. As mentioned previously, the above
form of the governing equations for the fractional-order beam do not admit
closed-form analytical solutions for the most general loading and boundary
conditions. Therefore, in the following, we develop a fractional-order finite
element method (f-FEM) to obtain the numerical solution of the nonlinear
fractional-order governing equations. We will then use the f-FEM to
investigate the characteristic response of a nonlinear fractional-order beam.

##  3 Fractional Finite Element Method (f-FEM)

The fractional-order nonlocal and nonlinear boundary value problem described
by Eqs. (18a-18e) are numerically solved using finite element method. The
fractional-order boundary value problem builds upon the FE methods developed
in the literature for integral models of nonlocal elasticity [54, 55].
However, several modifications are necessary due to both the specific form of
the attenuation functions used in the definition of the fractional-order
derivatives, and to the nonlocal continuum model adopted in this study. Note
that, although a f-FEM formulation was developed in [47], the study dealt
exclusively with linear fractional-order boundary value problems, while here
we present the approach in the context of a nonlinear BVP. We also highlight
that, although the f-FEM is presented here for nonlinear fractional-order
beams, its formulation can be easily extended to model the nonlocal and
nonlinear behavior of plates and shells.

###  3.1 f-FEM formulation

We formulate the f-FEM starting from a discretized form of the total potential
energy functional $\Pi[\textbf{u}(\textbf{x})]$ given in Eq. (17). For this
purpose, the beam domain $\Omega=~{}[0,L]$ along the $x_{1}$ direction, is
uniformly discretized into disjoint two-noded elements
$\Omega^{e}_{i}=(x_{1}^{i},x_{1}^{i+1})$, such that
$\cup_{i=1}^{N_{e}}\Omega^{e}_{i}=\Omega$, and
$\Omega^{e}_{j}\cap\Omega^{e}_{k}=\emptyset~{}\forall~{}j\neq k$. $N_{e}$ is
the total number of discretized elements and the length of each element
$\Omega^{e}\in\Omega$ is $l_{e}$. The unknown displacement field variables
$u_{0}(x_{1})$ and $w_{0}(x_{1})$ at any point $x_{1}\in\Omega^{e}_{i}$ are
evaluated by interpolating the corresponding nodal degrees of freedom for
$\Omega^{e}_{i}$. From the definitions given in Eq. (13), it is clear that the
interpolation of the axial and transverse displacement fields following Euler-
Bernoulli beam theory would require at least Lagrangian ($\mathcal{L}_{i}$,
$i=1,2$) and Hermitian ($\mathcal{H}_{j}$, $j=1,2,3,4$) interpolation
functions, respectively. Consequently, the axial $u_{0}(x_{1})$ and transverse
$w_{0}(x_{1})$ displacement fields can be written as:

$\\{u_{0}(x_{1})\\}=[S_{u}(x_{1})]\\{U_{l}(x_{1})\\}$ (21a)
$\\{w_{0}(x_{1})\\}=[S_{w}(x_{1})]\\{W_{l}(x_{1})\\}$ (21b) where the
localized displacement coordinate vectors $\\{U_{l}(x_{1})\\}$ and
$\\{W_{l}(x_{1})\\}$ corresponding to the two-noded 1D element enclosing the
point $x_{1}$ are given as:
$\\{U_{l}(x_{1})\\}^{T}=\begin{bmatrix}u_{0}^{(1)}&u_{0}^{(2)}\end{bmatrix}$
(21c)
$\\{W_{l}(x_{1})\\}^{T}=\begin{bmatrix}w_{0}^{(1)}&\frac{\mathrm{d}w_{0}^{(1)}}{\mathrm{d}x_{1}}&w_{0}^{(2)}&\frac{\mathrm{d}w_{0}^{(2)}}{\mathrm{d}x_{1}}\end{bmatrix}$
(21d) and the associated shape function matrices are:
$[S_{u}(x_{1})]=\begin{bmatrix}\mathcal{L}_{1}^{e}(x_{1})&\mathcal{L}_{2}^{e}(x_{1})\end{bmatrix}$
(21e)
$[S_{w}(x_{1})]=\begin{bmatrix}\mathcal{H}_{1}^{e}(x_{1})&\mathcal{H}_{2}^{e}(x_{1})&\mathcal{H}_{3}^{e}(x_{1})&\mathcal{H}_{4}^{e}(x_{1})\end{bmatrix}$
(21f)

The superscripts $(\cdot)^{(1)}$ and $(\cdot)^{(2)}$ in the above equation
denote the local node numbers for the two-noded element. Using the definition
of the RC derivative in Eq. (4), the fractional derivative
$D_{x_{1}}^{\alpha}\left[u_{0}(x_{1})\right]$ is expressed as:

$D_{x_{1}}^{\alpha}\left[u_{0}(x_{1})\right]=\frac{1}{2}(1-\alpha)\left[l_{A}^{\alpha-1}\int_{x_{1}-l_{A}}^{x_{1}}\frac{D^{1}[u_{0}(s_{1})]}{(x_{1}-s_{1})^{\alpha}}~{}\mathrm{d}s_{1}+l_{B}^{\alpha-1}\int_{x_{1}}^{x_{1}+l_{B}}\frac{D^{1}[u_{0}(s_{1})]}{(s_{1}-x_{1})^{\alpha}}~{}\mathrm{d}s_{1}\right]$
(22)

where $D^{1}(\cdot)$ is the first integer-order derivative and $s_{1}$ is a
dummy variable along the direction $\textbf{x}_{1}$ used within the definition
of the fractional-order derivative. The above expression can be recast as:

$D_{x_{1}}^{\alpha}\left[u_{0}(x_{1})\right]=\int_{x_{1}-l_{A}}^{x_{1}}A(x_{1},s_{1},l_{A},\alpha)~{}D^{1}\left[u_{0}(s_{1})\right]~{}\mathrm{d}s_{1}+\int_{x_{1}}^{x_{1}+l_{B}}A(x_{1},s_{1},l_{B},\alpha)~{}D^{1}\left[u_{0}(s_{1})\right]~{}\mathrm{d}s_{1}$
(23a) where
$A(x,y,l,\alpha)=\frac{(1-\alpha)}{2}~{}l^{\alpha-1}~{}\frac{1}{|x-y|^{\alpha}}$
(23b)

is the kernel of the fractional derivative. Note that kernel $A(x,y,l,\alpha)$
is a function of the relative distance between the points $x$ and $y$, and can
be interpreted as an analog form of the attenuation functions typically used
in integral models of nonlocal elasticity. Clearly, the attenuation decays as
a power-law in the distance with an exponent equal to the order $\alpha$ of
the fractional derivative. Similar expressions for
$D^{\alpha}_{x_{1}}\left[w_{0}(x_{1})\right]$ and
$D^{\alpha}_{x_{1}}\left[\mathrm{d}w_{0}(x_{1})/\mathrm{d}x_{1}\right]$ can
also be derived using an equivalent approach.

Note that $D^{\alpha}_{x_{1}}\left[u_{0}(x_{1})\right]$,
$D^{\alpha}_{x_{1}}\left[w_{0}(x_{1})\right]$ and
$D^{\alpha}_{x_{1}}\left[\mathrm{d}w_{0}(x_{1})/\mathrm{d}x_{1}\right]$
contain the integer-order derivatives $D^{1}_{s_{1}}[u_{0}(s_{1})]$,
$D^{1}_{s_{1}}[w_{0}(s_{1})]$ and
$D^{1}_{s_{1}}[\mathrm{d}w_{0}(s_{1})/\mathrm{d}s_{1}]$ (equivalent to
$D^{2}_{s_{1}}[w_{0}(s_{1})]$). These integer-order derivatives at $s_{1}$ are
evaluated in terms of the nodal variables corresponding to the element
$\Omega_{p}^{e}$, such that $s_{1}\in\Omega_{p}^{e}$. Using Eq. (21), the
integer order derivatives can be expressed as:

$D^{1}[u_{0}(s_{1})]=[B_{u}(s_{1})]\\{U_{l}(s_{1})\\},~{}~{}~{}D^{1}[w_{0}(s_{1})]=[B_{w}(s_{1})]\\{W_{l}(s_{1})\\},~{}~{}~{}D^{2}[w_{0}(s_{1})]=[B_{\theta}(s_{1})]\\{W_{l}(s_{1})\\}$
(24)

where the $[B_{\square}]$ ($\square\in\\{u,w,\theta\\}$) matrices are
expressed in terms of the associated shape functions for the two-noded element
as:

$[B_{u}(s_{1})]=\begin{bmatrix}\frac{\mathrm{d}\mathcal{L}^{e}_{1}}{\mathrm{d}s_{1}}&\frac{\mathrm{d}\mathcal{L}^{e}_{2}}{\mathrm{d}s_{1}}\end{bmatrix}$
(25a)
$[B_{w}(s_{1})]=\begin{bmatrix}\frac{\mathrm{d}\mathcal{H}_{1}^{e}}{\mathrm{d}s_{1}}&\frac{\mathrm{d}\mathcal{H}_{2}^{e}}{\mathrm{d}s_{1}}&\frac{\mathrm{d}\mathcal{H}_{3}^{e}}{\mathrm{d}s_{1}}&\frac{\mathrm{d}\mathcal{H}_{4}^{e}}{\mathrm{d}s_{1}}\end{bmatrix}$
(25b)
$[B_{\theta}(s_{1})]=\begin{bmatrix}\frac{\mathrm{d}^{2}\mathcal{H}_{1}^{e}}{\mathrm{d}s_{1}^{2}}&\frac{\mathrm{d}^{2}\mathcal{H}_{2}^{e}}{\mathrm{d}^{2}s_{1}}&\frac{\mathrm{d}^{2}\mathcal{H}_{3}^{e}}{\mathrm{d}^{2}s_{1}}&\frac{\mathrm{d}^{2}\mathcal{H}_{4}^{e}}{\mathrm{d}^{2}s_{1}}\end{bmatrix}$
(25c)

Using the above expression for the integer-order derivative
$D^{1}[u_{0}(x_{1})]$, the fractional-order derivative in Eq. (23) can be
written as:

$D_{x_{1}}^{\alpha}\left[u_{0}(x_{1})\right]=\int_{x_{1}-l_{A}}^{x_{1}+l_{B}}A_{p}(x_{1},s_{1},l_{A},l_{B},\alpha)[B_{u}(s_{1})]\\{U_{l}(s_{1})\\}\mathrm{d}s_{1}$
(26a) where the following definition for the attenuation function
$A_{p}(x_{1},s_{1},l_{A},l_{B},\alpha)$ is employed:
$A_{p}(x_{1},s_{1},l_{A},l_{B},\alpha)=\begin{cases}A(x_{1},s_{1},l_{A},\alpha)&~{}~{}s_{1}\in{(x_{1}-l_{A},x_{1})},\\\
A(x_{1},s_{1},l_{B},\alpha)&~{}~{}s_{1}\in{(x_{1},x_{1}+l_{B})}.\end{cases}$
(26b)

Similar expressions are derived for $D^{\alpha}_{x_{1}}[w_{0}(x_{1})]$ and
$D^{\alpha}_{x_{1}}[\mathrm{d}w_{0}(x_{1})/\mathrm{d}x_{1}]$:

$D_{x_{1}}^{\alpha}\left[w_{0}(x_{1})\right]=\int_{x_{1}-l_{A}}^{x_{1}+l_{B}}A_{p}(x_{1},s_{1},l_{A},l_{B},\alpha)[B_{w}(s_{1})]\\{W_{l}(s_{1})\\}\mathrm{d}s_{1}$
(27a)
$D_{x_{1}}^{\alpha}\left[\frac{\mathrm{d}w_{0}(x_{1})}{\mathrm{d}x_{1}}\right]=\int_{x_{1}-l_{A}}^{x_{1}+l_{B}}A_{p}(x_{1},s_{1},l_{A},l_{B},\alpha)[B_{\theta}(s_{1})]\\{W_{l}(s_{1})\\}\mathrm{d}s_{1}$
(27b)

Note that the evaluation of the fractional derivative requires a convolution
of the integer-order derivatives across the horizon of nonlocality
$(x_{1}-l_{A},x_{1}+l_{B})$. Thus, while obtaining the FE approximation in Eq.
(26a), the nonlocal contributions from the different finite elements in the
horizon must be properly attributed to the corresponding nodes. In order
perform this mapping of the nonlocal contributions from elements in the
horizon, we transform the nodal values $\\{U_{l}(s_{1})\\}$ and
$\\{W_{l}(s_{1})\\}$ into the respective global variable vectors $\\{U_{g}\\}$
and $\\{W_{g}\\}$ with the help of connectivity matrices in the following
manner:

$\\{U_{l}(s_{1})\\}=[\tilde{\mathcal{C}_{u}}(x_{1},s_{1})]\\{U_{g}\\}$ (28a)
$\\{W_{l}(s_{1})\\}=[\tilde{\mathcal{C}_{w}}(x_{1},s_{1})]\\{W_{g}\\}$ (28b)

The connectivity matrices $[\tilde{\mathcal{C}_{u}}(x_{1},s_{1})]$ and
$[\tilde{\mathcal{C}_{w}}(x_{1},s_{1})]$ are defined for the axial and
transverse displacement vectors, respectively. These matrices are designed
such that they are non-zero only if $s_{1}$ lies in the domain
$(x_{1}-l_{A},x_{1}+l_{B})$, the horizon of nonlocality for $x_{1}$. It is
evident that these matrices activate the contribution of the nodes enclosing
$s_{1}$ for the numerical evaluation of the convolution integral given in Eq.
(26a) [47]. Using the above formalism, Eq. (26a) can be recast as:

$D^{\alpha}_{x_{1}}\left[u_{0}(x_{1})\right]=[\tilde{B}_{u}(x_{1})]\\{U_{g}\\}$
(29a) where
$[\tilde{B}_{u}(x_{1})]=\int_{x_{1}-l_{A}}^{x_{1}+l_{B}}A_{p}(x_{1},s_{1},l_{A},l_{B},\alpha)[B_{u}(s_{1})][\tilde{\mathcal{C}_{u}}(x_{1},s_{1})]~{}\mathrm{d}s_{1}$
(29b)

Similarly, we express $D^{\alpha}_{x_{1}}\left[w_{0}(x_{1})\right]$ and
$D^{\alpha}_{x_{1}}\left[{\mathrm{d}w_{0}(x_{1})}/{\mathrm{d}x_{1}}\right]$
as:

$D^{\alpha}_{x_{1}}\left[w_{0}(x_{1})\right]=[\tilde{B}_{w}(x_{1})]\\{W_{g}\\}$
(30a)
$D^{\alpha}_{x_{1}}\left[\frac{\mathrm{d}w_{0}(x_{1})}{\mathrm{d}x_{1}}\right]=[\tilde{B}_{\theta}(x_{1})]\\{W_{g}\\}$
(30b) where
$[\tilde{B}_{w}(x_{1})]=\int_{x_{1}-l_{A}}^{x_{1}+l_{B}}A_{p}(x_{1},s_{1},l_{A},l_{B},\alpha)[B_{w}(s_{1})][\tilde{\mathcal{C}_{w}}(x_{1},s_{1})]~{}\mathrm{d}s_{1}$
(30c)
$[\tilde{B}_{\theta}(x_{1})]=\int_{x_{1}-l_{A}}^{x_{1}+l_{B}}A_{p}(x_{1},s_{1},l_{A},l_{B},\alpha)[B_{\theta}(s_{1})][\tilde{\mathcal{C}_{w}}(x_{1},s_{1})]~{}\mathrm{d}s_{1}$
(30d)

Note that $[B_{u}(s_{1})]$, $[B_{w}(s_{1})]$ and $[B_{\theta}(s_{1})]$ have
been given in Eq. (25).

###  3.2 Nonlinear f-FEM model

In the following, we will use the above formalism to obtain the algebraic
equations of motion for the geometrically nonlinear response of the
fractional-order nonlocal beam by minimizing the total potential energy of the
beam. The first variation of the deformation energy $\mathcal{U}$ defined in
Eq. (16) is given by:

$\delta\mathcal{U}=b\int_{0}^{L}\int_{-h/2}^{h/2}\delta\tilde{\epsilon}_{11}(x_{1},x_{3})~{}\tilde{\sigma}_{11}(x_{1},x_{3})\mathrm{d}x_{3}\mathrm{d}x_{1}$
(31)

Using the above expression for $\delta\mathcal{U}$ along with the nonlocal
axial strain derived in Eq. (13) for the Euler-Bernoulli beam theory and the
definitions of the stress resultants given in Eq. (19), we obtain the
following expression for the statement of the minimum potential energy
$\delta\Pi=0$:

$\begin{split}\int_{0}^{L}\left\\{\left(D_{x_{1}}^{\alpha}\left[\delta
u_{0}(x_{1})\right]+D_{x_{1}}^{\alpha}\left[w_{0}(x_{1})\right]~{}D_{x_{1}}^{\alpha}\left[\delta
w_{0}(x_{1})\right]\right)N(x_{1})-D_{x_{1}}^{\alpha}\left[\frac{\mathrm{d}\delta
w_{0}(x_{1})}{\mathrm{d}x_{1}}\right]M(x_{1})\right\\}\mathrm{d}x_{1}&\\\
-\int_{0}^{L}F_{t}(x_{1})\delta
w_{0}(x_{1})\mathrm{d}x_{1}-\int_{0}^{L}F_{a}(x_{1})\delta
u_{0}(x_{1})\mathrm{d}x_{1}&=0\end{split}$ (32)

Note that the variations $\delta u_{0}(x_{1})$ and $\delta w_{0}(x_{1})$ are
independent of each other. Now, by collecting the terms containing $\delta
u_{0}(x_{1})$ and $\delta w_{0}(x_{1})$, the above weak form is expressed in
the form of the following equations:

$\int_{0}^{L}\left\\{D_{x_{1}}^{\alpha}\left[\delta
u_{0}(x_{1})\right]N(x_{1})-F_{a}(x_{1})\delta
u_{0}(x_{1})\right\\}\mathrm{d}x_{1}=0$ (33a)
$\int_{0}^{L}\left\\{\left(D_{x_{1}}^{\alpha}\left[w_{0}(x_{1})\right]~{}D_{x_{1}}^{\alpha}\left[\delta
w_{0}(x_{1})\right]\right)N(x_{1})-D_{x_{1}}^{\alpha}\left[\frac{\mathrm{d}\delta
w_{0}(x_{1})}{\mathrm{d}x_{1}}\right]M(x_{1})-F_{t}(x_{1})\delta
w_{0}(x_{1})\right\\}\mathrm{d}x_{1}=0$ (33b)

We recast the above equations in terms of of the displacement field variables
by using the expressions for the stress resultants given in Eqs. (10) and
(19), and obtain the following:

$\int_{0}^{L}\Bigg{\\{}A_{11}\left(D_{x_{1}}^{\alpha}\left[\delta
u_{0}(x_{1})\right]\right)\big{(}D_{x_{1}}^{\alpha}\left[u_{0}(x_{1})\right]+\frac{1}{2}\left(D_{x_{1}}^{\alpha}\left[w_{0}(x_{1})\right]\right)^{2}\big{)}-F_{a}(x_{1})\delta
u_{0}(x_{1})\Bigg{\\}}\mathrm{d}x_{1}=0$ (34a)
$\begin{split}\int_{0}^{L}&\Bigg{\\{}A_{11}~{}\bigg{(}D_{x_{1}}^{\alpha}\left[w_{0}(x_{1})\right]~{}D_{x_{1}}^{\alpha}\left[\delta
w_{0}(x_{1})\right]\bigg{)}~{}\big{(}D_{x_{1}}^{\alpha}\left[u_{0}(x_{1})\right]+\frac{1}{2}\left(D_{x_{1}}^{\alpha}\left[w_{0}(x_{1})\right]\right)^{2}\big{)}\\\
&+D_{11}~{}\left(D_{x_{1}}^{\alpha}\left[\frac{\mathrm{d}\delta
w_{0}(x_{1})}{\mathrm{d}x_{1}}\right]\right)~{}\left(D_{x_{1}}^{\alpha}\left[\frac{\mathrm{d}w_{0}(x_{1})}{\mathrm{d}x_{1}}\right]\right)-F_{t}(x_{1})~{}\delta
w_{0}(x_{1})\Bigg{\\}}~{}\mathrm{d}x_{1}=0\end{split}$ (34b)

where $A_{11}=Ebh$ and $D_{11}=Ebh^{3}/12$ are the axial and bending stiffness
of the beam, respectively. The above equations in the weak form are converted
into algebraic equations using the f-FEM model discussed in §3.1. Now, by
using Eqs. (29-30) and enforcing minimization, the above integrals reduce to
the following algebraic equations:

$\begin{bmatrix}[\tilde{K}_{11}]&[\tilde{K}_{12}]\\\
[\tilde{K}_{21}]&[\tilde{K}_{22}]\end{bmatrix}\begin{Bmatrix}\\{U_{g}\\}\\\
\\{W_{g}\\}\end{Bmatrix}=\begin{Bmatrix}\\{F_{A}\\}\\\
\\{F_{T}\\}\end{Bmatrix}$ (35)

In the above system of equations, the different stiffness matrices are
obtained as:

$[\tilde{K}_{11}]=\int_{0}^{L}A_{11}[\tilde{B}_{u}(x_{1})]^{T}[\tilde{B}_{u}(x_{1})]~{}\mathrm{d}x_{1}$
(36a)
$[\tilde{K}_{12}]=\frac{1}{2}\int_{0}^{L}A_{11}\left(D_{x_{1}}^{\alpha}[w_{0}(x_{1})]\right)[\tilde{B}_{u}(x_{1})]^{T}[\tilde{B}_{w}(x_{1})]~{}\mathrm{d}x_{1}$
(36b)
$[\tilde{K}_{21}]=\int_{0}^{L}A_{11}\left(D_{x_{1}}^{\alpha}[w_{0}(x_{1})]\right)[\tilde{B}_{w}(x_{1})]^{T}[\tilde{B}_{u}(x_{1})]~{}\mathrm{d}x_{1}$
(36c)
$[\tilde{K}_{22}]=\int_{0}^{L}D_{11}[\tilde{B}_{\theta}(x_{1})]^{T}[\tilde{B}_{\theta}(x_{1})]~{}\mathrm{d}x_{1}+\frac{1}{2}\int_{0}^{L}\left[A_{11}\left(D_{x_{1}}^{\alpha}[w_{0}(x_{1})]\right)^{2}\right][\tilde{B}_{w}(x_{1})]^{T}[\tilde{B}_{w}(x_{1})]~{}\mathrm{d}x_{1}$
(36d) and the axial and transverse nodal force vectors are obtained as:
$\\{F_{A}\\}^{T}=\int_{0}^{L}F_{a}(x_{1})~{}[S_{u}(x_{1})]~{}\mathrm{d}x_{1}$
(36e)
$\\{F_{T}\\}^{T}=\int_{0}^{L}F_{t}(x_{1})~{}[S_{w}(x_{1})]~{}\mathrm{d}x_{1}$
(36f)

The solution of the system of algebraic equations in Eq. (35) provides the
nodal generalized displacement coordinates, which can then be used along with
the Euler-Bernoulli relations in Eq. (11) to determine the displacement field
at any point of the beam. The effect of the geometric nonlinearity is evident
from the expressions of the stiffness matrices given in Eq. (36), where it
appears that the stiffness matrices are a function of the deformed system
configuration. Note that the linear f-FEM model can be obtained by ignoring
the contribution of $D_{x_{1}}^{\alpha}[w_{0}(x_{1})]$ in Eq. (36). Further,
similar to classical beam models accounting for geometric nonlinearity, a
coupling between the axial and bending responses is obtained through the non-
zero $[\tilde{K}_{12}]$ and $[\tilde{K}_{21}]$ matrices.

In this study, we solve the nonlinear algebraic equations using an
incremental-iterative Newton-Raphson (NR) numerical scheme. Similar to
classical nonlinear models, the NR procedure for the fractional-order
nonlinear equations also requires the evaluation of the tangent stiffness
matrix. The procedure to derive the tangent stiffness matrix
$[\tilde{K}_{T}^{i}]$ corresponding to the $i$-th iteration of the NR method
follows from the classical approach:

$[\tilde{K}_{T}^{i}]=\left[\frac{\partial\\{R\\}}{\partial\\{X\\}}\right]^{i}$
(37)

where $\\{X\\}$ is a vector containing all the displacement coordinates of the
beam and $\\{R\\}$ is the residual vector evaluated at the current iteration.
The residual $\\{R\\}$ at the $i$-th iteration step is given by:

$\\{R(\\{X\\})\\}^{i}=[\tilde{K}_{S}^{i}]\\{X\\}^{i}-\\{F\\}$ (38)

where $\tilde{K}_{S}^{i}$ denotes the stiffness matrix of the entire system.
In order to obtain the solution of the nonlinear equations, we also adopt a
load increment procedure wherein the nonlocal response of the beam at each
load step is evaluated using the NR iteration:

$\\{X\\}^{i+1}=\\{X\\}^{i}-\left[\tilde{K}_{T}^{i}\right]^{-1}\\{R(\\{X\\})\\}^{i}$
(39)

The iterations at each load level are continued until the residual becomes
less than a chosen tolerance.

###  3.3 Numerical integration for nonlocal matrices

In this following, we provide the details of the numerical scheme used to
integrate both the stiffness and the tangent stiffness matrices. Note that the
numerical procedure to evaluate the force vector follows directly from
classical FE formulations. Hence, we do not provide a complete account of all
the steps. Evaluation of the stiffness matrices for the nonlocal system given
in Eq. (36) requires the evaluation of the nonlocal matrices
$[\tilde{B}_{\square}]$ $(\square\in\\{u,w,\theta\\})$. As seen in Eqs.
(29-30), this involves a convolution of the integer-order derivatives with the
fractional-order attenuation function $A_{p}(x_{1},s_{1},l_{A},l_{B},\alpha)$
over the nonlocal horizon $(x_{1}-l_{A},x_{1}+l_{B})$. Clearly, the numerical
model for fractional-order derivative involves additional integration over the
nonlocal horizon to account for nonlocal interactions. These nonlocal
interactions across the horizon of nonlocality have already been accounted for
numerically in [15, 54]. However, differently from these studies, the
attenuation function $A_{p}(x_{1},s_{1},l_{A},l_{B},\alpha)$ in the
fractional-order model involves an end-point singularity (more specifically,
at $x_{1}=s_{1}$) due to the nature of the kernel (see Eq. (26b)). The
fractional-order nonlocal interactions as well as the end-point singularity
were addressed in [47] and are briefly reviewed here.

In the following, we describe the procedure for the numerical evaluation of
the nonlinear stiffness matrix $[\tilde{K}_{12}]$. The same procedure can be
extended to evaluate the other stiffness matrices, hence their discussion is
not provided here. In order to perform the numerical integration of the
nonlocal and nonlinear stiffness matrix $[\tilde{K}_{12}]$, we adopt an
isoparametric formulation and introduce a natural coordinate system $\xi_{1}$.
The Jacobian of the transformation $x_{1}\rightarrow\xi_{1}$ is given as
$J(\xi_{1})$. Now, by using the Gauss-Legendre quadrature rule, the matrix
$[\tilde{K}_{12}]$ is approximated as:

$[\tilde{K}_{12}]\approx\sum\limits_{i=1}^{N_{e}}\sum\limits_{j=1}^{N_{GP}}\hat{w}_{j}J^{i}\left\\{A_{11}\left(D_{x_{1}}^{\alpha}\left[w_{0}(\xi_{1}^{i,j})\right]\right)\right\\}[\tilde{B}_{u}(\xi_{1}^{i,j})]^{T}[\tilde{B}_{w}(\xi_{1}^{i,j})]$
(40)

where $\xi_{1}^{i,j}$ is the $j-$th Gauss integration point in the $i-$th
element, $\hat{w}_{j}$ is the corresponding weight for numerical integration,
and $N_{GP}$ is the total number of Gauss points chosen for the numerical
integration, such that $j\in\\{1,...N_{GP}\\}$. $J^{i}$ is the Jacobian of the
coordinate transformation for the $i-$th element. The hat symbol on the weight
is used to distinguish it from the transverse displacement. As previously
highlighted, the matrices $[\tilde{B}_{\square}]$
($\square\in\\{u,w,\theta\\}$) (equivalently,
$[\tilde{B}_{\square}(\xi_{1}^{i,j})]$ in the discretized form) involve an
additional integration (see Eq. (29b)) due to the fractional-order nonlocality
and are evaluated in the following manner:

$\left[\tilde{B}_{\square}(x_{1}^{i,j})\right]=\int_{x_{1}^{i,j}-l_{A}}^{x_{1}^{i,j}+l_{B}}A_{p}(x_{1}^{i,j},s_{1},l_{A},l_{B},\alpha)[B_{\square}(s_{1})][\tilde{\mathcal{C}_{\square}}(x_{1}^{i,j},s_{1})]~{}\mathrm{d}s_{1}$
(41)

where $x_{1}^{i,j}$ is the Cartesian coordinate of the Gauss point
$\xi_{1}^{i,j}$ and $[B_{\square}(s_{1})]$ is given in Eq. (25). Using Eq.
(26b), we obtain the following expression for
$\left[\tilde{B}_{\square}(x_{1}^{i,j})\right]$:

$\left[\tilde{B}_{\square}(x_{1}^{i,j})\right]=\int_{x_{1}^{i,j}-l_{A}}^{x_{1}^{i,j}}\mathcal{I}_{L}\mathrm{d}s_{1}+\int_{x_{1}^{i,j}}^{x_{1}^{i,j}+l_{B}}\mathcal{I}_{R}\mathrm{d}s_{1}$
(42a) where, the integrands $\mathcal{I}_{L}$ and $\mathcal{I}_{R}$ are given
as:
$\mathcal{I}_{L}=A(x_{1}^{i,j},s_{1},l_{A},\alpha)[B_{\square}(s_{1})][\tilde{\mathcal{C}_{u}}(x_{1}^{i,j},s_{1})]$
(42b)
$\mathcal{I}_{R}=A(x_{1}^{i,j},s_{1},l_{B},\alpha)[B_{\square}(s_{1})][\tilde{\mathcal{C}_{u}}(x_{1}^{i,j},s_{1})]$
(42c)

Note that the limits of the integrals in Eq. (42a) span over the elements that
constitute the nonlocal horizon $(x_{1}^{i,j}-l_{A},x_{1}^{i,j}+l_{B})$ at
$x_{1}^{i,j}$. These limits for the convolution integral follow from the
terminals of the RC fractional derivative defined for the study of nonlocal
continuum in Eq. (4). These integrals are evaluated numerically in the
following manner:

$\int_{x_{1}^{i,j}-l_{A}}^{x_{1}^{i,j}}\mathcal{I}_{L}\mathrm{d}s_{1}\approx\underbrace{\int_{x_{1}^{i-N_{A}^{inf}}}^{x_{1}^{i-(N_{A}^{inf}-1)}}\mathcal{I}_{L}\mathrm{d}s_{1}+...\int_{x_{1}^{i-1}}^{x_{1}^{i}}\mathcal{I}_{L}\mathrm{d}s_{1}}_{\text{Gauss-
Legendre
Quadrature}}+\underbrace{\int_{x_{1}^{i}}^{x_{1}^{i,j}}\mathcal{I}_{L}\mathrm{d}s_{1}}_{\text{Singularity
at }x_{1}^{i,j}}$ (43a)
$\int_{x_{1}^{i,j}}^{x_{1}^{i,j}+l_{B}}\mathcal{I}_{R}\mathrm{d}s_{1}\approx\underbrace{\int_{x_{1}^{i,j}}^{x_{1}^{i+1}}\mathcal{I}_{R}\mathrm{d}s_{1}}_{\text{Singularity
at
}x_{1}^{i,j}}+\underbrace{\int_{x_{1}^{i+1}}^{x_{1}^{i+2}}\mathcal{I}_{R}\mathrm{d}s_{1}...+\int_{x_{1}^{i+(N_{B}^{inf}-1)}}^{x_{1}^{i+N_{B}^{inf}}}\mathcal{I}_{R}\mathrm{d}s_{1}}_{\text{Gauss-
Legendre Quadrature}}$ (43b)

In the above expressions, $N_{A}^{inf}$ and $N_{B}^{inf}$ are the number of
(complete) elements in the nonlocal horizon to the left and right side of
$x_{1}$, respectively. More specifically, $N_{A}^{inf}=\lceil
l_{A}/l_{e}\rceil$ and ${N_{B}^{inf}}=\lfloor{l_{B}/l_{e}}\rfloor$ where
$l_{e}$ is the length of the discretized element. The ceil
($\lceil\cdot\rceil$) and floor ($\lfloor{\cdot}\rfloor$) functions are used
to round the number of elements to the greater integer on the left side and
the lower integer on the right side. For $x_{1}^{i,j}$ close to the boundaries
of the beam ($x_{1}=\\{0,L\\}$), $N_{A}^{inf}$ and $N_{B}^{inf}$ are truncated
in order to account for asymmetric horizon lengths. This is essential to
satisfy frame-invariance as discussed in §2.

As discussed previously, due to the nature of the kernel of the fractional-
order derivative, an end-point singularity occurs in the integrals at the
Gauss point $x_{1}^{i,j}$ in the element $\Omega_{i}^{e}$. This is evident
from the definitions of the left and right integrals given in Eq. (43).
Following [47], this end-point singularity is circumvented by an analytical
evaluation of these integrals over the elements containing the singularities.
This analytical evaluation can be carried out by using the expressions for
$[\tilde{B}_{\square}]$ ($\square\in u,w,\theta$) given in Eq. (29b). The
integrals over the remaining elements (i.e. those without singularities) are
evaluated using the Gauss-Legendre quadrature method. The expression for this
quadrature based integration corresponding to the nonlocal contribution of the
$r-$th element in the horizon of nonlocality of the Gauss point $x_{1}^{i,j}$
is given by:

$\begin{split}\int_{x_{1}^{r}}^{x_{1}^{r+1}}A_{p}(x^{i,j}_{1},s_{1},l_{A},l_{B},\alpha)[B_{u}(s_{1})]&[\tilde{\mathcal{C}_{u}}(x_{1}^{i,j},s_{1})]~{}\mathrm{d}s_{1}=\\\
&\sum_{k=1}^{N_{GP}}\hat{w}_{k}J^{r}A_{p}(x^{i,j}_{1},{x}_{1}^{r,k},l_{A},l_{B},\alpha)[B_{u}({s}_{1}^{r,k})][\tilde{\mathcal{C}_{u}}(x_{1}^{i,j},{x}_{1}^{r,k})]\end{split}$
(44)

where $x_{1}^{r,k}$ is the Cartesian coordinate of the $k-$th Gauss point in
the $r-$th element along $\textbf{x}_{1}$, $\hat{w}_{k}$ is the corresponding
weight, and $J^{r}$ is the Jacobian of the transformation for $r-$th element.
The above numerical procedure is schematically illustrated in Fig. 4. We
highlight that $A_{p}$ is a function of the Cartesian coordinates, therefore
global coordinates should be used in its evaluation.

Figure 4: Illustration of the influence zone for the Gauss point $x_{1}^{i,2}$
and the evaluation of $[\tilde{B_{u}}(x_{1}^{i,2})]$.

As mentioned previously, the procedure followed for the evaluation of the
$[\tilde{K}_{12}]$ stiffness matrix extends directly to the other stiffness
matrices defined in Eq. (36), as well as the tangent stiffness matrix defined
in Eq. (37). Note that all these matrices are pre-assembled in the global form
due to the use of the connectivity matrices in Eq. (28). The resulting
matrices are then used in the NR iteration procedure outlined in Eq. (39) to
obtain the generalized displacement coordinates for the geometrically
nonlinear response of the fractional-order nonlocal beam.

##  4 Application and Numerical Results

We use the f-FEM developed in §3 to study the geometrically nonlinear and
fractional-order nonlocal elastic response of isotropic beams. More
specifically, the effect of the fractional order $\alpha$ and of the length
scales $l_{A}$, $l_{B}$ on the nonlocal elastic behavior of the nonlinear beam
are studied. For all the numerical studies reported in this section, the
aspect ratio of the beam is chosen to be $L/h=100$ in order to satisfy the
slender beam assumption of the Euler-Bernoulli beam theory. The length of the
beam is maintained at $L=1$ m, the width is kept constant at unity, and the
elastic modulus is set at $E=3$ GPa. It is assumed here that that the length
scales $l_{A}$ and $l_{B}$ for the isotropic beam are equal to a constant
$l_{f}$ for points sufficiently within the domain of the beam. These length
scales are truncated appropriately for points close to the beam boundaries as
discussed in §2 (see Fig. (3)). Note that, although the presented f-FEM is
capable of simulating the static response of the fractional-order nonlocal
beam under both axial and transverse loads, here we present only the results
for transverse loads. Both the case of a uniformly distributed load per unit
length (UDL) (in N/m) applied along the length of the beam, and of a point
load (PtL) (in N) applied at a specified point. The transverse displacement,
normal stress $\sigma_{11}$ and applied load are reported as:

$\overline{w}_{0}=\frac{w_{0}}{h};~{}~{}~{}~{}~{}\overline{\sigma}_{11}=\sigma_{11}\times\left(\frac{h}{L}\right)^{2}\frac{1}{q_{0}};~{}~{}~{}~{}~{}\overline{q}=q_{0}\times\frac{L}{h};~{}~{}~{}~{}~{}\overline{P}=P\times\frac{L}{h}$
(45)

In the following, we first present the results of the validation and
convergence study and then we present the static response of the fractional-
order nonlocal beam for the various loading conditions. For all the numerical
results, the mesh discretization was chosen such that $l_{f}/l_{e}=10$, where
$l_{e}$ is the length of the discretized element. This choice allows for
sufficient number of elements to be included in the horizon of nonlocality at
any point hence guaranteeing to accurately capture the nonlocal interactions.

###  4.1 Validation

We validated the performance and evaluated the accuracy of the f-FEM by using
the following strategies: validation #1: the f-FEM model was solved for
$\alpha=1.0$, corresponding to the classical elastic solid, and the results
were compared with those available in the literature[53]; validation #2: the
f-FEM model was solved for non-integer values of $\alpha$ while ignoring the
nonlinear effects in Eq. (36), and the results were compared against [47]; and
validation #3: the f-FEM results were compared with the exact solution for a
nonlinear and fractional-order elastic clamped-clamped beam subjected to a
distributed transverse load. The validations #1 and #2 can rather be
considered as checks on the robustness of the nonlinear f-FEM. Specific
details of the three validation strategies are provided here below.

Validation #1: in this study, the fractional order is taken equal to
$\alpha=1.0$. This assumption reduces the nonlocal constitutive model (strain-
displacement relations) in Eq. (8) to their integer-order counterparts for
classical elasticity. The nonlinear f-FEM is solved for this choice of
$\alpha$, $l_{f}$ being arbitrary for the local elastic solid, and the
numerical results are compared with [53] in the Fig. (5(a)). This comparison
is performed for beams corresponding to two different boundary conditions,
namely, pinned-pinned and clamped-clamped. As evident from Fig. (5(a)), the
match between the transverse displacement obtained numerically from the f-FEM
and the results in the literature is excellent. The error between them is less
than $1\%$ for all cases.

Validation #2: in this study, the nonlinear effects are ignored for the sake
of a comparison with the literature over fractional-order nonlocal and linear
elastic response [47]. This simplification is achieved by ignoring the term
$D_{x_{1}}^{\alpha}[w_{0}(x_{1})]$ in the definition for stiffness matrices
given in Eq. (36), which contributes to the nonlinearity. The resulting
linearized f-FEM model was solved for $\alpha=0.9$ and $l_{f}=L/10$ and the
numerical results are compared with [47] in the Fig. (5(b)). Similar to
validation #1, this comparison is also performed for two different boundary
conditions. As evident from Fig. (5(b)), the match between the transverse
displacement obtained numerically from the f-FEM and the results available in
the literature is excellent and the relative error is less than $1\%$ for all
cases.

(a) Comparison with the nonlinear local continuum.

(b) Comparison with the linear nonlocal continnum.

Figure 5: Numerical nonlinear f-FEM results for (a) nonlinear classical
continuum $(\alpha=1.0)$ compared with [53], and (b) linear nonlocal continuum
($\alpha=0.9$, $l_{f}=L/10$) compared with [47].

Validation #3: in this study, the following displacement field at the mid-
plane of the beam is assumed:

$u_{0}(x_{1})=U_{0}\left(1-\frac{x_{1}}{L}\right)\left(\frac{x_{1}}{L}\right)$
(46a)
$w_{0}(x_{1})=W_{0}\left(1-\frac{x_{1}}{L}\right)^{2}\left(\frac{x_{1}}{L}\right)^{2}$
(46b)

Note that the above displacement fields satisfy the clamped-clamped boundary
conditions for the beam given in Eq. (18). Although, the above displacement
fields are independent of the fractional constitutive parameters $\alpha$ and
$l_{f}$, the axial and the transverse forces required to generate the above
elastic response of the beam would be functions of these parameters. These
distributed loads are obtained from the strong form of the governing equations
given in Eq. (18) using the constitutive relations stated in Eqs. (10) and
(19). Note that these expressions derived from the governing equations contain
the appropriate information regarding the geometry. Moreover, the amplitude of
the load is dependent on the assumed displacement field distribution in Eq.
(46). The closed-form expressions of the distributed axial and transverse
loads derived for the exact solutions in Eq. (46) are used in the f-FEM model
given in Eq. (35). This strategy is similar to the one presented in [47]. The
numerical results obtained from the f-FEM for the mid-plane transverse
displacement $w_{0}(x_{1})$ are compared with the exact solution. The results,
shown in Fig. (6), are in excellent agreement. The error between the
numerically obtained f-FEM solutions and the exact solutions is less than
$3\%$ for all cases. Similar inverse approaches for validation have also been
adopted in [56, 57] for different numerical methods.

(a) Case 1.

(b) Case 2.

Figure 6: Numerical f-FEM results (a) Case 1: $W_{0}=10h$ and $U_{0}=0.05h$,
and (b) Case 2: $W_{0}=20h$ and $U_{0}=0.1h$, corresponding to the assumed
displacement field in Eq. (46) compared with the exact solutions for the
clamped-clamped boundary condition. Note that, in this validation procedure,
the transverse displacement field is assumed independent of $\alpha$ and
$l_{f}$.

###  4.2 Convergence

In the following, we present the results of the sensitivity analysis on the FE
mesh size. The convergence of the integer-order FEM with element
discretization, referred to as the $h-$refinement, is well established in the
literature [53]. In this study, we noted that in addition to the FE mesh size,
the convergence of the f-FEM also depends on the fractional order $\alpha$, on
the relative size of $l_{e}$, and on the nonlocal length scale $l_{f}$. It
appears that the convergence of the f-FEM depends on the strength of the
nonlocal interactions across the nonlocal horizon. Therefore, sufficient
number of elements $N_{inf}(=l_{f}/l_{e})$ should be available in the
influence zone of the fractional-order (nonlocal) interactions. Additionally,
mesh refinement across the length of the 1D beam would increase the spatial
resolution and hence, result in lower inconsistencies due to the truncation of
the nonlocal horizon caused by the ceil and floor operations connected to Eq.
(43). This approach would provide a better numerical approximation for the
nonlocal matrices, $[\tilde{B}_{\square}]$ in Eqs. (42a)-(43). Therefore,
convergence of the f-FEM is expected when the number of elements in the
influence zone $N_{inf}$, referred to as the “dynamic rate of convergence”
[55, 47], is sufficient to accurately capture the nonlocal interactions.

(a) Clamped-Clamped beam.

(b) Pinned-Pinned beam.

Figure 7: Numerical convergence of the f-FEM numerical model ($\alpha=0.8$ and
$l_{f}=L/10$).

In the context of the current study, we establish the convergence of the f-FEM
for two different boundary conditions. The fractional-order parameters are
chosen as $\alpha=0.8$ and $l_{f}=L/10$. The results are presented in Fig.
(7). In both cases, Fig. (7(a)) for the clamped-clamped and Fig. (7(b)) for
the pinned-pinned beams, the dynamic rate of convergence is noted to be
$N_{inf}\geq 5$. Additionally, a comprehensive comparison of the convergence
of the f-FEM for various values of the fractional parameters $\alpha$ and
$l_{f}$, for a clamped-clamped beam, is presented in Table 1. We have
discussed in §4.3, the reason behind choosing a restricted range of
$\alpha\in[0.5,1]$. In Table 1, the convergence may be assessed by comparing
the values along a prescribed column on the far right, corresponding to the
transverse displacement for a chosen set of fractional constitutive
parameters. As evident from Table 1, for higher values of $l_{f}$ and lower
values of $\alpha$, that is for stronger fractional-order nonlocal
interactions, larger number of elements are necessary for the convergence of
f-FEM. Note also that $N_{inf}=10$ provides a satisfactory numerical
convergence yielding a difference smaller than $1\%$ between successive
refinements. Therefore, the f-FEM simulations presented below will be carried
out using this mesh discretization.

Fractional Horizon Length, $l_{f}$ |  $N_{inf}$ | fractional-order, $\alpha$
---|---|---
$\alpha=1.0$ | $\alpha=0.9$ | $\alpha=0.8$ | $\alpha=0.7$ | $\alpha=0.6$ | $\alpha=0.5$
$l_{f}=L/5$ | 2 | 0.7352 | 0.7911 | 0.8454 | 0.9004 | 0.9612 | 1.0367
5 | 0.7394 | 0.7819 | 0.8202 | 0.8569 | 0.8967 | 0.9488
10 | 0.7426 | 0.7821 | 0.8168 | 0.8494 | 0.8844 | 0.9322
20 | 0.7428 | 0.7810 | 0.8140 | 0.8449 | 0.8775 | 0.9224
$l_{f}=L/10$ | 2 | 0.7410 | 0.7821 | 0.8223 | 0.8636 | 0.9095 | 0.9667
5 | 0.7426 | 0.7653 | 0.7856 | 0.8054 | 0.8273 | 0.8556
10 | 0.7428 | 0.7606 | 0.7756 | 0.7899 | 0.8057 | 0.8264
20 | 0.7429 | 0.7583 | 0.7708 | 0.7826 | 0.7956 | 0.8128
$l_{f}=L/20$ | 2 | 0.7424 | 0.7807 | 0.8187 | 0.8577 | 0.8994 | 0.9472
5 | 0.7428 | 0.7597 | 0.7748 | 0.7895 | 0.8053 | 0.8244
10 | 0.7426 | 0.7538 | 0.7627 | 0.7710 | 0.7802 | 0.7920
20 | 0.7429 | 0.7510 | 0.7568 | 0.7622 | 0.7684 | 0.7769
Table 1: Convergence of the f-FEM for clamped-clamped beam with different
fractional parameters.

###  4.3 Effect of the fractional-order constitutive parameters

In this section, we use the f-FEM to analyze the static response of the
fractional-order nonlocal and geometrically nonlinear Euler-Bernoulli beam.
More specifically, we analyze the effect of the fractional parameters $\alpha$
and $l_{f}$ on the response of the beam subject to different loading and
boundary conditions.

We start by considering a beam that is clamped at both ends and is subject to
a UDL of magnitude $q_{0}$ (in N/m). The nonlinear load-displacement curves of
the clamped-clamped beam are obtained and compared for different values of
$\alpha$ and $l_{f}$. The results of this comparison are provided in Fig. (8).
The effect of the fractional order $\alpha$, with $l_{f}$ being held constant,
is studied in Fig. (8(a)). Recall that for $\alpha=1.0$, the differential and
algebraic form of the governing equations and boundary conditions, reduce to
those of a geometrically nonlinear classical (local) elasticity study.
Additionally, we also compare the effect of the horizon of nonlocality over
the elastic response of the beam in Fig. (8(b)). As evident from Fig. (8), the
transverse displacement increases with increasing degree of nonlocality, that
is by reducing the value of $\alpha$ or by increasing the horizon length
$l_{f}/L$. These trends point to a reduction in the stiffness of the
fractional-order nonlocal structure. As expected, a convergence of the
fractional-order model to the local elastic response is noted for $\alpha$
close to $1.0$ and $l_{f}/L<<1$.

For the clamped-clamped beam described above, we also analyzed the effect of
the fractional constitutive parameters on the axial normal stress across the
thickness of the beam at its mid-length, that is at $x_{1}=L/2$. The results
are presented in Fig. (9) in terms of the non-dimensional stress following the
application of two different transverse loads. As expected of nonlinear
response, an increase in magnitude of the transverse load (for $\bar{q}=10)$
produces a non-zero axial stress at the mid-plane (see Fig. (9(b))) of the
beam. This axial stress caused by the transverse loading is a direct result of
the axial-transverse coupling induced by the geometric nonlinearity. For lower
magnitudes of the transverse load $\overline{q}=2$, a weak nonlinearity is
clear from the lower value of $\bar{\sigma}_{11}$ at the mid-plane (see Fig.
(9(a))). Note that, given the linear nature of the constitutive relations in
Eq. (10) used for the current analysis, the observations noted above for the
stress in the beam, directly extend to the fractional-order strain.

Before proceeding further, we briefly discuss how the degree of nonlocality
increases by reducing the order $\alpha$ and increasing the length-scale
$l_{f}$. Note that the convolution kernel used within the definition of the
fractional-order derivative is a power-law kernel which is monotonically
decreasing in nature. It appears that, by reducing the value of the order
$\alpha$, the magnitude of this kernel increases. More specifically, for a
fixed point $x$ interacting with a specific point $s$ in its horizon of
nonlocality,
$\left(\frac{1}{|x-s|}\right)^{\alpha_{1}}>\left(\frac{1}{|x-s|}\right)^{\alpha_{2}}\forall~{}\alpha_{1}<\alpha_{2}$.
Since the magnitude of the kernel increases, the strength of the corresponding
nonlocal interactions increases and so does the degree of nonlocality.
Similarly, by increasing the length scale $l_{f}$, the size of the horizon of
nonlocality increases. Hence, information corresponding to the nonlocal
interactions between a larger number of points within the solid are accounted
for in the formulation. Thus the degree of nonlocality increases with
increasing values of $l_{f}$.

(a) $\overline{w}$ vs $\alpha$.

(b) $\overline{w}$ vs $l_{f}$.

Figure 8: Nonlocal elastic transverse displacement of the clamped-clamped beam
compared for different values of the fractional-order constitutive properties.

(a) $\overline{q}=2$.

(b) $\overline{q}=10$.

Figure 9: Normal axial stress $\tilde{\sigma}_{11}$ in the nonlocal beam
subject to clamped-clamped boundary conditions and two cases of the transverse
load compared for different fractional-order constitutive properties.

We extend the above study to analyze a beam with both ends pinned. The effect
of the fractional model parameters over the nonlocal elastic response of the
geometrically nonlinear pinned-pinned beam are shown in Fig. (10). From these
figures, it may be noted that the transverse displacement increases while
reducing the value of $\alpha$ and increasing the nonlocal horizon length
$l_{f}$. These observations for the pinned-pinned beam agree very well with
those noted previously for the clamped-clamped beam. The consistency of these
observations regardless of the nature of the boundary conditions indicates the
robustness of the fractional order model and marks a net departure from the
paradoxical results obtained in the literature for more traditional integral
beam theories [27, 28, 58]. These inaccurate predictions have been traced back
to the ill-posed nature of the nonlocal constitutive relations as discussed in
detail in [29]. On the contrary, the fractional-order model of nonlocality is
free from mathematical ill-posedness and leads to consistent predictions of
softening for the nonlocal beam structures [47].

(a) $\overline{w}$ vs $\alpha$.

(b) $\overline{w}$ vs $l_{f}$.

Figure 10: Comparison of the nonlocal elastic transverse displacement of the
pinned-pinned beam for different values of the fractional-order constitutive
properties.

Finally, the nonlocal response of the clamped-clamped and pinned-pinned beams
are compared for different loading conditions and the results are shown in
Fig. (11). The linear and nonlinear load-displacement curves of the Euler-
Bernoulli beam subject to an UDL and with fractional constitutive parameters
fixed at $\alpha=0.8$ and $l_{f}=L/10$ are compared in Fig. (11(a)). A similar
comparison for the beam subject to a PtL applied at mid-length is shown in
Fig. (11(b)). It appears that the effect of the geometric nonlinearity is
significant for the nonlocal pinned-pinned beams, irrespective of the loading
conditions. Note that this is in agreement with the observations from the
classical theory of elasticity [53].

(a) UDL, $q_{0}$ per unit length.

(b) PtL, $P$.

Figure 11: Geometric non-linear effects compared for loading and boundary
conditions ($\alpha=0.8$, $l_{f}=L/10$).

An important remark should be made concerning the physically acceptable range
for the order $\alpha$. Although a consistent softening is observed for
decreasing value of $\alpha$, we highlight here that for very low values of
$\alpha$, non-physical results were noted in literature[40, 59]. Note that the
fractional-order kinematic relations of order $\alpha\in(0,1)$ define the
strain as a fractional-order interpolation of the zeroth and first-order
derivatives of the displacement. Hence, there exists a limit on the order of
the fractional derivative beyond which this fractional model of nonlocality
breaks down. In the nonlinear study, we note that this specific value is
$\alpha\approx 0.4$. It was also observed that this critical value for
$\alpha$ depends on other parameters including the relative ratio of
$l_{f}/L$, where $l_{f}$ is the nonlocal length scale and $L$ is the geometric
length of the structure. However, we emphasize that this breakdown is not a
characteristic of the f-FEM technique employed for the numerical analysis, but
instead a physical limit of the fractional-order modeling[47].

##  5 Conclusions

The current study laid the foundation for the analysis of geometrically
nonlinear and nonlocal solids via fractional-order operators while providing
accurate numerical schemes for their solution. More specifically, the
geometrically nonlinear formulation for fractional-order nonlocal Euler-
Bernoulli beams was developed and its solution was achieved via a dedicated
finite-element-based numerical method. The general fractional-order continuum
theory was extended to develop a frame-invariant and dimensionally consistent
nonlinear fractional-order strain-displacement theoretical framework. This
theory is well suited to capture scale effects, nonlocality, and heterogeneity
in both complex solids and interfaces undergoing large deformations. The
constitutive model for the nonlinear response of fractional-order beams was
derived following variational principles. The 1D beam model lays the
foundation for a more general nonlinear analysis of fractional-order nonlocal
solids. Owing to the complexity of the nonlocal models, a fractional-order FEM
(f-FEM) was also developed by means of potential energy minimization. This
nonlinear f-FEM provides an efficient and accurate computational approach to
obtain the numerical solution of nonlinear integro-differential governing
equations with a singular kernel. The robustness of this numerical approach
for the solution of nonlocal elasticity problems was demonstrated by showing
the consistency of the predicted nonlocal response subject to different loads
and boundary conditions. This latter aspect is in contrast with other
traditional approaches to nonlocal elasticity that, occasionally, provide
paradoxical results depending on the selection of both loading and boundary
conditions. Finally we note that, although the f-FEM was used in this study to
solve nonlinear beam models, the same computational framework can be easily
extended either to the simulation of systems with higher dimensions or to
other nonlinear physical processes described by fractional differ-integral
equations.

## Appendix 1

Frame invariance of the fractional-order kinematic relations:  To establish
frame invariance of the formulation, we begin with the definition given in Eq.
(3) for the deformation gradient tensor following fractional-order kinematics.
The frame-invariance of each of the individual terms in the definition would
establish the frame invariance of the complete tensor
$\tilde{\textbf{F}}_{X}$. For this purpose, we consider a rigid-body motion
superimposed on a general point X (see Fig. (1)) of the reference
configuration of the body as:

$\bm{\Psi}(\textbf{X})=\textbf{c}+\textbf{Q}\textbf{X},$ (47)

where Q, a proper orthogonal tensor, and c, a spatially constant term,
represent the rotation and translation of the reference frame, respectively.
Under this rigid-body motion, the fractional deformation gradient tensor
should transform as $\tilde{\textbf{F}}^{\Psi}_{X}=\textbf{Q}$ (similar to the
classical continuum elasticity) such that $\tilde{\textbf{F}}^{\Psi
T}_{X}\tilde{\textbf{F}}^{\Psi}_{X}=\textbf{I}$. This would result in trivial
strain measures for the rigid body transformation. The above rigid body motion
employed in the definition of $\tilde{\textbf{F}}_{X}$ given in Eq. (2a)
results in:

$\begin{split}\tilde{\textbf{F}}_{X_{ij}}^{\Psi}=\frac{1}{2}\Gamma(2-\alpha)\biggl{[}\frac{L_{A_{j}}^{\alpha-1}}{\Gamma(1-\alpha)}\int_{X_{A_{j}}}^{X_{j}}\frac{D^{1}_{R_{j}}\Psi_{i}(\textbf{R})}{(X_{j}-R_{j})^{\alpha}}\mathrm{d}R_{j}+\frac{L_{B_{j}}^{\alpha-1}}{\Gamma(1-\alpha)}\int_{X_{j}}^{X_{B_{j}}}\frac{D^{1}_{R_{j}}\Psi_{i}(\textbf{R})}{(R_{j}-X_{j})^{\alpha}}\mathrm{d}R_{j}\biggr{]}\end{split}$
(48)

where R is a dummy vector representing the spatial variable of a point within
the nonlocal horizon of X. For the rigid body motion given in Eq. (47), the
first integer-order derivative $D^{1}_{R_{j}}\Psi_{i}(\textbf{R})$ is

$D^{1}_{R_{j}}\Psi_{i}(\textbf{R})=Q_{ij}$ (49)

As discussed in §2.1, the length-scales $L_{A_{j}}$ and $L_{B_{j}}$ in Eq.
(48) are given by $L_{A_{j}}=X_{j}-X_{A_{j}}$ and $L_{B_{j}}=X_{B_{j}}-X_{j}$.
Using these relations and recurrence definition for the Gamma function
$\Gamma(2-\alpha)$, it immediately follows that the deformation gradient
tensor is $\tilde{\textbf{F}}_{X_{ij}}^{\Psi}=Q_{ij}$ identically. As
discussed previously, this results in null strain measures for rigid body
translation and rotations using the fractional-order kinematic relations.
Repeating the above procedure, similar arguments can be made for the frame
invariance of $\tilde{\textbf{F}}^{\Psi}_{x}$ and subsequently of
$\overset{\alpha}{{\textbf{F}}}$. This establishes the frame-invariance of the
model. Note that the above proof holds true independent of the values of
$L_{A_{j}}$ and $L_{B_{j}}$, and for a more general case $L_{A_{j}}\neq
L_{B_{j}}$, which occurs in the presence of asymmetric horizons noted at
points close to material boundaries and interfaces.

## Appendix 2

Derivation of governing equations: In the following, we provide the detailed
steps for obtaining the strong form of the governing differential equations
given in Eq. (18). The governing equations are obtained by minimizing the
potential energy given in Eq. (17). The minimization of the potential energy
results in two independent integrals which are given in Eq. (33) because of
the independent nature of the variations $\delta u_{0}(x_{1})$ and $\delta
w_{0}(x_{1})$. It follows from Eq. (33a) that:

$\int_{0}^{L}N(x_{1})D_{x_{1}}^{\alpha}\left[\delta
u_{0}(x_{1})\right]\mathrm{d}x_{1}=-\int_{0}^{L}{}^{R-RL}_{x_{1}-l_{B}}D_{x_{1}+l_{A}}^{\alpha}\left[N(x_{1})\right]\delta
u_{0}(x_{1})\mathrm{d}x_{1}+\left[{}_{x_{1}-l_{B}}I^{1-\alpha}_{x_{1}+l_{A}}\left[N(x_{1})\right]\delta
u_{0}\right]\big{|}_{0}^{L}$ (50a) In the above equation, the Riesz fractional
integral of $N(x_{1})$ is given as:
${}_{x_{1}-l_{B}}I^{1-\alpha}_{x_{1}+l_{A}}\left[N(x_{1})\right]=\frac{\Gamma(2-\alpha)}{2}\left(l_{B}^{\alpha-1}{}_{x_{1}-l_{B}}I_{x_{1}}^{1-\alpha}\left[N(x_{1})\right]-l_{A}^{\alpha-1}{}_{x_{1}}I_{x_{1}+l_{A}}^{1-\alpha}\left[N(x_{1})\right]\right)$
(50b)

where ${}_{x_{1}-l_{B}}I_{x_{1}}^{1-\alpha}$ and
${}_{x_{1}}I_{x_{1}+l_{A}}^{1-\alpha}$ are the left and right Riesz integrals
[60]. The simplification of Eq. (33) to Eq. (50a) is performed using
integration by parts and the definitions for the fractional derivatives. The
specific steps leading to the above simplification can be found in [47]. The
above discussed simplification extends directly to other integrals in Eq.
(33). Combining all the simplifications yields the first variation of the
potential energy as:

$\begin{split}\delta\Pi&=-\int_{0}^{L}\underbrace{\left[{}^{R-RL}_{x_{1}-l_{B}}D^{\alpha}_{x_{1}+l_{A}}[N(x_{1})]+F_{a}(x_{1})\right]}_{\text{Axial
Governing Equation}}\delta
u_{0}(x_{1})\mathrm{d}x_{1}+\underbrace{\left[{}_{x_{1}-l_{B}}I^{1-\alpha}_{x_{1}+l_{A}}[N(x_{1})]\delta
u_{0}(x_{1})\right]\Big{|}_{0}^{L}}_{\text{Axial Boundary Conditions}}\\\
&-\int_{0}^{L}\underbrace{\left[D^{1}_{x_{1}}\left({}^{R-RL}_{x_{1}-l_{B}}D^{\alpha}_{x_{1}+l_{A}}[M(x_{1})]\right)+{}^{R-RL}_{x_{1}-l_{B}}D^{\alpha}_{x_{1}+l_{A}}\left[N(x_{1})D_{x_{1}}^{\alpha}[w_{0}(x_{1})]\right]+F_{t}(x_{1})\right]}_{\text{Transverse
Governing Equation}}\delta w_{0}(x_{1})\mathrm{d}x_{1}\\\
&-\underbrace{\left[{}_{x_{1}-l_{B}}I^{1-\alpha}_{x_{1}+l_{A}}[M(x_{1})]\delta
D^{1}_{x_{1}}w_{0}(x_{1})-\left\\{{}_{x_{1}-l_{B}}^{R-RL}D^{\alpha}_{x_{1}+l_{A}}[M(x_{1})]+{}_{x_{1}-l_{B}}I^{1-\alpha}_{x_{1}+l_{A}}[N(x_{1})]\right\\}\delta
w_{0}(x_{1})\right]\Big{|}_{0}^{L}}_{\text{Transverse Boundary
Conditions}}\end{split}$ (51)

where $D^{1}_{x_{1}}(\cdot)$ denotes the first integer-order derivative with
respect to $x_{1}$. The governing differential equations and the associated
boundary conditions given in Eq. (18) are now obtained from the above Eq. (51)
by using first principle of variational calculus. Note that, for longitudinal
boundaries of the beam, i.e., $x_{1}=\\{0,L\\}$, $l_{A}$ and $l_{B}$ are equal
to zero in the Caputo and RL derivatives defined in Eqs. (4,20), respectively,
and the Riesz integrals given in Eq. (50b). These limit cases reduces the
boundary conditions to the integer-order counterparts as given in Eq. (18)
[47].

##  6 Acknowledgements

The following work was supported by the Defense Advanced Research Project
Agency (DARPA) under the grant #D19AP00052, and the National Science
Foundation (NSF) under the grant DCSD #1825837. The content and information
presented in this manuscript do not necessarily reflect the position or the
policy of the government. The material is approved for public release;
distribution is unlimited.

## References

  * [1] T. L. Szabo.  Time domain wave equations for lossy media obeying a frequency power law.  The Journal of the Acoustical Society of America, 96(1):491–500, 1994.
  * [2] K. A. Wear.  A stratified model to predict dispersion in trabecular bone.  IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 48(4):1079–1083, 2001.
  * [3] Z. E. A. Fellah, S. Berger, W. Lauriks, and C. Depollier.  Verification of kramers–kronig relationship in porous materials having a rigid frame.  Journal of sound and vibration, 270(4-5):865–885, 2004.
  * [4] A. Stulov and V. I. Erofeev.  Frequency-dependent attenuation and phase velocity dispersion of an acoustic wave propagating in the media with damages.  In Generalized Continua as Models for Classical and Advanced Materials, pages 413–423. Springer, 2016.
  * [5] E Kröner.  Elasticity theory of materials with long range cohesive forces.  International Journal of Solids and Structures, 3(5):731–742, 1967\.
  * [6] A. C. Eringen and D. G. B. Edelen.  On nonlocal elasticity.  International Journal of Engineering Science, 10(3):233–248, 1972\.
  * [7] W. Sumelka, T. Blaszczyk, and C. Liebold.  Fractional euler–bernoulli beams: Theory, numerical study and experimental validation.  European Journal of Mechanics-A/Solids, 54:243–251, 2015.
  * [8] Ö. Civalek and C. Demir.  A simple mathematical model of microtubules surrounded by an elastic matrix by nonlocal finite element method.  Applied Mathematics and Computation, 289:335–352, 2016.
  * [9] Z. Rahimi, W. Sumelka, and X-J. Yang.  Linear and non-linear free vibration of nano beams based on a new fractional non-local theory.  Engineering Computations, 34(5):1754–1770, 2017.
  * [10] R. H. J. Peerlings, M. G. D. Geers, R. De Borst, and W. A. M. Brekelmans.  A critical comparison of nonlocal and gradient-enhanced softening continua.  International Journal of solids and Structures, 38(44-45):7723–7746, 2001.
  * [11] E. C. Aifantis.  Update on a class of gradient theories.  Mechanics of materials, 35(3-6):259–280, 2003.
  * [12] Suman Guha, Sandeep Sangal, and Sumit Basu.  A review of higher order strain gradient theories of plasticity: Origins, thermodynamics and connections with dislocation mechanics.  Sadhana, 40(4):1205–1240, 2015.
  * [13] S. Sidhardh and M. C. Ray.  Element-free galerkin model of nano-beams considering strain gradient elasticity.  Acta Mechanica, 229(7):2765–2786, 2018.
  * [14] S. Sidhardh and M. C. Ray.  Size-dependent elastic response in functionally graded microbeams considering generalized first strain gradient elasticity.  The Quarterly Journal of Mechanics and Applied Mathematics, 2019\.
  * [15] C. Polizzotto.  Nonlocal elasticity and related variational principles.  International Journal of Solids and Structures, 38(42-43):7359–7380, 2001.
  * [16] Z. P. Bažant and M. Jirásek.  Nonlocal integral formulations of plasticity and damage: survey of progress.  Journal of Engineering Mechanics, 128(11):1119–1149, 2002.
  * [17] S. Sidhardh and M. C. Ray.  Effect of nonlocal elasticity on the performance of a flexoelectric layer as a distributed actuator of nanobeams.  International Journal of Mechanics and Materials in Design, 14(2):297–311, 2018.
  * [18] J. Yang, L. L. Ke, and S. Kitipornchai.  Nonlinear free vibration of single-walled carbon nanotubes using nonlocal timoshenko beam theory.  Physica E: Low-dimensional Systems and Nanostructures, 42(5):1727–1735, 2010.
  * [19] S. Srividhya, P. Raghu, A. Rajagopal, and J. N. Reddy.  Nonlocal nonlinear analysis of functionally graded plates using third-order shear deformation theory.  International Journal of Engineering Science, 125:1–22, 2018.
  * [20] J. N. Reddy.  Nonlocal nonlinear formulations for bending of classical and shear deformation theories of beams and plates.  International Journal of Engineering Science, 48(11):1507–1518, 2010.
  * [21] S. A. Emam.  A general nonlocal nonlinear model for buckling of nanobeams.  Applied Mathematical Modelling, 37(10-11):6929–6939, 2013.
  * [22] M. Lembo.  On nonlinear deformations of nonlocal elastic rods.  International Journal of Solids and Structures, 90:215–227, 2016\.
  * [23] A. R. Srinivasa and J. N. Reddy.  A model for a constrained, finitely deforming, elastic solid with rotation gradient dependent strain energy, and its specialization to von kármán plates and beams.  Journal of the Mechanics and Physics of Solids, 61(3):873–885, 2013\.
  * [24] J. N. Reddy and S. El-Borgi.  Eringen’s nonlocal theories of beams accounting for moderate rotations.  International Journal of Engineering Science, 82:159–177, 2014\.
  * [25] F. Najar, S. El-Borgi, J. N. Reddy, and K. Mrabet.  Nonlinear nonlocal analysis of electrostatic nanoactuators.  Composite Structures, 120:117–128, 2015.
  * [26] J. N. Reddy.  Nonlocal theories for bending, buckling and vibration of beams.  International Journal of Engineering Science, 45(2-8):288–307, 2007\.
  * [27] N. Challamel and C. M. Wang.  The small length scale effect for a non-local cantilever beam: a paradox solved.  Nanotechnology, 19(34):345703, 2008.
  * [28] P. Khodabakhshi and J. N. Reddy.  A unified integro-differential nonlocal model.  International Journal of Engineering Science, 95:60–75, 2015.
  * [29] Giovanni Romano, Raffaele Barretta, Marina Diaco, and Francesco Marotti de Sciarra.  Constitutive boundary conditions and paradoxes in nonlocal elastic nanobeams.  International Journal of Mechanical Sciences, 121:151–156, 2017\.
  * [30] R. Ansari, R. Gholami, and H. Rouhi.  Size-dependent nonlinear forced vibration analysis of magneto-electro-thermo-elastic timoshenko nanobeams based upon the nonlocal elasticity theory.  Composite Structures, 126:216–226, 2015.
  * [31] R. L. Bagley and P. J. Torvik.  A theoretical basis for the application of fractional calculus to viscoelasticity.  Journal of Rheology, 27(3):201–210, 1983.
  * [32] A. Chatterjee.  Statistical origins of fractional derivatives in viscoelasticity.  Journal of Sound and Vibration, 284(3-5):1239–1245, 2005.
  * [33] Sansit Patnaik and Fabio Semperlotti.  Application of variable-and distributed-order fractional operators to the dynamic analysis of nonlinear oscillators.  Nonlinear Dynamics, pages 1–20, 2020.
  * [34] F. Mainardi.  Fractional relaxation-oscillation and fractional diffusion-wave phenomena.  Chaos, Solitons & Fractals, 7(9):1461–1477, 1996.
  * [35] S. Buonocore, M. Sen, and F. Semperlotti.  Occurrence of anomalous diffusion and non-local response in highly-scattering acoustic periodic media.  New Journal of Physics, 2019.
  * [36] John P Hollkamp, Mihir Sen, and Fabio Semperlotti.  Analysis of dispersion and propagation properties in a periodic rod using a space-fractional wave equation.  Journal of Sound and Vibration, 441:204–220, 2019.
  * [37] J. P. Hollkamp, M. Sen, and F. Semperlotti.  Model-order reduction of lumped parameter systems via fractional calculus.  Journal of Sound and Vibration, 419:526–543, 2018.
  * [38] J. P. Hollkamp and F. Semperlotti.  Application of fractional order operators to the simulation of ducts with acoustic black hole terminations.  Journal of Sound and Vibration, 465:115035, 2020.
  * [39] A. Carpinteri, P. Cornetti, and A. Sapora.  A fractional calculus approach to nonlocal elasticity.  The European Physical Journal Special Topics, 193(1):193, 2011.
  * [40] W. Sumelka.  Thermoelasticity in the framework of the fractional continuum mechanics.  Journal of Thermal Stresses, 37(6):678–706, 2014.
  * [41] W. Sumelka.  Fractional calculus for continuum mechanics–anisotropic non-locality.  Bulletin of the Polish Academy of Sciences Technical Sciences, 64(2):361–372, 2016.
  * [42] Nick Laskin and G Zaslavsky.  Nonlinear fractional dynamics on a lattice with long range interactions.  Physica A: Statistical Mechanics and its Applications, 368(1):38–54, 2006.
  * [43] M. Di Paola and M. Zingales.  Long-range cohesive interactions of non-local continuum faced by fractional calculus.  International Journal of Solids and Structures, 45(21):5642–5659, 2008.
  * [44] G. Cottone, M. Di Paola, and M. Zingales.  Fractional mechanical model for the dynamics of non-local continuum.  In Advances in numerical methods, pages 389–423. Springer, 2009\.
  * [45] Gioacchino Alotta, Giuseppe Failla, and Massimiliano Zingales.  Finite element method for a nonlocal timoshenko beam model.  Finite Elements in Analysis and Design, 89:77–92, 2014.
  * [46] G. Alotta, G. Failla, and M. Zingales.  Finite-element formulation of a nonlocal hereditary fractional-order timoshenko beam.  Journal of Engineering Mechanics, 143(5):D4015001, 2017.
  * [47] Sansit Patnaik, Sai Sidhardh, and Fabio Semperlotti.  A Ritz-based finite element method for a fractional-order boundary value problem of nonlocal elasticity.  International Journal of Solids and Structures, 2020.
  * [48] Krzysztof Szajek, Wojciech Sumelka, Tomasz Blaszczyk, and Krzysztof Bekus.  On selected aspects of space-fractional continuum mechanics model approximation.  International Journal of Mechanical Sciences, 167:105287, 2020.
  * [49] Sansit Patnaik, Sai Sidhardh, and Fabio Semperlotti.  Fractional-order models for the static and dynamic analysis of nonlocal plates.  arXiv preprint arXiv:2002.10244, 2020.
  * [50] Sai Sidhardh, Sansit Patnaik, and Fabio Semperlotti.  Thermoelastic response of fractional-order nonlocal and geometrically nonlinear beams.  arXiv preprint arXiv:2003.10215, 2020.
  * [51] Y. Zheng, C. Li, and Z. Zhao.  A note on the finite element method for the space-fractional advection diffusion equation.  Computers & Mathematics with Applications, 59(5):1718–1726, 2010\.
  * [52] B. Jin, R. Lazarov, and Z. Zhou.  A petrov–galerkin finite element method for fractional convection-diffusion equations.  SIAM Journal on Numerical Analysis, 54(1):481–503, 2016.
  * [53] J. N. Reddy.  An Introduction to Nonlinear Finite Element Analysis: with applications to heat transfer, fluid mechanics, and solid mechanics.  OUP Oxford, 2014.
  * [54] A. A. Pisano, A. Sofi, and P. Fuschi.  Nonlocal integral elasticity: 2d finite element based solutions.  International Journal of Solids and Structures, 46(21):3836–3849, 2009.
  * [55] A. Norouzzadeh and R. Ansari.  Finite element analysis of nano-scale timoshenko beams using the integral model of nonlocal elasticity.  Physica E: Low-dimensional Systems and Nanostructures, 88:194–200, 2017.
  * [56] Pinghui Zhuang, Fawang Liu, Vo Anh, and Ian Turner.  Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term.  SIAM Journal on Numerical Analysis, 47(3):1760–1781, 2009.
  * [57] H Zhang, Fawang Liu, Mantha S Phanikumar, and Mark M Meerschaert.  A novel numerical method for the time variable fractional order mobile–immobile advection–dispersion model.  Computers & Mathematics with Applications, 66(5):693–701, 2013\.
  * [58] J. Fernández-Sáez, R. Zaera, J. A. Loya, and J. N. Reddy.  Bending of euler–bernoulli beams using eringen’s integral formulation: a paradox resolved.  International Journal of Engineering Science, 99:107–116, 2016\.
  * [59] W Sumelka.  Non-local kirchhoff–love plates in terms of fractional calculus.  Archives of Civil and Mechanical Engineering, 15(1):231–242, 2015\.
  * [60] I. Podlubny.  Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, volume 198.  Elsevier, 1998.