# Thermodynamics of fractional-order nonlocal continua and its application to
the thermoelastic response of beams

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

###### Abstract

This study presents a comprehensive framework for constitutive modeling of a
frame-invariant fractional-order approach to nonlocal thermoelasticity in
solids. For this purpose, thermodynamic and mechanical balance laws are
derived for nonlocal solids modeled using the fractional-order continuum
theory. This includes revisiting the Cauchy’s hypothesis for surface traction
vector in order to account for long-range interactions across the domain of
nonlocal solid. Remarkably, it is shown that the fractional-order model allows
the rigorous localized application of thermodynamic balance principles unlike
existing integral approaches to nonlocal elasticity. Further, the mechanical
governing equations of motion for the fractional-order solids obtained here
are consistent with existing results from variational principles. These
fractional-order governing equations involve self-adjoint operators and admit
unique solutions, in contrast to analogous studies following the local
Cauchy’s hypothesis. To illustrate the efficacy of this framework, case-
studies for the linear and the geometrically nonlinear responses of nonlocal
beams subject to combined thermomechanical loads are considered here.
Comparisons with existing integer-order integral nonlocal approaches highlight
a consistent softening response of nonlocal structures predicted by the
fractional-order framework, irrespective of the boundary and thermomechanical
loading conditions. This latter aspect addresses an important incongruence
often observed in strain-based integral approaches to nonlocal elasticity.

Keywords- Fractional calculus; Nonlocal elasticity; Constitutive Modeling;
Geometric nonlinearity; Thermoelasticity
† All correspondence should be addressed to<EMAIL_ADDRESS>
##  1 Introduction

Several theoretical and experimental studies have shown that size-dependent
effects, also referred to as nonlocal effects, are prominent in the response
of complex structures of great relevance for many real-world applications.
These size-dependent effects can be traced back to medium heterogeneity,
existence of surface stresses, presence of thermal loads, and even medium
geometry. More specifically, in the case of micro- and nano-structures, size-
dependent effects have been traced back to the existence of surface and
interface stresses due to nonlocal atomic interactions and Van der Waals
forces [1, 2, 3]. In the case of macroscale structures, nonlocal effects can
result from an ensemble of factors including material heterogeneity,
interactions between layers (e.g. in FGMs or composite media) or unit cells
(e.g. in periodic media), and geometric inhomogeneity [4, 5, 6]. In other
terms, nonlocal governing equations for macrostructures often result from a
process of homogenization of the initial inhomogeneous system. Further,
geometric effects such as changes in curvature have also been shown to induce
nonlocal size-dependent effects in nano-, micro-, and macro-structures [1, 3].

The above mentioned complex slender structures have important applications in
engineering and biotechnology. As an example, macroscale structures made from
functionally graded materials (FGM) or sandwiched designs have been largely
used in weight-critical applications such as aerospace, naval, and automotive
systems [7, 8]. Similarly, thin films, carbon nanotubes, monolayer graphene
sheets and micro tubules have far-reaching applications in atomic devices,
micro/nano-electromechanical devices, as well as sensors and biological
implants [1, 9]. Independently on the spatial scale, the key design
constraints in the above applications include restrictions on space and
weight. As a result, structural assemblies for lightweight applications are
typically made of a combination of slender components like beams, plates, and
shells. There are several applications where these structures are subject to
large and rapidly varying mechanical and thermal loads that drive the system
into a geometrically nonlinear regime. A practical example includes the
analysis of supersonic or hypersonic aerospace systems where the combination
of large and quickly varying aero-thermomechanical loads induces highly a
nonlinear response [10, 7, 8]. Similarly, the ability to account for coupled
thermomechanical nonlinear effects is also critical in applications involving
nano- and micro-structures such as, for example, in the design of biological
implants, measurement devices, and sensors [9, 11]. Despite the undeniable
need for proper theoretical frameworks and computational tools capable of
simulating the thermomechanical response of nonlinear and nonlocal structures,
only a limited amount of studies focusing on geometrically nonlinear
thermomechanical response of nonlocal slender structures are available in the
literature. In the following, we briefly review the main characteristics of
these studies and discuss key limitations.

Seminal works from Kron̈er [12] and Eringen [13] have explored the role of
nonlocality in elasticity and laid its theoretical foundation. The key
principle behind nonlocal theories relies on the idea that all the particles
located within a prescribed area, typically indicated as the horizon of
nonlocality, influence one another by means of long-range cohesive forces.
This interaction between particles is accounted for by using gradient or
integral relations for the strain field within the constitutive equations.
These approaches lead to so-called “weak” gradient methods or “strong”
integral methods, respectively. Integral methods [13, 14, 15] capture nonlocal
effects by re-defining the constitutive law in the form of a convolution
integral of either the strain or the stress field over the horizon of
nonlocality, whereas gradient elasticity theories [16, 17, 18, 19] account for
the nonlocal behavior by introducing strain or stress gradient dependent terms
in the stress-strain constitutive law. As emphasized earlier, several
applications involving nonlocal slender structures also experience thermal
loads and geometric nonlinearities. Although several studies are available on
the topics of geometrically nonlinear response of nonlocal slender structures
[20, 21, 9] and on the thermomechanical response of nonlocal structures [14,
22, 23, 11], theoretical and numerical methods capable of addressing the
combined geometrically nonlinear thermomechanical response of nonlocal
structures have not been specifically addressed.

As mentioned previously, the classical studies on nonlocal (either linear or
nonlinear) elasticity and nonlocal (linear) thermoelasticity encounter some
key shortcomings. As an example, gradient theories experience difficulties
when enforcing the boundary conditions associated with the strain gradient-
dependent terms [16, 17]. On the other side, the integral methods are better
suited to deal with boundary conditions but they lead to mathematically ill-
posed governing equations. This mathematical ill-posedness leads to erroneous
predictions such as the absence of nonlocal effects or the occurrence of
hardening behavior (not consistent with integral models) for certain
combinations of boundary conditions [24]. In this class of problems, the ill-
posedness stems from the fact that the constitutive relation between the
bending moment field and the curvature is a Fredholm integral of the first
kind, whose solution does not generally exists and, if it exists, it is not
necessarily unique [24]. Additionally, in both these classes of methods, there
are no available explicit relations to estimate the stress at a given point
given the strain at that particular point. This latter aspect prevents the
application of variational principles [25, 26] and has critical implications
on the development of a thermodynamic framework for the classical nonlocal
approaches. More specifically, the modeling of nonlocality through nonlocal
stress-strain constitutive relations allows only for a weak application (in a
domain integral sense) and prevents the localized (point-wise) application of
the thermodynamic balance laws. As discussed in [14], the weak application of
thermodynamic balance laws, particularly the second law, leads to
inconsistencies in the nonlocal continuum framework.

In recent years, fractional calculus has emerged as a powerful mathematical
tool to model a variety of nonlocal and multiscale phenomena. Fractional
derivatives, which are a differ-integral class of operators, are intrinsically
multiscale and provide a natural way to account for nonlocal effects [27].
Given the multiscale nature of fractional operators, fractional calculus has
found wide-spread applications in nonlocal elasticity [28, 29, 30, 31, 32, 5,
6, 33]. In a series of papers, Patnaik et al. [6, 34, 35, 36, 37] have shown
that a nonlocal continuum approach based on fractional-order kinematic
relations provides an effective way to address the previously mentioned
shortcomings of classical approaches to nonlocal elasticity. In addition, the
formulation developed in these works is fully frame-invariant. Note that,
unlike gradient elasticity methods, additional essential boundary conditions
are not required when using the Caputo definition of the fractional operator
[5, 6]. Further, the nonlocal model based on fractional-order kinematic
relations allows the application of variational principles and leads to well-
posed governing equations that admit unique solutions [34, 36]. Finally we
note that the fractional-order continuum theories are shown to be effective in
developing reduced-order models for inhomogeneous systems, particularly
periodic structures: periodic bar [5], acoustic black hole structures [33],
and grid-stiffened plates [6].

In this study, we build upon the fractional-order nonlocal continuum theory
and develop a frame-invariant thermodynamic framework for the nonlocal solids.
The overall goal of this study is two fold. First, to develop a thermodynamic
framework for the fractional-order continuum formulation. We will show that
the fractional-order continuum formulation allows for a rigorous application
of all thermodynamic principles. More specifically, the use of fractional-
order kinematic relations prevents the requirement of additional integral
constitutive stress-strain relations as seen in classical nonlocal approaches
(see, for example, [14]). In fact, these additional residuals due to nonlocal
interactions are assimilated within the corresponding fractional-order
analogues. The most direct result is that the formulation does not require
additional constraints associated with the thermodynamic balance laws. It
follows that the thermodynamic balance laws in the fractional-order theory are
free from nonlocal residual terms, hence greatly simplifying the constitutive
modeling of the nonlocal continuum and enabling a rigorous implementation of
the thermodynamic principles at each point in the solid. The latter
observation highlights an important benefit and a key motivation to pursue a
fractional-order formulation to nonlocal thermoelasticity as it sets the
fractional approach apart from other classical nonlocal theories (either
integral[38, 39] or differential[40]) which typically can only satisfy
thermodynamic equilibrium in a weak (global) sense. Then, by enforcing the
mechanical balance laws, we derive the governing equations of motion for the
fractional-order nonlocal continuum. In this process, we employ a modification
to Cauchy’s hypothesis for the traction vector in order to include the effect
of long-range interactions. The thermodynamic framework is also of critical
importance because it allows extending the energy-based methodologies typical
of classical elasticity theories to the analysis of nonlocal structures.
Examples of their applications include the linear buckling [41] and post-
buckling [42] response of nonlocal structures that could not be tackled by
standard integral models of nonlocal elasticity. The second objective of this
study is to highlight the performance of the framework via a specific case-
study focusing on a fractional-order Euler-Bernoulli beam subject to combined
thermomechanical loads. For this purpose, we extend the fractional-order
finite element model (f-FEM), originally developed in [35], to accurately
solve the geometrically nonlinear fractional-order thermomechanical equations.

We note that several fractional-order thermomechanical models have been
previously developed and presented in the literature [43, 44, 45, 46].
However, these studies have focused primarily on the use of fractional-order
operators to model complex thermal exchanges. More specifically, time-
fractional operators have been employed to rewrite the heat conduction
equation [43, 46] in order to model dissipative effects associated with the
thermal processes. Space-fractional operators have been used in [44, 47],
still within the heat transfer governing equation, to model anomalous forms of
spatial diffusion. We emphasize that, differently from the present study,
these previous works considered a local form of the stress-strain constitutive
relation. On the contrary, our work does not consider fractional heat
transfer, hence the heat conduction equation used in our study matches the
classical integer-order form. In this regard, we merely note that a recent
study [47] has shown that the use of a space-fractional heat conduction
equation leads to inconsistencies in the application of the second law of
thermodynamics. This latter observation motivated us to use the classical
(integer-order) heat conduction equation in our study.

The remainder of the paper is structured as follows: we begin with the
development of a constitutive model for the fractional-order approach to
nonlocal elasticity. This involves the thermodynamic and mechanical balance
laws applied to the fractional-order solids. Later, to highlight the
significance of the fractional-order models, we present the case-study of both
a linear and a geometrically nonlinear response of nonlocal beams subject to
combined thermomechanical loads and solved numerically via finite element
techniques.

##  2 Constitutive model for fractional-order thermoelasticity

In this section, we develop the constitutive model for fractional-order
nonlocal thermoelasticity. As discussed earlier, the nonlocal beam theory
presented in this study builds upon the formulation of a fractional-order
nonlocal continuum presented in [6]. This formulation is a generalization of
the seminal works on fractional-order continuum theories for nonlocal solids
developed in [48, 30, 31]. In the following, we review the key highlights of
the continuum theory, and proceed with the development of constitutive model
for nonlocal thermoelasticity.

Note the following notation used throughout the manuscript: $(\dot{\square})$
denotes the first integer-order derivative with respect to time, comma
notation in the subscript $\square_{,j}$ will be used to denote integer-order
spatial derivative with respect to the coordinate $x_{j}$, and Einstein
summation is implied for repeated indices.

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

Analogous to the classical approach to continuum mechanics, the response of a
nonlocal solid can be analyzed by introducing two configurations, namely, the
reference (undeformed) and the current (deformed) configurations. The motion
of the body from the reference configuration (denoted as $\bm{X}$) to the
current configuration (denoted as $\bm{x}$) is assumed as:

$\bm{x}=\bm{\Phi}(\bm{X},t)$ (1)

such that $\bm{\Phi}(\bm{X},t)$ is a bijective mapping operation. The above
mapping operation is used to model the differential line elements
$\mathrm{d}\tilde{\bm{X}}$ and $\mathrm{d}\tilde{\bm{x}}$ in the undeformed
and deformed configurations of the nonlocal solid using fractional-order
operators. The fractional-order deformation gradient tensor
$\overset{\alpha}{\textbf{F}}(\bm{x},\bm{X})$ defined with respect to nonlocal
line elements is given by:

