A Note on Gradient/Fractional One-Dimensional Elasticity and Viscoelasticity

Abstract

An introductory discussion on a (weakly non-local) gradient generalization of some one-dimensional elastic and viscoelastic models, and their fractional extension is provided. Emphasis is placed on the possible implications of micro- and nano- engineering problems, including small-scale structural mechanics and composite materials, as well as collagen biomechanics and nanomaterials.

1. Introduction

A gradient generalization of Hooke’s law of elasticity has been proposed in the past by the last author to model behavior at the nanoscale, where heterogeneity and size effects become important [1,2]. A recent review with applications to various branches of elasticity theory and an extension to include fractional derivatives can be found in [3,4].

Among the main results of the above works were the removal of singularities from dislocation lines and crack tips, the interpretation of size effects, and the description of experimentally observed nanoscale phenomena that classical elasticity could not capture. The basic mechanism in the mathematical formulation was to introduce the Laplacian of strain and/or stress into the standard elasticity constitutive equation, as a representative measure of the underlying heterogeneity at the micron-scale and nanoscale. This Laplacian-enhanced gradient generalization of classical elasticity theory was extended to include fractional and fractal considerations in cases where the underlying micro/nano-structure could not be described with the use of integer derivatives and Euclidean geometry [5,6]. It is noted, in relation to this, that the standard Laplacian operator of integer order has a direct straightforward extension to a fractional Laplacian operator of non-integer order with corresponding nice properties, and, thus, easily derived fractional results can be obtained.

The purpose of the present article is to make use of the aforementioned developments by focusing on simple one-dimensional elastic and viscoelastic models in order to better illustrate the potential of such generalizations and their implications for various problems in current engineering and biomedical technology, where micro/nanoscale effects need to be accounted for in a manner that is as simple as possible, yet effective. In particular, reference is made to micro/nano metallic and polymeric ultra-thin rods used in the fabrication of nanocomposite materials and miniature devices/frames for energy storage and medical implants, as well as to one-dimensional tendon and collagen fiber problems that biomechanics deal with.

In Section 2 we briefly discuss the one-dimensional extension of the classical Hooke’s law of elasticity to include spatial derivatives of the strain of integers and fractional orders. In Section 3 we discuss the corresponding generalizations of the Kelvin–Voigt model of viscoelasticity. Finally, in Section 4, a discussion about a number of related issues that need to be further considered and a number of future research directions that need to be further explored is provided.

2. Gradient and Fractional Hooke’s Elasticity Model

2.1. Gradient Elastic Bar

The classical Hooke’s law for an elastic bar in one dimension reads

$$\sigma =E\epsilon ,$$

(1)

where $\sigma$ is the stress, $\epsilon$ ($=\partial u/\partial x$, u denotes displacement and x the spatial dimension) is the strain, and E is the elastic modulus. The gradient generalization of Equation (1) amounts to replacing the (macroscopic) strain, $\epsilon$, with its local (microscopic) gradient counterpart, $\epsilon -{\ell }^2\left({\partial }^2\epsilon /\partial x^2\right)$, where ℓ is an internal length related to the underlying substructure (e.g., size/spacing of inclusions in a composite material, grain size in a polycrystal). If static equilibrium is assumed under the action of an applied stress ${\sigma }_0$, it follows that the governing differential equation for a semi-infinite bar ($0\le x<\infty$) is

$$\epsilon -{\ell }^2\frac{{\partial }^2\epsilon }{\partial x^2}=\frac{{\sigma }_0}{E}={\epsilon }_0,$$

(2)

with a general solution given by the expression

$$\epsilon \left(x\right)=\frac{{\sigma }_0}{E}+c{e}^{-x/\ell },$$

(3)

on the assumption that $\epsilon$ is finite as $x\to \infty$, or that the gradient solution converges to the classical solution at the free boundary of a very long bar. By also assuming that, at the finite boundary ($x=0$) the microstrain vanishes ($\epsilon =0$), we obtain $c=-\frac{{\sigma }_0}{E}$, and thus

$$\epsilon \left(x\right)=\frac{{\sigma }_0}{E}\left(1-e^{-x/\ell }\right),$$

(4)

confirming again the fact that, at large distances from the fixed end of the bar, the strain approaches the value predicted by classical elasticity. Other types of boundary conditions and; other one-dimensional configurations may also be considered: e.g., a bar of length L ($-L/2\le x\le L/2$) with symmetric loading through a stress ${\sigma }_0$ applied at its two ends.

On returning to Equation (4), it follows that the local (micro) strain at the end of a bar of length L is

$$\epsilon \left(L\right)=\frac{{\sigma }_0}{E}\left(1-e^{-L/\ell }\right),$$

(5)

