# Displacement-Driven Approach to Nonlocal Elasticity

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

###### Abstract

This study presents a physically consistent displacement-driven reformulation
of the concept of action-at-a-distance, which is at the foundation of nonlocal
elasticity. In contrast to existing approaches that adopts an integral stress-
strain constitutive relation, the displacement-driven approach is predicated
on an integral strain-displacement relation. The most remarkable consequence
of this reformulation is that the (total) strain energy is guaranteed to be
convex and positive-definite without imposing any constraint on the symmetry
of the kernels. This feature is critical to enable the application of nonlocal
formulations to general continua exhibiting asymmetric interactions;
ultimately a manifestation of material heterogeneity. Remarkably, the proposed
approach also enables a strong satisfaction of the locality recovery condition
and of the laws of thermodynamics, which are not foregone conclusions in most
classical nonlocal elasticity theories. Additionally, the formulation is
frame-invariant and the nonlocal operator remains physically consistent at
boundaries. The study is complemented by a detailed analysis of the dynamic
response of the nonlocal continuum and of its intrinsic dispersion leading to
the consideration that the choice of nonlocal kernels should depend on the
specific material. Examples of exponential or power-law kernels are presented
in order to demonstrate the applicability of the method to different classes
of nonlocal media. The ability to admit generalized kernels reinforces the
generalized nature of the displacement-driven approach over existing integral
methodologies, which typically lead to simplified differential models based on
exponential kernels. The theoretical formulation is also leveraged to simulate
the static response of nonlocal beams and plates illustrating the intrinsic
consistency of the approach, which is free from unwanted boundary effects.

Keywords— Nonlocal elasticity, Integral methods, Displacement-driven approach,
Energy methods, Nonlocal beams, Nonlocal plates
All correspondence should be addressed to<EMAIL_ADDRESS>or
<EMAIL_ADDRESS>

Highlights

  * •

Frame-invariant displacement-driven approach to nonlocal elasticity.

  * •

Locality recovery condition, boundary effects and thermodynamic consistency
analyzed.

  * •

Well-posed nonlocal governing equations derived from a convex positive-
definite system.

  * •

Dispersion relations and material-informed kernels.

  * •

Application to the static response of nonlocal beams and plates.

##  1 Introduction

Experimental and theoretical investigations have shown that, irrespective of
the spatial scale, size-dependent or nonlocal effects can become prominent in
the response of several structures in the most diverse fields of engineering,
nanotechnology, biotechnology, and even medicine [1, 2, 3, 4, 5, 6, 7, 8]. The
study and analysis of nonlocal effects has increasingly become an important
component of structural design and analysis. Specific examples of size-
dependent structures include sensors, biological implants, micro/nano-
electromechanical devices and even macroscale structures such as stiffened
panels in naval and aerospace applications [6, 5, 7, 8, 4]. In the case of
micro- and nano-structures, for example those involving carbon nanotubes or
graphene sheets, nonlocal effects can be traced back to interatomic
interactions, and even medium heterogeneity [8, 9, 10, 11, 12, 13]. These
nonlocal forces operate at the nanoscale and hence do not significantly
contribute to the overall macroscopic response. For these reasons, nonlocal
effects have long been considered negligible at the macro scales. However, at
the macro scales, nonlocal effects can originate from material heterogeneities
[2, 3, 14, 15, 5, 4, 16, 17] or even from intentionally nonlocal designs [7,
6]. While long-range interactions occur naturally in heterogeneous materials
(e.g. interactions between the layers of a functionally graded material), they
can also be induced artificially via carefully designed structural links in
the intentional nonlocal designs. The existence of size-dependent effects
leads to a softening or stiffening of the structure when compared to the
predictions made from classical (local) continuum mechanics relations.
Additionally, size-dependent structures exhibit anomalous wavenumber-frequency
dispersion characteristics leading to wavenumber/frequency dependent wave
speeds, contrary to constant wave speeds achieved via local descriptions [18,
19, 9, 20, 14, 21]. Consequently, the ability to accurately model size-
dependent effects has profound implications for the diverse applications in
the aforementioned examples. The inability of the classical continuum theory
in capturing size-dependent effects prevented its use in these types of
applications and led to the development of nonlocal continuum theories.

Seminal studies by Kroner [22], Edelen et al. [23], and Eringen et al. [24, 1]
laid the theoretical foundation of nonlocal elasticity, and explored its role
in the modeling of nonlocal size-dependent structures. The mathematical
description of the nonlocal continuum theory proposed in these seminal studies
relied on the introduction of additional contributions, resulting from long-
range nonlocal interactions, in terms of a convolution integral of the strain
field in the constitutive equations. Over the years, several researchers [16,
25, 26, 27, 28, 29, 30, 31, 32, 33] have proposed different modifications to
the formulation presented in the seminal studies [22, 23, 24]. From a
mathematical standpoint, these approaches belong to a class of the so-called
strong integral methods that capture nonlocal effects by re-defining the
stress-strain constitutive law in the form of a convolution integral of either
the strain or the stress field over a certain spatial domain (the so-called
horizon of nonlocality). Depending on whether the nonlocal contributions are
modeled using the strain or the stress fields, the integral methods can be
classified as strain-driven [24, 25, 26, 27] or stress-driven approaches [33].
Further discussions addressing the origin of nonlocal effects, existing
theories of nonlocal elasticity, and their applications can be found in this
recent review study [34].

Although the strain- and stress-driven integral approaches to nonlocal
elasticity have been able to address a multitude of aspects typical of the
response of size-dependent nonlocal structures, some important challenges
still remain open. In both these approaches, the nonlocal continuum
formulation is based on a nonlocal definition of the stress-strain
constitutive relation that adopts a local definition of the strain-
displacement kinematic relation. In other terms, the effect of nonlocality is
accounted only via the constitutive stress-strain relations [34]. Indeed, we
will show that it is this assumption that is at the root of several of the
challenges faced by existing approaches. These challenges are summarized in
the following:

  * •

The strain-driven approach leads to a Fredholm integral of the first-kind,
which results in ill-posed boundary-value problems. Consequently, its
application to model nonlocal structures of practical engineering interest
(e.g. beams) leads to inconsistent (also called ’paradoxical’) predictions
when the same analysis is performed under different loading conditions [35,
36, 33, 37, 31]. The issue of an ill-posed formulation was addressed by
adopting either stress-driven approaches or enhanced mixed-phase (local and
nonlocal) approaches, which led to nonlocal constitutive boundary conditions
[33, 37, 31]. Mathematically speaking, these nonlocal constitutive boundary
conditions are additional constraints required to guarantee well-posedness of
the integral equations. Notably, the stress-driven approach has yet to achieve
a practical numerical implementation (e.g. finite element method) and
currently relies on analytical methods to simulate the nonlocal response.
While the analytical approach does provide interesting insights, it prevents
the application of the stress-driven theory to more complex scenarios
involving, as an example, higher dimensional structures subject to general
loading conditions or nonlinear behavior.

  * •

Further, the most direct result of either the strain- or the stress-driven
constitutive formulations consists in the inability to yield a positive-
definite and convex strain energy density. This follows from the fact that the
strain energy density consists of the product between a function and its
convolution with a smoothing kernel, which does not yield a quadratic form.
Consequently, the only option to achieve a positive-definite strain energy is
to require a symmetric nonlocal kernel [25, 33]. While possible, such a
condition restricts the application of the resulting theory to structures
exhibiting asymmetric long-range interactions [16, 29, 15, 10, 6, 17]. This
aspect was also very recently discussed in [38] which highlighted the misuse
of Eringen’s stress-strain nonlocal constitutive relation to model structures
made from functionally graded materials.

  * •

In addition to restrictions on the symmetry of the convolution kernel, several
studies also restrict functional nature of the kernel. More specifically, most
strain- and stress-driven integral approaches reduce to simplified
differential models using the special exponential kernel [8, 39, 9, 33, 13].
The use of an exponential kernel while being a plausible option, prevents the
application of the resulting theory to a wide class of structures,
particularly those exhibiting anomalous attenuation and dispersion (as we will
show in §6). Examples include highly scattering and multifractal media [18,
19, 40, 41, 20, 17] and structures with intentionally designed nonlocal
connections [7, 6].

  * •

The use of a nonlocal stress-strain constitutive relation violates the
locality recovery condition which states that, for a physically consistent
nonlocal formulation, a uniform stress field and a length-scale independent
strain energy must be recovered from a uniform strain field [27]. Notably,
this limitation was addressed very recently in [31] via an enhanced mixed-
phase model, albeit in a weaker form that satisfies only the local stress
recovery condition.

  * •

As noted in [26], the use of a nonlocal stress-strain constitutive behavior is
not robust to boundary and surface effects in finite domains. This is a direct
result of an inconsistent (and incomplete) truncation of the nonlocal kernel
(often, exponential) at boundaries.

Some of the above mentioned challenges are also described in [31]. In this
study, we attempt to address these conceptual as well as practical limitations
of the existing approaches to nonlocal elasticity via fundamental changes in
the nonlocal constitutive modeling. More specifically, we reformulate the
concept of action-at-a-distance by using nonlocal kinematic relations with a
generalized kernel (not restricted to exponential functions). In this study,
we will analyze the specific advantages of the proposed approach and
illustrate how this method is well equipped to address the different
challenges faced by existing integral approaches.

The concept of nonlocal kinematics can be traced back to a seminal study of
Drapaca et al. on modeling nonlocal solids by using fractional calculus [42].
Fractional-order operators are a class of differ-integral operators, uniquely
equipped to model temporal and spatial memory effects in complex materials
[19, 28, 40, 43, 44, 14, 45]. The authors in [42] reformulated the classical
(local) deformation gradient tensor by using fractional-order derivatives and
used the fractional deformation tensor to build the nonlocal continuum theory.
This approach was later refined in [46] and applied to the analysis of slender
engineering structures. A key limitation of the above mentioned studies [42,
46] consists in an ad-hoc imposition of force and moment balance principles.
The governing equations in these studies were derived by a direct replacement
of the local stress tensor (in the classical continuum governing equations)
with the nonlocal stress tensor. However, as shown in [25], the application of
global balance principles to derive governing equations with nonlocal
operators leads to the presence of (additional) integral terms. Further, the
traction boundary conditions in these studies were identical to Cauchy’s
postulate for surface tractions. The presence of nonlocal interactions
however, leads to a modification of Cauchy’s postulate resulting in nonlocal
traction conditions [47, 48, 49]. A significant benefit of adopting
fractional-order (nonlocal) kinematics was exposed in [50, 51, 12], where the
authors could define a positive-definite nonlocal continuum formulation
ultimately enabling the use of variational principles. This latter approach
had both theoretical advantages including the formulation of well-posed
governing equations [50, 51, 12], the incorporation of nonlinear effects [52,
53, 54], and thermodynamic consistency [49], as well as practical implications
like the ability to develop finite element solutions [50]. These advantages
allowed the application of the fractional-order kinematic approach to model
the response of different nonlocal structures subject to mechanical and
thermomechanical loading conditions. In all the above mentioned studies,
consistent predictions free from boundary and loading effects were obtained.

In this study, we extend the fractional-order kinematic approach proposed in
[42, 46, 14, 12], to develop a well-posed and generalized nonlocal kinematic
approach. In the following, we will refer to this approach as the
displacement-driven approach to nonlocal elasticity. Broadly speaking, this
approach uses a general differ-integral (nonlocal) definition to capture the
strain-displacement relations of a nonlocal continuum. It is exactly this
feature that allows the displacement-driven approach to overcome the
limitations of existing integral approaches. The process of formulating the
displacement-driven approach, analyzing its characteristics, and exploring its
application, into three main tasks:

  * •

First, the displacement-driven approach to nonlocal elasticity is developed
based on generalized nonlocal differ-integral operators introduced within the
classical strain-displacement relations. This approach will lead to convex and
positive-definite formulation without requiring a symmetric kernel. In this
process, we will also analyze the satisfaction of locality recovery condition,
the frame-invariance of the formulation, and the consistent behavior of the
nonlocal operator at boundaries. We will complete the analysis of the
displacement-driven continuum approach by casting the formulation within a
thermodynamic framework and demonstrating its thermodynamic consistency. In
this regard, we will show how this approach allows a strong (point-wise)
imposition of thermodynamic balance laws, consistent with the principles of
thermodynamics.

  * •

Following the development of the displacement-driven constitutive relations,
we will derive the equilibrium equations describing the response of the
nonlocal solid by using variational principles. Further, we will derive
dispersion relations describing the transient dynamics of the nonlocal bulk
solid. By taking the example of either an exponential or a power-law kernel,
we will demonstrate how the specific choice of the kernel should not be
universal but rather determined based on experimental dispersion
characteristics carrying the signature of the specific material.

  * •