$\overset{\alpha}{\textbf{F}}=\frac{\mathrm{d}\tilde{\bm{x}}}{\mathrm{d}\tilde{\bm{X}}}$
(2)

In analogy with classical strain measures, the nonlocal strain can be defined
using the fractional-order differential line elements as
$\mathrm{d}\tilde{\bm{x}}\mathrm{d}\tilde{\bm{x}}-\mathrm{d}\tilde{\bm{X}}\mathrm{d}\tilde{\bm{X}}$.
Following the above definition for fractional-order deformation-gradient
tensor, the strain in the nonlocal solid is expressed as:

$\overset{\alpha}{\textbf{E}}=\frac{1}{2}\left(\overset{\alpha}{\textbf{F}}^{T}\overset{\alpha}{\textbf{F}}-\textbf{I}\right)$
(3)

Extending the above formalism, the Lagrangian strain tensor in the nonlocal
medium is given by[6, 35]:

$\mathop{\textbf{E}}^{\alpha}=\frac{1}{2}\bigr{(}\underbrace{\nabla^{\alpha}{\textbf{U}}_{X}+\nabla^{\alpha}{\textbf{U}}_{X}^{T}}_{\tilde{\bm{\epsilon}}:~{}\text{Linear
strain}}+\nabla^{\alpha}{\textbf{U}}_{X}^{T}\nabla^{\alpha}{\textbf{U}}_{X}\bigl{)}$
(4)

where $\textbf{U}(\bm{X})=\bm{x}(\bm{X})-\bm{X}$ denotes the displacement
field. The fractional gradient denoted by $\nabla^{\alpha}\textbf{U}_{X}$ is
given as $\nabla^{\alpha}\textbf{U}_{X_{ij}}=D^{\alpha}_{X_{j}}U_{i}$ and
consists of space-fractional derivatives. The space-fractional derivative
$D^{\alpha}_{\bm{X}}\textbf{U}(\bm{X},t)$ is taken according to a Riesz-Caputo
(RC) definition with order $\alpha\in(0,1)$ and it is defined on the interval
$\bm{X}\in(\bm{X}_{A},\bm{X}_{B})\subseteq\mathbb{R}^{3}$ in the following
manner:

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

where $\Gamma(\cdot)$ is the Gamma function, and
${}^{C}_{\bm{X}_{A}}D^{\alpha}_{\bm{X}}\textbf{U}$ and
${}^{C}_{\bm{X}}D^{\alpha}_{\bm{X}_{B}}\textbf{U}$ are the left- and right-
handed Caputo derivatives of U, respectively. The complete expression for the
nonlocal strain in Eq. (3) includes nonlinear terms required when accounting
for large deformations. In this expression, we also highlight the linear
component $\bm{\tilde{\epsilon}}$ of the fractional-order strain tensor that
will be employed in studies of infinitesimal deformations in nonlocal solids.
Unless otherwise specified, fractional-order strains in nonlocal solids will
refer to the linear component $\bm{\tilde{\epsilon}}$. Before proceeding, it
is worth discussing certain implications of this definition of the fractional-
order derivative. The interval of the fractional derivative
$(\bm{X}_{A},\bm{X}_{B})$ defines the horizon of nonlocality (also called
attenuation range in classical nonlocal elasticity). The length scale
parameters $\textbf{L}_{A}^{\alpha-1}$ and $\textbf{L}_{B}^{\alpha-1}$ ensure
the dimensional consistency of the deformation gradient tensor, and along with
the term $\frac{1}{2}\Gamma(2-\alpha)$ ensure the frame invariance of the
strain-displacement relations [6, 31]. These length scales are independent
parameters, which may or may not be equal to one another, hence resulting in
an asymmetric horizon of nonlocality $(\bm{X}_{A},\bm{X}_{B})$ at the point X.
The unequal values for the length scales ensure a truncation of the nonlocal
region of influence for points close to the external boundaries and
discontinuities within the solid. For complete details the reader should refer
to [6, 34]. The choice of $\textbf{L}_{A}=\textbf{L}_{B}=l_{f}$ everywhere in
the solid, where $l_{f}$ is a constant, identically reduces the above
formulation to the fractional-order kinematic relations proposed in [31].

In this formulation, nonlocality was introduced by using fractional-order
kinematic relations. The fractional-order definition of the strain has
critical implications on the thermodynamic framework for the fractional-order
model of a nonlocal continuum. In the following section, we will show that the
above approach to nonlocality allows the first and second law of
thermodynamics to be enforced in a strong (localized) sense. In other terms,
the fundamental laws of thermodynamics can be applied in a strict sense at
each point in the nonlocal continuum as opposed to what happens in classical
nonlocal approaches.

###  2.2 Thermomechanical balance laws for fractional-order thermoelasticity

####  2.2.1 Thermodynamic balance laws

In this section, we cast the fractional-order nonlocal model presented above
within a thermodynamic framework. Consider a nonlocal solid $\Omega$ that
undergoes the arbitrary motion $\bm{\Phi}$ which places a particle
$p\in\Omega$ at $\bm{y}=\bm{\Phi}(p,t)$ at time $t$. At this instant, the
solid occupies a domain $\mathcal{B}_{t}$ bounded by a surface
$\partial\mathcal{B}_{t}$. Further, consider a part of the solid
$\mathcal{P}\in\Omega$ that occupies the domain
$\mathcal{D}_{t}=\bm{\Phi}(\mathcal{P},t)$ and such that point
$p\in\mathcal{P}$.

In order to enforce the first law of thermodynamics (i.e. the conservation of
energy) over the arbitrary domain $\mathcal{D}_{t}$, we consider the energy
balance applied to the domain $\mathcal{D}_{t}$. Later, we will show that the
global balance laws hold true over any arbitrary domain $\mathcal{D}_{t}$ for
the fractional-order continuum theory. This is unlike classical (integer-
order) theories for nonlocal elasticity that allow the balance law to be
applied only over the entire domain $\mathcal{B}_{t}$, that is in a global
(weak) sense. The arbitrary choice of the domain, made possible by the
fractional-order framework, also allows $\mathcal{D}_{t}$ to be reduced to an
infinitesimal domain surrounding the point $p$ so to give the localized
version of the thermodynamic law. The application of the first law of
thermodynamics to the domain $\mathcal{D}_{t}$ allows writing the total energy
$\mathcal{E}(\mathcal{D}_{t},t)$ as the sum of the total internal energy
$\mathcal{U}(\mathcal{D}_{t},t)$ and the kinetic energy
$\mathcal{T}(\mathcal{D}_{t},t)$:

$\mathcal{E}(\mathcal{D}_{t},t)=\underbrace{\int_{\mathcal{D}_{t}}\frac{1}{2}\rho(\bm{x},t)\dot{u}_{i}(\bm{x},t)\dot{u}_{i}(\bm{x},t)~{}\mathrm{d}V}_{\mathcal{T}(\mathcal{D}_{t},t)}+\underbrace{\int_{\mathcal{D}_{t}}\rho(\bm{x},t)e(\bm{x},t)~{}\mathrm{d}V}_{\mathcal{U}(\mathcal{D}_{t},t)}$
(6)

where $\rho$ is the mass density in the current configuration,
$\dot{\textbf{u}}(=\dot{\bm{y}})$ is the spatial velocity of the particle at a
generic point $p$ defined as the time derivative of the displacement field. In
the above equation, $e(\bm{x},t)$ is the internal energy per unit mass, which
is postulated to be a function of the mechanical and thermal state variables.

Recall that, for the classical (integer-order) local elastic solid, the
internal energy of a point is a function of the integer-order strain and
entropy defined at that point. For the nonlocal solid, in addition to the
local strain energy, the internal energy at a point must also include the
energy contributions from long-range cohesive forces by other points within
the solid. Given the fractional-order kinematic relations described in §2.1,
the contribution of the additional energy from nonlocal interactions is fully
captured in the fractional-order strain $\tilde{\bm{\epsilon}}(\bm{x},t)$. It
is immediate to see that the internal energy density evaluated at a point
$\bm{x}$ is a function of the fractional-order strain
($\tilde{\bm{\epsilon}}(\bm{x},t)$). Consequently, we have the internal energy
$e=e(\tilde{\bm{\epsilon}},\tilde{\eta})$ corresponding to the thermoelastic
response defined entirely in terms of the state variables, that are the
fractional-order strain $\tilde{\bm{\epsilon}}(\bm{x},t)$ and the entropy per
unit mass $\tilde{\bm{\eta}}(\bm{x},t)$. This functional relationship is in
sharp contrast with the thermodynamic framework for classical nonlocal
approaches. More specifically, the thermodynamic framework for classical
(integer-order) nonlocality leads to
$e=e({\bm{\epsilon}},\mathcal{R}(\bm{\epsilon}),\tilde{\eta})$ where
$\bm{\epsilon}(\bm{x},t)$ is the local strain field, and
$\mathcal{R}(\bm{\epsilon})$ is a linear integral operator which models
nonlocality in the solid [39, 49, 14]. The additional functional relationship
via $\mathcal{R}(\bm{\epsilon})$ is necessary in these approaches to account
for energy contributions due to nonlocal interactions. To this regard, note
that the fractional-order strain $\tilde{\bm{\epsilon}}$ combines the local
integer-order strain $\bm{\epsilon}$ and its integral
$\mathcal{R}(\bm{\epsilon})$ into a single term (see Eq. (5)) [6]. Indeed,
this is precisely the reason that allows expressing the internal energy
density as $e=e(\tilde{\bm{\epsilon}},\tilde{\eta})$. It will be shown that
the latter observation is significant as it allows the first law of
thermodynamics to be applied in a strict sense at every point in the domain
without additional constraints. For this, we consider the first law of
thermodynamics over the arbitrary domain $\mathcal{D}_{t}$:

$\dot{\mathcal{E}}(\mathcal{D}_{t},t)=\int_{\mathcal{D}_{t}}\rho(\bm{x},t)b_{i}(\bm{x},t)\dot{u}_{i}(\bm{x},t)~{}\mathrm{d}V+\int_{\partial\mathcal{D}_{t}}\tilde{\text{t}}_{i}(\bm{x},t,\bm{n})\dot{u}_{i}(\bm{x},t)~{}\mathrm{d}A+\int_{\mathcal{D}_{t}}\rho(\bm{x},t)r(\bm{x},t)~{}\mathrm{d}V-\int_{\partial\mathcal{D}_{t}}h(\bm{x},t,\bm{n})~{}\mathrm{d}A$
(7)

where the first two terms on the right side correspond to the mechanical power
supplied to domain $\mathcal{D}_{t}$ by the body forces $\textbf{b}(\bm{x},t)$
and the surfaces forces applied per unit area
$\tilde{\textbf{t}}(\bm{x},t,\bm{n})$, respectively. Here, $\bm{n}(\bm{x},t)$
denotes the outward normal to the surface $\partial\mathcal{D}_{t}$. Note that
the surface forces on $\partial\mathcal{D}_{t}$ are applied by points external
to the domain of interest. These surface forces include interactions within
the infinitesimal neighborhood, as in (classical) local elasticity, and also
the long-range interactions of the nonlocal solid. On the contrary, the body
forces $\textbf{b}(\bm{x},t)$ applied at every point in the domain by external
sources are local as they are independent of the nonlocal interactions within
the solid. Additional terms in the above equation correspond to the rate of
change in thermal energy of the body due to internal heat generation at the
rate of $r(\bm{x},t)$ per unit volume, and the heat flux out of the body at
the rate of $h(\bm{x},t,\bm{n})$ through the surface
$\partial\mathcal{D}_{t}$. Using Eq. (6), the first law of thermodynamics can
be expressed as:

$\begin{split}\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathcal{D}_{t}}\rho(\bm{x},t)\left(\frac{1}{2}\dot{u}_{i}(\bm{x},t)\dot{u}_{i}(\bm{x},t)+e(\bm{x},t)\right)~{}\mathrm{d}V&=\int_{\mathcal{D}_{t}}\rho(\bm{x},t)b_{i}(\bm{x},t)\dot{u}_{i}(\bm{x},t)~{}\mathrm{d}V+\int_{\partial\mathcal{D}_{t}}\tilde{\text{t}}_{i}(\bm{x},t,\bm{n})\dot{u}_{i}(\bm{x},t)~{}\mathrm{d}A\\\
&+\int_{\mathcal{D}_{t}}\rho(\bm{x},t)r(\bm{x},t)~{}\mathrm{d}V-\int_{\partial\mathcal{D}_{t}}h(\bm{x},t,\bm{n})~{}\mathrm{d}A\end{split}$
(8)

The above equation corresponding to the first thermodynamic balance law has
important implications. Unlike analogous results from classical models of
nonlocal elasticity, the above energy balance law holds true for arbitrary
domain $\mathcal{D}_{t}$ within the nonlocal solid. This integral form over
the arbitrary domain $\mathcal{D}_{t}$ can be extended to develop the
expressions for a localized imposition of the first law of thermodynamics at
any point $p\in\mathcal{D}_{t}$, at any arbitrary time $t$. Recall that the
internal energy density of a nonlocal solid must include contributions from
long-range forces. For this purpose, modifications are introduced in classical
integer-order models of nonlocal elasticity via the integral operator
$\mathcal{R}(\mathbf{\epsilon})$ which restricts a localized imposition of the
first thermodynamic law [14].