which, as already mentioned earlier, for a very long bar ($L\gg \ell$) is the same as the strain predicted by the classical macroscopic Hooke’s law represented in Equation (1). In addition, for a bar for which the ratio of the internal length ℓ over its length L cannot be neglected, Equation (5) can be written in the form of Equation (1), i.e.,

$$\sigma =E_{\mathrm{eff}}\epsilon ;E_{\mathrm{eff}}=E{\left[1-e^{-L/\ell }\right]}^{-1},$$

(6)

where $E_{\mathrm{eff}}$ is an effective elastic modulus. Possible implications of Equation (6) or its counterparts for other types of boundary conditions and configurations can be explored in conjunction with experiments and some efforts along these lines for standard problems of structural mechanics (trusses, frames, beams, shear-moment diagrams) are currently under investigation in [7,8].

2.2. A Fractional Extension

The fractional counterpart of Equation (2) reads

$$\epsilon +{\ell }^{\alpha }{\left(-\mathrm{\Delta }\right)}^{\alpha /2}\epsilon =\frac{{\sigma }_0}{E},$$

(7)

where ${\ell }^{\alpha }$ is a fractional internal length (in units of ${\left[L\right]}^{\alpha }$), and ${\left(-\mathrm{\Delta }\right)}^{\alpha /2}$ denotes the Riesz fractional Laplacian [9] of order $\alpha$ ($1<\alpha <2$), defined in terms of the Fourier transform as

$$\mathcal{F}\left\{{\left(-\mathrm{\Delta }\right)}^{\alpha /2}f\left(x\right)\right\}\left(k\right)={\left|k\right|}^{\alpha }\mathcal{F}\left\{f\left(x\right)\right\}\left(k\right),$$

(8)

where $f\left(x\right)$ denotes a Fourier integrable function. Other definitions of a fractional Laplacian are possible, and a detailed discussion is provided in [10]. In any case, it is beneficial to express the fractional Laplacian operator as the convolution of the integer-order Laplacian with an appropriate power-law kernel [11] (see also [12]). In the one-dimensional case, the above process provides

$$\begin{array}{cc}\hfill {(-\mathrm{\Delta })}^{\alpha /2}f\left(x\right)& ={\mathcal{F}}^{-1}\left\{{\left|k\right|}^{\alpha }\mathcal{F}\left\{f\right\}\left(k\right)\right\}\left(x\right)\hfill \\ \hfill & =-{\mathcal{F}}^{-1}\left\{{\left|k\right|}^{\alpha -2}\left[-{\left|k\right|}^2\mathcal{F}\left\{f\right\}\left(k\right)\right]\right\}\left(x\right)\hfill \\ \hfill & =-{\mathcal{F}}^{-1}\left\{{\left|k\right|}^{\alpha -2}\mathcal{F}\left\{\mathrm{\Delta }f\right\}\left(k\right)\right\}\left(x\right).\hfill \end{array}$$

(9)

In Equation (9), the Laplacian $\mathrm{\Delta }$ coincides with the second order spatial derivative ${\partial }_x^2$, for which the property $\mathcal{F}\left\{\mathrm{\Delta }f\right\}\left(k\right)=-{\left|k\right|}^2\mathcal{F}\left\{f\right\}\left(k\right)$ holds. In the case of unbounded one-dimensional space, $-\infty <x<\infty$ with $1<\alpha <2$, Equation (9) can be expressed in terms of the Riesz fractional integral operator [9] as ${(-\mathrm{\Delta })}^{\alpha /2}f\left(x\right)=-\left({\mathcal{I}}^{2-\alpha }{\partial }_x^{2}f\right)\left(x\right)$, or, equivalently, in convolution form through the Riesz potential $-{\int }_{\mathbb{R}}{R}_{\alpha }(x-\xi ){\partial }_{\xi }^{2}f\left(\xi \right)d\xi$, where $R_{\alpha }\left(x\right):{\mathcal{F}}^{-1}\left\{{\left|k\right|}^{-(2-\alpha )}\right\}\left(x\right)={\tilde{\gamma }}_{\alpha }^{-1}{\left|x\right|}^{1-\alpha }$ with ${\tilde{\gamma }}_{\alpha }=2^{2-\alpha }\sqrt{\pi }\mathrm{\Gamma }(1-\frac{\alpha }{2})/\mathrm{\Gamma }\left(\frac{\alpha -1}{2}\right)$. After some simplifications we obtain

$${(-\mathrm{\Delta })}^{\alpha /2}f\left(x\right)=\frac{1}{2cos(\pi \alpha /2)\mathrm{\Gamma }(2-\alpha )}{\int }_{-\infty }^{\infty }\frac{1}{{|x-\xi |}^{\alpha -1}}\frac{{\partial }^{2}f}{\partial {\xi }^2}d\xi .$$