Finally, we will apply the displacement-driven approach to model the response
of nonlocal beams and plates. For this purpose, we will employ the shear-
deformable Timoshenko beam and Mindlin plate models. Since the Euler-Bernoulli
beam and the Kirchhoff plate formulations can be recovered as special cases
from the Timoshenko and Mindlin formulations, the latter choice allows
obtaining more general results. Using a finite element procedure, we will
numerically simulate the static response of the nonlocal structural elements
for various combinations of nonlocal kernels, kernel parameters, and loading
conditions. The choice of using a finite element based solver follows from the
flexibility afforded by this method in dealing with general loading conditions
as well as integral boundary conditions, typical of nonlocal models [25, 33,
12]. The numerical results will show that the displacement-driven approach
predict a consistent (softening) behavior irrespective of the nature of the
kernel and of the loading conditions.

The remainder of the paper is structured as follows: first we discuss the
formulation and the relevance of nonlocal elasticity via the displacement
driven approach in §2. The analysis of the frame-invariance of the formulation
and of the behavior of the nonlocal operator at boundaries follows in §3. The
thermodynamic consistency of the formulation is discussed in §4. In §5, we
derive the governing equations for the displacement-driven approach. The
dispersion relations for a 1D unbounded nonlocal solid are derived in §6 and
are further used to discuss the importance of adopting material-informed
kernels. Finally, the proposed approach is applied to model the response of
nonlocal beams and plates in §7.

##  2 Reformulation of nonlocal elasticity: Displacement-driven approach

As mentioned in the introduction, the motivation to reformulate the classical
approaches to nonlocal elasticity follows from the non-convex nature of the
strain energy obtained via these approaches. To better highlight this aspect,
we first provide a brief review of the classical integral approaches to
nonlocal elasticity.

In the classical approaches (either strain- or stress-driven) to nonlocal
elasticity, the stress and strain at a given point within the nonlocal
continuum are related by an integral (convolution) operator. This integral
operator enriches the response of a point with the information of the response
of a collection of points within a characteristic distance known as the
horizon of nonlocality. Depending on the specific approach, the stress-strain
constitutive relation consists of the convolution of a spatially
(monotonically) decaying kernel with either the stress or the strain field in
the following fashion [24, 25, 33]:

$\text{Strain-driven approach:
}\sigma_{ij}(\bm{x})=\int_{\overline{\Omega}}\mathbb{C}_{ijkl}\mathcal{K}_{\varepsilon}(\bm{x},\bm{x}^{\prime})\varepsilon_{kl}(\bm{x}^{\prime})\mathrm{d}\bm{x}^{\prime}$
(1a) $\text{Stress-driven approach:
}\varepsilon_{ij}(\bm{x})=\int_{\overline{\Omega}}\mathbb{S}_{ijkl}\mathcal{K}_{\sigma}(\bm{x},\bm{x}^{\prime})\sigma_{kl}(\bm{x}^{\prime})\mathrm{d}\bm{x}^{\prime}$
(1b)

where ${\sigma}_{ij}(\bm{x})$ and ${\varepsilon}_{ij}(\bm{x})$ denote the
stress and strain tensor at a point $\bm{x}$ within the domain of the solid
$\Omega_{\infty}$, and $\bm{x}^{\prime}$ is a point within the horizon of
nonlocality $\overline{\Omega}\subseteq\Omega_{\infty}$. $\bm{\mathbb{C}}$ and
$\bm{\mathbb{S}}$ are the fourth-order material elasticity and compliance
tensors, respectively. $\mathcal{K}_{\varepsilon}(\bm{x},\bm{x}^{\prime})$ and
$\mathcal{K}_{\sigma}(\bm{x},\bm{x}^{\prime})$ are the nonlocal kernels that
determine the strength of the long range interactions between the points
$\bm{x}$ and $\bm{x}^{\prime}$, in the strain- and stress-driven approaches,
respectively.

The constitutive relations in Eq. (1) result in the following expressions for
the total strain energy ($\Pi$) and the strain energy density
($\mathbb{U}(\bm{x})$) of the nonlocal solid [25, 33]:

$\text{Strain-driven approach:
}\Pi=\frac{1}{2}\int_{\Omega_{\infty}}\underbrace{\left[\varepsilon_{kl}(\bm{x})\int_{\overline{\Omega}}\mathbb{C}_{ijkl}\mathcal{K}_{\varepsilon}(\bm{x},\bm{x}^{\prime})\varepsilon_{kl}(\bm{x}^{\prime})\mathrm{d}\bm{x}^{\prime}\right]}_{\mathbb{U}(\bm{x}):\text{
Strain energy density}}\mathrm{d}\bm{x}$ (2a) $\text{Stress-driven approach:
}\Pi=\frac{1}{2}\int_{\Omega_{\infty}}\underbrace{\left[\sigma_{ij}(\bm{x})\int_{\overline{\Omega}}\mathbb{S}_{ijkl}\mathcal{K}_{\sigma}(\bm{x},\bm{x}^{\prime})\sigma_{kl}(\bm{x}^{\prime})\mathrm{d}\bm{x}^{\prime}\right]}_{\mathbb{U}(\bm{x}):\text{
Strain energy density}}\mathrm{d}\bm{x}$ (2b)

It appears from the above expressions that, in either the strain- or stress-
driven approach, the strain energy density cannot be reduced to a quadratic
form, and hence it is not always guaranteed to be positive-definite. While it
is possible to achieve a positive value for the total strain energy, this
imposes restrictions on the nature of the kernel, that is, the kernel must be
positive and symmetric in nature [25, 33]. It was argued, in the introduction,
that this latter requirement on symmetry can prevent the application of the
resulting theory to several cases of practical relevance [38]. As an example,
consider a collection of particles that is either heterogeneous or subject to
thermal gradients. In such case, the interactions on either side of a given
particle are not symmetric [16, 29, 15, 10, 17, 6, 38], hence violating the
assumption of a symmetric kernel. In other terms, from a physical perspective,
the assumption of a symmetric kernel can be considered valid, only for
isotropic structures with an underlying symmetry in the arrangement of
particles [11]. Mathematically speaking, the restriction of a symmetric kernel
is a direct result of the fact that the strain energy density functional is
not guaranteed to be positive or convex in nature (simply because it consists
of the product of the integral of a function and the function itself, which is
not convex in nature). As noted in classical (local) continuum mechanics, a
well-posed elastostatic formulation is guaranteed if and only if the strain
energy density is both convex and positive-definite in nature.

A possible route to achieve a convex and positive-definite formulation
consists in maintaining the quadratic form of the strain energy density,
analogous to local elastic continuum mechanics, such that the stress is a
work-conjugate of the strain:

$\mathbb{U}=\frac{1}{2}\sigma_{ij}\varepsilon_{ij}\equiv\frac{1}{2}\mathbb{C}_{ijkl}\varepsilon_{ij}\varepsilon_{kl}\equiv\frac{1}{2}\mathbb{S}_{ijkl}\sigma_{ij}\sigma_{kl}$
(3)

The use of the above quadratic form for the strain energy requires nonlocality
to be introduced directly into the strain tensor or equivalently into the
stress tensor while using a non-integral (local) stress-strain constitutive
relation. Remarkably, this can be achieved by insisting that the kinematic
strain-displacement relations are nonlocal (differ-integral) in nature while
the stress-strain constitutive relations follow the local continuum mechanics
formulation. By following the latter argument, the infinitesimal strain tensor
and the Cauchy stress tensor are defined as:

${\bm{\varepsilon}}=\frac{1}{2}\left[\overline{\bm{\nabla}}{\bm{u}}+\overline{\bm{\nabla}}{\bm{u}}^{T}\right]$
(4a) $\bm{\sigma}=\bm{\mathbb{C}}:\bm{\varepsilon}$ (4b)

where $\bm{u}(\bm{x})$ is the displacement field. $\overline{\bm{\nabla}}$ is
a differ-integral gradient operator defined as:

$\overline{\bm{\nabla}}(\cdot)=\overline{D}_{x}(\cdot)\hat{x}+\overline{D}_{y}(\cdot)\hat{y}+\overline{D}_{z}(\cdot)\hat{z}$
(5)

where $\\{\hat{x},\hat{y},\hat{z}\\}$ are the Cartesian basis vectors.
$\overline{D}_{x_{j}}(\cdot)$ is a differ-integral operator similar to the
convolution operations in Eq. (1) defined in the following manner:

$\overline{D}_{x_{j}}\phi=\int_{\overline{\Omega}_{j}}c^{*}_{j}\mathcal{K}(\bm{x},\bm{x}^{\prime})\left[D^{1}_{x^{\prime}_{j}}\phi\right]\mathrm{d}x^{\prime}_{j}$
(6)

where $\phi$ is an arbitrary function. $c^{*}_{j}$ is a scalar multiplier for
the differ-integral operator in the $\hat{x}_{j}$ direction, which will be
used to enforce frame-invariance and dimensional consistency of the
formulation in §3. Further, $D^{1}_{x^{\prime}_{j}}\phi$ indicates the
classical first-order derivative $\mathrm{d}\phi/\mathrm{d}x^{\prime}_{j}$.
For a nonlocal horizon $\overline{\Omega}_{j}=[x_{-_{j}},x_{+_{j}}]$ in the
$\hat{x}_{j}$ direction about the point $\bm{x}$, the nonlocal operator can be
expanded as:

$\overline{D}_{x_{j}}\phi=c_{-_{j}}^{*}\int_{x_{-_{j}}}^{x_{j}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\left[D^{1}_{x^{\prime}_{j}}\phi\right]\mathrm{d}x^{\prime}_{j}+c_{+_{j}}^{*}\int_{x_{j}}^{x_{+_{j}}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\left[D^{1}_{x^{\prime}_{j}}\phi\right]\mathrm{d}{x}^{\prime}_{j}$
(7)

where $c_{-_{j}}^{*}$ and $c_{+_{j}}^{*}$ are independent scalar multipliers.
The length of the horizon of nonlocality on either side of the point $\bm{x}$
in the $\hat{x}_{j}$ direction are denoted as $l_{-}$ and $l_{+}$, such that
$l_{-}=x_{j}-x_{-_{j}}$ and $l_{+}=x_{+_{j}}-x_{j}$. In general,
$l_{-_{j}}\neq l_{+_{j}}$ implying that the horizon of nonlocality is not
necessarily symmetric like in cases involving material interfaces or
boundaries falling within the horizon. This aspect is schematically
illustrated in Fig. (1).

Figure 1:  Schematic illustration of the horizon of nonlocality at different
points in a 2D domain. In case of a domain characterized by symmetric nonlocal
interactions the horizon of nonlocality at the point $\bm{x}_{2}$ is equal to
$l_{f}$ on either sides in the $\hat{x}$ direction. However, the horizon of
nonlocality at $\bm{x}_{1}$ must be truncated in the $\hat{x}$, such that
$l_{-_{x}}^{*}<l_{+_{x}}$, due to the presence of the boundary. In a similar
way, it follows that for a point located on the boundary $x=0$ it must be
$l_{-_{x}}=0$.

Revisiting the strain energy of the solid obtained via the above constitutive
formulation, we obtain:

$\mathbb{U}=\frac{1}{2}\underbrace{C_{ijkl}\left[\overline{D}_{x_{j}}u_{i}+\overline{D}_{x_{i}}u_{j}\right]}_{\text{Nonlocal
stress}}\underbrace{\left[\overline{D}_{x_{l}}u_{k}+\overline{D}_{x_{k}}u_{l}\right]}_{\text{Nonlocal
strain}}$ (8)

which immediately demonstrates its convex nature, irrespective of whether the
nonlocal kernel is symmetric or asymmetric in nature. We denominate this
nonlocal approach as the “displacement-driven approach to nonlocal elasticity”
since nonlocality is embedded via the nonlocal strain-displacement kinematic
relations Eq. (4a). Note that the stress-strain constitutive relations in Eq.
(4b) are non-integral (local) in nature contrary to the classical stress- and
strain-driven approaches to nonlocal elasticity [25, 33]. Remarkably, this
functional form for the stress-strain constitutive relation trivially
satisfies the locality recovery condition [27]. More specifically, it ensures
that a uniform strain field induces a uniform stress field, and consequently,
we obtain a length-scale independent strain energy density via Eq. (3).

The displacement-driven formulation deserves some additional remarks. While
the approach was presented starting based on strain energy arguments, the same
formulation can be derived via standard continuum mechanics arguments. In this
study, we favored the former approach to clearly highlight the connection to a
convex and positive-definite nonlocal formulation. When following a continuum
mechanics approach, the displacement-driven formulation can be derived by
starting from a nonlocal deformation gradient tensor which relates the
differential line elements in the deformed and undeformed configurations by
using the differ-integral gradient operator defined in Eq. (5). More
specifically, the nonlocal deformation gradient tensor which maps a
differential element in the undeformed configuration to the deformed
configuration can be defined as:

$\overline{F}_{ij}=\overline{D}_{X_{j}}x_{i}$ (9)

where $\bm{X}$ and $\bm{x}$ denote points within the deformed and undeformed
coordinates, respectively. Following the above definition for the nonlocal
deformation gradient tensor, a fully nonlinear Lagrangian and Eulerian
description of the strain tensor can be derived by using the scalar difference
of the differential line elements in the deformed and undeformed
configurations, similar to classical continuum formulation. The same approach,
albeit by using specialized power-law kernels which transforms the gradient
operator in Eq. (5) into a nonlocal gradient operator containing the well
established fractional-order derivatives, can be found in [42, 46, 14].

Second, note that the model relies on the important hypothesis of nonlocal
kinematics. As highlighted above, in the corresponding continuum description,
the differential line elements are defined using differ-integral deformation
gradients similar to [42, 46, 14]. While this hypothesis might appear
unsettling in the context of the more classical continuum approaches, we
emphasize that such definition is simply a reformulation of the concept of
action-at-a-distance. This hypothesis results in assuming that the response of
a selected point within the solid is affected directly by the response of a
collection of points within the horizon of nonlocality. Given that the differ-
integral operator is applied directly to the displacement field (see Eq.
(4a)), from a physical standpoint, the formulation accounts for long-range
interactions that are proportional to the relative displacement of distant
points within the horizon. It follows that, a change in length of an
infinitesimal line at the point $\bm{x}$ between the reference and the current
configurations would be affected directly by the response of the points within
the nonlocal horizon of $\bm{x}$. This is indeed a possible reformulation of
action-at-a-distance that is often implemented in terms of long-range forces
[14, 50].

##  3 Frame-invariance: Enforcement and Implications

The reformulation of the kinematic relations via differ-integral operators
requires a thorough investigation of the frame-invariance, which states that a
rigid-body motion should not induce a change in the strain energy of the body.
In this section, we will derive the necessary constraints for the scalar
multipliers $c_{-}^{*}$ and $c_{+}^{*}$ such that frame-invariance is
satisfied at all points within the nonlocal continuum. Remarkably, we will
show that the satisfaction of frame-invariance via the scalar multipliers also
ensures a consistent behavior of the nonlocal operator at material boundaries.

###  3.1 Enforcing frame-invariance

Consider a rigid-body motion superimposed on the reference configuration of
the body as:

$\bm{\chi}(\bm{X},t)=\bm{C}(t)+\bm{Q}(t)\bm{X},$ (10)

where $\bm{C}(t)$ is a spatially constant term representing a translation and
$\bm{Q}(t)$ is a proper orthogonal tensor denoting a rotation. Under this
rigid-body motion, the deformation gradient tensor
$\overline{\bm{F}}^{\bm{\chi}}$ should be an orthogonal tensor such that
${\overline{\bm{F}}^{\bm{\chi}}}^{T}\overline{\bm{F}}^{\bm{\chi}}=\textbf{I}$.
More specifically, the nonlocal deformation gradient tensor should transform
as $\overline{\bm{F}}^{\bm{\chi}}=\bm{Q}$ in order to ensure that the strain
measures corresponding to the rigid-body motion are null. From the definition
of $\overline{\bm{F}}^{\bm{\chi}}$ given in Eq. (9) it follows that:

$\begin{split}\overline{D}_{x_{j}}\chi_{i}=c_{-_{j}}^{*}\int_{x_{-_{j}}}^{x_{j}}\mathcal{K}(\bm{x},\bm{x}^{\prime}){D^{1}_{x^{\prime}_{j}}\chi_{i}(\bm{x}^{\prime},t)}\mathrm{d}x^{\prime}_{j}+c_{+_{j}}^{*}\int_{x_{j}}^{x_{+_{j}}}\mathcal{K}(\bm{x},\bm{x}^{\prime}){D^{1}_{x^{\prime}_{j}}\chi_{i}(\bm{x}^{\prime},t)}\mathrm{d}{x}^{\prime}_{j}\end{split}$
(11)

Further, noting that $D^{1}_{x^{\prime}_{j}}C_{i}(t)=0$ and $\bm{Q}=\bm{Q}(t)$
we obtain:

$D^{1}_{x^{\prime}_{j}}\chi_{i}(\bm{x}^{\prime},t)=Q_{ik}x_{k,j}=Q_{ik}\delta_{kj}=Q_{ij}$
(12)

Thus, under the rigid body motion $\bm{\chi}$:

$\overline{D}_{x_{j}}\chi_{i}=\left[c_{-_{j}}^{*}\int_{x_{-_{j}}}^{x_{j}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\mathrm{d}x^{\prime}_{j}+c_{+_{j}}^{*}\int_{x_{j}}^{x_{+_{j}}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\mathrm{d}\bm{x}^{\prime}_{j}\right]Q_{ij}$
(13)

To achieve frame-invariance, that is enforcing
$\overline{D}_{x_{j}}\chi_{i}=Q_{ij}$ at all points within the continuum and
at all time instants, we require that:

$c_{-_{j}}^{*}=\frac{1}{2\int_{x_{-_{j}}}^{x_{j}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\mathrm{d}x^{\prime}_{j}}$
(14a)
$c_{+_{j}}^{*}=\frac{1}{2\int_{x_{j}}^{x_{+_{j}}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\mathrm{d}\bm{x}^{\prime}_{j}}$
(14b)

Notably, the above expressions for $c_{-_{j}}^{*}$ and $c_{+_{j}}^{*}$ also
ensure the dimensional consistency of the formulation. Here below we provide
the expressions of $\\{c_{-_{j}}^{*},c_{+_{j}}^{*}\\}$ for some nonlocal
kernels commonly adopted in literature:

  * •

Exponential kernel:
$\mathcal{K}(\bm{x},\bm{x}^{\prime})=\exp\left(-\frac{|\bm{x}-\bm{x}^{\prime}|}{l_{0}}\right)$,
where $l_{0}$ is a material parameter:

$c_{-_{j}}^{*}=\frac{1}{2l_{0}}\left[1-\exp^{\left(-\frac{l_{-}}{l_{0}}\right)}\right]^{-1}$
(15a)
$c_{+_{j}}^{*}=\frac{1}{2l_{0}}\left[1-\exp^{\left(-\frac{l_{+}}{l_{0}}\right)}\right]^{-1}$
(15b)

  * •

Power-law kernel:
$\mathcal{K}(\bm{x},\bm{x}^{\prime})=|\bm{x}-\bm{x}^{\prime}|^{-\alpha}/\Gamma(1-\alpha)$,
where $\Gamma(\cdot)$ is the well known Gamma function and $\alpha\in(0,1]$ is
a material parameter:

$c_{-_{j}}^{*}=\frac{1}{2}\Gamma(2-\alpha)l_{-}^{\alpha-1}$ (16a)
$c_{+_{j}}^{*}=\frac{1}{2}\Gamma(2-\alpha)l_{+}^{\alpha-1}$ (16b)

While the above results are limited to the exponential and the power-law
kernels, the approach is general and can be extended to other nonlocal kernels
such as, for example, bell and conical curves employed in classical strain-
and stress-driven models of nonlocal elasticity [25]. In each case, the
constants $c_{-}^{*}$ and $c_{+}^{*}$ must be derived by using Eq. (14). We
will analyze the specific impact of different nonlocal kernels in more detail
in §6.

###  3.2 Behavior of the nonlocal operator at boundaries

A detailed analysis of the expressions for the scalar multipliers
$c_{-_{j}}^{*}$ and $c_{+_{j}}^{*}$ in Eq. (14) presents some interesting
insights on the nature of the kernel at material boundaries that will be of
practical relevance. For the most general definition of the kernel
$\mathcal{K}(\bm{x},\bm{x}^{\prime})$ and the location of the point $\bm{x}$
recall that $l_{-}\neq l_{+}$. It immediately follows from Eq. (14)
$c_{-_{j}}^{*}\neq c_{+_{j}}^{*}$. This latter condition, apart from ensuring
an exact satisfaction of frame-invariance, also ensures the completeness of
the nonlocal kernel at points close to the material boundaries or at points
located on the boundary itself. In order to demonstrate this latter statement,
we investigate the behavior of the nonlocal operator
$\overline{D}_{x_{j}}\phi$ at a point located on the boundary of the nonlocal
solid given in Fig. (1).

For a point (say $\bm{x}_{0}$) located on one of the boundaries (identified by
the normal in the $j^{th}$ direction), either $l_{-_{j}}\rightarrow 0$ or
$l_{+_{j}}\rightarrow 0$ (see Fig. (1)). Here below we present the expression
for the nonlocal operator when $l_{-_{j}}\rightarrow 0$ (similar expressions
hold when $l_{+_{j}}\rightarrow 0$). This limiting case, upon using the
expression for $c_{-_{j}}^{*}$ in Eq. (14), gives:

$\lim_{l_{-_{j}}\to 0}\overline{D}_{x_{j}}\phi=\lim_{l_{-_{j}}\to
0}\left[\frac{1}{2\int_{x_{0_{j}}-l_{-_{j}}}^{x_{0_{j}}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\mathrm{d}x^{\prime}_{j}}\int_{x_{0_{j}}-l_{-_{j}}}^{x_{0_{j}}}\mathcal{K}(\bm{x},\bm{x}^{\prime})D^{1}_{x^{\prime}_{j}}\phi\mathrm{d}x^{\prime}_{j}+c^{*}_{+_{j}}\int_{x_{0_{j}}}^{x_{+_{j}}}\mathcal{K}(\bm{x},\bm{x}^{\prime})D^{1}_{x^{\prime}_{j}}\phi\mathrm{d}x^{\prime}_{j}\right]$
(17)

Since the interval length ($=l_{-_{j}}$) of the left integral within the
right-hand side of the above equation is very small,
$D^{1}_{x^{\prime}_{j}}\phi$ can be assumed to be constant and equal to the
boundary condition:

$D^{1}_{x^{\prime}_{j}}\phi=D^{1}_{x_{j}}\phi\Big{|}_{\bm{x}_{0}}$ (18)

Substituting Eq. (18) in Eq. (17) leads to the following:

$\lim_{l_{-_{j}}\to 0}\overline{D}_{x_{j}}\phi=\lim_{l_{-_{j}}\to
0}\left[\frac{D^{1}_{x_{j}}\phi\Big{|}_{\bm{x}_{0}}}{2\cancelto{}{\int_{x_{0_{j}}-l_{-_{j}}}^{x_{0_{j}}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\mathrm{d}x^{\prime}_{j}}}\cancelto{}{\int_{x_{0_{j}}-l_{-_{j}}}^{x_{0_{j}}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\mathrm{d}x^{\prime}_{j}}+c^{*}_{+_{j}}\int_{x_{0_{j}}}^{x_{+_{j}}}\mathcal{K}(\bm{x},\bm{x}^{\prime})D^{1}_{x^{\prime}_{j}}\phi\mathrm{d}x^{\prime}_{j}\right]$
(19)

which immediately yields:

$\lim_{l_{-_{j}}\to
0}\overline{D}_{x_{j}}\phi=\frac{1}{2}\underbrace{D^{1}_{x_{j}}\phi\Big{|}_{\bm{x}_{0}}}_{\begin{subarray}{c}\text{Local
effect due}\\\ \text{to truncation of}\\\ \text{nonlocal
horizon}\end{subarray}}+\underbrace{c^{*}_{+_{j}}\int_{x_{0_{j}}}^{x_{+_{j}}}\mathcal{K}(\bm{x},\bm{x}^{\prime})D^{1}_{x^{\prime}_{j}}\phi\mathrm{d}x^{\prime}_{j}}_{\text{Remaining
nonlocal interaction}}$ (20)

From Eq. (20), it is immediate to observe that while the right-handed operator
captures nonlocality ahead of the point $\bm{x}_{0}$ (in the $j^{th}$
direction), the left-handed operator is reduced to the classical first-order
derivative. From a mathematical standpoint, the definition of the nonlocal
operator ensures that, independent of the nature of the kernel, we recover a
delta-function (and a corresponding local behavior) every time the size of the
nonlocal horizon on any given side of a point approaches (or is equal to)
zero. This suggests that the truncation of the nonlocal horizon and the
corresponding convolution at the boundary has been accounted for in a
consistent manner. From an atomistic or molecular dynamic perspective, the
reduction of the left-handed operator to the classical first-order derivative
is physically analogous to the annihilation (or removal) of long-range
interactions between atoms or molecules close to material boundaries [11].

Before proceeding further, we merely note that the equality
$c_{-_{j}}^{*}=c_{+_{j}}^{*}$ holds if and only if the nonlocal kernel is
symmetric and $l_{-_{j}}=l_{+_{j}}$. The above conditions require that the
body has uniform material properties (such that the nonlocal kernel exhibits
uniform spatial characteristics) and the point under consideration has a
symmetric horizon of nonlocality (possible for points located sufficiently far
from the material boundaries).

##  4 Considerations on the thermodynamics

In this section, we explore the implications of a displacement-driven approach
on the thermodynamic formulation of the nonlocal solid and we will show how
this approach is uniquely equipped to enforce the laws of thermodynamics in a
strong (point-wise) sense.

First law of thermodynamics: For a nonlocal solid, the statement of the
conservation of energy for a point $\bm{x}$ includes, in addition to the local
strain energy, the energy associated with the long-range interactions with
points contained within its horizon of nonlocality. The functional
relationship between the energy density and the different energy exchanges
(local and nonlocal) is typically expressed as
$e=e(\bm{\varepsilon}_{l},\mathcal{R}(\bm{\varepsilon}_{l}),{\eta})$ where
$\bm{\varepsilon}_{l}$ is the local strain field, $\mathcal{R}(\bm{\epsilon})$
denotes a linear integral operator which models nonlocality in the solid, and
$\eta$ is the entropy of the solid. Broadly speaking, the nonlocal strain in
the displacement driven approach allows combining the local strain
$\bm{\varepsilon}_{l}$ and its integral $\mathcal{R}(\bm{\varepsilon}_{l})$.
More specifically, given the nonlocal kinematic relations described in §2, the
contribution of the energy contained in the long-range interactions are fully
captured in the nonlocal strain. In this context, recall the argument that the
nonlocal kinematic relations were a reformulation of the concept of action-at-
a-distance, and hence carry information of nonlocal interactions within the
solid. It immediately follows that the internal energy density of the solid
can be expressed explicitly as a function of the nonlocal strain.
Mathematically, we have $e=e(\bm{\varepsilon},\eta)$ $\bm{\varepsilon}$ is the
nonlocal strain and $\eta$ is the entropy of the solid. This latter conclusion
is remarkable as it allows the first law of thermodynamics to be applied in a
strict sense at every point in the domain:

$\dot{e}={\sigma}_{ij}{\varepsilon}_{ij}+h-q_{i,i}~{}\forall\bm{x}\in\Omega$
(21)

where the notation $(\dot{\square})$ denotes the first-order derivative with
respect to time, $h$ is the heat generated internally per unit volume, and
$\bm{q}$ is the heat flux density.

Second law of thermodynamics: The second law of thermodynamics states that the
internal entropy production rate within a solid is non-negative for all points
inside the solid, that is $\dot{\eta}_{0}\geq 0~{}\forall~{}\bm{x}\in\Omega$.
In order to apply the second law of thermodynamics to the displacement-driven
formulation, we consider the internal entropy production rate:

$\dot{{\eta}}_{\text{0}}=\dot{{\eta}}-\left[\frac{h}{T}-\left(\frac{q_{i}}{T}\right)_{,i}\right]$
(22)

where $\nabla(\cdot)$ denotes the classical divergence operator and $T$
denotes the absolute temperature of the solid. In analogy with the classical
approach, we introduce the Legendre transformation $\psi=e-T\eta$, where
$\psi$ denotes the Helmholtz free energy. By following the rationale presented
in defining $e=e(\varepsilon,\eta)$, it immediately follows that
$\psi=\psi(\varepsilon,T)$. By using the Legendre transformation along with
Eq. (22) we obtain:

$T\dot{{\eta}}_{\text{0}}={\sigma}_{ij}{\varepsilon}_{ij}-\dot{\psi}-\eta\dot{T}-T_{,i}\frac{q_{i}}{T}\geq
0$ (23)

Note that the statement of the second law of thermodynamics for the
displacement-driven approach is analogous to the classical Clausius-Duhem
inequality in its function form, hence presenting a clear departure from
several classical nonlocal approaches characterized by the presence of
additional terms. These additional terms, resulting from a functional
dependence of $\psi$ on $\mathcal{R}(\bm{\varepsilon}_{l})$, casting the
second law of thermodynamics in strong-form. Indeed, these terms disappear
only when a weak-form is considered (that is
$\int_{\Omega}\dot{{\eta}}_{\text{0}}\geq 0$). The ability to satisfy the
second law of thermodynamics only in a weak sense is the main source of
physical inconsistencies of classical nonlocal formulation, as highlighted in
[25].

The Clausius-Duhem inequality in Eq. (23) can be used to derive
thermodynamically consistent constitutive equations which should be equal to
the ones presented in §2. By substituting the functional relationship
$\psi=\psi(\varepsilon,T)$, the inequality in Eq. (23) is simplified as:

$T\dot{{\eta}}_{\text{0}}=\left({\sigma}_{ij}-\frac{\partial\psi}{\partial{\varepsilon}_{ij}}\right)\dot{{\varepsilon}}_{ij}-\left(\eta+\frac{\partial\psi}{\partial
T}\right)\dot{T}-T_{,i}\frac{q_{i}}{T}\geq 0$ (24)

Since the above inequality must hold for arbitrary choices of the independent
fields ${\bm{\varepsilon}}$ and $T$ at all times, we obtain the following
constitutive laws:

${\sigma}_{ij}=\frac{\partial\psi}{\partial{\varepsilon}_{ij}}~{}\forall\bm{x}\in\Omega$
(25a) ${\eta}=-\frac{\partial\psi}{\partial T}~{}\forall\bm{x}\in\Omega$ (25b)

Finally, by using the above constitutive relations within Eq. (23), we obtain:

$T\dot{{\eta}}_{\text{int}}=-T_{,i}\frac{q_{i}}{T}\geq
0~{}\forall\bm{x}\in\Omega$ (26)

which establishes the strong-form of the second law of thermodynamics.

##  5 Governing equations

In this section, we derive the strong-form of the equilibrium equations
governing the response of the nonlocal solid by using variational principles.
Note that the governing equations can also be derived by using a more
classical Newtonian approach based on force and moment balance over a
representative volume element of the nonlocal domain. Although both approaches
yield identical results, the variational approach is deemed more direct given
that the explicit use of forces and moments balance would require additional
care in accounting for the effect of the nonlocal interactions, as described
in detail in [25, 49].

We derive the strong-form of the governing equations by using the Hamilton’s
principle:

$\int_{t_{1}}^{t_{2}}\delta\left(\Pi-V-T\right)\mathrm{d}t=0$ (27)

In the above equation, $V$ and $T$ denote the work done by externally applied
forces and the kinetic energy of the system. Recall that $\Pi$ denotes the
total strain energy of the formulation. Substituting the expressions for the
different physical variables, we obtain:

$\int_{t_{1}}^{t_{2}}\delta\Bigg{[}\underbrace{\frac{1}{2}\int_{\Omega}(\bm{\sigma}:\bm{\varepsilon})\mathrm{d}\mathbb{V}}_{\text{Strain
Energy}}-\underbrace{\int_{\Omega}(\bm{\overline{b}}\cdot\bm{u})\mathrm{d}\mathbb{V}-\int_{\partial\Omega}(\bm{\overline{t}}\cdot\bm{u})\mathrm{d}\mathbb{A}}_{\text{External
Work}}-\underbrace{\frac{1}{2}\int_{\Omega}\rho(\bm{\dot{u}}\cdot{\bm{\dot{u}}})\mathrm{d}\mathbb{V}}_{\text{Kinetic
Energy}}\Bigg{]}\mathrm{d}t=0$ (28)

where $\mathrm{d}\mathbb{V}$ and $\mathrm{d}\mathbb{A}$ indicate volume and
area elements of the nonlocal solid, respectively. $\rho$ denotes the density
of the solid, and $\bm{\overline{b}}$ and $\bm{\overline{t}}$ are prescribed
values of the body force per unit volume and the surface traction per unit
area, respectively.

By simplifying Eq. (28) using principles of variational calculus, we obtain
the governing equation as:

$\bm{\widetilde{\nabla}}\cdot\bm{\sigma}+\bm{\overline{b}}=\rho\ddot{\bm{u}}$
(29)

The associated boundary conditions are obtained as:

$\bm{\tilde{I}}_{\bm{\hat{n}}}\cdot\bm{\sigma}=\bm{\overline{t}}~{}~{}~{}~{}\text{or}~{}~{}~{}~{}{\bm{u}}=\bm{\overline{u}}~{}~{}~{}~{}\forall~{}\bm{x}\in\partial\Omega$
(30)

The operator $\bm{\tilde{I}}_{\bm{\hat{n}}}(\cdot)$ is defined as:

$\bm{\tilde{I}}_{\bm{\hat{n}}}(\cdot)=n_{x}\tilde{I}_{x}(\cdot)\hat{x}+n_{y}\tilde{I}_{y}(\cdot)\hat{y}+n_{z}\tilde{I}_{z}(\cdot)\hat{z}$
(31)

where $\bm{\hat{n}}=n_{x}\hat{x}+n_{y}\hat{y}+n_{z}\hat{z}$ is the unit normal
vector at the surface of the solid. In the above equation,
$\tilde{I}_{x_{j}}(\cdot)$ is an integral operator defined in the following
manner:

$\tilde{I}_{x_{j}}\phi=c_{-_{j}}^{*}\int_{x_{j}-l_{+_{j}}}^{x_{j}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\phi\mathrm{d}x_{j}^{\prime}+c_{+_{j}}^{*}\int^{x_{j}+l_{-_{j}}}_{x_{j}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\phi\mathrm{d}x_{j}^{\prime}$
(32)

Further, the gradient operator denoted by $\bm{\widetilde{\nabla}}(\cdot)$ is
a differ-integral operator (analogous to $\overline{\nabla}$ in Eq. (6))
containing the above defined integral operators:

$\bm{\widetilde{\nabla}}^{\alpha_{m}}(\cdot)=D^{1}_{x}\left[\tilde{I}_{x}(\cdot)\right]\hat{x}+D^{1}_{y}\left[\tilde{I}_{y}(\cdot)\right]\hat{y}+D^{1}_{z}\left[\tilde{I}_{z}(\cdot)\right]\hat{z}\equiv\widetilde{D}_{x}(\cdot)\hat{x}+\widetilde{D}_{y}(\cdot)\hat{y}+\widetilde{D}_{z}(\cdot)\hat{z}$
(33)

The detailed derivation of the governing equations is provided in the
Appendix. Note that the natural boundary conditions are integral (and hence,
nonlocal) in nature, consistent with existing integral approaches to nonlocal
elasticity [25, 33]. From a physical perspective, the nonlocal natural
boundary conditions account for the combined effect of long-range forces
exerted on the boundary points by points located within the nonlocal solid.
While the use of variational principles allowed obtaining the nonlocal natural
boundary conditions in a relatively straight-forward manner, the use of global
conservation principles (mechanical balance laws) requires a careful
modification of the surface tractions to account for the nonlocal interactions
of the boundary points (mentioned previously). In other terms, Cauchy’s
postulate for surface tractions requires a modification to account for the
nonlocal interactions. This latter aspect has been addressed in detail in [47,
48, 49].

##  6 Dispersion relations and choice of kernels

The nonlocal elastodynamic framework presented in §5 allows for the derivation
of wavenumber-frequency dispersion relations describing the transient response
of the nonlocal bulk solid. It is anticipated that the dispersion relations
functionally depend on the nature of the nonlocal kernel used to model
nonlocality within the solid. In this section, we will derive the dispersion
relations for an unbounded 1D nonlocal solid and use them to gain insight into
the specific effects of the nature of the nonlocal kernel.

The governing equation describing a 1D nonlocal solid can be extracted from
Eq. (29) as:

$E\widetilde{D}_{x}\overline{D}_{x}u=\rho\ddot{u}$ (34)

where $E$ is the modulus of elasticity of the isotropic solid. To obtain the
dispersion relation, we substitute in the elastodynamic equation the following
ansatz:

$u(x,t)=u_{0}e^{i(kx-\omega t)}$ (35)

where $u_{0}$ is the amplitude of the longitudinal wave, $k$ denotes the wave-
number, $\omega$ denotes the angular frequency of free longitudinal
vibrations, and $i=\sqrt{-1}$. Using the expressions for the differ-integral
operators in Eq. (34), we obtain the following expression for the dispersion
relation in an unbounded nonlocal solid:

$\omega^{2}+ike^{-ikx}\frac{E}{\rho}D^{1}_{x}\left[\int^{\infty}_{-\infty}\mathcal{K}(x,x^{\prime})\int^{\infty}_{-\infty}\mathcal{K}(x^{\prime},s^{\prime})e^{iks^{\prime}}\mathrm{d}s^{\prime}\mathrm{d}x^{\prime}\right]=0$
(36)

The above expression provides important insights into the nature of the
formulation. First, the expression suggests that for a well-posed formulation
resulting in bounded wave speeds, it is essential that the convolution
integrals in Eq. (36) exist and are bounded. This latter condition is always
satisfied if the nonlocal kernel: 1) decays spatially with the increasing
inter-particle distance, that is,
$D^{1}_{x-x^{\prime}}\mathcal{K}(x,x^{\prime})<0~{}\forall\\{x,x^{\prime}\\}$,
and 2) is positive-definite (more concretely, the Fourier transform of the
kernel is positive everywhere [55]). Note also that, a decaying kernel, also
guarantees the existence of the integrals in Eq. (14) which are essential to
enable a frame-invariant formulation. This condition is also consistent with
existing approaches to nonlocal elasticity, which invariably use decaying
kernels [24, 25, 28, 39, 33]. At this stage, following the analysis of frame-
invariance and of the dispersion relations, we are well-equipped to summarize
all the sufficient requirements for the nonlocal kernel in order to ensure a
physically consistent formulation: the kernel should be positive-definite and
must decay monotonically with the interparticle distance.

Note that further simplifications of the expression in Eq. (36) requires the
specific functional definition of the nonlocal kernel. From a practical
perspective, this suggests that the applicability of a given functional form
of the nonlocal kernel, to model a certain class of nonlocal structures,
depends on whether the dispersion characteristics obtained by simplification
of Eq. (36) match the experimentally obtained dispersion characteristics of
the class of materials under consideration. This latter aspect becomes evident
by considering the example of the exponential and the power-law kernels
introduced in §2. The dispersion relations for the exponential and the power-
law kernels are obtained by using Eq. (36) as:

$\text{Exponential kernel:
}\left(\frac{\omega}{k}\right)^{2}={\frac{E}{\rho}}\left(\frac{1}{1+k^{2}l_{0}^{2}}\right)^{2}$
(37a) $\text{Power-law kernel:
}\left(\frac{\omega}{k}\right)^{2}={\frac{E}{\rho}}\big{[}\underbrace{\cos\left({\pi+\alpha\pi}\right)}_{\text{Propagation}}+i\underbrace{\sin\left({\pi+\alpha\pi}\right)}_{\text{Attenuation}}\big{]}\left(kl^{*}\right)^{2(\alpha-1)}$
(37b)

While the derivation of the dispersion relation is straightforward for the
exponential kernel (given the exponential nature of both the kernel and the
assumed wave solution), some special relations are required to derive the
dispersion relations for the power-law kernel. These relations (that concern
the fractional-order derivatives of an exponential function) as well as the
derivation of the dispersion relations for the power-law kernel can be found
in [12]. In the power-law dispersion, $l^{*}=1$m is a dimensional constant
which is used to enforce dimensional consistency for the 1D unbounded solid
[41, 12]. This constant is required for the 1D unbounded solid because frame-
invariance is trivially satisfied for the only possible rigid-body motion that
is, translation, in this case [12]. The nature of the dispersion relations
obtained for the exponential and power-law kernels presented above lead to the
following observations and conclusions:

  * •

The use of an exponential kernel is suitable for modeling nonlocal structures
which are dispersive but not attenuating in nature. Examples of such
structures can include micro- and nano-structures such as carbon nanotubes,
thin films, monolayer graphene sheets, and nanobeams made of metals [1, 8, 9,
39, 21]. Additionally, macrostructures made of functionally graded materials
and sandwiched cores, also exhibit similar dispersive characteristics [3, 4].

  * •

The power-law kernel is suitable for materials which exhibit anomalous
attenuation and dispersion. Note that the real and imaginary parts of the
dispersion relation in Eq. (37b) correspond to the propagating and attenuating
component of the wave. Examples of such structures can include macrostructures
such as highly scattering media (e.g. fractal, porous, and layered media) and
even animal tissues [18, 19, 40, 20]. Periodic media [41, 17, 14] and
structures with intentionally designed nonlocal geometric features (such as,
for example, acoustic black hole metamaterials [7] and metasurfaces with
intentionally designed nonlocality [6]) also exhibit similar characteristics.
Remarkably, the dispersion as well as the attenuation in the wave speeds in
several classes of these structures, exhibit a power-law dependence on the
wave-number as also captured by the use of the power-law kernel [18, 19, 40,
14].

  * •

The dispersion characteristics obtained using the exponential kernel matches
exactly the dispersion relation obtained via classical differential models
recovered from Eringen’s approach by using exponential kernels [39, 9]. The
motivation behind adopting the differential model follows from a good match of
the predicted dispersion against experimental observations, particularly for
different nanostructures such as, for example, graphene tubules and carbon
nanotubes. As established in [37], this procedure of obtaining a simplified
differential model (based on the classical strain-driven approach) leads to
ill-posed mathematical formulations with inconsistent predictions for
different static loading conditions. On the contrary, the displacement-driven
approach enables a well-posed static framework while recovering the same
dynamic characteristics.

  * •

Since the differential approaches to nonlocal elasticity are recovered by
assuming (only) the exponential kernel within Eringen’s approach [8, 9, 39,
13], these approaches are not suitable to model structures which exhibit
anomalous attenuation-dispersion characteristics. This observation suggests
that the displacement-driven approach is more general, since the kernel can be
modified to capture the behaviour of a wider class of size-dependent
structures.

The above discussion highlights the relevance of the generalized procedure,
involving analysis of both singular (e.g. power-law) and non-singular (e.g.
exponential) kernels, adopted in this study. The discussion also suggests that
the application of any nonlocal approach to modeling of practical structures
must involve the a priori selection of a suitable kernel, guided by
experimentally observed dispersion characteristics.

##  7 Application to nonlocal beams and plates

In this section, we will use the displacement-driven approach to analyze the
static response of slender nonlocal structures, including a Timoshenko beam
and a Mindlin plate. Similar to previous analyses, we will consider separately
the effect of either the exponential or the power-law kernel, on the static
response of the nonlocal structures subject to different loading conditions.
For each structure, we will first present the nonlocal constitutive
formulation, followed by the numerical results generated using the fractional-
order finite element method (f-FEM). The f-FEM is a highly general formulation
that is well-suited for both singular and non-singular kernels [50]. The
algorithm of the f-FEM is discussed in detail in [51] and, for the sake of
brevity, details are not repeated here. For completeness, we note that the
nonlocal formulation involving the exponential kernel can also be simulated by
adapting the classical nonlocal finite element procedure [25, 56, 57], from a
strain-driven approach to a displacement-driven approach. Since we have
considered shear-deformable structures in the following analyses, we adopted a
reduced-order Gauss quadrature integration for the evaluation of the shear
stiffness matrices in order to avoid shear-locking effects [58]. More
specifically, for the nonlocal Timoshenko beam we used $2$ and $1$ Gauss
quadrature elements to numerically integrate the bending and shear stiffness
matrices, respectively. Analogously, for the nonlocal Mindlin plate, we used
$2\times 2$ and $1\times 1$ Gauss quadrature points for the bending and shear
stiffness matrices, respectively.

Analogous to classical finite element procedures, the numerical solution in
the f-FEM is obtained by discretization of the Hamiltonian of the nonlocal
structure by using an isoparametric formulation. Consequently, we only provide
the weak-form of the governing equations for both the nonlocal Timoshenko beam
and the nonlocal Mindlin plate. The corresponding strong-form of the governing
equations can be easily obtained by following the derivation of the 3D
governing equations in §5. In fact, the strong-form nonlocal governing
equations do not generally admit closed-from analytical expressions for the
most general loading conditions. Hence, for the sake of brevity, we do not
provide again the strong-form of the governing equations.

###  7.1 Application to nonlocal beams

In this section, we apply the displacement-driven approach to model and
analyze the static response of nonlocal shear-deformable beams. We consider a
nonlocal beam of uniform cross-section with length $L_{b}$, width $b_{b}$ and
thickness $h_{b}$. A schematic of the beam, along with the chosen Cartesian
coordinate axes, is provided in Fig. (2). As evident from Fig. (2), the
coordinate axes are chosen such that the mid-plane of the beam is identified
as $z=0$. The displacement field components at a spatial point $\bm{x}(x,y,z)$
of the beam can be expressed by following the shear-deformable Timoshenko
formulation as [58]:

$u(\bm{x})=u_{0}(x)-z\theta_{x}(x)$ (38a) $w(\bm{x})=w_{0}(x)$ (38b)

where $u_{0}$ and $w_{0}$ denote the axial and transverse displacement of the
mid-plane of the beam. $\theta_{x}$ denotes the rotation of the normal to the
mid-plane of the beam, about the $\hat{x}$ axis.

Figure 2: Schematic of the rectangular beam illustrating the different
geometric parameters.

Employing the nonlocal strain-displacement relation given in Eq. (4a), we
obtain the following expressions for the strains developed in the nonlocal
beam:

${\varepsilon}_{xx}(\bm{x})=\overline{D}_{x}u_{0}(x)-z\overline{D}_{x}\theta_{x}(x)$
(39a) ${\gamma}_{xz}=\overline{D}_{x}w_{0}(x)-\theta_{x}(x)$ (39b)

The stress developed in the beam can be be derived from the stress-strain
constitutive relation in Eq. (4b).

Using the nonlocal strains and stresses generated in the beam, the total
strain energy can be expressed as:

$\Pi=\int_{0}^{L_{b}}\int_{\frac{-b_{b}}{2}}^{\frac{b_{b}}{2}}\int_{\frac{-h_{b}}{2}}^{\frac{h_{b}}{2}}\left[\sigma_{xx}\varepsilon_{xx}+\kappa_{s}\sigma_{xz}\gamma_{xz}\right]\mathrm{d}z\mathrm{d}y\mathrm{d}x$
(40)

where $k_{s}$ is the shear correction factor chosen as $\kappa_{s}=5/6$
throughout this study [58]. Finally the work done by external loads on the
beam can be expressed as:

$V=\int_{0}^{L_{b}}\int_{\frac{-b_{b}}{2}}^{\frac{b_{b}}{2}}\int_{\frac{-h_{b}}{2}}^{\frac{h_{b}}{2}}\left[F_{x}u_{0}+F_{z}w_{0}+M_{\theta_{x}}\theta_{x}\right]\mathrm{d}z\mathrm{d}y\mathrm{d}x$
(41)

where $\\{F_{x},F_{z}\\}$ are the loads applied externally in the $\hat{x}$
and $\hat{z}$ directions, respectively. Further, $M_{\theta_{x}}$ is the
moment applied about the $\hat{y}$ axis. The nonlocal equilibrium differential
equations describing the static response of the nonlocal beam can now be
derived by applying variational principles ($\delta(\Pi-V)=0$).

####  7.1.1 Results and discussion

In this section, we analyze the effect of different kernels and different
loading conditions on the static bending response of the nonlocal beam. In all
simulations, we considered an isotropic beam having Young’s modulus $E=30$GPa
and Poisson’s ratio $\nu=0.3$. The geometric dimensions of the beam were taken
as $L_{b}=1$m and $b_{b}=h_{b}=L_{b}/10$. As mentioned previously, we analyzed
the effect of two different kernels: the exponential kernel and the power-law
kernel, on the response of the nonlocal beam. In each case, we have assumed an
isotropic and symmetric horizon of nonlocality such that
$l_{-_{x}}=l_{+_{x}}=l_{f}$ for points sufficiently far from the boundaries.
Recall from §2 and §3 that this symmetry is broken ($l_{-_{x}}\neq l_{+_{x}}$)
for points close to the boundaries (see Fig. (1)). Note that while the choice
of the different material properties and geometric parameters was somewhat
arbitrary, their numerical values do not affect the applicability of the model
and the generality of the results presented here below.

Using the above described numerical scheme, we analyzed two different loading
conditions: (a) a cantilever beam subject to a transverse force at its end-
point, and (b) a simply-supported beam subject to a uniformly distributed
transverse load (UDTL). For each boundary condition, we obtained the response
of the beam for the following different kernel parameters:

  * •

Exponential kernel: the kernel parameter $l_{0}$ was varied in
$(0,0.005]\text{m}~{}(=(0,L_{b}/200])$ and the nonlocal horizon parameter
$l_{f}$ was varied in $[0.5,1]\text{m}~{}(=[L_{b}/2,L_{b}])$. The lowest value
for $l_{0}$ was chosen as $l_{0}=10^{-6}$ instead of zero.

  * •

Power-law kernel: the nonlocal horizon $l_{f}$ was varied in
$[0.5,1]\text{m}~{}(=L_{b}/2,L_{b}])$ and the kernel parameter $\alpha$ was
varied in $[0.7,1]$.

The numerical results, expressed in terms of the maximum transverse
displacement of the beam, are presented in Fig. (3) and Fig. (4) for the
exponential and the power-law kernel, respectively. Note that the maximum
transverse displacement of the nonlocal beam is non-dimensionalized against
the analogous value obtained for a local beam. The non-dimensionalized maximum
transverse displacement is indicated by $\overline{w}$. Results show that the
proposed displacement-driven approach predicts a consistent softening
response, that is $\overline{w}>1$, for different combinations of the kernel
parameters and loading conditions. Given this general observation, we proceed
to critically analyze the effects of the two different kernels and of their
parameters.

Figure 3: Bending response measured in terms of the maximum non-
dimensionalized displacement of the nonlocal beam, modeled using the
exponential kernel. The beam is (a) a cantilever beam subject to an transverse
load at its end-point, and (b) simply-supported at all its edges and subject
to a UDTL. The response is parameterized for different values of the kernel
parameter $l_{0}$ and the nonlocal horizon $l_{f}$.

The results presented for the exponential kernel in Fig. (3), lead to the
following observations:

  * •

Effect of $\bm{l_{0}}$: In the limiting case $l_{0}\to 0$, the classical
(local) solution $\overline{w}\to 1$ is recovered [33, 31]. This observation
is consistent with the fact that $l_{0}=0$ reduces the exponential kernel to
the Dirac-delta operator, which results in a reduction of the nonlocal model
to its local elastic counterpart. Note also that an increase in the value of
$l_{0}$ first leads to an increase in the degree of softening followed by an
asymptotic decline, as a function of $l_{0}$. Nonetheless, irrespective of the
specific value of $l_{0}$, the nonlocal beam exhibits a softening effect, that
is, $\overline{w}>1$. This behavior is observed independently of the loading
conditions and is a direct result of the nature of the exponential kernel.

Note that the strength of the nonlocal interactions between two material
points, separated by a certain distance (say $d_{0}$) within the nonlocal
beam, is directly related to the value of the nonlocal kernel. In this case,
the strength of the nonlocal interactions can be given as
$\propto[\exp(-d/l_{0})]/l_{0}$ (the $l_{0}$ in the denominator appears from
the multipliers $c^{*}_{-_{x}}$ and $c^{*}_{-_{x}}$). Consider the following
function for a fixed value of $d_{0}$:

$g(l_{0})=\frac{1}{l_{0}}\exp\left(-\frac{d_{0}}{l_{0}}\right)$ (42)

A quick analysis of the above function reveals that $g(l_{0})$ increases
initially with an increase in $l_{0}$, to reach a global maximum at
$l_{0}=d_{0}$. Following the maximum, it asymptotically approaches zero with
increase in $l_{0}$, that is, $\lim_{l_{0}\to\infty}g(l_{0})=0$. It is exactly
this behavior that is reflected in the static response curves of the nonlocal
beam in Fig. (3). Notably, this behavior is also observed in the response of
atomic lattice models [59, 21, 31]. This observation is also consistent with
our discussion in §6, where we indicated that the dispersion relation obtained
from the exponential kernel is well suited to model the dynamic response of
nano- and micro-structures, typically presented in the form of atomic lattice
chains. Further, the result presented above is also consistent with
predictions made for nonlocal beams in recently developed well-posed integral
approaches [31].

  * •

Effect of $\bm{l_{f}}$: It appears that the length of the horizon of
nonlocality $l_{f}$ does not significantly affect the degree of softening when
$\mathcal{O}(l_{f}/L_{b})=1$. Recall that the parameter $l_{f}$ physically
represents the size of the horizon of nonlocality, that is it determines the
distance beyond which two particles no longer interact via long-range forces.
In this context, the apparent insensitivity of $\overline{w}$ to changes in
$l_{f}$ (in the chosen range) is a direct result of the fact that
$\mathcal{K}(x,x^{\prime})=e^{-|x-x^{\prime}|/l_{0}}\approx 0$ for
$|x-x^{\prime}|>5l_{0}$. Thus, the effect of the nonlocal kernel is always
accounted completely (in the chosen range for $l_{f}$) and does not change
appreciably with a further increase in $l_{f}$. However, note that $l_{f}$
plays a critical role (along with $c^{*}_{-_{x}}$ and $c^{*}_{+_{x}}$) in
achieving frame-invariance and consistency to different boundary conditions.

Figure 4: Bending response measured in terms of the maximum non-
dimensionalized displacement of the nonlocal beam, modeled using the power-law
kernel. The beam is (a) a cantilever beam subject to an transverse load at its
end-point, and (b) simply-supported at all its edges and subject to a UDTL.
The response is parameterized for different values of the kernel parameter
$\alpha$ and the nonlocal horizon $l_{f}$.

We now proceed to analyze the results obtained for the nonlocal beam modeled
via the power-law kernel. A detailed analysis of the results presented in Fig.
(4) leads to the following conclusions:

  * •

Effect of $\bm{\alpha}$: For $\alpha=1$ we recover the classical (local)
solution $\overline{w}=1$, consistent with the nature of the power-law kernel
[14, 46]. Next, observe that a decrease in the value of $\alpha$ leads to a
consistent increase in the degree of softening, irrespective of the nature of
the loading conditions. A decrease in the value of $\alpha$ leads to an
increase in the value of the power-law kernel $|x-x^{\prime}|^{-\alpha}$ for a
fixed value of $|x-x^{\prime}|$. This suggests that a decrease in the value of
$\alpha$ leads to an increase in the degree of nonlocality and a higher degree
of softening. We merely note that the above findings are also consistent with
results presented in existing fractional-order approaches to nonlocal
elasticity (which employ the power-law kernel) [50, 51]. We highlight that
there exists a critical value of $\alpha$ ($=0.4$ for the nonlocal beam),
beyond which the beam undergoes excessive softening and the model breaks down.
More detailed discussions on this aspect can be found in [50, 52, 46].

  * •

Effect of $\bm{l_{f}}$: An increase in the value of $l_{f}$ leads to
significant softening in the case of a power-law kernel. The enhanced
softening is a direct result of the fact that, an increase in the value of
$l_{f}$ leads to an increase in the size of the horizon of nonlocality.
Consequently, a larger number of points within the solid contribute to the
nonlocal interactions leading to an increase in the degree of nonlocality.
Again, this observation is supported by results presented in fractional-order
approaches to nonlocal elasticity [50, 51]. Note that, unlike constraints on
the lower limit for $\alpha$, the value of $l_{f}$ can approach zero. In fact,
in the limiting case $l_{f}\to 0$, we recover the classical local solution
[49, 12].

###  7.2 Application to nonlocal plates

We conclude this study by analyzing the response of nonlocal shear-deformable
plate, in a manner similar to the analysis of nonlocal beams in §7.1. For this
purpose, we consider a rectangular plate of uniform thickness, illustrated in
Fig. (5). The length, width, and thickness of the plate are denoted by
$L_{p}$, $B_{p}$ and $h_{p}$, respectively. The Cartesian reference frame, as
indicated in Fig. (5), is chosen such that $z=0$ denotes the mid-plane of the
plate. In the chosen reference frame, the displacement field components at a
spatial point $\bm{x}(x,y,z)$ can be expressed by following the shear-
deformable Mindlin formulation as [58]:

$u(\bm{x})=u_{0}(x,y)-z\theta_{x}(x,y)$ (43a)
$v(\bm{x})=v_{0}(x,y)-z\theta_{y}(x,y)$ (43b) $w(\bm{x})=w_{0}(x,y)$ (43c)

where $u_{0}$, $v_{0}$, and $w_{0}$ denote the displacements at the mid-plane
of the plate along the $\hat{x}$, $\hat{y}$, and $\hat{z}$ directions.
$\theta_{x}$ and $\theta_{y}$ denote the rotations of the transverse normal at
the mid-plane of the plate, about the $\hat{y}$ and $\hat{x}$ axes.

Figure 5: Schematic of the rectangular plate illustrating the different
geometric parameters.

The nonzero infinitesimal nonlocal strains, corresponding to the deformation
field in Eq. (43), are obtained by using Eq. (4a) as:

$\varepsilon_{xx}(\bm{x})=\overline{D}_{x}u_{0}(x,y)-z\overline{D}_{x}\theta_{x}(x,y)$
(44a)
$\varepsilon_{yy}(\bm{x})=\overline{D}_{y}v_{0}(x,y)-z\overline{D}_{y}\theta_{y}(x,y)$
(44b)
$\gamma_{xy}(\bm{x})=\overline{D}_{y}u_{0}(x,y)+\overline{D}_{x}v_{0}(x,y)-z\left[\overline{D}_{y}\theta_{x}(x,y)+\overline{D}_{y}\theta_{y}(x,y)\right]$
(44c) $\gamma_{xz}(\bm{x})=\overline{D}_{x}w_{0}(x,y)-\theta_{x}(x,y)$ (44d)
$\gamma_{yz}(\bm{x})=\overline{D}_{y}w_{0}(x,y)-\theta_{y}(x,y)$ (44e)

The nonlocal stresses developed in the nonlocal plate can now be evaluated
using the stress-strain constitutive relation in Eq. (4b).

Using the nonlocal strains and stresses generated in the plate, the total
strain energy can be expressed as:

$\Pi=\int_{0}^{L_{p}}\int_{0}^{B_{p}}\int_{\frac{-h_{p}}{2}}^{\frac{h_{p}}{2}}\left[\sigma_{xx}\varepsilon_{xx}+\sigma_{yy}\varepsilon_{yy}+\sigma_{xy}\gamma_{xy}+\kappa_{s}\sigma_{xz}\gamma_{xz}+\kappa_{s}\sigma_{yz}\gamma_{yz}\right]\mathrm{d}z\mathrm{d}y\mathrm{d}x$
(45)

Finally the work done by external loads on the plate can be expressed as:

$V=\int_{0}^{L_{p}}\int_{0}^{B_{p}}\int_{\frac{-h_{p}}{2}}^{\frac{h_{p}}{2}}\left[F_{x}u_{0}+F_{y}v_{0}+F_{z}w_{0}+M_{\theta_{x}}\theta_{x}+M_{\theta_{y}}\theta_{y}\right]\mathrm{d}z\mathrm{d}y\mathrm{d}x$
(46)

where $\\{F_{x},F_{y},F_{z}\\}$ are the loads applied externally in the
$\hat{x}$, $\hat{y}$, and $\hat{z}$ directions, respectively. Further,
$\\{M_{\theta_{x}},M_{\theta_{y}}\\}$ are the moments applied about the
$\hat{y}$ and $\hat{x}$ axes, respectively. The nonlocal equilibrium
differential equations describing the static response of the nonlocal plate
can now be derived by applying variational principles. The response of the
nonlocal plate is numerically simulated using the f-FEM and the results are
presented in the following section.

####  7.2.1 Results and discussion

In this section, we analyze the effect of different kernels and different
boundary conditions on the static bending response of the nonlocal plate. In
all the simulations, we considered an isotropic plate with material properties
identical to the nonlocal beam in §7.1. The in-plane dimensions of the plate
were chosen as $L_{p}=B_{p}=1$m, and thickness was chosen as $h_{p}=L_{p}/10$.
We considered the effect of the exponential kernel and the power-law kernel,
similar to the analysis conducted for the nonlocal beam. In each case, we have
assumed an isotropic and symmetric horizon of nonlocality such that
$l_{-_{j}}=l_{+_{j}}=l_{f}$ for points sufficiently far from the boundaries.
We analyzed the static response of the plate subject to a uniformly
distributed transverse load (UDTL) for two different kinds of boundary
conditions: (a) the plate clamped at all the edges, and (b) the plate simply-
supported at all its edges, for different combinations of the kernel
parameters. The constraints on the generalized displacement coordinates
corresponding to these boundary conditions are as follows [58]:

$\text{Clamped
plate}:\left\\{\begin{matrix}&\hat{x}=\\{0,L_{p}\\}:u_{0}=v_{0}=w_{0}=\theta_{x}=\theta_{y}=0\\\
&\hat{y}=\\{0,L_{p}\\}:u_{0}=v_{0}=w_{0}=\theta_{x}=\theta_{y}=0\end{matrix}\right.$
(47a) $\text{Simply-supported
plate}:\left\\{\begin{matrix}&\hat{x}=\\{0,L_{p}\\}:v_{0}=w_{0}=\theta_{y}=0\\\
&\hat{y}=\\{0,B_{p}\\}:u_{0}=w_{0}=\theta_{x}=0\end{matrix}\right.$ (47b)

For each boundary condition, we obtained the response of the plate for the
following different kernel parameters:

  * •

Exponential kernel: the kernel parameter $l_{0}$ was varied in
$(0,0.005]\text{m}~{}(=(0,L_{p}/200])$ and the nonlocal horizon parameter
$l_{f}$ was varied in $[0.5,1]\text{m}~{}(=[L_{p}/2,L_{p}])$. The lowest value
for $l_{0}$ was chosen as $l_{0}=10^{-6}$ instead of zero. This is because the
choice of $l_{0}=0$ leads to a singularity within the exponential kernel,
wherein the exponential kernel reduces to a Dirac-delta operator [37].

  * •

Power-law kernel: the nonlocal horizon parameter $l_{f}$ was varied in
$[0.5,1]\text{m}~{}(=[L_{p}/2,L_{p}])$ and the kernel parameter $\alpha$, also
called as the fractional-order [14], was varied in $[0.7,1]$. Similar to the
nonlocal beam, there exists a constraint on the lower limit of $\alpha$
(=$0.4$ for the nonlocal plate), while the value of $l_{f}$ can approach to
zero [51].

The numerical results, in terms of the maximum transverse displacement
obtained at the mid-point of the plate, are presented in Fig. (6) and Fig. (7)
for the exponential and the power-law kernel, respectively. Analogous to the
nonlocal beam, the transverse displacement of the mid-point of the nonlocal
plate is non-dimensionalized against the analogous value obtained for a local
Mindlin plate. As evident from the results, the same considerations noted for
the nonlocal Timoshenko beam directly extend to the nonlocal Mindlin plate.

Figure 6: Bending response measured in terms of the non-dimensionalized
displacement at the center point of the nonlocal plate subject to a UDTL and
modeled using the exponential kernel. The plate is (a) clamped at all its
edges and (b) simply-supported at all its edges. The response is parameterized
for different values of the kernel parameter $l_{0}$ and the nonlocal horizon
$l_{f}$. Figure 7: Bending response measured in terms of the non-
dimensionalized displacement at the center point of the nonlocal plate subject
to a UDTL and modeled using the power-law kernel. The plate is (a) clamped at
all its edges and (b) simply-supported at all its edges. The response is
parameterized for different values of the kernel parameter $\alpha$ and the
nonlocal horizon $l_{f}$.

###  7.3 Few additional remarks

The displacement-driven approach predicts a consistent softening response for
both nonlocal beams and plates, irrespective of the nature of the kernel and
the loading conditions. This is a significant outcome when compared to the
often incoherent predictions of classical (integer-order) strain-driven
integral approaches [35, 36, 37]. The results obtained from the displacement-
driven approach are free from stiffening effects or even the absence of
nonlocal effects noted in some strain-driven methods, under selected loading
conditions. This latter statement does not suggest that stiffening effects are
not expected in nonlocal structures; it only implies that stiffening effects
should not be predicted by strain-driven approaches [37, 31].

Indeed, several nonlocal structures do exhibit stiffening effects [37, 12,
31]. In this context, we highlight here that the displacement-driven approach
can be extended in a straight-forward and physically consistent manner to
capture both stiffening and softening effects. This has been established
exclusively for power-law kernels in [12]. Further, the displacement-driven
approach is also well-suited for the geometrically nonlinear analyses of
nonlocal structures. The extension to geometrically nonlinear analyses has
been performed in [52, 53, 54] to analyse the nonlinear bending and
postbuckling response of nonlocal beams and plates described by power-law
kernels. The above discussion suggests that the displacement-driven approach
presents a promising alternative to the analysis of nonlocal structures.

##  8 Conclusions

We presented an approach to model the nonlocal behavior in an elastic solid
that leverages strain-displacement relations formulated in integral form. This
method was dubbed displacement-driven approach to clearly differentiate it
from the existing stress- and strain-driven techniques. The mathematical,
physical, and thermodynamic consistency of the method were carefully
demonstrated and shown to address some critical deficiencies of existing
integral-type models of nonlocal elasticity. This result may also be seen as a
generalization of fractional order formulations of nonlocal elasticity that
restricts the nonlocal kernel to power-law type. The well-posedness of the
elastic governing equations, which is a requisite for modeling well-posed
elastic behavior, is achieved by the present approach. This observation is a
direct consequence of the positive-definite and convex deformation energy
density achievable by using the displacement-driven approach. The formulation
also presents some other important benefits including the generalization to
include asymmetric kernels arising from material heterogeneities, and the
proper truncation of the nonlocal horizon at physical boundaries to maintain
the completeness of the kernel. These features will be fundamental as nonlocal
theories start percolating into a broad spectrum of practical applications. We
also derived the dispersion relations for a 1D infinite solid and showed the
practical influence of some commonly used nonlocal kernels. It was found that
different kernels capture different features of the dynamic behavior, hence
their selection for dynamic analyses should be properly pondered. Finally,
numerical simulations were performed to investigate the static response of
both nonlocal Timoshenko beams and nonlocal Mindlin plates formulated based on
the displacement-driven approach. Numerical results highlighted the very
consistent physical behavior enabled by the displacement-driven approach. This
consistency was observable based on the monotonic softening of the elastic
response following an increase in the degree of nonlocality; a behavior
observed irrespective of the choice of the nonlocal kernel, and of the loading
and boundary conditions. In conclusion, the displacement-driven formulation
provides a powerful framework to model nonlocal effects in elastic solids
while enabling the use of different kernels that are representative of
different dynamic material behavior.

## Acknowledgements

The authors gratefully acknowledge the financial support of the National
Science Foundation (NSF) under grants MOMS #1761423 and DCSD #1825837, and the
Defense Advanced Research Project Agency (DARPA) under grant #D19AP00052. S.P.
acknowledges the partial support received through the R. H. Kohr Graduate
Fellowship. 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.

## Competing Interests

The authors declare that they have no conflict of interest.

## Appendix: Derivation of the strong-form of equilibrium equations

In the following, we present the details of the steps involved in the
derivation of differ-integral (nonlocal) equilibrium equation for the nonlocal
solid. For this purpose, we recall the Hamilton’s principle in Eq. (28).

$\int_{t_{1}}^{t_{2}}\delta\left[\underbrace{\frac{1}{2}\int_{\Omega}(\bm{\sigma}:\bm{\varepsilon})\mathrm{d}\mathbb{V}}_{\text{Strain
Energy}}-\underbrace{\int_{\Omega}(\bm{\overline{b}}\cdot\bm{u})\mathrm{d}\mathbb{V}-\int_{\partial\Omega}(\bm{\overline{t}}\cdot\bm{u})\mathrm{d}\mathbb{A}}_{\text{External
Work}}-\underbrace{\frac{1}{2}\int_{\Omega}\rho(\bm{\dot{u}}\cdot{\bm{\dot{u}}})\mathrm{d}\mathbb{V}}_{\text{Kinetic
Energy}}\right]\mathrm{d}t=0$ (28)

The variations of terms corresponding to the external work and kinetic energy
are identical to their counterparts in local elasticity and hence, in the
following, we will derive only the variation of the strain energy:

$\mathcal{I}_{1}=\int_{t_{1}}^{t_{2}}\delta\left[\frac{1}{2}\int_{\Omega}(\bm{\sigma}:\bm{\varepsilon})\mathrm{d}\mathbb{V}\right]\mathrm{d}t$
(48)

The above expression can be simplified by using the symmetries in the strain
and stress tensors as:

$\mathcal{I}_{1}=\int_{t_{1}}^{t_{2}}\underbrace{\int_{\Omega}{\sigma}_{ij}(\bm{x})~{}\delta\overline{D}_{x_{j}}\left[u_{i}(\bm{x})\right]\mathrm{d}\mathbb{V}}_{\mathcal{I}_{2}}\mathrm{d}t$
(49)

Now, by using the definition of the nonlocal operator
$\overline{D}_{x_{j}}(\cdot)$ from Eq. (7) we obtain:

$\begin{split}\mathcal{I}_{2}=\int_{\Omega_{x_{i}}}\int_{\Omega_{x_{k}}}\int_{\Omega_{x_{j}}}\sigma_{ij}(\bm{x})\left[c_{-_{j}}^{*}\int_{x_{-j}}^{x_{j}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\left[\frac{\mathrm{d}\delta
u_{i}(x^{\prime}_{j})}{\mathrm{d}x_{j}^{\prime}}\right]\mathrm{d}x_{j}^{\prime}+c_{+_{j}}^{*}\int_{x_{j}}^{x_{+j}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\left[\frac{\mathrm{d}\delta
u_{i}(x^{\prime}_{j})}{\mathrm{d}x_{j}^{\prime}}\right]\mathrm{d}x_{j}^{\prime}\right]\mathrm{d}x_{j}\mathrm{d}x_{k}\mathrm{d}x_{i}\end{split}$
(50)

Each of the two terms within the inner integral in the above expression, can
be re-written by changing the order of integration over the nonlocal spatial
variable $x_{j}^{\prime}$ and local spatial variable along the same direction
$x_{j}$. Recall also that $x_{-j}=x_{j}-l_{-}$ and $x_{+j}=x_{j}+l_{+}$. This
gives:

$\mathcal{I}_{21}=c_{-_{j}}^{*}\int_{\Omega_{x_{j}}}\sigma_{ij}(\bm{x})\int_{x_{j}-l_{-}}^{x_{j}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\left[\frac{\mathrm{d}\delta
u_{i}(x^{\prime}_{j})}{\mathrm{d}x_{j}^{\prime}}\right]\mathrm{d}x_{j}^{\prime}\mathrm{d}x_{j}=c_{-_{j}}^{*}\int_{\Omega_{x_{j}}}\left[\frac{\mathrm{d}\delta
u_{i}(x^{\prime}_{j})}{\mathrm{d}x_{j}^{\prime}}\right]\underbrace{\int_{x_{j}}^{x_{j}+l_{-}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\sigma_{ij}(\bm{x})\mathrm{d}x_{j}}_{{I}^{-}_{x_{j}^{\prime}}\left[\sigma_{ij}(x_{j}^{\prime})\right]}\mathrm{d}x_{j}^{\prime}$
(51a)
$\mathcal{I}_{22}=c_{+_{j}}^{*}\int_{\Omega_{x_{j}}}\sigma_{ij}(\bm{x})\int_{x_{j}}^{x_{j}+l_{+}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\left[\frac{\mathrm{d}\delta
u_{i}(x^{\prime}_{j})}{\mathrm{d}x_{j}^{\prime}}\right]\mathrm{d}x_{j}^{\prime}\mathrm{d}x_{j}=c_{+_{j}}^{*}\int_{\Omega_{x_{j}}}\left[\frac{\mathrm{d}\delta
u_{i}(x^{\prime}_{j})}{\mathrm{d}x_{j}^{\prime}}\right]\underbrace{\int_{x_{j}-l_{+}}^{x_{j}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\sigma_{ij}(\textbf{x})\mathrm{d}x_{j}}_{{I}^{+}_{x_{j}^{\prime}}\left[\sigma_{ij}(x_{j}^{\prime})\right]}\mathrm{d}x_{j}^{\prime}$
(51b)

The indicated expressions in the above equation, corresponding to the
convolution of the nonlocal stress tensor, are denoted by
${I}^{-}_{x_{j}}\left[\sigma_{ij}(x_{j}^{\prime})\right]$ and
${I}^{+}_{x_{j}}\left[\sigma_{ij}(x_{j}^{\prime})\right]$. We again use
integration-by-parts to separate the variation of displacement field in the
above expressions which results in the following expressions:

$\begin{split}\mathcal{I}_{21}=c_{-_{j}}^{*}\left[{I}^{-}_{x_{j}}\left[\sigma_{ij}(x_{j}^{\prime})\right]\right]\delta
u_{i}(x^{\prime}_{j})\bigg{|}_{\partial\Omega_{x_{j}}}-\int_{\Omega_{x_{j}}}\delta
u_{i}(x_{j}^{\prime})\left[\frac{\mathrm{d}}{\mathrm{d}x_{j}^{\prime}}\left(c_{-_{j}}^{*}{I}^{-}_{x_{j}}\left[\sigma_{ij}(x_{j}^{\prime})\right]\right)\right]\mathrm{d}x_{j}^{\prime}\end{split}$
(52a)
$\mathcal{I}_{22}=c_{+_{j}}^{*}\left[{I}^{+}_{x_{j}}\left[\sigma_{ij}(x_{j}^{\prime})\right]\right]\delta
u_{i}(x^{\prime}_{j})\bigg{|}_{\partial\Omega_{x_{j}}}-\int_{\Omega_{x_{j}}}\delta
u_{i}(x_{j}^{\prime})\left[\frac{\mathrm{d}}{\mathrm{d}x_{j}^{\prime}}\left(c_{+_{j}}^{*}{I}^{+}_{x_{j}}\left[\sigma_{ij}(x_{j}^{\prime})\right]\right)\right]\mathrm{d}x_{j}^{\prime}$
(52b)

Combining the above results yields the following simplified expression:

$\mathcal{I}_{2}=-\int_{\Omega_{x_{i}}}\int_{\Omega_{x_{k}}}\int_{\Omega_{x_{j}}}\left[\frac{\mathrm{d}\tilde{I}_{x_{j}}\sigma_{ij}(x_{j})}{\mathrm{d}x_{j}}\right]\delta
u_{i}(x_{j})\mathrm{d}x_{j}\mathrm{d}x_{k}\mathrm{d}x_{i}+\int_{\Omega_{x_{i}}}\int_{\Omega_{x_{k}}}\left\\{\tilde{I}_{x_{j}}\sigma_{ij}(x_{j})\right\\}\bigg{|}_{\partial\Omega_{x_{j}}}\mathrm{d}x_{k}\mathrm{d}x_{i}$
(53)

where, $\tilde{I}_{x_{j}}\sigma_{ij}(x_{j})$ follows the definition of the
integral operator defined in Eq. (32). In obtaining the above result, we
interchanged the dummy variables $x_{j}^{\prime}$ and $x_{j}$ are interchanges
below for the sake of brevity. Finally, by substituting the above result in
Eq. (49), we obtain:

$\mathcal{I}_{1}=\int_{t_{1}}^{t_{2}}\left[-\int_{\Omega}\left[\left(\bm{\widetilde{\nabla}}^{\alpha}\cdot\bm{\sigma}\right)\cdot\delta\bm{u}\right]\mathrm{d}\mathbb{V}+\\\
\int_{\partial\Omega}\left[\left(\bm{I}^{1-\alpha_{1}}_{\bm{\hat{n}}}\cdot\bm{\sigma}\right)\cdot\delta\bm{u}\right]\mathrm{d}\mathbb{A}\right]\mathrm{d}t$
(54)

Applying the fundamental principle of variational calculus, we recover the
governing equations in Eq. (29) and boundary conditions in Eq. (30).

## References

  * [1] Eringen, A. C.  Linear theory of nonlocal elasticity and dispersion of plane waves.  _International Journal of Engineering Science_ 10, 425–435 (1972).
  * [2] Bazant, Z. P.  Instability, ductility, and size effect in strain-softening concrete.  _ASCE J Eng Mech Div_ 102, 331–344 (1976).
  * [3] Romanoff, J. & Reddy, J.  Experimental validation of the modified couple stress timoshenko beam theory for web-core sandwich panels.  _Composite Structures_ 111, 130–137 (2014).
  * [4] Romanoff, J., Karttunen, A. T., Varsta, P., Remes, H. & Reinaldo Goncalves, B.  A review on non-classical continuum mechanics with applications in marine engineering.  _Mechanics of Advanced Materials and Structures_ 27, 1065–1075 (2020).
  * [5] Alotta, G., Di Paola, M., Pinnola, F. P. & Zingales, M.  A fractional nonlocal approach to nonlinear blood flow in small-lumen arterial vessels.  _Meccanica_ 1–16 (2020).
  * [6] Zhu, H., Patnaik, S., Walsh, T. F., Jared, B. H. & Semperlotti, F.  Nonlocal elastic metasurfaces: Enabling broadband wave control via intentional nonlocality.  _Proceedings of the National Academy of Sciences_ 117, 26099–26108 (2020).
  * [7] Nair, S.  _Nonlocal Acoustic Black Hole Metastructures: Achieving Ultralow Frequency and Broadband Vibration Attenuation_.  Ph.D. thesis, Purdue University Graduate School (2019).
  * [8] Pradhan, S. & Murmu, T.  Small scale effect on the buckling of single-layered graphene sheets under biaxial compression via nonlocal continuum mechanics.  _Computational materials science_ 47, 268–274 (2009).
  * [9] Wang, C. Y., Murmu, T. & Adhikari, S.  Mechanisms of nonlocal effect on the vibration of nanoplates.  _Applied Physics Letters_ 98, 153101 (2011).
  * [10] Trovalusci, P., De Bellis, M. L. & Masiani, R.  A multiscale description of particle composites: From lattice microstructures to micropolar continua.  _Composites Part B: Engineering_ 128, 164–173 (2017).
  * [11] Tuna, M., Kirca, M. & Trovalusci, P.  Deformation of atomic models and their equivalent continuum counterparts using eringen’s two-phase local/nonlocal model.  _Mechanics Research Communications_ 97, 26–32 (2019).
  * [12] Patnaik, S., Sidhardh, S. & Semperlotti, F.  Towards a unified approach to nonlocal elasticity via fractional-order mechanics.  _International Journal of Mechanical Sciences_ 189, 105992 (2020).
  * [13] Lal, R. & Dangi, C.  Thermomechanical vibration of bi-directional functionally graded non-uniform timoshenko nanobeam using nonlocal elasticity theory.  _Composites Part B: Engineering_ 172, 724–742 (2019).
  * [14] Patnaik, S. & Semperlotti, F.  A generalized fractional-order elastodynamic theory for non-local attenuating media.  _Proceedings of the Royal Society A_ 476, 20200200 (2020).
  * [15] Trovalusci, P., De Bellis, M. L., Ostoja-Starzewski, M. & Murrali, A.  Particulate random composites homogenized as micropolar materials.  _Meccanica_ 49, 2719–2727 (2014).
  * [16] Silling, S. A.  Reformulation of elasticity theory for discontinuities and long-range forces.  _Journal of the Mechanics and Physics of Solids_ 48, 175–209 (2000).
  * [17] Hollkamp, J. P. & Semperlotti, F.  Application of fractional order operators to the simulation of ducts with acoustic black hole terminations.  _Journal of Sound and Vibration_ 465, 115035 (2020).
  * [18] Fellah, Z. E. A., Chapelon, J. Y., Berger, S., Lauriks, W. & Depollier, C.  Ultrasonic wave propagation in human cancellous bone: Application of biot theory.  _The Journal of the acoustical Society of america_ 116, 61–73 (2004).
  * [19] Fellah, Z. E. A., Berger, S., Lauriks, W. & Depollier, C.  Verification of kramers–kronig relationship in porous materials having a rigid frame.  _Journal of sound and vibration_ 270, 865–885 (2004).
  * [20] Buonocore, S., Sen, M. & Semperlotti, F.  Occurrence of anomalous diffusion and non-local response in highly-scattering acoustic periodic media.  _New Journal of Physics_ 21, 033011 (2019).
  * [21] Lim, C., Zhang, G. & Reddy, J.  A higher-order nonlocal elasticity and strain gradient theory and its applications in wave propagation.  _Journal of the Mechanics and Physics of Solids_ 78, 298–313 (2015).
  * [22] Kröner, E.  Elasticity theory of materials with long range cohesive forces.  _International Journal of Solids and Structures_ 3, 731–742 (1967).
  * [23] Edelen, D. G. & Laws, N.  On the thermodynamics of systems with nonlocality.  _Archive for Rational Mechanics and Analysis_ 43, 24–35 (1971).
  * [24] Eringen, A. C. & Edelen, D. G. B.  On nonlocal elasticity.  _International Journal of Engineering Science_ 10, 233–248 (1972).
  * [25] Polizzotto, C.  Nonlocal elasticity and related variational principles.  _International Journal of Solids and Structures_ 38, 7359–7380 (2001).
  * [26] Polizzotto, C., Fuschi, P. & Pisano, A.  A strain-difference-based nonlocal elasticity model.  _International Journal of Solids and Structures_ 41, 2383–2401 (2004).
  * [27] Polizzotto, C., Fuschi, P. & Pisano, A.  A nonhomogeneous nonlocal elasticity model.  _European Journal of Mechanics-A/Solids_ 25, 308–333 (2006).
  * [28] Lazopoulos, K. A.  Non-local continuum mechanics and fractional calculus.  _Mechanics research communications_ 33, 753–757 (2006).
  * [29] Askari, E. _et al._ Peridynamics for multiscale materials modeling.  In _Journal of Physics: Conference Series_ , vol. 125, 012078 (IOP Publishing, 2008).
  * [30] Di Paola, M., Failla, G., Pirrotta, A., Sofi, A. & Zingales, M.  The mechanically based non-local elasticity: An overview of main results and future challenges.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 371, 20120433 (2013).
  * [31] Pisano, A. A., Fuschi, P. & Polizzotto, C.  Integral and differential approaches to eringen’s nonlocal elasticity models accounting for boundary effects with applications to beams in bending.  _ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik_ e202000152 (2021).
  * [32] Lazopoulos, K. & Lazopoulos, A.  On plane $\lambda$-fractional linear elasticity theory.  _Theoretical and Applied Mechanics Letters_ 10, 270–275 (2020).
  * [33] Romano, G. & Barretta, R.  Nonlocal elasticity in nanobeams: the stress-driven integral model.  _International Journal of Engineering Science_ 115, 14–27 (2017).
  * [34] Shaat, M., Ghavanloo, E. & Fazelzadeh, S. A.  Review on nonlocal continuum mechanics: Physics, material applicability, and mathematics.  _Mechanics of Materials_ 103587 (2020).
  * [35] Challamel, N. & Wang, C. M.  The small length scale effect for a non-local cantilever beam: A paradox solved.  _Nanotechnology_ 19, 345703 (2008).
  * [36] Fernández-Sáez, J., Zaera, R., Loya, J. A. & Reddy, J. N.  Bending of Euler–Bernoulli beams using Eringen’s integral formulation: A paradox resolved.  _International Journal of Engineering Science_ 99, 107–116 (2016).
  * [37] Romano, G., Barretta, R., Diaco, M. & de Sciarra, F. M.  Constitutive boundary conditions and paradoxes in nonlocal elastic nanobeams.  _International Journal of Mechanical Sciences_ 121, 151–156 (2017).
  * [38] Batra, R.  Misuse of eringen’s nonlocal elasticity theory for functionally graded materials.  _International Journal of Engineering Science_ 159, 103425.
  * [39] Askes, H. & Aifantis, E. C.  Gradient elasticity in statics and dynamics: an overview of formulations, length scale identification procedures, finite element implementations and new results.  _International Journal of Solids and Structures_ 48, 1962–1990 (2011).
  * [40] Magin, R. L.  Fractional calculus models of complex dynamics in biological tissues.  _Computers & Mathematics with Applications_ 59, 1586–1593 (2010).
  * [41] Hollkamp, J. P., Sen, M. & Semperlotti, F.  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).
  * [42] Drapaca, C. & Sivaloganathan, S.  A fractional model of continuum mechanics.  _Journal of Elasticity_ 107, 105–123 (2012).
  * [43] Carpinteri, A., Cornetti, P. & Sapora, A.  Nonlocal elasticity: An approach based on fractional calculus.  _Meccanica_ 49, 2551–2569 (2014).
  * [44] Mashayekhi, S., Miles, P., Hussaini, M. Y. & Oates, W. S.  Fractional viscoelasticity in fractal and non-fractal media: Theory, experimental validation, and uncertainty analysis.  _Journal of the Mechanics and Physics of Solids_ 111, 134–156 (2018).
  * [45] Failla, G. & Zingales, M.  Advanced materials modelling via fractional calculus: challenges and perspectives (2020).
  * [46] Sumelka, W.  Thermoelasticity in the framework of the fractional continuum mechanics.  _Journal of Thermal stresses_ 37, 678–706 (2014).
  * [47] Dell’Isola, F., Seppecher, P. & Madeo, A.  How contact interactions may depend on the shape of cauchy cuts in nth gradient continua: approach “à la d’alembert”.  _Zeitschrift für angewandte Mathematik und Physik_ 63, 1119–1141 (2012).
  * [48] Dell’Isola, F., Seppecher, P. & Corte, A. D.  The postulations á la d’alembert and á la cauchy for higher gradient continuum theories are equivalent: a review of existing results.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 471, 20150415 (2015).
  * [49] Sidhardh, S., Patnaik, S. & Semperlotti, F.  Thermodynamics of fractional-order nonlocal continua and its application to the thermoelastic response of beams.  _European Journal of Mechanics-A/Solids_ 104238 (2021).
  * [50] Patnaik, S., Sidhardh, S. & Semperlotti, F.  A Ritz-based finite element method for a fractional-order boundary value problem of nonlocal elasticity.  _International Journal of Solids and Structures_ 202, 398–417 (2020).
  * [51] Patnaik, S., Sidhardh, S. & Semperlotti, F.  Fractional-order models for the static and dynamic analysis of nonlocal plates.  _Communications in Nonlinear Science and Numerical Simulation_ 95, 105601 (2020).
  * [52] Sidhardh, S., Patnaik, S. & Semperlotti, F.  Geometrically nonlinear response of a fractional-order nonlocal model of elasticity.  _International Journal of Nonlinear Mechanics_ 125, 103529 (2020).
  * [53] Patnaik, S., Sidhardh, S. & Semperlotti, F.  Geometrically nonlinear analysis of nonlocal plates using fractional calculus.  _International Journal of Mechanical Sciences_ 179, 105710 (2020).
  * [54] Sidhardh, S., Patnaik, S. & Semperlotti, F.  Analysis of the post-buckling response of nonlocal plates via fractional order continuum theory.  _Journal of Applied Mechanics_ 88, 1–22 (2020).
  * [55] Bažant, Z. P. & Chang, T. P.  Instability of nonlocal continuum and strain averaging.  _Journal of Engineering Mechanics_ 110, 1441–1450 (1984).
  * [56] Ansari, R., Torabi, J. & Norouzzadeh, A.  Bending analysis of embedded nanoplates based on the integral formulation of eringen’s nonlocal theory using the finite element method.  _Physica B: Condensed Matter_ 534, 90–97 (2018).
  * [57] Sidhardh, S. & Ray, M. C.  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, 297–311 (2018).
  * [58] Reddy, J. N.  _Mechanics of laminated composite plates and shells: theory and analysis_ (CRC press, 2003).
  * [59] Eringen, A. C.  On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves.  _Journal of applied physics_ 54, 4703–4710 (1983).