Before proceeding further, we briefly discuss key characteristics of the
traction vector $\tilde{\textbf{t}}(\bm{x},t,\bm{n})$. The surface traction at
a point $\bm{x}$ is defined over an imaginary surface normal to the vector
$\bm{n}(\bm{x},t)$ and passing through the point $\bm{x}$. This vector
captures the forces acting on the surface due to the interaction of point
$\bm{x}$ with the rest of the solid. Local elasticity which considers
interactions restricted within an infinitesimal domain surrounding the point,
employs the classical Cauchy’s hypothesis to define the traction vector, that
is $\textbf{t}=\bm{n}\cdot\bm{\sigma}$, where $\bm{\sigma}$ is the classical
stress-tensor defined at the point of interest. In the fractional-order
formulation, suitable modifications to the above definition are required to
account for additional forces acting on the surface due to long-range
interactions. For the fractional-order approach, we present the following
generalization of the Cauchy’s postulate:

The surface traction $\tilde{\textbf{t}}(\bm{x},t,\bm{n})$ acting on an
imaginary surface $\partial\mathcal{D}_{t}$ perpendicular to the normal vector
$\bm{n}(\bm{x},t)$ and passing through the point $\bm{x}$ is given as:

$\tilde{\textbf{t}}(\bm{x},t,\bm{n})={\bm{I}}_{\bm{n}}^{1-\alpha}\cdot\tilde{\bm{\sigma}}(\bm{x},t)$
(9)

where $\tilde{\bm{\sigma}}$ is the nonlocal stress evaluated at the point
$\bm{x}$. The integral operator ${\bm{I}}_{\bm{n}}^{1-\alpha}$ defined over
the horizon of influence of point $\bm{x}$ is given as:

$\bm{I}^{1-\alpha}_{\bm{n}}(\tilde{\bm{\sigma}})={I}^{1-\alpha}_{x_{1}}(n_{1}\cdot\tilde{\bm{\sigma}})\hat{e}_{1}+{I}^{1-\alpha}_{x_{2}}(n_{2}\cdot\tilde{\bm{\sigma}})\hat{e}_{2}+{I}^{1-\alpha}_{x_{3}}(n_{3}\cdot\tilde{\bm{\sigma}})\hat{e}_{3}$
(10)

such that $\bm{n}=n_{i}\hat{e}_{i}$ ($i=1,2,3$), $\hat{e}_{i}$ are the
orthonormal basis vectors. Further, ${I}^{1-\alpha}_{i}(\cdot)$ is a Riesz-
type fractional integral defined as:

$I^{1-\alpha}_{x_{i}}\chi=\frac{1}{2}\Gamma(2-\alpha)\left[l_{B_{i}}^{\alpha-1}\left({}_{{x_{i}}-l_{B_{i}}}I_{{x_{i}}}^{1-\alpha}\chi\right)+l_{A_{i}}^{\alpha-1}\left({}_{{x_{i}}}I_{{x_{i}}+l_{A_{i}}}^{1-\alpha}\chi\right)\right]$
(11)

where ${}_{x_{i}-l_{B_{i}}}I_{{x_{i}}}^{1-\alpha}\chi$ and
${}_{{x_{i}}}I_{{x_{i}}+l_{A_{i}}}^{1-\alpha}\chi$ are the left and right
fractional integrals to the order $\alpha\in(0,1)$ of an arbitrary function
$\chi$, in the $\hat{e}_{i}$ direction. In this study, the definition of the
normal vector is restricted to orthonormal triads of Cartesian coordinates.
This is due to the limited developments in fractional vector calculus, as
discussed in [50, 51]. A detailed discussion on the above expression for the
traction vector in a nonlocal solid is outlined in the Appendix.

Analogous to classical local elasticity, the traction vector defined in Eq.
(9) is used within the global balance laws (Eq. (8)) to obtain the
corresponding strong form. Additionally, we define the rate of heat transfer
across the surface as $h(\bm{x},t,\bm{n})=\bm{q}(\bm{x},t)\cdot\bm{n}$, where
the vector field $\bm{q}$ is referred to as the heat flux out of the body. By
substituting the definitions of the nonlocal traction and the heat flux into
Eq. (8), and then applying the divergence theorem gives:

$\int_{\mathcal{D}_{t}}\left[\rho\frac{De}{Dt}+\rho(\ddot{{u}}_{i}-b_{i})\dot{u}_{i}-\left(\bm{I}^{1-\alpha}_{n_{i}}\tilde{\sigma}_{ij}\dot{u}_{j}\right)_{,i}-\left(\rho
r-q_{i,i}\right)\right]\mathrm{d}V=0$ (12)

Here, $D\cdot/Dt$ is the material derivative evaluated in the Eulerian
configuration. By assuming continuous field variables and using the
localization lemma we obtain:

$\rho\frac{De}{Dt}+\rho(\ddot{{u}}_{i}-b_{i})\dot{u}_{i}-\left(\bm{I}^{1-\alpha}_{n_{i}}\tilde{\sigma}_{ij}\dot{u}_{j}\right)_{,i}=\rho
r-q_{i,i}~{}~{}\forall\bm{x}\in\mathcal{D}_{t}$ (13)

The above result clearly illustrates the possibility of enforcing the first
law of thermodynamics (i.e. the energy balance) at any arbitrary point
$p\in\mathcal{D}_{t}$ of the nonlocal solid. Note that, according to the
fractional-order continuum theory, we do not require additional constraints to
be imposed over the state variables. This is unlike the previous studies based
on integer-order theories of nonlocal elasticity which required the energy
balance to be complemented by an additional condition on the nonlocal
residuals. Thus, a localized imposition of the energy balance law for nonlocal
solids is derived here following the fractional-order models. In order to
simplify the above expression to a more recognizable form, we use the
mechanical balance laws for linear and angular momentum of the nonlocal solid,
derived later in §2.2.2 (refer Eqs. (20) and (23)). Employing the mechanical
balance laws for the nonlocal solid, and following standard algebraic
operations, we write:

$\rho\frac{De}{Dt}=\tilde{{\sigma}}_{ij}\dot{\tilde{\epsilon}}_{ij}+\rho
r-q_{i,i}~{}~{}\forall\bm{x}\in\mathcal{B}_{t}$ (14)

Note that the strain tensor $\tilde{\epsilon}_{ij}$ in the above expression is
evaluated using fractional-order derivatives defined in Eq. (2). We emphasize
again that the above localized form obtained for the fractional-order approach
is in net contrast with classical nonlocal approaches where the conservation
of the first law can only be applied in a weak sense (see, for example, [14]).

Next, we apply the second law of thermodynamics to the fractional-order
continuum model. Recall the Clausius-Duhem inequality applied over the entire
solid $\mathcal{B}_{t}$ states:

$\frac{d}{dt}\int_{\mathcal{B}_{t}}\rho(\bm{x},t)\tilde{\eta}_{\text{int}}(\bm{x},t)\mathrm{d}V=\frac{d}{dt}\int_{\mathcal{B}_{t}}\rho(\bm{x},t)\tilde{\eta}(\bm{x},t)\mathrm{d}V-\int_{\mathcal{B}_{t}}\frac{\rho(\bm{x},t)r(\bm{x},t)}{T(\bm{x},t)}+\int_{\partial\mathcal{B}_{t}}\frac{h(\bm{x},t,\bm{n}(\bm{x},t))}{T(\bm{x},t)}\mathrm{d}A\geq
0$ (15)

where $T$ denotes the temperature of the solid, and
$\tilde{\eta}_{\text{int}}$ is the internal entropy production density. Recall
that the localized form of the second law of thermodynamics states that the
internal entropy production rate is non-negative for all points inside the
solid, that is $D{\tilde{\eta}}_{\text{int}}/Dt\geq
0~{}\forall~{}\bm{x}\in\mathcal{B}_{t}$. Classical approaches to nonlocal
thermoelasticity satisfy this inequality only in a weak sense, that is the
integral form given in the above equation [13, 39, 49]. A detailed discussion
of the consequent physical inconsistencies can be found in [14].

In analogy with the classical approach, we introduce the Legendre
transformation $\psi=e-T\tilde{\eta}$, where $\psi$ denotes the Helmholtz free
energy per unit mass. It follows that $\psi=\psi(\tilde{\bm{\epsilon}},T)$
which is different from classical nonlocal approaches wherein
$\psi=\psi(\bm{\epsilon},\mathcal{R}(\bm{\epsilon}),T)$ [14]. By using the
Legendre transformation along with Eqs. (14) and (15), we obtain the following
local form of the second law of thermodynamics for the nonlocal solid:

$\rho
T\frac{D{\tilde{\eta}}_{\text{int}}}{Dt}=\tilde{\sigma}_{ij}\dot{\tilde{\epsilon}}_{ij}-\rho\frac{D\psi}{Dt}-\rho\eta\dot{T}-T_{,i}\frac{q_{i}}{T}\geq
0$ (16)

Remarkably, the above fractional-order inequality matches, in its functional
form, the classical Clausius-Duhem inequality. Equation (16) also highlights a
clear difference compared with classical nonlocal formulations that require
additional terms within the inequality as a result of the functional
dependence of $\psi$ on $\mathcal{R}(\bm{\epsilon})$. As discussed in [14],
these additional terms within the inequality disappear only when a weak form
is considered. However, as mentioned previously, satisfying the second law of
thermodynamics only in a weak sense leads to nonphysical results.

Thus, it appears that the fractional-order continuum theory for nonlocal
solids allows the thermodynamic balance laws to be applied in a localized
form; a key observation in order to establish the thermodynamic consistency of
the fractional-order continuum theory and to derive rigorous constitutive
models for the nonlocal solid. Note that the thermodynamic balance principles
for the nonlocal solid have been presented in the Eulerian setting over the
current configuration $\mathcal{B}_{t}$ of the solid at time instant $t$.
However, assuming small deformation, the above results can be extended for the
domain in reference configuration $\mathcal{B}_{0}$. For a generalized study
of large deformations, resulting from both geometric and material
nonlinearities simultaneously, this assumption will not be valid. In such
cases, the Lagrangian analogues for the governing equations can be derived
from the above results[52].

####  2.2.2 Mechanical balance laws

In the above discussion on thermodynamic balance laws, we have obtained the
energy balance and entropy inequality laws for a fractional-order nonlocal
solid. For the sake of completeness, in the following we also derive the
mechanical balance laws corresponding to conservation of linear and angular
momentum. For this purpose, we continue the discussion taking the point of
view of the current configuration $\mathcal{B}_{t}$. Similar to the procedure
outlined in §2.2.1, we begin with the integral balance laws over an arbitrary
domain $\mathcal{D}_{t}$ and derive the localized forms of the governing
equations for a point $p\in\mathcal{D}_{t}$. For the domain under
consideration, the statement for the balance of linear momentum is:

$\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathcal{D}_{t}}\rho(\bm{x},t)\dot{\textbf{u}}(\bm{x},t)~{}\mathrm{d}V=\int_{\mathcal{D}_{t}}~{}\rho(\bm{x},t)\textbf{b}(\bm{x},t)\mathrm{d}V+\int_{\partial\mathcal{D}_{t}}~{}\tilde{\textbf{t}}(\bm{x},t,\bm{n})~{}\mathrm{d}A$
(17a) and the balance of angular momentum in the absence of external couples
is expressed as:
$\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathcal{D}_{t}}\bm{y}\times\rho(\bm{x},t)\dot{\textbf{u}}(\bm{x},t)~{}\mathrm{d}V=\int_{\mathcal{D}_{t}}~{}\bm{y}\times\rho(\bm{x},t)\textbf{b}(\bm{x},t)\mathrm{d}V+\int_{\partial\mathcal{D}_{t}}~{}\bm{y}\times\tilde{\textbf{t}}(\bm{x},t,\bm{n})~{}\mathrm{d}A$
(17b)

where the operator ${}^{\prime}\times^{\prime}$ denotes the exterior product.
First, the above balance laws are simplified by imposing the classical
(integer-order) result for the conservation of mass [52]. Thereafter, by
substituting the expression for surface traction $\tilde{\textbf{t}}$ given in
Eq. (9) for the fractional-order solid, we obtain:

$\int_{\mathcal{D}_{t}}\rho(\bm{x},t)\ddot{\textbf{u}}(\bm{x},t)~{}\mathrm{d}V=\int_{\mathcal{D}_{t}}~{}\rho(\bm{x},t)\textbf{b}(\bm{x},t)\mathrm{d}V+\int_{\partial\mathcal{D}_{t}}~{}{\bm{I}}_{\bm{n}}^{1-\alpha}\cdot\tilde{\bm{\sigma}}(\bm{x},t)~{}\mathrm{d}A$
(18a)
$\int_{\mathcal{D}_{t}}\bm{y}\times\rho(\bm{x},t)\ddot{\textbf{u}}(\bm{x},t)~{}\mathrm{d}V=\int_{\mathcal{D}_{t}}~{}\bm{y}\times\rho(\bm{x},t)\textbf{b}(\bm{x},t)\mathrm{d}V+\int_{\partial\mathcal{D}_{t}}~{}\bm{y}\times{\bm{I}}_{\bm{n}}^{1-\alpha}\cdot\tilde{\bm{\sigma}}(\bm{x},t)~{}\mathrm{d}A$
(18b)

where we use: $\dot{\bm{y}}=\dot{\textbf{u}}$. Using the definition of the
Reisz-integral within the definition of the surface traction and applying the
divergence theorem we obtain [51]:

$\int_{\mathcal{D}_{t}}\rho(\bm{x},t)\ddot{{u}}_{j}(\bm{x},t)~{}\mathrm{d}V=\int_{\mathcal{D}_{t}}~{}\rho(\bm{x},t){b}_{j}(\bm{x},t)\mathrm{d}V+\int_{\partial\mathcal{D}_{t}}~{}\left({{I}}_{n_{i}}^{1-\alpha}\tilde{{\sigma}}_{ij}(\bm{x},t)\right)_{,i}~{}\mathrm{d}V$
(19a)
$\int_{\mathcal{D}_{t}}\varepsilon_{ijk}{y}_{j}~{}\rho(\bm{x},t)\ddot{{u}}_{k}(\bm{x},t)~{}\mathrm{d}V=\int_{\mathcal{D}_{t}}~{}\varepsilon_{ijk}~{}{y}_{j}\rho(\bm{x},t){b}_{k}(\bm{x},t)\mathrm{d}V+\int_{\partial\mathcal{D}_{t}}~{}\varepsilon_{ijk}~{}\left({y}_{j}{\bm{I}}_{n_{m}}^{1-\alpha}\tilde{{\sigma}}_{mk}(\bm{x},t)\right)_{,m}~{}\mathrm{d}V$
(19b)

Applying the localization lemma to Eq. (19a), we arrive at the following
result:

$\rho(\bm{x},t)\ddot{{u}}_{j}(\bm{x},t)=\rho(\bm{x},t){b}_{j}(\bm{x},t)\mathrm{d}V+\mathfrak{D}_{x_{i}}^{\alpha}\tilde{{\sigma}}_{ij}(\bm{x},t),~{}~{}~{}~{}\forall\bm{x}\in\mathcal{B}_{t}$
(20)

which is the localized linear momentum balance law for fractional-order
solids. Here, $\mathfrak{D}_{x_{i}}^{\alpha}(\cdot)$ is the Riesz Riemann-
Liouville derivative of order $\alpha$ which is defined as [34]:

$\mathfrak{D}^{\alpha}_{x_{i}}\chi=\frac{1}{2}\Gamma(2-\alpha)\left[l_{B}^{\alpha-1}\left({}^{~{}~{}~{}RL}_{x_{i}-l_{B}}D^{\alpha}_{x_{i}}\chi\right)-l_{A}^{\alpha-1}\left({}^{RL}_{x_{i}}D^{\alpha}_{x_{i}+l_{A}}\chi\right)\right]$
(21)

where $\chi$ is an arbitrary function and
${}^{~{}~{}~{}RL}_{x_{i}-l_{B}}D^{\alpha}_{x_{i}}\chi$ and
${}^{RL}_{x_{i}}D^{\alpha}_{x_{i}+l_{A}}\chi$ are the left- and right-handed
Riemann Liouville derivatives of $\chi$ to the order $\alpha$, respectively.
In the derivation of the localized form of the linear momentum balance law,
given in Eq. (20) from Eq. (19a), we used the following relation:

$\mathfrak{D}^{\alpha}_{x_{i}}\chi=\frac{\mathrm{d}}{\mathrm{d}x_{i}}\left[I^{1-\alpha}_{i}\chi\right]$
(22)

The governing equations of motion derived here using the balance principles
agree very well with the results obtained via variational principles [51].
Using the above result for linear momentum balance and $y_{i,j}=\delta_{ij}$
($\delta_{ij}$ is the Krönecker delta), the strong form of the angular
momentum balance for fractional-order solids reduces to the symmetry condition
of the nonlocal stress tensor:

$\tilde{\sigma}_{ij}=\tilde{\sigma}_{ji},~{}~{}~{}~{}\forall\bm{x}\in\mathcal{B}_{t}$
(23)

It is interesting to note the fractional-order Riesz-type Riemann-Liouville
fractional-order derivative (divergence) of the nonlocal stress tensor in the
elastodynamic equation (Eq. 20). This result differ from existing studies on
fractional-order continuum theories that employ a first-order derivative of
the nonlocal stress ($D^{1}_{x_{i}}\tilde{\sigma}_{ij}$) tensor in the same
equation (see, for example, [31]). This difference arises due to the
consideration of long-range interactions within the expression for traction
given in Eq. (9). The additional contributions to the surface traction at a
point in the nonlocal solid follow from the long-range interactions (see
Appendix). Note that the Riesz-type Riemann-Liouville operator in the
mechanical governing equations given in Eq. (20) is self-adjoint and the
system is positive-definite. This result was established and proved in [34].
These observations are clearly in contrast with either classical integer-order
approaches to nonlocal elasticity, that have shown it is not possible to
define a self-adjoint quadratic potential energy [53, 54], or some fractional-
order models [31, 32]. Given the self-adjoint and positive-definite nature of
our formulation, the resulting system of equations is well-posed and admit a
unique solution [34]. We will show in §5 that this well-posedness results in a
consistent softening behavior of the structure with increasing degree of
nonlocality irrespective of the thermomechanical load distributions and
boundary conditions. This result is significant because it bypasses a key
inconsistency observed in classical nonlocal models and associated with the
non self-adjointness of the operators [53, 54].

###  2.3 Constitutive framework for fractional-order thermoelasticity

####  2.3.1 Constitutive modeling

The inequality in Eq. (16) is used to derive the thermodynamically-consistent
constitutive equations for fractional-order nonlocal elasticity. By
substituting the expression for the time derivative of the Helmholtz free
energy, the inequality in Eq. (16) is expressed as:

$T\dot{\tilde{\eta}}_{\text{int}}=\left(\tilde{\sigma}_{ij}-\rho\frac{\partial\psi}{\partial\tilde{\epsilon}_{ij}}\right)\dot{\tilde{\epsilon}}_{ij}-\rho\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 all thermoelastic processes as well
as for arbitrary choices of the independent fields
$\dot{\tilde{\epsilon}}_{ij}$ and $\dot{{T}}$, we obtain the following
constitutive laws:

$\tilde{\sigma}_{ij}={\rho}\frac{\partial\psi}{\partial\dot{\tilde{\epsilon}}_{ij}},~{}~{}\forall~{}\bm{x}\in\mathcal{B}_{t}$
(25a) $\tilde{\eta}=-\frac{\partial\psi}{\partial
T},~{}~{}\forall~{}\bm{x}\in\mathcal{B}_{t}$ (25b)

Under assumptions of linear elasticity, the above equations are the fractional
analogues for the Duhamel-Neumann’s laws for classical thermoelasticity.
Further, by using the above constitutive relations within Eq. (16), the
inequality reduces to:

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

which establishes the second law of thermodynamics for a fractional order
nonlocal solid. The relations in Eqs. (25,26) can be expressed as:

Theorem: The constitutive relations for fractional-order nonlocal
thermoelasticity do not violate the Clausius-Duhem inequality if they are of
the form given in Eq. (25) and subject to Eq. (26).

A few additional comments on this thermodynamic framework are needed. At a
first glance, the form of the stress-strain constitutive relation in Eq. (25a)
might be deceiving as it appears to lead to a classical constitutive relation.
Although this is, formally, a correct statement it does not entirely capture
the nature of Eq. (25a). As highlighted earlier, nonlocality was modeled using
fractional-order kinematic relations given in Eq. (4). Therefore, the stress
defined through the Eq. (25a) is also nonlocal in nature. In addition, this
construction of nonlocality (i.e. based on fractional-order kinematic
relations) allows the application of variational principles, ensures a
quadratic form of the potential energy of the system, and leads to well-posed
nonlocal governing equations [34, 36].

Note also that, in the above study, we assumed an integer-order Fourier heat
conduction law, that is $q_{i}=-kT_{,i}$, where $k$ is the material
conductivity constant such that $k>0$. It is immediate that the inequality in
Eq. (26) obtained from the second law is trivially satisfied for the integer-
order heat conduction law. Finally, following from the latter remark, a space-
fractional thermal conduction law defined as
$q_{i}=-k{D}_{X_{i}}^{\alpha}\theta$ was proposed in [44]. However, as shown
in [47], the space-fractional heat conduction law violates the second law of
thermodynamics. Thus, we limit the scope of the current study to fractional-
order constitutive modeling for the mechanical field but integer-order models
for the thermal fields.

####  2.3.2 Linear fractional-order thermoelasticity

In this section, we derive the constitutive relations for linear thermoelastic
response of fractional-order nonlocal solids. For this purpose, recalling that
the fractional-order nonlocal formulation allows a localized implementation of
the thermodynamic principles, we write the Helmholtz free energy density for
the thermoelastic response following the typical approach for local elasticity
(albeit by using the fractional-order strain). Further, we use the Helmholtz
free energy density to construct the material constitutive relations for
nonlocal thermoelasticity. In this study on linear elastic behavior of
nonlocal solids using fractional-order theories of thermoelasticity, we extend
the linear material constitutive relations to the geometrically nonlinear
response of the fractional-order nonlocal solids (Saint Venant-Kirchhoff
material model). Therefore, assuming small deformations, the distinction
between reference and current configurations for the domain $\Omega$ vanishes.
This implies the constitutive relations given in Eq. (25) for domain $\Omega$
hold true in all configurations ($\mathcal{B}_{0}$ & $\mathcal{B}_{t}$).

We can cast the constitutive relations in Eq. (25) in the form:

$\tilde{\sigma}_{ij}=\frac{\partial\mathcal{W}}{\partial\tilde{\epsilon}_{ij}},~{}~{}~{}~{}~{}~{}\tilde{\eta}=-\rho_{0}^{-1}\frac{\partial\mathcal{W}}{\partial
T},~{}~{}~{}\forall\bm{X}\in\Omega$ (27)

where the total free energy $\mathcal{W}=\rho_{0}\psi$ is expressed in terms
of the mass density $\rho_{0}$ in reference configuration. We write the free
energy for an isotropic material as a series expansion of the fractional-order
strain $\tilde{\epsilon}_{ij}$ and the temperature difference
$\theta=T-T_{0}$, which is the difference between the temperature $T$ at any
point within the continuum and the uniform ambient temperature $T_{0}$ at the
reference state. The free energy, expanded in power series with respect to the
strain invariants and the temperature difference, is given as [55]:

$\begin{split}\mathcal{W}=a_{0}+a_{1}\tilde{J}_{1}+a_{2}\tilde{J}_{2}+a_{3}\tilde{J}_{3}+a_{4}\theta+a_{5}\tilde{J}_{1}^{2}+a_{6}\tilde{J}_{2}^{2}+a_{7}\tilde{J}_{3}^{2}+a_{8}\tilde{J}_{1}\tilde{J}_{2}&\\\
+a_{9}\tilde{J}_{1}\tilde{J}_{3}+a_{10}\tilde{J}_{2}\tilde{J}_{3}+a_{11}\theta^{2}+a_{12}\tilde{J}_{1}\theta+a_{13}\tilde{J}_{2}\theta+a_{14}\tilde{J}_{3}\theta+\text{h.o.t}&\end{split}$
(28a) where $a_{k}$ are material constants and
$\tilde{J}_{1}=\tilde{\epsilon}_{ii};~{}~{}~{}\tilde{J}_{2}=\frac{1}{2}\left(\tilde{\epsilon}_{ii}\tilde{\epsilon}_{jj}-\tilde{\epsilon}_{ij}\tilde{\epsilon}_{ij}\right);~{}~{}~{}\tilde{J}_{3}=\text{det}(\tilde{\epsilon}_{ij})$
(28b)

are the invariants of the nonlocal strain tensor $\tilde{\bm{\epsilon}}$.
Assuming that the solid is stress free in the undeformed state and the free
energy is restricted to linear isotropic thermoelasticity (i.e. ignoring the
higher order terms in Eq. (28a)), we obtain the following expression for
$\mathcal{W}$:

$\mathcal{W}=a_{2}\tilde{J}_{2}+a_{5}\tilde{J}_{1}^{2}+a_{11}\theta^{2}+a_{12}\tilde{J}_{1}\theta$
(29)

where the material constants are given as [56]:

$a_{2}=-2\mu;~{}~{}~{}a_{5}=\frac{1}{2}\left(\lambda+2\mu\right);~{}~{}~{}a_{11}=-\frac{\rho_{0}C_{v}^{0}}{2T_{0}};~{}~{}~{}a_{12}=(3\lambda+2\mu)\alpha_{0}$
(30)