(10)

Equation (10) provides us a Caputo-type expression for the one-dimensional fractional Laplacian in infinite space. However, in the case where the domain under consideration has finite boundaries, or occupies a semi-infinite space, additional terms arising from the boundaries should be included. In particular, for the problem of a semi-infinite bar with a fixed boundary at $x=0$, the fractional Laplacian of the strain $\epsilon$ also includes two additional boundary terms proportional to the value of the strain $\epsilon \left(0\right)$ and its derivative ${\partial }_x\epsilon \left(0\right)$ at $x=0$, respectively. This aspect, which does not arise in the non-fractional treatment of the problem at hand, will be discussed elsewhere for more general configurations and gradient models. For the present case and configuration, however, it turns out that Equation (10) leads to

$$\begin{array}{c}{(-\mathrm{\Delta })}^{\alpha /2}\epsilon \left(x\right)=\frac{1}{2cos(\pi \alpha /2)}\left\{\frac{1}{\mathrm{\Gamma }(2-\alpha )}{\int }_{-\infty }^{\infty }\frac{\theta \left(\xi \right)}{{|x-\xi |}^{\alpha -1}}\frac{{\partial }^2\epsilon }{\partial {\xi }^2}d\xi +\right.\hfill \\ \hfill \left.\frac{1}{\mathrm{\Gamma }(2-\alpha )}\frac{{\partial }_x\epsilon \left(0\right)}{{\left|x\right|}^{\alpha -1}}+\frac{sgn\left(x\right)}{\mathrm{\Gamma }(1-\alpha )}\frac{\epsilon \left(0\right)}{{\left|x\right|}^{\alpha }}\right\},\end{array}$$

(11)

where $\theta \left(\xi \right)$ denotes the Heaviside function. It then follows [10] that we can conclude that ${(-\mathrm{\Delta })}^{\alpha /2}=-{\mathcal{D}}_{x_{+}}^{\alpha }$, where ${\mathcal{D}}_{x_{+}}^{\alpha }$ denotes the Riesz fractional derivative in $[0,\infty )$, defined as

$${\mathcal{D}}_{x_{+}}^{\alpha }\epsilon \left(x\right)=-\frac{1}{2cos(\pi \alpha /2)}\left({}_0^{RL}{D}_x^{\alpha }\epsilon \left(x\right)+{}_x^{RL}{D}_{\infty }^{\alpha }\epsilon \left(x\right)\right),$$

(12)

with ${}_0^{RL}{D}_x^{\alpha }$ and ${}_x^{RL}{D}_{\infty }^{\alpha }$ denoting the left- and right-sided Riemann–Liouville fractional derivatives given, respectively, by

$${}_0^{RL}{D}_x^{\alpha }\epsilon \left(x\right)=\frac{1}{\mathrm{\Gamma }(2-\alpha )}\frac{{\partial }^2}{\partial x^2}{\int }_0^x\frac{\epsilon \left(\xi \right)}{{\left(x-\xi \right)}^{\alpha -1}}d\xi ,$$

(13)

and

$${}_x^{RL}{D}_{\infty }^{\alpha }\epsilon \left(x\right)=\frac{1}{\mathrm{\Gamma }(2-\alpha )}\frac{{\partial }^2}{\partial x^2}{\int }_x^{\infty }\frac{\epsilon \left(\xi \right)}{{\left(\xi -x\right)}^{\alpha -1}}d\xi .$$

(14)

The equivalence between the fractional Laplacian ${(-\mathrm{\Delta })}^{\alpha /2}$ and the Riesz fractional derivative ${\mathcal{D}}_{x_{+}}^{\alpha }$ allows us to rewrite Equation (7) as $\epsilon -{\ell }^{\alpha }{\mathcal{D}}_{x_{+}}^{\alpha }\epsilon ={\sigma }_0/E$. Then, by applying the Fourier sine transform,

$$\begin{array}{c}\hfill {\mathcal{F}}_s\left\{\epsilon \right\}\left(k\right)=\tilde{\epsilon }\left(k\right)=\sqrt{\frac{2}{\pi }}{\int }_0^{\infty }\epsilon \left(x\right)sin\left(kx\right)dx,\end{array}$$

(15)

$$\begin{array}{c}\hfill {\mathcal{F}}_s^{-1}\left\{\tilde{\epsilon }\right\}\left(x\right)=\epsilon \left(x\right)=\sqrt{\frac{2}{\pi }}{\int }_0^{\infty }\tilde{\epsilon }\left(k\right)sin\left(kx\right)dk,\end{array}$$

(16)

and taking into account the boundary condition $\epsilon \left(0\right)=0$, we are led to the following integral expression