The material constants $\lambda$ and $\mu$ are the isothermal Lamé parameters
for the isotropic solid, $\alpha_{0}$ is the coefficient of volumetric thermal
expansion and $C_{v}^{0}$ is the specific heat at constant strain. The
parameters $\alpha_{0}$ and $C_{v}^{0}$ are all defined in the reference state
at $T_{0}$. Thus, the Helmholtz free energy density (per unit volume) for the
thermoelastic response of a nonlocal isotropic solid is given by [57, 56]:

$\mathcal{W}=\frac{1}{2}\lambda\tilde{\epsilon}_{kk}\tilde{\epsilon}_{ll}+\mu\tilde{\epsilon}_{ij}\tilde{\epsilon}_{ij}-(3\lambda+2\mu)\alpha_{0}\tilde{\epsilon}_{kk}\theta-\frac{\rho_{0}C_{v}^{0}}{2T_{0}}\theta^{2}$
(31)

By using the above expression for $\psi$ together with Eq. (25), the
thermoelastic constitutive relations relating the different physical
quantities for the isotropic solid are obtained as:

$\tilde{\sigma}_{ij}=\lambda\delta_{ij}\tilde{\epsilon}_{kk}+2\mu\tilde{\epsilon}_{ij}-(3\lambda+2\mu)\alpha_{0}\delta_{ij}\theta$
(32a)
$\tilde{\eta}=\rho_{0}^{-1}(3\lambda+2\mu)\alpha_{0}\tilde{\epsilon}_{kk}+\frac{C_{v}^{0}}{T_{0}}\theta$
(32b)

Using the above results, the Helmholtz free energy in Eq. (31) is recast in
the following manner:

$\mathcal{W}=\frac{1}{2}~{}\tilde{\sigma}_{ij}~{}\tilde{\epsilon}_{ij}-\frac{1}{2}~{}\rho_{0}~{}\tilde{\eta}~{}\theta$
(33)

##  3 Thermoelastic Euler-Bernoulli nonlocal beam model

In this section, we use the thermoelastic constitutive relations developed for
the nonlocal solid to analyze the thermoelastic response of a fractional-order
Euler-Bernoulli beam. Building on [34, 35], we derive the geometrically
nonlinear governing equations and the corresponding boundary conditions for
the thermoelastic boundary value problem (BVP) using variational principles.

###  3.1 Geometrically nonlinear constitutive relations

Consider a nonlocal beam subject to distributed transverse mechanical and
thermal loads as illustrated in Fig. (1). As indicated in the schematic, the
Cartesian coordinates for the current study are chosen such that $x_{3}=\pm
h/2$ coincides with the top and bottom surfaces of the beam, and $x_{1}=0$ and
$x_{1}=L$ are the ends of the beam along the longitudinal direction. The
surface $x_{3}=0$ coincides with the mid-plane of the beam and the origin of
the reference frame is chosen at the intersection of the mid-plane with the
left-end of the beam.

Figure 1: Schematic of an elastic beam subject to distributed transverse
mechanical load $F_{t}(x_{1})$ and thermal load $T(x_{1})$.

The axial and transverse components of the displacement field
$\textbf{u}(x_{1},x_{3})$ are denoted by $u_{1}(x_{1},x_{3})$ and
$u_{3}(x_{1},x_{3})$, respectively. These displacement fields are given by the
Euler-Bernoulli theory as:

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

where $u_{0}(x_{1})$ and $w_{0}(x_{1})$ are the mid-plane axial and transverse
displacements, respectively. For a geometrically nonlinear analysis, assuming
moderate rotations ($10^{\circ}-15^{\circ}$) but small strains, the
fractional-order Lagrangian strain tensor in Eq. (4) can be further simplified
using von-Kármán relations. The resulting fractional-order von-Kármán strain-
displacement relations are given as [35]:

$\tilde{\epsilon}_{11}(x_{1},x_{3})=D_{x_{1}}^{\alpha}u_{1}(x_{1})+\frac{1}{2}\left[D_{x_{1}}^{\alpha}u_{3}(x_{1})\right]^{2}$
(35)

where $D^{\alpha}_{x_{1}}(\cdot)$ denotes the fractional-order RC derivative
defined in Eq. (5). Note that $\tilde{\epsilon}_{11}$ here is nonlinear in
transverse displacements. As discussed previously, for the RC derivative used
above, the 1D domain $(\bm{x}_{A},\bm{x}_{B})$ along the mid-plane of the beam
is the horizon of nonlocal interaction at $\bm{x}(x_{1},0)$. The end-points of
the nonlocal horizon $\bm{x}_{A}(x_{A_{1}},0)$ and $\bm{x}_{B}(x_{B_{1}},0)$
are the terminals of the left- and right-handed Caputo derivatives within the
RC derivative. It follows from Eq. (5) that $l_{A}=x_{1}-x_{A_{1}}$ and
$l_{B}=x_{B_{1}}-x_{1}$ are the length scales along $\hat{\bm{x}}_{1}$ to the
left and right hand sides of the point $\bm{x}(x_{1},0)$, respectively. For
the Euler-Bernoulli beam displacement field given in Eq. (34), a non-zero
expression for the transverse shear strain would be obtained. However, for the
slender beam assumed here the rigidity to transverse shear deformation is much
higher when compared to its bending rigidity. Therefore, we neglect the
contribution of the transverse shear deformation towards the deformation
energy of the fractional-order nonlocal solid in the subsequent analysis.

By combining the fractional-order nonlinear axial strain in the above equation
along with the Euler-Bernoulli displacement field given in Eq. (34), the axial
strain can be recast as:

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

In the above equation, $\tilde{\epsilon}_{0}(x_{1})$ and
$\tilde{\kappa}(x_{1})$ denote the fractional-order axial and bending strains,
respectively. They are expressed in terms of the mid-plane field variables as:

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

Note that we make use of the von-Kármán definition for geometrically nonlinear
strains. Following this approach, large deformations are considered only in
the transverse direction. For studies based on this definition of
geometrically nonlinear strains, the linear elastic constitutive relations
developed in §2 can still be employed [58]. Therefore, the axial stress in the
nonlocal isotropic solid subject to thermoelastic loads may be written from
Eq. (32a) to be:

$\tilde{\sigma}_{11}(x_{1},x_{3})=E\big{(}\tilde{\epsilon}_{11}(x_{1},x_{3})-\alpha_{0}\theta(x_{1},x_{3})\big{)}$
(38)

where $E$ is the Young’s modulus of the isotropic solid and $\alpha_{0}$ is
the coefficient of thermal expansion for the isotropic solid, as defined
earlier. Using the above defined fractional-order strain and stress fields,
the deformation energy $\mathcal{U}$ of the nonlocal beam is obtained as:

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

where $\Omega$ denotes the volume occupied by the beam. The total potential
energy functional of the beam subject to distributed axial ($F_{a}(x_{1})$)
and transverse forces ($F_{t}(x_{1})$) acting on the mid-plane, assuming no
body forces, is given by:

$\Pi[\textbf{u}(\bm{x})]=\mathcal{U}-\underbrace{\int_{0}^{L}F_{a}(x_{1})u_{0}(x_{1})\mathrm{d}x_{1}}_{\text{Work
done by axial
loads}}-\underbrace{\int_{0}^{L}F_{t}(x_{1})w_{0}(x_{1})\mathrm{d}x_{1}}_{\text{Work
done by transverse loads}}$ (40)

We now derive the governing equations and the associated boundary conditions
for the thermoelastic response of the nonlocal beam in the strong form by
imposing optimality conditions on the above functional. While we derived the
governing equations for the fractional-order solid in §2.2.2 following the
balance principles for intuitive purposes, here we employ the variational
methods due to its ability to treat geometrically non-linear systems. Before
presenting the governing equations, we highlight that the objective of this
study is to evaluate the elastic response of a 1D beam when subject to
combined thermal and mechanical loads. The thermal load consists of a steady-
state temperature distribution applied along the length of the beam on the
face at $x_{3}=\pm h/2$. Thus, the independent variation of temperature field
$\delta T$, and thereby $\delta\theta$, is identically zero.

###  3.2 Governing equations

The fractional-order governing equations for the thermoelastic response of
geometrically nonlinear and nonlocal beams are obtained using variational
principles (i.e. by minimizing the total potential energy given in Eq. (40)).
They are given as follows:

$\mathfrak{D}^{\alpha}_{x_{1}}\mathcal{N}(x_{1})+F_{a}(x_{1})=0~{}~{}\forall~{}x_{1}\in(0,L)$
(41a)
$D^{1}_{x_{1}}\left[\mathfrak{D}^{\alpha}_{x_{1}}\mathcal{M}(x_{1})\right]+\mathfrak{D}^{\alpha}_{x_{1}}\left[\mathcal{N}(x_{1})D_{x_{1}}^{\alpha}\left[w_{0}(x_{1})\right]\right]+F_{t}(x_{1})=0~{}~{}\forall~{}x_{1}\in(0,L)$
(41b)

The corresponding essential and natural boundary conditions are obtained as:

$\mathcal{N}(x_{1})=0~{}~{}\text{or}~{}~{}\delta
u_{0}(x_{1})=0~{}~{}\forall~{}x_{1}\in\\{0,L\\}$ (42a)
$\mathcal{M}(x_{1})=0~{}~{}\text{or}~{}~{}\delta\left[D^{1}_{x_{1}}w_{0}(x_{1})\right]=0~{}~{}\forall~{}~{}x_{1}\in\\{0,L\\}$
(42b)
$D^{1}_{x_{1}}\mathcal{M}(x_{1})+\mathcal{N}(x_{1})D^{1}_{x_{1}}\left[w_{0}(x_{1})\right]=0~{}~{}\text{or}~{}~{}\delta
w_{0}(x_{1})=0~{}~{}\forall~{}x_{1}\in\\{0,L\\}$ (42c)

Note that the detailed steps leading to the above fractional-order nonlinear
governing equations extend directly from the geometrically nonlinear analysis
of fractional-order beams presented in [35], hence they are not provided here.
In the above Eqs. (41,42), $D^{1}_{x_{1}}(\cdot)$ denotes the first integer-
order derivative with respect to the axial variable $x_{1}$. Note that the
fractional derivative $\mathfrak{D}^{\alpha}_{x_{1}}(\cdot)$ is defined over
the interval $(x_{1}-l_{B},x_{1}+l_{A})$ unlike the fractional derivative
$D^{\alpha}_{x}(\cdot)$ which is defined over the interval
$(x_{1}-l_{A},x_{1}+l_{B})$. This change in the terminals of the interval of
the Riesz Riemann-Liouville fractional derivative follows from the standard
integration by parts technique used to simplify the variational integrals (see
[34]). Further, $\mathcal{N}(x_{1})$ and $\mathcal{M}(x_{1})$ are axial and
bending stress resultants defined in the following manner:

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

By using the constitutive relations for a homogeneous isotropic solid given in
Eq. (38) along with the above definitions, the stress resultants are obtained
as:

$\mathcal{N}(x_{1})=A_{11}\tilde{\epsilon}_{0}(x_{1})-N_{\theta}(x_{1})$ (44a)
$\mathcal{M}(x_{1})=-D_{11}\tilde{\kappa}(x_{1})-M_{\theta}(x_{1})$ (44b)

where $A_{11}=Ebh$ and $D_{11}=Ebh^{3}/12$ are the axial and bending stiffness
coefficients of the beam, respectively. The thermal resultants
$N_{\theta}(x_{1})$ and $M_{\theta}(x_{1})$ for the isotropic beam are given
as:

$\left\\{N_{\theta}(x_{1}),M_{\theta}(x_{1})\right\\}=Eb\alpha_{0}\int_{-h/2}^{h/2}~{}\\{1,~{}x_{3}\\}~{}\theta(x_{1},x_{3})~{}\mathrm{d}x_{3}$
(45)

Note that for a general distribution of material properties across the
thickness of the beam, additional terms due to the bending-extension coupling
would be noted in the Eqs. (43,45). Upon ignoring the nonlinear terms in the
governing equations given above, striking similarity may be noted to the
governing equations in Eq. (20) derived using linear momentum balance law.

In the following, we discuss a few characteristics of the thermomechanical
governing equations given in Eq. (41). First, observe that the stress
resultants given in Eq. (44) introduce the thermoelastic variables into the
governing equations in Eq. (41). In the absence of thermal loads
($\theta(x_{1},x_{3})=0$), the constitutive models reduce to the expressions
derived in [35] for the geometrically nonlinear fractional-order nonlocal
beams. Owing to the nonlinear nature of the structural response, the axial and
transverse displacement fields are coupled unlike what seen in the linear
elastic case [34]. Second, we emphasize that the fractional-order model for
the linear thermoelastic response of a nonlocal beam can be obtained by
ignoring the nonlinear terms in the constitutive relations developed above.
The linear thermoelastic model of the fractional-order nonlocal beam will be
discussed further in §5. Finally, the classical thermoelastic models are
recovered for $\alpha=1$.

##  4 Nonlinear fractional finite element model (f-FEM)

Given the nonlinear and integro-differential nature of the governing
equations, it is unlikely to obtain closed form solutions for the most general
loading and boundary conditions. Therefore, we employ a fractional-order
finite element method to obtain the numerical solution of the nonlinear
governing equations. The f-FEM developed to solve the thermomechanical
fractional-order BVP builds upon the numerical solvers developed for
fractional-order models of nonlocal elasticity [34, 35]. Note that, although
the f-FEM is developed and applied for thermoelastic response of fractional-
order beams, it can be easily extended to higher dimensional structures like
plates and shells.

Analogously to traditional FEM, the f-FEM is formulated starting from a
discretized form of the total potential energy functional
$\Pi[\textbf{u}(\bm{x})]$ given in Eq. (40). For this purpose, the 1D domain
$\Omega=[0,L]$ of the beam indicated in Fig. (1) is uniformly discretized into
disjoint two noded elements $\Omega^{e}_{i}=(x_{1}^{i},x_{1}^{i+1})$ of length
$l_{e}$ such that $\cup_{i=1}^{N_{e}}\Omega^{e}_{i}=\Omega$, $N_{e}$ being the
total number of discretized elements. It is immediate that
$\Omega^{e}_{j}\cap\Omega^{e}_{k}=\emptyset~{}\forall~{}j\neq k$. The unknown
field variables $u_{0}(x_{1})$ and $w_{0}(x_{1})$ in Eq. (41) can now be
evaluated at any point $x_{1}\in\Omega^{e}_{i}$ by interpolating the
corresponding nodal values for $\Omega^{e}_{i}$ as:

$\\{u_{0}(x_{1})\\}=[\mathcal{L}(x_{1})]\\{U_{e}(x_{1})\\};~{}~{}~{}\\{w_{0}(x_{1})\\}=[\mathcal{H}(x_{1})]\\{W_{e}(x_{1})\\}$
(46)

where $\\{U_{e}(x_{1})\\}$ and $\\{W_{e}(x_{1})\\}$ are the axial and
transverse displacement degrees of freedom of the two-noded element
$\Omega^{e}_{i}$. $[\mathcal{L}(x_{1})]$ and $[\mathcal{H}(x_{1})]$ are the
Lagrangian and Hermitian interpolation functions, respectively, chosen to
enforce the continuity of the axial and transverse displacement fields for the
Euler-Bernoulli beam theory. In terms of the above discussed numerical
discretization of the nonlocal domain, the fractional derivatives can be
expressed as:

$D^{\alpha}_{x_{1}}\left[u_{0}(x_{1})\right]=[\tilde{B}_{u}(x_{1})]\\{U_{g}\\};~{}D^{\alpha}_{x_{1}}\left[w_{0}(x_{1})\right]=[\tilde{B}_{w}(x_{1})]\\{W_{g}\\};~{}D^{\alpha}_{x_{1}}\left[D^{1}_{x_{1}}w_{0}(x_{1})\right]=[\tilde{B}_{\theta}(x_{1})]\\{W_{g}\\}$
(47a) where the matrices $[\tilde{B}_{\square}(x_{1})]$ corresponding to
nonlocal strain-displacement matrices are given as:
$[\tilde{B}_{\square}(x_{1})]=\int_{x_{1}-l_{A}}^{x_{1}+l_{B}}\mathcal{A}(x_{1},s_{1},l_{A},l_{B},\alpha)[B_{\square}(s_{1})][\tilde{\mathcal{C}_{\square}}(x_{1},s_{1})]\mathrm{d}s_{1}$
(47b) The kernel $\mathcal{A}(x_{1},s_{1},l_{A},l_{B},\alpha)$ in the above
equation is:
$\mathcal{A}(x_{1},s_{1},l_{A},l_{B},\alpha)=\begin{cases}\frac{1}{2}(1-\alpha)l_{A}^{\alpha-1}{(x_{1}-s_{1})^{-\alpha}}&~{}~{}s_{1}\in{(x_{1}-l_{A},x_{1})}\\\
\frac{1}{2}(1-\alpha)l_{B}^{\alpha-1}{(s_{1}-x_{1})^{-\alpha}}&~{}~{}s_{1}\in{(x_{1},x_{1}+l_{B})}\end{cases}$
(47c) and the matrices $[B_{\square}(x_{1})]$ $\square\in\\{u,w,\theta\\}$ are
the integer-order strain-displacement matrices given by:
$[B_{u}(s_{1})]=\frac{\mathrm{d}[\mathcal{L}(s_{1})]}{\mathrm{d}s_{1}};~{}[B_{w}(s_{1})]=\frac{\mathrm{d}[\mathcal{H}(s_{1})]}{\mathrm{d}s_{1}};~{}[B_{\theta}(s_{1})]=\frac{\mathrm{d}^{2}[\mathcal{H}(s_{1})]}{\mathrm{d}s_{1}^{2}}$
(47d)

Complete details of the steps involved in the numerical evaluation of system
matrices following an additional convolution integral in Eq. (47b) for
nonlocal matrices and resolving the singular kernel in Eq. (47c) for
fractional-order derivatives are provided in [34, 35].

We use the FE approximations of the different fractional-order derivatives to
obtain the algebraic governing equations corresponding to the geometrically
nonlinear thermoelastic response of the fractional-order beam. In the interest
of a more compact notation, the functional dependence of the different
physical quantities on the spatial variables will be implied, unless stated to
be constant. The first variation of the potential energy function
$\Pi[\textbf{u}]$ defined in Eq. (40) is obtained as:

$\delta\Pi=b\int_{0}^{L}\int_{-h/2}^{h/2}\delta\tilde{\epsilon}_{11}~{}\tilde{\sigma}_{11}\mathrm{d}x_{3}\mathrm{d}x_{1}-\int_{0}^{L}F_{t}\delta
w_{0}\mathrm{d}x_{1}-\int_{0}^{L}F_{a}\delta u_{0}\mathrm{d}x_{1}$ (48)

By using the strain-displacement relations in Eq. (36) and the stress-
resultants in Eq. (43) we obtain:

$\delta\Pi=\int_{0}^{L}\left\\{\mathcal{N}\left[D_{x_{1}}^{\alpha}\left(\delta
u_{0}\right)\right]+\mathcal{N}\left[D_{x_{1}}^{\alpha}w_{0}\right]\left[D_{x_{1}}^{\alpha}\left(\delta
w_{0}\right)\right]-\mathcal{M}\left[D_{x_{1}}^{\alpha}\left[D^{1}_{x_{1}}\left(\delta
w_{0}\right)\right]\right]-F_{a}\delta u_{0}-F_{t}\delta
w_{0}\right\\}\mathrm{d}x_{1}$ (49)

The minimum potential energy principle, $\delta\Pi=0$, is enforced to obtain
the algebraic equations of equilibrium. More specifically, by using the
numerical approximations developed for the different fractional derivatives
(see Eq. (47)) and then enforcing the minimization of the total potential
energy, we obtain the following system of nonlinear algebraic equations in
globally assembled vectors of nodal displacement degrees of freedom
$\\{U_{g}\\}$ and $\\{W_{g}\\}$:

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

where the different stiffness matrices are given by:

$[\tilde{K}_{11}]=\int_{0}^{L}A_{11}[\tilde{B}_{u}(x_{1})]^{T}[\tilde{B}_{u}(x_{1})]\mathrm{d}x_{1}$
(51a)
$[\tilde{K}_{12}]=\frac{1}{2}\int_{0}^{L}\underbrace{A_{11}\left(D_{x_{1}}^{\alpha}[w_{0}(x_{1})]\right)[\tilde{B}_{u}(x_{1})]^{T}[\tilde{B}_{w}(x_{1})]}_{\text{Nonlinear
matrix}}\mathrm{d}x_{1}$ (51b)
$[\tilde{K}_{21}]=\int_{0}^{L}\underbrace{A_{11}\left(D_{x_{1}}^{\alpha}[w_{0}(x_{1})]\right)[\tilde{B}_{w}(x_{1})]^{T}[\tilde{B}_{u}(x_{1})]}_{\text{Nonlinear
matrix}}\mathrm{d}x_{1}$ (51c)
$[\tilde{K}_{22}]=\int_{0}^{L}D_{11}[\tilde{B}_{\theta}(x_{1})]^{T}[\tilde{B}_{\theta}(x_{1})]~{}\mathrm{d}x_{1}+\frac{1}{2}\int_{0}^{L}\underbrace{\left[A_{11}\left(D_{x_{1}}^{\alpha}[w_{0}(x_{1})]\right)^{2}\right][\tilde{B}_{w}(x_{1})]^{T}[\tilde{B}_{w}(x_{1})]}_{\text{Nonlinear
matrix}}\mathrm{d}x_{1}$ (51d)

The axial and transverse force vectors due to the mechanical and thermal loads
are given as:

$\\{F_{A}\\}^{T}=\int_{0}^{L}F_{a}(x_{1})[S_{u}(x_{1})]\mathrm{d}x_{1}$ (52a)
$\\{F_{T}\\}^{T}=\int_{0}^{L}F_{t}(x_{1})[S_{w}(x_{1})]\mathrm{d}x_{1}$ (52b)
$\\{F_{A_{\theta}}\\}^{T}=\int_{0}^{L}N_{\theta}(x_{1})[B_{u}(x_{1})]\mathrm{d}x_{1}$
(52c)
$\\{F_{T_{\theta}}\\}^{T}=\int_{0}^{L}N_{\theta}(x_{1})\left(D_{x_{1}}^{\alpha}[w_{0}(x_{1})]\right)~{}[B_{w}(x_{1})]\mathrm{d}x_{1}-\int_{0}^{L}M_{\theta}(x_{1})[B_{\theta}(x_{1})]\mathrm{d}x_{1}$
(52d)

Note that the geometric nonlinearity introduces additional nonlinear
thermomechanical coupled terms. This nonlinear behavior is dependent on the
thermal properties of the beam as evident from the expressions of the thermal
stress resultants given in Eq. (45). In fact, these nonlinear effects are
expected to be significant at high temperature. These additional nonlinear
thermomechanical terms can be accounted for in two ways: approach #1: the
terms are treated as an external nonlinear thermal force; and approach #2: the
contribution of these terms is accounted via the stiffness matrix of the
system. The equivalence of the results obtained through both these approaches
and a comparison of their accuracy and stability is presented in [59]. In this
study, we follow the approach #1 so that the linear analysis of the nonlinear
system, for small displacements, becomes straightforward without requiring
changes to the stiffness matrix.

The algebraic equations (50) are solved for the nodal values of the
generalized displacement coordinates for an isotropic beam subject to
distributed thermal and mechanical loads. The solution to these equations
along with Eq. (34) gives the displacement field at any point within the beam.
The geometric nonlinearity in the system is highlighted by the deformation
dependent stiffness terms in Eqs. (51b)-(51d). Further, as previously
discussed, the additional nonlinear thermomechanical terms are introduced into
the model as a nonlinear transverse force as evident from Eq. (52d). Given the
nonlinear nature of the FE algebraic equations, a Newton-Raphson (NR)
iterative numerical scheme was adopted to obtain the solution of the Eq. (50).
Similar to classical nonlinear models, the NR procedure for the fractional-
order nonlinear equations also requires the evaluation of the tangent
stiffness matrix. The procedure to evaluate the tangent stiffness matrix as
well as the NR scheme can be found in [35].

The linear f-FEM for the thermoelastic response of the nonlocal isotropic beam
can be obtained from the above model by ignoring the contribution of the
nonlinear coupling term, that is $(D_{x_{1}}^{\alpha}[w_{0}(x_{1})])^{2}$ in
the system matrices as well as in the force vectors. Note that the axial and
transverse displacement fields for the linear elastic response due to
thermomechanical loads are decoupled. Finally, the f-FEM reduces to a local
thermoelastic study of beams when the fractional-order is set to $\alpha=1$.

##  5 Numerical results and discussion

We use the numerical model developed in §4 to analyze both the linear and the
geometrically nonlinear thermoelastic response of fractional-order nonlocal
isotropic beams. In order to satisfy the underlying assumptions of the Euler-
Bernoulli beam theory, the beam is assumed to be slender with an aspect ratio
of $L/h=100$. In the following studies, the length of the beam is maintained
at $L=1$m and the width of the beam is considered to be unity. The beam is
assumed made out of aluminum that is $E=70$GPa and $\alpha_{0}=23\times
10^{-6}~{}{\text{K}}^{-1}$. The constitutive parameters of the fractional-
order continuum model, order $\alpha$ and length scales $l_{A}$ and $l_{B}$
are provided wherever necessary. The length scales $l_{A}$ and $l_{B}$ at a
point within the domain of the isotropic beam are considered equal, that is
$l_{A}=l_{B}=l_{f}$. However, following the discussion in [34], these length
scales are truncated for points close to geometric boundaries of the beam. The
variable nature of the nonlocal length scales is demonstrated in the schematic
given in Fig. 2.