$$\epsilon \left(x\right)=\frac{2{\sigma }_0}{\pi E}{\int }_0^{\infty }\frac{sin\left(kx\right)}{k\left(1+{\ell }^{\alpha }{k}^{\alpha }\right)}dk.$$

(17)

Equation (17) is the required fractional generalization of Equation (4). To confirm this, we observe that, in the limit $\alpha \to 2$, we recover Equation (4) through the integral (3.725–1, p. 429 of [13])

$$\epsilon \left(x\right){|}_{\alpha \to 2}=\frac{2{\sigma }_0}{\pi E}{\int }_0^{\infty }\frac{sin\left(kx\right)}{k\left(1+{\ell }^2{k}^2\right)}dk=\frac{{\sigma }_0}{E}\left(1-e^{-x/\ell }\right).$$

(18)

By noting that $sin\left(kx\right)/k={\int }_0^{x}cos\left(k\xi \right)d\xi$, we can express Equation (17) as $\epsilon \left(x\right)=\frac{2{\sigma }_0}{E}{\int }_0^x\varphi \left(\xi \right)d\xi$, where

$$\varphi \left(x\right)=\frac{1}{\pi }{\int }_0^{\infty }\frac{cos\left(kx\right)}{1+{\ell }^{\alpha }{k}^{\alpha }}dk.$$

(19)

The function $\varphi \left(x\right)$ is identified as the fundamental solution of the fractional Helmholtz equation, i.e., $\varphi +{\ell }^{\alpha }{(-\mathrm{\Delta })}^{\alpha /2}\varphi =\delta \left(x\right)$. In the limit $\alpha \to 2$, we have ${\varphi |}_{\alpha \to 2}=(\pi /2\ell )exp(-x/\ell )$ (see 3.723–2, p. 428 of [13]). By using a number of standard results in fractional calculus [4,14,15], we can obtain the following Mellin–Barnes integral representation for the function $\varphi \left(x\right)$, expressed in terms of the Fox H-function $H_{p,q}^{m,n}$ [9,16] as

$$\varphi \left(x\right)=\frac{1}{2\alpha \ell \sqrt{\pi }}{H}_{1,3}^{2,1}\left[\frac{x}{2\ell }|\begin{array}{c}\left(1-\frac{1}{\alpha },\frac{1}{\alpha }\right)\hfill \\ \left(1-\frac{1}{\alpha },\frac{1}{\alpha }\right),\left(0,\frac{1}{2}\right),\left(\frac{1}{2},\frac{1}{2}\right)\hfill \end{array}\right].$$

(20)

From Equation (20) it can be verified that $\varphi \left(x\right)$ is consistent with ${\varphi |}_{\alpha \to 2}=(\pi /2\ell )$$exp(-x/\ell )$ from Equation (19), via the application of the Legendre duplication formula. It then follows that $\epsilon \left(x\right)$ can be expressed in terms of the H-function through Equation (20) as

$$\epsilon \left(x\right)=\frac{{\sigma }_0}{\alpha \ell \sqrt{\pi }E}{\int }_0^x{H}_{1,3}^{2,1}\left[\frac{\xi }{2\ell }|\begin{array}{c}\left(1-\frac{1}{\alpha },\frac{1}{\alpha }\right)\hfill \\ \left(1-\frac{1}{\alpha },\frac{1}{\alpha }\right),\left(0,\frac{1}{2}\right),\left(\frac{1}{2},\frac{1}{2}\right)\hfill \end{array}\right]d\xi .$$

(21)

Equation (21) clearly satisfies the boundary condition $\epsilon \left(0\right)=0$ and can be utilized to express the local strain at the end of a bar of length L through the integral

$$\epsilon \left(L\right)=\frac{{\sigma }_0}{\alpha \ell \sqrt{\pi }E}{\int }_0^L{H}_{1,3}^{2,1}\left[\frac{\xi }{2\ell }|\begin{array}{c}\left(1-\frac{1}{\alpha },\frac{1}{\alpha }\right)\hfill \\ \left(1-\frac{1}{\alpha },\frac{1}{\alpha }\right),\left(0,\frac{1}{2}\right),\left(\frac{1}{2},\frac{1}{2}\right)\hfill \end{array}\right]d\xi .$$

(22)

In analogy with Equation (6), we define the fractional effective elastic modulus $E_{eff}^{\alpha }$ as

$$E_{eff}^{\alpha }=E\alpha \ell \sqrt{\pi }{\left[{\int }_0^L{H}_{1,3}^{2,1}\left[\frac{\xi }{2\ell }|\begin{array}{c}\left(1-\frac{1}{\alpha },\frac{1}{\alpha }\right)\hfill \\ \left(1-\frac{1}{\alpha },\frac{1}{\alpha }\right),\left(0,\frac{1}{2}\right),\left(\frac{1}{2},\frac{1}{2}\right)\hfill \end{array}\right]d\xi \right]}^{-1}.$$

(23)

For a long bar $L\gg \ell$, the strain $\epsilon$ given by Equation (21) recovers the classical Hookean expression ${\epsilon \left(L\right)|}_{L\to \infty }={\sigma }_0/E$, since, in the limit $L\to \infty$, the following relation holds

$${\epsilon |}_{L\to \infty }=\frac{{\sigma }_0}{\alpha \ell \sqrt{\pi }E}{\int }_0^{\infty }{H}_{1,3}^{2,1}\left[\frac{\xi }{2\ell }|\begin{array}{c}\left(1-\frac{1}{\alpha },\frac{1}{\alpha }\right)\hfill \\ \left(1-\frac{1}{\alpha },\frac{1}{\alpha }\right),\left(0,\frac{1}{2}\right),\left(\frac{1}{2},\frac{1}{2}\right)\hfill \end{array}\right]d\xi =\frac{{\sigma }_0}{E}.$$

(24)

Furthermore, in the limit $\ell \to 0$, the classical Hooke’s law of Equation (1) is recovered. An alternative expression for the strain can be obtained through the direct application of the Mellin transform in Equation (17), yielding the following equivalent H-function representation

$$\epsilon \left(x\right)=\frac{{\sigma }_0}{\alpha \ell \sqrt{\pi }E}{H}_{1,3}^{2,1}\left[\frac{x}{2\ell }|\begin{array}{c}\left(1,\frac{1}{\alpha }\right)\hfill \\ \left(1,\frac{1}{\alpha }\right),\left(\frac{1}{2},\frac{1}{2}\right),\left(0,\frac{1}{2}\right)\hfill \end{array}\right].$$

(25)

It is noted that Equation (25) can also be obtained by the direct integration of Equation (21). In the limit $\alpha \to 2$, Equation (25) reduces to its integer-order counterpart of Equation (4).

3. Gradient and Fractional Kelvin–Voigt Viscoelasticity Model

3.1. Gradient Kelvin–Voigt Viscoelastic Bar

The Kelvin–Voigt solid is a representative example of linear viscoelastic behavior, exhibiting creep and damping. It reads

$$\sigma =E\epsilon +\mu \dot{\epsilon },$$

(26)

indicating that the stress is made up by an elastic component ($E\epsilon$) and a viscous one ($\mu \dot{\epsilon }$) where $\dot{\epsilon }$ is the strain rate and $\mu$ denotes the viscosity coefficient taken as constant.

For a constant applied stress $\sigma ={\sigma }_0$ (the same as in Section 2), the integration of Equation (26) gives

$$\epsilon =\frac{{\sigma }_0}{E}\left(1-e^{-t/\tau }\right),$$

(27)

where $\tau =\mu /E$, and the initial condition $\epsilon \left(0\right)=0$ was used. For long times ($t\to \infty$), Equation (27) shows that a one-dimensional bar, made by this solid, extends (or “creeps”) under a constant applied stress from its initial length where $\epsilon \left(0\right)=0$ (i.e., before ${\sigma }_0$ is applied) to a final length where $\epsilon (\infty )={\sigma }_0/E$, (i.e. the value attained instantaneously for an elastic bar). In addition to this “creeping” behavior, the model exhibits damping, i.e., the loading and unloading paths/curves do not coincide (as for elastic behavior) and the area enclosed by the two curves is the energy dissipated during the cycle. This type of anelastic behavior is attributed to microscopic processes, such as production/annihilation, and the internal motion of structural defects (vacancy clusters in metals, fragments of molecular chains in polymers, microcracks in rocks).

In fact, general viscoelastic behavior may be described by the following set of equations (e.g., [17,18])

$$\sigma =f(\epsilon ,\alpha );\dot{\alpha }=g(\epsilon ,\alpha ),$$

(28)

where $\left(f,g\right)$ are constitutive functions, and $\alpha$ denotes a set of internal variables. In this connection, it should be pointed out that the last author suggested [19,20] to replace the “evolution law” for the internal variables as given by Equation (28)${}_2$ with a “complete balance law” containing both a rate and a flux term.

For an internal variable $\alpha$ with a diffusive transport (and linear behavior as assumed for all models discussed herein), Equations (28)${}_1$ and (28)${}_2$ are then replaced by

$$\sigma =E\epsilon -\mathrm{\Lambda }\alpha ;\dot{\alpha }=D\frac{{\partial }^2\alpha }{\partial x^2}-M\alpha -N\epsilon ,$$

(29)