Figure 2: Schematic illustration of the variable nature of the nonlocal length
scales for different points along the length of the beam. The length scale on
the left side of the point P, $l_{A}^{P}$ is truncated such that
$l^{P}_{A}<l^{Q}_{A}$. Similarly, $l^{R}_{B}<l^{Q}_{B}$.

We analyzed numerically the effect of the fractional parameters $\alpha$ and
$l_{f}$ on the response of the beam subject to different loading and boundary
conditions. Both the linear and nonlinear cases are considered.

Before presenting the results, we make a few remarks concerning the validation
and convergence of the f-FEM. In this regard, the f-FEM procedure has already
been validated for linear BVPs in [6] and for nonlinear BVPs in [35]. Further,
as discussed in [6, 35], the convergence of the f-FEM with finer element
discretization is controlled by the “dynamic rate of convergence” defined as:
$N_{inf}(=l_{f}/l_{e})$, where $l_{e}$ is the length of the discretized
element. This parameter was shown to be dependent on the fractional model
parameters $\alpha$ and $l_{f}$ that determine the strength of the nonlocal
interactions between distant elements [6, 35]. Following these convergence
studies, the mesh discretization $N_{inf}=10$ was chosen. This choice allows
for sufficient number of elements to be included in the horizon of nonlocality
at any point, in order to accurately capture the nonlocal interactions [6,
35].

Linear thermoelastic response: we considered a simply supported beam subject
to a uniformly distributed transverse load (UDTL) of magnitude $q_{0}$ (in
N/m) and to the following thermal load :

$\theta(x_{1},x_{3})=\theta_{1}\left(1+\frac{2x_{3}}{h}\right)$ (53)

Note that the above distribution is obtained from solving the Fourier’s
conduction law corresponding to uniform temperatures being applied at the top
and bottom suyrfaces ($x_{3}=\pm h/2$) of the isotropic beam. Here, the bottom
surface of the beam is maintained at the ambient temperature $T_{0}$. It
follows, from Eq. (53), that the temperature of the top surface is
$T_{1}(=T_{0}+2\theta_{1})$. The nonlocal elastic response to these
thermomechanical loads for different values of fractional constitutive
parameters $\alpha$ and $l_{f}$ are compared in Fig. (3). The transverse
displacement along the length of the simply supported beam was explored for
different values of $\alpha$ while maintaining $l_{f}$ constant (Fig. (3(a))).
Similarly, the transverse displacement was evaluated for different values of
$l_{f}$ while maintaining $\alpha$ constant (Fig. (3(b))). The increase in
transverse displacement with the increasing degree of nonlocality, achieved
either by reducing $\alpha$ (see Fig. (3(a))) or by increasing $l_{f}$ (see
Fig. (3(b))), points towards the reduction of the stiffness of the fractional-
order beam. We emphasize that the consistent softening of the structure with
increasing degree of nonlocality was also observed for beams subject to
different loading and boundary conditions. In order to facilitate the analysis
of the thermoelastic response of the nonlocal beam, we have provided, as a
reference, the transverse displacement of the local beam ($\alpha=1$) for two
different cases: (a) absence of thermal load, i.e. $\theta_{1}=0$; and (b)
linear thermal load given in Eq. (53). From Fig. (3), note that the nonlocal
results converge to the local results for $\alpha$ approaching $1$ and
$l_{f}/L<<1$.

(a) $w_{0}$ vs $\alpha$.

(b) $w_{0}$ vs $l_{f}$.

Figure 3: Transverse displacement corresponding to the linear response of a
simply supported beam for $q_{0}=10^{4}$N/m and $\theta_{1}=10$K. The plot is
parameterized for different values of the fractional model parameters.

(a) $\overline{\sigma}_{11}$ vs $\alpha$.

(b) $\overline{\sigma}_{11}$ vs $l_{f}$.

Figure 4: Normalized values of normal axial stress $\overline{\sigma}_{11}$
across the thickness corresponding to the linear response of a simply
supported beam for $q_{0}=10^{4}$N/m and $\theta_{1}=10$K. The plot is
parameterized for different values of the fractional model parameters.

(a) $\epsilon_{11}$ vs $\alpha$.

(b) $\epsilon_{11}$ vs $l_{f}$.

Figure 5: Normal axial strain $\overline{\sigma}_{11}$ corresponding to the
linear response of a simply supported beam for $q_{0}=10^{4}$N/m and
$\theta_{1}=10$K. The plot is parameterized for different values of the
fractional model parameters.

Extending these studies, the normal axial stress $\sigma_{11}$ and the
fractional-order axial normal strain $\epsilon_{11}$ are studied in Figs. (4)
and (5). Here, the normal axial stress $\sigma_{11}$ is normalized as
follows[60]:

$\overline{\sigma}_{11}(L/2,x_{3})=\frac{1}{q_{0}}\left(\frac{h}{L}\right)^{2}\tilde{\sigma}_{11}(L/2,x_{3})$
(54)

In the Fig. 4, (marginally) higher values of $\bar{\sigma}_{11}$ are noted
corresponding to an increasing degree of nonlocality. More specifically, an
increase in the values of the normalized stress are noted for lower values of
the fractional-order $\alpha$ in Fig. (4(a)), and higher value of the nonlocal
length scale $l_{f}$ in Fig. (4(b)). These observations point to a consistent
reduction in the stiffness of the fractional-order nonlocal structures, which
agrees with the previously noted higher values for the transverse displacement
$\bar{w}$ when compared against the response of a local beam (see Fig. (3)).
Finally, note the constant value of shift in stress profiles at mid-surface
caused by the application of the thermal load. Additional results
corresponding to the nonlocal axial strain $\tilde{\epsilon}_{11}$ along the
length of the beam are presented in Fig. 5.

In the following, we compare the predictions of the fractional-order approach
to nonlocal thermoelasticity with classical integer-order nonlocal theories
available in the literature. More specifically, we compare the predictions of
the fractional-order model with the stress-driven integral model presented in
[15]. For this purpose, we consider the axial displacement of a doubly-clamped
beam subject to the following thermal load in the absence of any mechanical
load:

$\theta(x_{1},x_{3})=\theta_{1}\left(1-\frac{x_{1}}{L}\right)\frac{x_{1}}{L}$
(55)

The axial displacement of the doubly-clamped beam, obtained via the
fractional-order model, is presented in Fig. (6(a)) for different values of
the fractional-order $\alpha$. As evident from Fig. (6(a)), a consistent
softening behavior is obtained when considering increasing degree of
nonlocality in the fractional model and a constant thermal load. This
observation also complements the results obtained in previous studies on
fractional-order nonlocal elasticity involving only mechanical loads [34]. On
the contrary, the stress-driven integral approach does not predict either a
consistent stiffening or softening response of the same beam for different
values of the nonlocal constitutive parameters. In fact, by increasing the
degree of nonlocality the beam initially softens and then stiffens, hence
leading the authors to label the system as being ’unpredictable’. This result
in [15] obtained via the stress-driven integral models for nonlocal
thermoelasticity also contrasts with the consistent stiffening predicted by
the same model for purely mechanical loads in [24]. Recall that the
fractional-order continuum theories are successful in capturing the softening
effects of nonlocal interactions similar to the Eringen’s integral theories.
The above discussion also highlights an important advantage of the fractional-
order approaches to nonlocal elasticity, that is the consistency of the
predicted response even in a multi-physics scenario.

Further, as evident from Figs. (3,6), the fractional-order model predicts a
reduction in the stiffness of the structure irrespective of the nature of the
boundary conditions. This is evident from the consistent softening response
shown by the cantilever beams subject to thermal loads given in Eq. (55) as
depicted in Fig. (6(b)). This observation of consistent stiffening agrees with
the previously studied cases of the doubly clamped and simply supported beams.
This is contrary to differential models for nonlocal thermoelasticity, which
predict a stiffening response for the cantilever beam and a softening response
for other boundary conditions [61]. As discussed in the introduction, the
absence of such paradoxical results in the fractional-order approach to
nonlocal elasticity adopted here follows from the positive-definite
deformation energy density and self-adjoint nature of the fractional-order
nonlocal governing equations [34, 36]. Finally, we make an important remark
concerning the physically acceptable range for the order $\alpha$. As
demonstrated above, the degree of nonlocality increases with decreasing
$\alpha$ leading to a consistent softening of the structure. However, as shown
in [31, 32], results for very low values of $\alpha$ ($\approx 0.2$, which
indicates a very strong nonlocality in the fractional sense) lead to non-
physical solutions. This is also illustrated in the response of fully clamped
and cantilever beams subject tot thermal loads in Fig. 6. Hence, there exists
a limit on the order of the RC fractional derivative [32]. In other terms, it
can be concluded that the fractional-calculus based continuum models break
down for values of alpha close to the lower integer limit. We emphasize that
this breakdown is not a characteristic of the f-FEM technique as the same
observation is also noted when using finite difference methods to obtain the
numerical solutions (see, for example, [31, 32]). A detailed discussion
supported by a possible physical explanation for this loss of consistency in
elastic response with reducing fractional order is available in [34].

(a) Clamped-Clamped Beam.

(b) Cantilever Beam.

Figure 6: Axial displacement corresponding to the linear response of beams
subject to (a) clamped-clamped (b) clamped-free (cantilever), boundary
conditions for $q_{0}=0$N/m and a thermal load $\theta_{1}=100$K. The plot is
parameterized for different values of the fractional order with a fixed length
scale $l_{f}=L/5$.

Nonlinear thermoelastic response: in this case, the beam was subjected to a
UDTL of magnitude $q_{0}$ (in N/m) and a uniform thermal field
$\theta(x_{1},x_{3})=\theta_{0}$ (in K) above reference temperature. First, we
considered a beam clamped at both ends and subject to the thermomechanical
loads described above. The transverse displacement of the beam for a fixed
UDTL and varying magnitude of the thermal load was obtained and compared for
different values of $\alpha$ and $l_{f}$. The magnitude of the UDTL was fixed
at $q_{0}=5\times 10^{4}$N/m and the value of the uniform thermal field
$\theta_{0}$ was varied in order to analyze the effect of the thermal load on
the response of the beam. Additionally, in order to analyze the effect of the
fractional model parameters on the response of the beam, the transverse
displacement of the beam was compared for different values of $\alpha$ and
$l_{f}$. The results of this study are presented in Fig. (7) in terms of the
thermal load versus displacement. The displacement values presented in Fig.
(7) correspond to the maximum displacement of the mid-plane of the beams,
obtained at $x_{1}=L/2$. The effect of the fractional-order $\alpha$ with
$l_{f}$ being held constant is compared in Fig. (7(a)). The effect of $l_{f}$
for fixed $\alpha$ is presented in Fig. (7(b)). The result for the local case
($\alpha=1$) is provided in both cases as a reference. As evident from Fig.
(7), the nonlocal beam exhibits consistent softening with increasing thermal
loads and increasing degree of nonlocality. As observed earlier, the
thermoelastic response of the nonlocal beam converges to the corresponding
local elastic response for $\alpha$ approaching $1$ and $l_{f}/L\ll 1$.

(a) max($w_{0}$) vs $\alpha$.

(b) max($w_{0}$) vs $l_{f}$.

Figure 7: Transverse displacement at the mid-point of a clamped-clamped beam
subject to $q_{0}=5\times 10^{4}$N/m and thermal load $\theta_{0}$. The curves
are parameterized for different values of the fractional-order nonlocal
parameters.

Finally, the study was repeated for a beam pinned at both ends. The effect of
the fractional model parameters over the geometrically nonlinear response of
the pinned-pinned beam subject to thermomechanical loads is presented in Fig.
(8). Observations analogous to those drawn for the clamped-clamped beam can be
noted for this case. Remarkably, the fractional-order approach to the modeling
of nonlocal elasticity exhibits good consistency across a variety of boundary
and loading conditions for both the linear and geometrically nonlinear
responses. This behavior differs sharply from the paradoxical results reported
in the literature for either gradient or integral based approaches to nonlocal
elasticity [62, 24].

(a) max($w_{0}$) vs $\alpha$.

(b) max($w_{0}$) vs $l_{f}$.

Figure 8: Transverse displacement at the mid-point of a pinned-pinned beam
subject to $q_{0}=5\times 10^{4}$N/m and thermal load $\theta_{0}$. The curves
are parameterized for different values of the fractional-order nonlocal
parameters.

##  6 Conclusions

This study established the thermodynamic framework for fractional-order models
of nonlocal thermoelasticity. One of the most significant outcome of this
formulation is the ability to rigorously enforce, in a point-wise manner, the
thermodynamic balance laws in a nonlocal medium. This result stands in stark
contrast with respect to traditional integer-order methods that can satisfy
the thermodynamic laws only in a (weak) integral sense. An important
consequence of this property of the fractional-order framework is the
substantial simplification of the formulation of the free energy density and
of the resulting constitutive models of nonlocal thermoelasticity. The
thermodynamically-consistent fractional-order continuum theory is well suited
to develop accurate models to capture nonlocal interactions, heterogeneity,
and scale effects in complex elastic solids operating in a thermomechanical
environment. Additionally, the mechanical balance laws are also derived from
balance principles. This approach to the development of the governing
equations involved a modified Cauchy’s lemma for surface tractions in order to
include additional forces due to nonlocal effects. The resulting mechanical
governing equations are consistent with self-adjoint linear operators
admitting unique solutions. This is also an important difference of the
present framework in comparison to existing integer- and fractional-order
approaches for nonlocal solids available in the literature.