where all the coefficients are constants. For a constant applied stress $\sigma ={\sigma }_0$, the time differentiation of Equation (29)${}_1$ in Fourier space gives ${\dot{\alpha }}_k=(E/\mathrm{\Lambda }){\dot{\epsilon }}_k$, while, by Fourier transforming Equation (29)${}_2$, we obtain ${\dot{\alpha }}_k=-\left(D{k}^2+M\right){\alpha }_k-N{\epsilon }_k$, where the subscript k denotes a Fourier transformed quantity. We see that by neglecting the diffusion-like term we obtain $-\mathrm{\Lambda }{\alpha }_k=\left(N\mathrm{\Lambda }/M\right){\epsilon }_k+\left(E/M\right){\dot{\epsilon }}_k$, or $-\mathrm{\Lambda }\alpha =\left(N\mathrm{\Lambda }/M\right)\epsilon +\left(E/M\right)\dot{\epsilon }$, and, thus, the expression of the Kelvin–Voigt model given by Equation (26) and the creep formula given by Equation (27) result from a proper identification of the coefficients.

If the diffusion term is retained in Equation (29)${}_2$ and we adiabatically eliminate the term ${\dot{\alpha }}_k$ (${\dot{\alpha }}_k\approx 0$) in the Fourier space, we obtain ${\alpha }_k=-\left[N/\left(M+D{k}^2\right)\right]{\epsilon }_k$ which upon Taylor expansion, gives ${\alpha }_k=-\left(\frac{N}{M}-\frac{DN}{M^2}{k}^2\right){\epsilon }_k$, or $\alpha =-\frac{N}{M}\epsilon -\frac{DN}{M^2}\frac{{\partial }^2\epsilon }{\partial x^2}$. Thus, $\alpha$ can be eliminated from Equation (29)${}_1$ in favor of the second spatial derivative of $\epsilon$, resulting in the gradient elastic model of Section 2. It follows that, by assuming a pair of internal variables, one behaving as in the previous paragraph and another as in the present paragraph, we can arrive at a gradient extension of the Kelvin–Voigt model of the form

$$\sigma =E\left(\epsilon -{\ell }^2\frac{{\partial }^2\epsilon }{\partial x^2}\right)+\mu \dot{\epsilon }.$$

(30)

It is now possible to provide an analysis of Equation (30) and its fractional counterpart to interpret, among other things, size effects in creep. For a constant applied stress $(\sigma ={\sigma }_0)$, Equation (30) can be written in the form

$$\dot{\epsilon }=A\frac{{\partial }^2\epsilon }{\partial x^2}+B\epsilon +\mathrm{\Gamma },$$

(31)

where the constants $(A,B,\mathrm{\Gamma })$ relate in a straightforward manner to those in Equation (30), and a simple change of variables can reduce Equation (31) to the classical Fick diffusion equation. Other type of models, including a fourth-order derivative in $\epsilon$, as well as spatial derivatives of $\dot{\epsilon }$ can also be derived and some related spatio-temporal viscoelastic models will be discussed in the future, as this is not within the scope of the present paper.

3.2. A Fractional Extension

For the remainder of this Section, we discuss a fractional extension of the classical Kelvin–Voigt model given by Equation (26). By replacing the integer order derivative $\dot{\epsilon }=\partial \epsilon /\partial t$ with its time fractional counterpart, ${\partial }^{\beta }\epsilon /\partial t^{\beta }$, where $\beta$ is a non-integer number $0<\beta <1$, we first derive the corresponding creep equation associated with the fractional index $\beta$, i.e., the fractional counterpart of Equation (27).

One of the most commonly utilized choices for ${\partial }^{\beta }\epsilon /\partial t^{\beta }$ is the time-fractional derivative of Caputo type ${}_0^C{D}_t^{\beta }\epsilon$, defined as

$${}_0^C{D}_t^{\beta }\epsilon (x,t)=\frac{1}{\mathrm{\Gamma }(1-\beta )}{\int }_0^t\frac{1}{{(t-\tau )}^{\beta }}\frac{\partial \epsilon (x,\tau )}{\partial \tau }d\tau ,0<\beta <1,$$

(32)

which, for $\beta =1$, reduces to the usual first-order integer derivative, $\partial \epsilon /\partial t$. It then follows that the time-fractional counterpart of Equation (26) reads

$$\sigma =E\epsilon +{\mu }_{\beta }{}_0^C{D}_t^{\beta }\epsilon ,$$

(33)

where ${\mu }_{\beta }$ denotes the fractional viscosity coefficient. For a constant applied stress ${\sigma }_0$, Equation (33) is solved by applying the Laplace transform and taking into account the initial condition $\epsilon \left(0\right)=0$, resulting in the expression

$$\epsilon =\frac{{\sigma }_0}{E}\left[1-E_{\beta }\left(-\frac{t^{\beta }}{{\tau }_{\beta }}\right)\right],$$

(34)

where ${\tau }_{\beta }={\mu }_{\beta }/E$ (in units of ${\left[T\right]}^{\beta }$), and $E_{\beta }$ denotes the Mittag–Leffler function of fractional calculus [21,22], defined as $E_{\beta }\left(\zeta \right)={\sum }_{\nu =0}^{\infty }{\zeta }^{\nu }/\mathrm{\Gamma }(\beta \nu +1)$. The Mittag–Leffler function is a proper generalization of the exponential function, since $E_1\left(\zeta \right)=exp\left(\zeta \right);\zeta \in \mathbb{R}$. Therefore, the above expression reduces to the integer-order Kelvin–Voigt expression of Equation (27) for $\beta \to 1$. It follows that, due to the presence of $E_{\beta }$, Equation (34) exhibits a much slower than exponential decay $\sim O\left(t^{-\beta }\right)$ as $t\to \infty$.

Next, we extend the time-fractional Kelvin–Voigt model of Equation (33) by also including a spatial gradient term proportional to the second derivative of the strain ${\partial }_x^2\epsilon$ in its elastic component. In relation to this, it is noted that a Laplacian term proportional to $\dot{\epsilon }$ (i.e., ${\partial }_x^2\dot{\epsilon }$ in one dimension) could also be considered as suggested by the last author [2] and partially implemented by others [23,24], but this is outside the scope of the present article. This is also the case for other possible models including the fractional Laplacians of both $\epsilon$ and $\dot{\epsilon }$, as well as $\sigma$ and $\dot{\sigma }$.

On returning to a particular illustrative example, we consider the following time-fractional counterpart of Equation (30),

$$\sigma =E\left(\epsilon -{\ell }^2\frac{{\partial }^2\epsilon }{\partial x^2}\right)+{\mu }_{\beta }{}_0^C{D}_t^{\beta }\epsilon .$$

(35)

We further consider the problem of a bar under constant stress ${\sigma }_0$. In this case, Equation (35) is solved through the combined application of Fourier sine (in space) and Laplace (in time) transforms, by also taking into account the initial condition ${\epsilon |}_{t\to 0}=0$, as well as the boundary condition ${\epsilon |}_{x\to 0}=0$. The resulting expression reads

$$\epsilon (x,t)=\frac{{\sigma }_0}{E}{\int }_0^{\infty }{e}^{-\zeta }{\mathcal{W}}_{-\beta ,1}\left(-\zeta {\tau }_{\beta }{t}^{-\beta }\right)erf\left(\frac{x}{2\ell \sqrt{\zeta }}\right)d\zeta ,$$

(36)

where ${\mathcal{W}}_{-\beta ,1}$ is a particular case of the Wright function ${\mathcal{W}}_{\lambda ,\mu }\left(\zeta \right)={\sum }_{\nu =0}^{\infty }{\zeta }^{\nu }/\mathrm{\Gamma }(\lambda \nu +\mu )\nu !$, and $erf\left(\zeta \right)$ denotes the error function, defined as $erf\left(\zeta \right)=2{\pi }^{-1/2}{\int }_0^{\zeta }{e}^{-{\xi }^2}d\xi$. The behavior of Equation (36) in the limit $t\to \infty$ is given by the formula of Equation (4), since ${\epsilon |}_{t\to \infty }\approx {\sigma }_0/E{\int }_0^{\infty }{e}^{-\zeta }erf\left(\frac{x}{2\ell \sqrt{\zeta }}\right)d\zeta$. Similarly, in the limit $\ell \to 0$, it turns out that Equation (36) reduces to the time-fractional solution of Equation (34).

As a final illustration of the various possible models, we consider the fractional extension of the gradient Kelvin–Voigt model of Equation (33) by replacing the first-order time derivative ${\partial }_t$ with the Caputo fractional time derivative ${}_0^C{D}_t^{\beta }$, as well as the second spatial derivative ${\partial }_x^2$ with the Riesz fractional Laplacian $-{(-\mathrm{\Delta })}^{\alpha /2}$. The resulting space–time-fractional gradient Kelvin–Voigt constitutive equation reads

$$\sigma =E\left[\epsilon +{\ell }^{\alpha }{\left(-\mathrm{\Delta }\right)}^{\alpha /2}\epsilon \right]+{\mu }_{\beta }{}_0^C{D}_t^{\beta }\epsilon ,$$

(37)

where ${\ell }^{\alpha }$ is a fractional internal length, and ${\mu }_{\beta }$ is the corresponding fractional viscosity coefficient. For a constant applied stress ${\sigma }_0$, Equation (37) is solved by employing the Fourier sine transform (see Section 2) combined with the Laplace transform. Inversion of both transforms yields the following expression for $\epsilon (x,t)$