The efficacy of the fractional-order modeling approach was illustrated by
applying the framework to the analysis of the static response of a nonlocal
Euler-Bernoulli beam subject to combined thermomechanical loads. Numerical
results, obtained using the fractional finite element method, highlighted the
extremely robust nature of the fractional models by illustrating the
consistency of the predicted nonlocal response across different
thermomechanical loads and boundary conditions. A comparison with the existing
integer-order nonlocal models was also provided in order to illustrate the
advantages offered by the fractional-order theory of nonlocal
thermoelasticity.

##  7 Acknowledgements

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

## References

  * [1] L. J. Sudak.  Column buckling of multiwalled carbon nanotubes using nonlocal continuum mechanics.  Journal of Applied Physics, 94(11):7281–7287, 2003.
  * [2] S. C. Pradhan and T. Murmu.  Small scale effect on the buckling of single-layered graphene sheets under biaxial compression via nonlocal continuum mechanics.  Computational Materials Science, 47(1):268–274, 2009.
  * [3] C. Y. Wang, T. Murmu, and S. Adhikari.  Mechanisms of nonlocal effect on the vibration of nanoplates.  Applied Physics Letters, 98(15):153101, 2011.
  * [4] J. Romanoff, J. N. Reddy, and J. Jelovica.  Using non-local timoshenko beam theories for prediction of micro-and macro-structural responses.  Composite Structures, 156:410–420, 2016.
  * [5] J. P. Hollkamp, M. Sen, and F. Semperlotti.  Analysis of dispersion and propagation properties in a periodic rod using a space-fractional wave equation.  Journal of Sound and Vibration, 441:204–220, 2019.
  * [6] S. Patnaik and F. Semperlotti.  A generalized fractional-order elastodynamic theory for non-local attenuating media.  Proceedings of the Royal Society A, 476(2238):20200200, 2020.
  * [7] M. A. Kouchakzadeh, M. Rasekh, and H. Haddadpour.  Panel flutter analysis of general laminated composite plates.  Composite Structures, 92(12):2906–2915, 2010.
  * [8] P. Marzocca, S. A. Fazelzadeh, and M. Hosseini.  A review of nonlinear aero-thermo-elasticity of functionally graded panels.  Journal of Thermal Stresses, 34(5-6):536–568, 2011.
  * [9] S. A. Emam.  A general nonlocal nonlinear model for buckling of nanobeams.  Applied Mathematical Modelling, 37(10-11):6929–6939, 2013.
  * [10] L. Librescu, P. Marzocca, and W. A. Silva.  Supersonic/hypersonic flutter and postflutter of geometrically imperfect circular cylindrical panels.  Journal of Spacecraft and Rockets, 39(5):802–812, 2002.
  * [11] F. Ebrahimi and E. Salari.  Nonlocal thermo-mechanical vibration analysis of functionally graded nanobeams in thermal environment.  Acta Astronautica, 113:29–50, 2015.
  * [12] E. Kröner.  Elasticity theory of materials with long range cohesive forces.  International Journal of Solids and Structures, 3(5):731–742, 1967\.
  * [13] A. C. Eringen and D. G. B. Edelen.  On nonlocal elasticity.  International Journal of Engineering Science, 10(3):233–248, 1972\.
  * [14] C. Polizzotto.  Nonlocal elasticity and related variational principles.  International Journal of Solids and Structures, 38(42-43):7359–7380, 2001.
  * [15] R. Barretta, M. Čanadija, R. Luciano, and F. M. de Sciarra.  Stress-driven modeling of nonlocal thermoelastic behavior of nanobeams.  International Journal of Engineering Science, 126:53–67, 2018.
  * [16] R. H. J. Peerlings, M. G. D. Geers, R. De Borst, and W. A. M. Brekelmans.  A critical comparison of nonlocal and gradient-enhanced softening continua.  International Journal of Solids and Structures, 38(44-45):7723–7746, 2001.
  * [17] E. C. Aifantis.  Update on a class of gradient theories.  Mechanics of materials, 35(3-6):259–280, 2003.
  * [18] S. Sidhardh and M. C. Ray.  Inclusion problem for a generalized strain gradient elastic continuum.  Acta Mechanica, 229(9):3813–3831, 2018.
  * [19] S. Sidhardh and M. C. Ray.  Exact solutions for elastic response in micro-and nano-beams considering strain gradient elasticity.  Mathematics and Mechanics of Solids, page 1081286518761182.
  * [20] J. Yang, L. L. Ke, and S. Kitipornchai.  Nonlinear free vibration of single-walled carbon nanotubes using nonlocal timoshenko beam theory.  Physica E: Low-dimensional Systems and Nanostructures, 42(5):1727–1735, 2010.
  * [21] A. R. Srinivasa and J. N. Reddy.  A model for a constrained, finitely deforming, elastic solid with rotation gradient dependent strain energy, and its specialization to von kármán plates and beams.  Journal of the Mechanics and Physics of Solids, 61(3):873–885, 2013\.
  * [22] L. E. Shen, H-S. Shen, and C-L. Zhang.  Nonlocal plate model for nonlinear vibration of single layer graphene sheets in thermal environments.  Computational Materials Science, 48(3):680–685, 2010.
  * [23] A. Tounsi, S. Benguediab, A. Semmah, M. Zidour, et al.  Nonlocal effects on thermal buckling properties of double-walled carbon nanotubes.  Advances in Nano Research, 1(1):1, 2013.
  * [24] G. Romano, R. Barretta, M. Diaco, and F. M. de Sciarra.  Constitutive boundary conditions and paradoxes in nonlocal elastic nanobeams.  International Journal of Mechanical Sciences, 121:151–156, 2017\.
  * [25] J. K. Phadikar and S. C. Pradhan.  Variational formulation and finite element analysis for nonlocal elastic nanobeams and nanoplates.  Computational materials science, 49(3):492–499, 2010.
  * [26] A. Anjomshoa.  Application of ritz functions in buckling analysis of embedded orthotropic circular and elliptical micro/nano-plates based on nonlocal elasticity theory.  Meccanica, 48(6):1337–1353, 2013.
  * [27] I. Podlubny.  Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, volume 198.  Elsevier, 1998.
  * [28] G. Cottone, M. Di Paola, and M. Zingales.  Elastic waves propagation in 1d fractional non-local continuum.  Physica E: Low-dimensional Systems and Nanostructures, 42(2):95–103, 2009.
  * [29] M. Di Paola and M. Zingales.  Long-range cohesive interactions of non-local continuum faced by fractional calculus.  International Journal of Solids and Structures, 45(21):5642–5659, 2008.
  * [30] A. Carpinteri, P. Cornetti, and A. Sapora.  Nonlocal elasticity: an approach based on fractional calculus.  Meccanica, 49(11):2551–2569, 2014.
  * [31] W. Sumelka and T. Blaszczyk.  Fractional continua for linear elasticity.  Archives of Mechanics, 66(3):147–172, 2014.
  * [32] W. Sumelka, T. Blaszczyk, and C. Liebold.  Fractional euler–bernoulli beams: Theory, numerical study and experimental validation.  European Journal of Mechanics-A/Solids, 54:243–251, 2015.
  * [33] J. P. Hollkamp and F. Semperlotti.  Application of fractional order operators to the simulation of ducts with acoustic black hole terminations.  Journal of Sound and Vibration, 465:115035, 2020.
  * [34] S. Patnaik, S. Sidhardh, and F. Semperlotti.  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\.
  * [35] S. Sidhardh, S. Patnaik, and F. Semperlotti.  Geometrically nonlinear response of a fractional-order nonlocal model of elasticity.  International Journal of Nonlinear Mechanics, 125:103529, 2020.
  * [36] S. Patnaik, S. Sidhardh, and F. Semperlotti.  Fractional-order models for the static and dynamic analysis of nonlocal plates.  Communications in Nonlinear Science and Numerical Simulation, page 105601, 2020.
  * [37] S. Patnaik, S. Sidhardh, and F. Semperlotti.  Geometrically nonlinear analysis of nonlocal plates using fractional calculus.  International Journal of Mechanical Sciences, 179:105710, 2020.
  * [38] A. C. Eringen.  Linear theory of nonlocal elasticity and dispersion of plane waves.  International Journal of Engineering Science, 10(5):425–435, 1972\.
  * [39] A. C. Eringen.  Theory of nonlocal thermoelasticity.  International Journal of Engineering Science, 12(12):1063–1077, 1974.
  * [40] A. C. Eringen.  On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves.  Journal of applied physics, 54(9):4703–4710, 1983.
  * [41] S. Sidhardh, S. Patnaik, and F. Semperlotti.  Fractional-order structural stability: Formulation and application to the critical load of slender structures.  arXiv preprint arXiv:2008.11528, 2020.
  * [42] S. Sidhardh, S. Patnaik, and F. Semperlotti.  Analysis of the post-buckling response of nonlocal plates via fractional order continuum theory.  Journal of Applied Mechanics, https://doi.org/10.1115/1.4049224, 2020.
  * [43] Y. Z. Povstenko.  Fractional heat conduction equation and associated thermal stress.  Journal of Thermal Stresses, 28(1):83–102, 2004.
  * [44] Y. Z. Povstenko.  Thermoelasticity that uses fractional heat conduction equation.  Journal of Mathematical Sciences, 162(2):296–305, 2009.
  * [45] Y. Povstenko.  Fractional thermoelasticity, volume 219.  Springer, 2015.
  * [46] M. Žecová and J. Terpák.  Heat conduction modeling by using fractional-order derivatives.  Applied Mathematics and Computation, 257:365–373, 2015.
  * [47] L. Vázquez, J. Trujillo, and M. P. Velasco.  Fractional heat equation and the second law of thermodynamics.  Fractional Calculus and Applied Analysis, 14(3):334–342, 2011.
  * [48] C. S. Drapaca and S. Sivaloganathan.  A fractional model of continuum mechanics.  Journal of Elasticity, 107(2):105–123, 2012.
  * [49] F. Balta and E. S. Şuhubi.  Theory of nonlocal generalised thermoelasticity.  International Journal of Engineering Science, 15(9-10):579–588, 1977.
  * [50] V. E. Tarasov.  Fractional vector calculus and fractional maxwell’s equations.  Annals of Physics, 323(11):2756–2778, 2008.
  * [51] S. Patnaik, S. Sidhardh, and F. Semperlotti.  Towards a unified approach to nonlocal elasticity via fractional-order mechanics.  International Journal of Mechanical Sciences, 189:105992, 2020.
  * [52] P. Germain, P. Suquet, and Q. S. Nguyen.  Continuum thermodynamics.  ATJAM, 50:1010–1020, 1983.
  * [53] JN Reddy.  Nonlocal nonlinear formulations for bending of classical and shear deformation theories of beams and plates.  International Journal of Engineering Science, 48(11):1507–1518, 2010.
  * [54] N. Challamel, Z. Zhang, C. M. Wang, J. N. Reddy, Q. Wang, T. Michelitsch, and B. Collet.  On nonconservativeness of eringen’s nonlocal elasticity in beam mechanics: correction from a discrete-based approach.  Archive of Applied Mechanics, 84(9-11):1275–1292, 2014.
  * [55] J. T. Oden.  Finite element analysis of nonlinear problems in the dynamical theory of coupled thermoelasticity.  Nuclear Engineering and Design, 10(4):465–475, 1969.
  * [56] P. A. Kelly.  Lecture notes on material models in continuum solid mechanics.  2020\.
  * [57] K. Washizu.  Variational methods in elasticity and plasticity, volume 3.  Pergamon press Oxford, 1975.
  * [58] P. G. Ciarlet.  A justification of the von kármán equations.  Archive for Rational Mechanics and Analysis, 73(4):349–389, 1980\.
  * [59] G. N. Praveen and J. N. Reddy.  Nonlinear transient thermoelastic analysis of functionally graded ceramic-metal plates.  International Journal of Solids and Structures, 35(33):4457–4476, 1998.
  * [60] P. Khodabakhshi and J. N. Reddy.  A unified integro-differential nonlocal model.  International Journal of Engineering Science, 95:60–75, 2015.
  * [61] A. M. Zenkour.  Nonlocal thermoelasticity theory without energy dissipation for nano-machined beam resonators subjected to various boundary conditions.  Microsystem Technologies, 23(1):55–65, 2017.
  * [62] N. Challamel and C. M. Wang.  The small length scale effect for a non-local cantilever beam: a paradox solved.  Nanotechnology, 19(34):345703, 2008.