$$\epsilon (x,t)=\frac{{\sigma }_0}{E}{\int }_0^{\infty }{e}^{-\zeta }{\mathcal{W}}_{-\beta ,1}\left(-\zeta {\tau }_{\beta }{t}^{-\beta }\right){\Phi }_{\alpha }\left(\frac{x}{2\ell {\zeta }^{\frac{1}{\alpha }}}\right)d\zeta ,$$

(38)

where ${\Phi }_{\alpha }\left(\zeta \right)$ is the fractional generalization of the error function, defined through the expression

$${\Phi }_{\alpha }\left(\zeta \right)\frac{4}{\alpha \sqrt{\pi }}{\int }_0^{\zeta }{}_1{\Psi }_1\left[\begin{array}{c}\left(\frac{1}{\alpha },\frac{2}{\alpha }\right)\hfill \\ \left(\frac{1}{2},1\right)\hfill \end{array};-{\xi }^2\right]d\xi .$$

(39)

In Equation (39), we introduced the generalized Wright function, ${}_1{\Psi }_1$, which is an extension of the generalized Hypergeometric function ${}_p{F}_q$ and defined as

$${}_p{\Psi }_q\left(\zeta \right)={}_p{\Psi }_q\left[\begin{array}{ccc}\left(a_1,A_1\right)\hfill & \dots \hfill & \left(a_p,A_p\right)\hfill \\ \left(b_1,B_1\right)\hfill & \dots \hfill & \left(b_q,B_q\right)\hfill \end{array};\zeta \right]=\sum _{\nu =0}^{\infty }\frac{{\prod }_{j=1}^p\mathrm{\Gamma }(a_j+A_j\nu )}{{\prod }_{j=1}^q\mathrm{\Gamma }(b_j+B_j\nu )}\frac{{\zeta }^{\nu }}{\nu !}.$$

(40)

By considering the long time asymptotics (as $t\to \infty )$ of Equation (38) or equivalently, when the viscosity effects are negligible (in the limit ${\mu }_{\beta }\to 0$), it can be confirmed that the solution of Equation (38) reduces to the expression for the strain of the one-dimensional bar given by Equation (21) which, as shown earlier, is also equivalent to Equation (25). Furthermore, if spatial gradients are neglected, i.e., when $\ell \to 0$, the integral expression of Equation (38) is readily evaluated, yielding the time-fractional relaxation solution given by Equation (34). Finally, in the limit $\alpha \to 2$, Equation (38) reduces to the integer-order solution of Equation (36), due to the fact that ${\Phi }_{\alpha }\left(\zeta \right)$ coincides with $erf\left(\zeta \right)$ in that limit. Exact analytical expressions for Equation (38) and its asymptotics, along with further generalizations and specific applications, will be provided in a future publication.

4. Conclusions and Discussion

A concise account on certain 1D gradient/fractional models for elasticity and viscoelasticity was provided to illustrate the many various possibilities that exist for generalizing their classical counterparts. The choice of the most suitable generalized models depends on their simplicity and their ability to capture new phenomena and experimental observations that their classical counterparts are not able to describe.

One-dimensional strain-hardening models and other rheological (e.g. Maxwell’s viscoelastic and Brinkman’s plastic) models can also be enhanced with gradient terms and their implications can be discussed in a similar manner. Their 3D counterparts can also be addressed, since the fractional Laplacian has been well studied, and a substantial number of results are available [10].

It follows that this type of gradient and fractional generalizations of classical deformation models will only deserve further attention when they can describe phenomena (pattern formation, size effects) not captured by standard theories. In fact, these theories have been developed on the basis of observations and experiments for macroscopic bodies where the underlying micro/nano structural heterogeneity, and spatio-temporal pattern formation could be averaged out in the corresponding macroscopic model development. With the current growth of experimental probes to detect material behavior at the micron scale and nanoscale, the relevance of generalized non-local gradient and fractional/fractal models may be important. The present article is a modest contribution in this direction.

A final remark concerns the generalization of these one-dimensional (1D) models to three dimensions (3D). In addition to questions related to 3D equilibrium and compatibility requirements, which do not arise in a non-trivial manner in the present 1D case, one also needs to address the issue of frame indifference, i.e., “the invariance requirements under translations and rotations”, as a reviewer indicated by pointing, for example, to reference [25]. For a more general discussion on fractional aspects and invariance considerations in continuum mechanics formulations, the reader can consult a most interesting recent reference [26]. Another pertinent issue is the question of identifying the length scale ℓ in relation to three-dimensional micro/nano structures in various materials, which still remains an open problem to collectively elaborate upon.