Calculation of the Relaxation Modulus in the Andrade Model by Using the Laplace Transform

Abstract

In the framework of the theory of linear viscoelasticity, we derive an analytical expression of the relaxation modulus in the Andrade model $G_{\alpha }\left(t\right)$ for the case of rational parameter $\alpha =m/n\in (0,1)$ in terms of Mittag–Leffler functions from its Laplace transform ${\tilde{G}}_{\alpha }\left(s\right)$. It turns out that the expression obtained can be rewritten in terms of Rabotnov functions. Moreover, for the original parameter $\alpha =1/3$ in the Andrade model, we obtain an expression in terms of Miller-Ross functions. The asymptotic behaviours of $G_{\alpha }\left(t\right)$ for $t\to 0^{+}$ and $t\to +\infty$ are also derived applying the Tauberian theorem. The analytical results obtained have been numerically checked by solving the Volterra integral equation satisfied by $G_{\alpha }\left(t\right)$ by using a successive approximation approach, as well as computing the inverse Laplace transform of ${\tilde{G}}_{\alpha }\left(s\right)$ by using Talbot’s method.

1. The Andrade Model in Linear Viscoelasticity

In the framework of linear viscoelasticity theory, a “transient” phase of deformation occurs right after the elastic response in creep phenomena and is marked by a strain rate that changes over time [1]. Among the rheological laws that exhibit a transient phase, the Andrade model has been effectively used to describe the behavior of various materials. This model was initially introduced by Andrade in 1910 to describe the elongation of metallic wires under constant tensile stress [2]. Its main feature is a transient that exhibits a fractional power function time dependence ∼$t^{\alpha }$. In their empirical stress–strain relationship, Andrade proposed the exponent $\alpha =1/3$ [3]. Nevertheless, later laboratory investigations have shown that values within the range of $0<\alpha <1$ are indeed possible for certain materials [4]. It is worth noting that during the last dozen years, there has been an increasing interest in the Andrade model [5,6]. Using the modern formulation of the Andrade model [7], the creep compliance $J_{\alpha }\left(t\right)$ (i.e., the strain per unit stress) is given by

(1)$J_{\alpha }\left(t\right)=J_U+\beta t^{\alpha }+{\displaystyle \frac{t}{\eta }},t\ge 0,$

where $J_U$ is the unrelaxed compliance, $\eta$ is the steady state Newtonian viscosity, $\beta$ is the magnitude of the inelastic contribution, and $\alpha$ represents the frequency of the compliance. The number of free parameters that appear in (1) can be reduced by adopting a useful parametrization given in [8], whose validity is discussed in [4]. Essentially, this parametrization performs the following change of variables:

(2)$J_U={\displaystyle \frac{1}{\mu }},\beta ={\displaystyle \frac{1}{\mu {\tau }^{\alpha }}},\eta =\mu \tau ,$

thus, we obtain

(3)$J_{\alpha }\left(t\right)={\displaystyle \frac{1}{\mu }}\left[1+{\left({\displaystyle \frac{t}{\tau }}\right)}^{\alpha }+{\displaystyle \frac{t}{\tau }}\right],t\ge 0,$

where the parameters $\mu ,\tau >0$ have a clear and important physical meaning. The aim of this paper is to analytically calculate in closed form the relaxation modulus $G_{\alpha }\left(t\right)$ (i.e., the stress per unit strain) from the creep compliance $J_{\alpha }\left(t\right)$ given in (3) by using the inverse Laplace transform.

It is worth noting that the computation of $G_{\alpha }\left(t\right)$ from $J_{\alpha }\left(t\right)$ is also possible by numerically solving the following Volterra integral equation of the second kind ([9] [Eqn. 2.87]):

(4)$G_{\alpha }\left(t\right)=\mu \left[1-{\int }_0^t{\displaystyle \frac{d{J}_{\alpha }\left(t^{\prime }\right)}{d{t}^{\prime }}}{G}_{\alpha }\left(t-t^{\prime }\right)d{t}^{\prime }\right],$

just as has been done for other models such as that of Jeffrey-Lomnitz rheological law [10]. Another numerical approach is to obtain the relaxation modulus in the Laplace domain ${\tilde{G}}_{\alpha }\left(s\right)=\mathcal{L}\left[G_{\alpha }\left(t\right);s\right]$, and then numerically evaluate the inverse Laplace transform to obtain $G_{\alpha }\left(t\right)$. However, an analytical solution is more desirable since it shows the role and weight of the model parameters explicitly. In addition, analytical approaches are generally more computationally efficient. In the present approach, we just analytically calculate the inverse Laplace transform of ${\tilde{G}}_{\alpha }\left(s\right)$. This result has been applied by some of the authors in [11], and here we provide the mathematical details of the calculation.

Next, we derive ${\tilde{G}}_{\alpha }\left(s\right)$ in the Andrade model. For this purpose, apply the Laplace transform to (3) in order to obtain

(5)$\mathcal{L}\left[J_{\alpha }\left(t\right);s\right]={\tilde{J}}_{\alpha }\left(s\right)={\displaystyle \frac{1}{\mu s^2}}\left[s+\mathrm{\Gamma }\left(1+\alpha \right){\tau }^{-\alpha }{s}^{1-\alpha }+{\displaystyle \frac{1}{\tau }}\right].$

However, since the following relation is satisfied in any linear viscoelasticity rheology ([9] [Eqn. 2.8])

(6)${\tilde{J}}_{\alpha }\left(s\right){\tilde{G}}_{\alpha }\left(s\right)={\displaystyle \frac{1}{s^2}},$

we conclude

(7)${\tilde{G}}_{\alpha }\left(s\right)={\displaystyle \frac{\mu \tau }{s\tau +\mathrm{\Gamma }\left(\alpha +1\right){\left(s\tau \right)}^{1-\alpha }+1}}.$

As mentioned above, the range of values that are interesting for $\alpha$ are between 0 and 1, so we restrict our study to $\alpha \in (0,1)$. It is worth noting that (7) can also be obtained rewriting the Volterra integral Equation (4) in terms of the Laplace convolution product [12], i.e.,

(8)$G_{\alpha }\left(t\right)=\mu \left(1-G_{\alpha }\left(t\right)\ast {\displaystyle \frac{d{J}_{\alpha }\left(t\right)}{dt}}\right),$

thus taking the Laplace transform in (8) and solving for ${\tilde{G}}_{\alpha }\left(s\right)$, we arrive at (7).

The manuscript is organized as follows. In Section 2, we perform the calculation of $G_{\alpha }\left(t\right)$ for $\alpha =1/3$ (as first suggested by Andrade) in terms of Miller-Ross functions. Section 3 generalizes this result for rational $\alpha$ in terms of a finite sum of Mittag–Leffler functions, which in turn can be expressed as a linear combination of Rabotnov functions. It is worth highlighting that although in principle $\alpha$ can be a real number, the value it actually acquires in rheological models is a positive fractional number less than unity. Section 4 calculates the asymptotic behaviour of $G_{\alpha }\left(t\right)$ for $t\to 0^{+}$ and $t\to +\infty$ by using the Tauberian theorem. Section 5 shows some numerical verifications on the expressions of $G_{\alpha }\left(t\right)$ derived in the previous sections, as well as on their asymptotic behaviors. Finally, we collect our conclusions in Section 6.

2. Laplace Inversion in Terms of Miller-Ross Functions for $\mathbf{\alpha }=\mathbf{1}/\mathbf{3}$

Consider in (7) the following change of variables:

(9)$b={\displaystyle \frac{1}{\tau }},c={\displaystyle \frac{\mathrm{\Gamma }\left(\alpha +1\right)}{{\tau }^{\alpha }}},$

to rewrite ${\tilde{G}}_{\alpha }\left(s\right)$ as

(10)${\tilde{G}}_{\alpha }\left(s\right)={\displaystyle \frac{\mu }{b+s+c{s}^{1-\alpha }}}.$

In order calculate the inverse Laplace transform of ${\tilde{G}}_{\alpha }\left(s\right)$ for the case $\alpha =1/3$, i.e.,

(11)${\tilde{G}}_{1/3}\left(s\right)={\displaystyle \frac{\mu }{b+s+c{s}^{2/3}}},$

apply the identity

(12)$\left(x+y\right)\left(x^2-xy+y^2\right)=x^3+y^3,$

taking $x=b+s$, and $y=c{s}^{2/3}$, to arrive at

(13)${\displaystyle \frac{1}{\mu }}{\tilde{G}}_{1/3}\left(s\right)=b^2\underset{\tilde{R}\left(s\right)}{\underbrace{{\displaystyle \frac{1}{p\left(s\right)}}}}+2b\underset{\tilde{Q}\left(s\right)=s\tilde{R}\left(s\right)}{\underbrace{s{\displaystyle \frac{1}{p\left(s\right)}}}}+\underset{s\tilde{Q}\left(s\right)}{\underbrace{{\displaystyle \frac{s^2}{p\left(s\right)}}}}-bc\underset{{\tilde{U}}_{2/3}\left(s\right)}{\underbrace{{\displaystyle \frac{s^{2/3}}{p\left(s\right)}}}}-c\underset{s{\tilde{U}}_{2/3}\left(s\right)}{\underbrace{s{\displaystyle \frac{s^{2/3}}{p\left(s\right)}}}}+c^2\underset{{\tilde{U}}_{4/3}\left(s\right)}{\underbrace{{\displaystyle \frac{s^{4/3}}{p\left(s\right)}}}},$

where

(14)$p\left(s\right)=s^3+\left(3b+c^3\right)s^2+3b^{2}s+b^3=\prod _{k=1}^3\left(s-s_k\right),$

and $s_k$ are the roots of the cubic equation $p\left(s\right)=0$, i.e.,

(15)$\begin{array}{c}{\displaystyle s_1=-\frac{1}{3}\left(M-N-L\right),}\hfill \\ {\displaystyle s_2=-\frac{1}{3}\left(M+e^{i\pi /3}N+e^{-i\pi /3}L\right),}\hfill \\ {\displaystyle s_3=-\frac{1}{3}\left(M+e^{-i\pi /3}N+e^{i\pi /3}L\right),}\hfill \end{array}$

being

(16)$\begin{array}{c}{\displaystyle M=3b+c^3,}\hfill \\ {\displaystyle L=\sqrt[3]{\frac{3}{2}\sqrt{3b^3{c}^6\left(27b+4c^3\right)}-\frac{27}{2}{b}^2{c}^3-9b{c}^6-c^9},}\hfill \\ {\displaystyle N=\frac{\left(6b+c^3\right)c^3}{L}.}\hfill \end{array}$

First, rewrite $\tilde{R}\left(s\right)$ as

(17)$\tilde{R}\left(s\right)={\displaystyle \frac{1}{p\left(s\right)}}=\sum _{k=1}^3{\displaystyle \frac{{\alpha }_k}{s-s_k}},{\alpha }_k=\prod _{j\ne k}^3{\displaystyle \frac{1}{s_k-s_j}}.$

Note that a simple algebraic calculation shows that

(18)$\begin{array}{c}{\displaystyle \sum _{k=1}^3{\alpha }_k=0,}\hfill \\ {\displaystyle \sum _{k=1}^3{\alpha }_k{s}_k=0.}\hfill \end{array}$

Now, define the following functions

(19)$R\left(t\right)={\mathcal{L}}^{-1}\left[\tilde{R}\left(s\right);t\right],$

thus,

(20)$Q\left(t\right)={\mathcal{L}}^{-1}\left[\tilde{Q}\left(s\right);t\right]={\mathcal{L}}^{-1}\left[s\tilde{R}\left(s\right);t\right]=R^{\prime }\left(t\right)+R\left(0\right)\delta \left(t\right),$

and

(21)${\mathcal{L}}^{-1}\left[s\tilde{Q}\left(s\right);t\right]=Q^{\prime }\left(t\right)+Q\left(0\right)\delta \left(t\right)=R^{\prime \prime }\left(t\right)+R^{\prime }\left(0\right){\delta }^{\prime }\left(t\right).$

Furthermore,

(22)$U_{\nu }\left(t\right)={\mathcal{L}}^{-1}\left[{\tilde{U}}_{\nu }\left(s\right);t\right]={\mathcal{L}}^{-1}\left[{\displaystyle \frac{s^{\nu }}{p\left(s\right)}};t\right],$

thus

(23)${\mathcal{L}}^{-1}\left[s{\tilde{U}}_{\nu }\left(s\right);t\right]=U_{\nu }^{\prime }\left(t\right)+U_{\nu }\left(0\right)\delta \left(t\right).$

Therefore, the inverse Laplace transform of ${\tilde{G}}_{1/3}\left(s\right)$ is given by

(24)$\begin{array}{ccc}\hfill G_{1/3}\left(t\right)& =& {\mathcal{L}}^{-1}\left[{\tilde{G}}_{1/3}\left(s\right);t\right]\hfill \\ & =& \mu \left\{b^{2}R\left(t\right)+2b\left[R^{\prime }\left(t\right)+R\left(0\right)\delta \left(t\right)\right]+R^{\prime \prime }\left(t\right)+R^{\prime }\left(0\right){\delta }^{\prime }\left(t\right)\right.\hfill \\ & & +\left.bc{U}_{2/3}\left(t\right)-c\left[U_{2/3}^{\prime }\left(t\right)+U_{2/3}\left(0\right)\delta \left(t\right)\right]+c^2{U}_{4/3}\left(t\right)\right\}.\hfill \end{array}$

2.1. Calculation of $R\left(T\right)$

According to (17) and (19), we have

(25)$R\left(t\right)={\mathcal{L}}^{-1}\left[\tilde{R}\left(s\right);t\right]=\sum _{k=1}^3{\alpha }_k{\mathcal{L}}^{-1}\left[{\displaystyle \frac{1}{s-s_k}};t\right]=\sum _{k=1}^3{\alpha }_{k}exp\left(s_{k}t\right),$

so, that

(26)$\begin{array}{c}{\displaystyle R^{\prime }\left(t\right)=\sum _{k=1}^3{\alpha }_k{s}_{k}exp\left(s_{k}t\right),}\\ {\displaystyle R^{″}\left(t\right)=\sum _{k=1}^3{\alpha }_k{s}_k^{2}exp\left(s_{k}t\right).}\end{array}$

Therefore, from (18), we obtain

(27)$\begin{array}{c}{\displaystyle R\left(0\right)=\sum _{k=1}^3{\alpha }_k=0,}\hfill \\ {\displaystyle R^{\prime }\left(0\right)=\sum _{k=1}^3{\alpha }_k{s}_k=0.}\hfill \end{array}$

2.2. Calculation of $U_{\nu }\left(T\right)$

According to (22), (17), and the inverse Laplace transform ([13] [Eqn. 2.1.2(9)])

(28)${\mathcal{L}}^{-1}\left[{\displaystyle \frac{s^{\nu }}{s-a}};t\right]=a^{\nu }{e}^{at}P\left(-\nu ,at\right),$

where

(29)$P\left(\nu ,x\right)={\displaystyle \frac{\gamma \left(\nu ,x\right)}{\mathrm{\Gamma }\left(\nu \right)}}={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\nu \right)}}{\int }_0^x{z}^{\nu -1}{e}^{-z}dz,$

denotes the normalized lower incomplete gamma function ([14] [Eqn. 8.2.4]), we have

(30)$\begin{array}{ccc}\hfill U_{\nu }\left(t\right)& =& {\mathcal{L}}^{-1}\left[{\displaystyle \frac{s^{\nu }}{p\left(s\right)}};t\right]\hfill \\ & =& \sum _{k=1}^3{\alpha }_k{\mathcal{L}}^{-1}\left[{\displaystyle \frac{s^{\nu }}{s-s_k}};t\right]\hfill \\ & =& \sum _{k=1}^3{\alpha }_k{s}_k^{\nu }exp\left(s_{k}t\right)P\left(-\nu ,s_{k}t\right).\hfill \end{array}$

Therefore, apply the property ([14] [Eqn. 8.4.2])

(31)$\underset{x\to 0}{lim}{x}^{\nu }P\left(-\nu ,x\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(1-\nu \right)}},$

to calculate, according to (18), that

(32)$\begin{array}{ccc}\hfill U_{\nu }\left(0\right)& =& \underset{t\to 0}{lim}{t}^{-\nu }\sum _{k=1}^3{\alpha }_k{\left(s_{k}t\right)}^{\nu }P\left(-\nu ,s_{k}t\right)\hfill \\ & =& \underset{t\to 0}{lim}{\displaystyle \frac{t^{-\nu }}{\mathrm{\Gamma }\left(1-\nu \right)}}\sum _{k=1}^3{\alpha }_k=0.\hfill \end{array}$

Furthermore, from the derivative formula

(33)${\displaystyle \frac{d}{dt}}\left[e^{at}P\left(-\nu ,at\right)\right]=a\left[e^{at}P\left(-\nu ,at\right)+{\displaystyle \frac{{\left(at\right)}^{-\nu -1}}{\mathrm{\Gamma }\left(-\nu \right)}}\right],$

we calculate, according also to (18), that

(34)$\begin{array}{ccc}\hfill U_{\nu }^{\prime }\left(t\right)& =& \sum _{k=1}^3{\alpha }_k{s}_k^{\nu +1}\left[exp\left(s_{k}t\right)P\left(-\nu ,s_{k}t\right)+{\displaystyle \frac{{\left(s_{k}t\right)}^{-\nu -1}}{\mathrm{\Gamma }\left(-\nu \right)}}\right]\hfill \\ & =& \sum _{k=1}^3{\alpha }_k{s}_k^{\nu +1}exp\left(s_{k}t\right)P\left(-\nu ,s_{k}t\right)+{\displaystyle \frac{t^{-\nu -1}}{\mathrm{\Gamma }\left(-\nu \right)}}\sum _{k=1}^3{\alpha }_k\hfill \\ & =& \sum _{k=1}^3{\alpha }_k{s}_k^{\nu +1}exp\left(s_{k}t\right)P\left(-\nu ,s_{k}t\right).\hfill \end{array}$

2.3. Calculation of $G_{1/3}\left(T\right)$

Insert (25)–(34) into (24) to arrive at the following result, after simplification,

(35)$\begin{array}{ccc}& & G_{1/3}\left(t\right)=\mu \sum _{k=1}^3{\alpha }_{k}exp\left(s_{k}t\right)\hfill \\ & & \left\{\left(b+s_k\right)\left[b+s_k-c{s}_k^{2/3}P\left(-{\displaystyle \frac{2}{3}},s_{k}t\right)\right]+c^2{s}_k^{4/3}P\left(-{\displaystyle \frac{4}{3}},s_{k}t\right)\right\},t\ge 0.\hfill \end{array}$

It is worth noting that we can rewrite (35) in terms of the Miller-Ross functions ([9] [Eqn. E.37]), defined as

(36)${\mathrm{E}}_t\left(\nu ,a\right)={\displaystyle \frac{a^{-\nu }{e}^{at}}{\mathrm{\Gamma }\left(\nu \right)}}\gamma \left(\nu ,at\right)=a^{-\nu }{e}^{at}P\left(\nu ,at\right),$

thus, after simplification, we arrive at the following result:

(37)$G_{1/3}^{\mathrm{MR}}\left(t\right)=\mu \sum _{k=1}^3{\alpha }_k\left\{\left(b+s_k\right)\left[\left(b+s_k\right)e^{s_{k}t}-c{\mathrm{E}}_t\left(-{\displaystyle \frac{2}{3}},s_k\right)\right]+c^2{\mathrm{E}}_t\left(-{\displaystyle \frac{4}{3}},s_{k}t\right)\right\},$

where the superscript $\mathrm{MR}$ takes into account that $G_{1/3}\left(t\right)$ is given in terms of Miller-Ross functions. Furthermore, according to (9), the parameters $b=1/\tau$ and $c=\mathrm{\Gamma }\left({\displaystyle \frac{4}{3}}\right){\tau }^{-1/3}$. Moreover, $s_k$ and ${\alpha }_k$ are given in (15)–(17), respectively.

3. Laplace Inversion in Terms of Rabotnov Functions

Consider (7) for the case $\alpha ={\displaystyle \frac{m}{n}}\in \mathbb{Q}$,

(38)$0<m<n,n,m\in \mathbb{N},$

and perform the change of variables $r={\left(s\tau \right)}^{1/n}$,

(39)${\tilde{G}}_{m/n}\left(r\right)={\displaystyle \frac{\mu \tau }{\underset{p_{n,m}\left(r\right)}{\underbrace{r^n+\mathrm{\Gamma }\left(1+\frac{m}{n}\right)r^{n-m}+1}}}},$

where $p_{n,m}\left(r\right)$ is a polynomial of n-th order. If $p_{n,m}\left(r\right)$ has non-repeated roots $r_k$, $\left(k=1,\dots ,n\right)$, then, according to ([15] [Eqn. 17:13:10]), we have

(40)${\displaystyle \frac{1}{p_{n,m}\left(r\right)}}=\sum _{k=1}^n{\displaystyle \frac{1}{p_{n,m}^{\prime }\left(r_k\right)\left(r-r_k\right)}},$

thus

(41)${\tilde{G}}_{m/n}\left(r\right)=\mu \tau \sum _{k=1}^n{\displaystyle \frac{1}{p_{n,m}^{\prime }\left(r_k\right)\left(r-r_k\right)}},$

and

(42)${\tilde{G}}_{m/n}\left(s\right)=\mu \tau \sum _{k=1}^n{\displaystyle \frac{1}{p_{n,m}^{\prime }\left(r_k\right)\left(s^{1/n}{\tau }^{1/n}-r_k\right)}}.$

Define,

(43)$\begin{array}{ccc}\hfill f_{m/n}\left(t\right)& =& {\mathcal{L}}^{-1}\left[{\tilde{G}}_{m/n}\left(s/\tau \right);t\right]\hfill \\ & =& \mu \tau \sum _{k=1}^n{\displaystyle \frac{1}{p_{n,m}^{\prime }\left(r_k\right)}}{\mathcal{L}}^{-1}\left[{\displaystyle \frac{1}{s^{1/n}-r_k}}\right],\hfill \end{array}$

and apply ([15] [Eqn. 45:14:4])

(44)${\mathcal{L}}^{-1}\left[{\displaystyle \frac{s^{\mu -\nu }}{s^{\mu }-a}};t\right]=t^{\nu -1}{\mathrm{E}}_{\mu ,\nu }\left(a{t}^{\mu }\right),$

where ${\mathrm{E}}_{\alpha ,\beta }\left(z\right)$ denotes the two-parameter Mittag–Leffler function ([15] [Eqn. 45:14:2]),

(45)${\mathrm{E}}_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }{\displaystyle \frac{z^k}{\mathrm{\Gamma }\left(\beta +\alpha k\right)}}.$

Therefore, for $\mu =\nu ={\displaystyle \frac{1}{n}}$, we have

(46)${\mathcal{L}}^{-1}\left[{\displaystyle \frac{1}{s^{1/n}-a}};t\right]=t^{1/n-1}{\mathrm{E}}_{{\displaystyle \frac{1}{n}},{\displaystyle \frac{1}{n}}}\left(a{t}^{1/n}\right).$

Insert (46) in (43) to arrive at

(47)$f_{m/n}\left(t\right)=\mu \tau t^{1/n-1}\sum _{k=1}^n{\displaystyle \frac{{\mathrm{E}}_{\frac{1}{n},\frac{1}{n}}\left(r_k{t}^{1/n}\right)}{p_{n,m}^{\prime }\left(r_k\right)}}.$

Finally, apply the property ([13] [Eqn. 1.1.1(3)]),

(48)${\mathcal{L}}^{-1}\left[\tilde{F}\left(s\right);t\right]=F\left(t\right)\leftrightarrow {\mathcal{L}}^{-1}\left[\tilde{F}\left(as\right);t\right]={\displaystyle \frac{1}{a}}F\left({\displaystyle \frac{t}{a}}\right),$

to obtain

(49)$G_{m/n}\left(t\right)={\mathcal{L}}^{-1}\left[{\tilde{G}}_{m/n}\left(s\right);t\right]={\displaystyle \frac{1}{\tau }}{f}_{m/n}\left({\displaystyle \frac{t}{\tau }}\right),$

i.e.,

(50)$G_{m/n}\left(t\right)=\mu {\left({\displaystyle \frac{t}{\tau }}\right)}^{1/n-1}\sum _{k=1}^n{\displaystyle \frac{{\mathrm{E}}_{\frac{1}{n},\frac{1}{n}}\left(r_k{\left(t/\tau \right)}^{1/n}\right)}{p_{n,m}^{\prime }\left(r_k\right)}},t\ge 0,$

where remember that $r_k$ are the n non-repeated roots of the polynomial:

(51)$p_{n,m}\left(r\right)=r^n+\mathrm{\Gamma }\left(1+{\displaystyle \frac{m}{n}}\right)r^{n-m}+1.$

Note that the solution is a linear combination of Rabotnov functions ([9] [Eqn. E.46]):

(52)${\mathrm{R}}_{\nu }\left(\mu ,t\right)=t^{\nu }{\mathrm{E}}_{\nu +1,\nu +1}\left(\mu t^{\nu +1}\right),$

thereby,

(53)$G_{m/n}\left(t\right)=\mu \sum _{k=1}^n{\displaystyle \frac{{\mathrm{R}}_{1/n-1}\left(r_k,t/\tau \right)}{p_{n,m}^{\prime }\left(r_k\right)}}.$

For the particular case $\alpha =m/n=1/3$, we have

(54)$G_{1/3}^{\mathrm{R}}\left(t\right)=\mu \sum _{k=1}^3{\displaystyle \frac{{\mathrm{R}}_{-2/3}\left(r_k,t/\tau \right)}{2\mathrm{\Gamma }\left(\frac{4}{3}\right)r_k+3r_k^2}},t\ge 0,$

where the superscript $\mathrm{R}$ takes into account that $G_{1/3}\left(t\right)$ is given in terms of Rabotnov functions, and $r_k$ are the three distinct roots of the cubic equation,

(55)$r^3+\mathrm{\Gamma }\left({\displaystyle \frac{4}{3}}\right)r^2+1=0,$

i.e.,

(56)$\begin{array}{c}{r}_1\approx -1.40184,\hfill \\ r_2\approx 0.254432-0.805364i,\hfill \\ r_3\approx 0.254432+0.805364i.\hfill \end{array}$

4. Asymptotic Behaviour via Tauberian Theorem

Next, we will obtain the asymptotic behaviour of the relaxation modulus $G_{\alpha }\left(t\right)$ as $t\to 0^{+}$ and as $t\to +\infty$ from its Laplace transform ${\tilde{G}}_{\alpha }\left(s\right)$ by using the following version of the Tauberian theorem, (for other version of the Tauberian theorem, see [16]).

Theorem 1. 

Consider that the Laplace transform of a function $f\left(t\right)$ is given by $\tilde{f}\left(s\right)=\mathcal{L}\left[f\left(t\right);s\right]$. The asymptotic behaviour of $\tilde{f}\left(s\right)$ as $s\to +\infty$ is given by

(57) $\tilde{f}\left(s\right)\approx \mathcal{L}\left[g\left(t\right);s\right],s\to +\infty ,$

where $g\left(t\right)$ is the asymptotic behaviour of $f\left(t\right)$ as $t\to 0^{+}$. Furthermore, the asymptotic behaviour of $\tilde{f}\left(s\right)$ as $s\to 0^{+}$ is given by

(58) $\tilde{f}\left(s\right)\approx \mathcal{L}\left[h\left(t\right);s\right],s\to 0^{+},$

where $h\left(t\right)$ is the asymptotic behaviour of $f\left(t\right)$ as $t\to +\infty$.

Proof. 

On the one hand, consider that the asymptotic behaviour of $f\left(t\right)$ as $t\to 0^{+}$ is given by

(59)$f\left(t\right)\approx g\left(t\right)=\sum _{k=1}^N{a}_k{t}^{b_k},t\to 0^{+},$

with $N=0,1,2,\dots$, and $0\le b_1<b_2<\dots <b_N$. Apply the Laplace transform to (59) in order to obtain

$\tilde{f}\left(s\right)=\mathcal{L}\left[f\left(t\right);s\right]\approx \mathcal{L}\left[g\left(t\right);s\right]=\sum _{k=1}^N{a}_k\mathcal{L}\left[t^{b_k};s\right]=\sum _{k=1}^N{\displaystyle \frac{a_k\mathrm{\Gamma }\left(b_k+1\right)}{s^{b_k+1}}}.$

Therefore, we obtain the asymptotic behaviour of $\tilde{f}\left(s\right)$ as $s\to +\infty$,

$\tilde{f}\left(s\right)\approx \sum _{k=0}^N{\displaystyle \frac{a_k\mathrm{\Gamma }\left(b_k+1\right)}{s^{b_k+1}}},s\to +\infty ,$

as we wanted to prove.

On the other hand, consider that the asymptotic behaviour of $\tilde{f}\left(s\right)$ as $s\to 0^{+}$ is given by

(60)$\tilde{f}\left(s\right)\approx \sum _{k=1}^N{c}_k{s}^{d_k},s\to 0^{+},$

with $N=0,1,2,\dots$, and $0\le d_1<d_2<\dots <d_N$. Apply the inverse Laplace transform to (60) in order to obtain

$f\left(t\right)={\mathcal{L}}^{-1}\left[\tilde{f}\left(s\right);t\right]=\sum _{k=1}^N{c}_k{\mathcal{L}}^{-1}\left[s^{d_k};t\right]=\sum _{k=1}^N{\displaystyle \frac{c_k}{\mathrm{\Gamma }\left(-d_k\right)t^{d_k+1}}},$

Therefore, we obtain the asymptotic behaviour of $f\left(t\right)$ as $t\to +\infty$,

$f\left(t\right)\approx h\left(t\right)=\sum _{k=1}^N{\displaystyle \frac{c_k}{\mathrm{\Gamma }\left(-d_k\right)t^{d_k+1}}},t\to +\infty ,$

as we wanted to prove. □

4.1. Asymptotic Behaviour for $T\to +\infty$

We know that the Laplace transform of the relaxation modulus in the Andrade model is (7)

(61)${\tilde{G}}_{\alpha }\left(s\right)={\displaystyle \frac{\mu \tau }{s\tau +\mathrm{\Gamma }\left(\alpha +1\right){\left(s\tau \right)}^{1-\alpha }+1}},\alpha \in \left(0,1\right).$

Since

(62)${\displaystyle \frac{1}{1-x}}=1+x+x^2+\dots ,\left|x\right|<1,$

we have that

(63)${\displaystyle \frac{1}{1+x}}\approx 1-x,x\to 0^{+},$

thus, taking $x=s\tau +\mathrm{\Gamma }\left(\alpha +1\right){\left(s\tau \right)}^{1-\alpha }$, we get

(64)${\tilde{G}}_{\alpha }\left(s\right)\approx \mu \tau \left[1-s\tau -\mathrm{\Gamma }\left(\alpha +1\right){\left(s\tau \right)}^{1-\alpha }\right],s\to 0^{+}.$

Since $\alpha \in \left(0,1\right)$,

(65)${\tilde{G}}_{\alpha }\left(s\right)\approx \mu \tau \left[1-\mathrm{\Gamma }\left(\alpha +1\right){\left(s\tau \right)}^{1-\alpha }\right],s\to 0^{+}.$

According to (58), the asymptotic behaviour of $G_{\alpha }\left(t\right)$ as $t\to +\infty$ is calculated as

(66)$\begin{array}{ccc}\hfill G_{\alpha }\left(t\right)& \approx & \mu \tau {\mathcal{L}}^{-1}\left[1-\mathrm{\Gamma }\left(\alpha +1\right){\left(s\tau \right)}^{1-\alpha };t\right]\hfill \\ \hfill & =& \mu \tau \delta \left(t\right)-\mu {\displaystyle \frac{\mathrm{\Gamma }\left(\alpha +1\right)}{\mathrm{\Gamma }\left(\alpha -1\right)}}{\left({\displaystyle \frac{t}{\tau }}\right)}^{\alpha -2}.\hfill \end{array}$

Applying the factorial property of the Gamma function, $\mathrm{\Gamma }\left(z+1\right)=z\mathrm{\Gamma }\left(z\right)$, and taking into account that $\delta \left(t\right)=0$ as $t\to +\infty$, we conclude that

(67)$G_{\alpha }\left(t\right)\approx \mu \alpha \left(1-\alpha \right){\left({\displaystyle \frac{t}{\tau }}\right)}^{\alpha -2},t\to +\infty .$

Taking more terms in the expansion of ${\tilde{G}}_{\alpha }\left(s\right)$ as $s\to 0^{+}$, we can calculate more terms of $G_{\alpha }\left(t\right)$ as $t\to +\infty$ by using (58). Thereby, we obtain

(68)$G_{\alpha }\left(t\right)\approx \mu \alpha \left[\left(1-\alpha \right){\left({\displaystyle \frac{t}{\tau }}\right)}^{\alpha -2}-2{\displaystyle \frac{{\mathrm{\Gamma }}^2\left(1+\alpha \right)}{\mathrm{\Gamma }\left(2\alpha -1\right)}}{\left({\displaystyle \frac{t}{\tau }}\right)}^{2\alpha -3}+\dots \right],t\to +\infty .$

For the particular case $\alpha ={\displaystyle \frac{1}{3}}$ in (67), we have

(69)$G_{1/3}\left(t\right)\approx {\displaystyle \frac{2}{9}}\mu {\left({\displaystyle \frac{t}{\tau }}\right)}^{-5/3},t\to +\infty .$

As a consistency test, we can obtain (69) from the the expression given in (50) for $G_{m/n}\left(t\right)$ with $m=1$ and $n=3$, and the asymptotic formula ([17] [Eqn. 18.1(22)]),

(70)${\mathrm{E}}_{\alpha ,\beta }\left(z\right)\approx -\sum _{k=1}^N{\displaystyle \frac{z^{-k}}{\mathrm{\Gamma }\left(\beta -\alpha k\right)}}+O\left({\left|z\right|}^{-N}\right),z\to \infty ,\left|arg\left(-z\right)\right|<\left(1-{\displaystyle \frac{\alpha }{2}}\right)\pi .$

4.2. Asymptotic Behaviour for $T\to 0^{+}$

Rewrite the Laplace transform of the relaxation modulus (61), as follows:

(71)${\tilde{G}}_{\alpha }\left(s\right)={\displaystyle \frac{\mu }{s}}\left[{\displaystyle \frac{1}{1+\mathrm{\Gamma }\left(\alpha +1\right){\left(s\tau \right)}^{-\alpha }+{\left(s\tau \right)}^{-1}}}\right].$

Note that, for $a\in \left(0,1\right)$ (thus $1-a\in \left(0,1\right)$), and $A\in \mathbb{R}$, we have that

(72)$\underset{y\to +\infty }{lim}{\displaystyle \frac{1+A{y}^{-a}+y^{-1}}{1+A{y}^{-a}}}=\underset{y\to +\infty }{lim}{\displaystyle \frac{y+A{y}^{1-a}+1}{y+A{y}^{1-a}}}=1,$

thus, taking $A=\mathrm{\Gamma }\left(\alpha +1\right)$, $a=\alpha \in \left(0,1\right)$, and $y=s\tau \to +\infty$, (i.e., $s\to +\infty$, since $\tau >0$), we get

(73)$1+\mathrm{\Gamma }\left(\alpha +1\right){\left(s\tau \right)}^{-\alpha }+{\left(s\tau \right)}^{-1}\approx 1+\mathrm{\Gamma }\left(\alpha +1\right){\left(s\tau \right)}^{-\alpha },s\to +\infty .$

Apply (73) to (71), in order to obtain

(74)${\tilde{G}}_{\alpha }\left(s\right)\approx {\displaystyle \frac{\mu }{s}}\left[{\displaystyle \frac{1}{1+\mathrm{\Gamma }\left(\alpha +1\right){\left(s\tau \right)}^{-\alpha }}}\right],s\to +\infty .$

Now, perform the change of variables $x=1/z$ in (63),

(75)${\displaystyle \frac{1}{1+\frac{1}{z}}}\approx 1-{\displaystyle \frac{1}{z}},z\to +\infty ,$

and take $z={\displaystyle \frac{{\left(s\tau \right)}^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}}\to +\infty$ (i.e., $s\to +\infty$), to arrive at

(76)${\tilde{G}}_{\alpha }\left(s\right)\approx \mu \left[{\displaystyle \frac{1}{s}}-{\displaystyle \frac{\mathrm{\Gamma }\left(\alpha +1\right)}{{\tau }^{\alpha }{s}^{\alpha +1}}}\right],s\to +\infty .$

According to (57), the asymptotic behaviour of $G_{\alpha }\left(t\right)$ as $t\to 0^{+}$ is calculated as

(77)$G_{\alpha }\left(t\right)\approx \mu {\mathcal{L}}^{-1}\left[{\displaystyle \frac{1}{s}}-{\displaystyle \frac{\mathrm{\Gamma }\left(\alpha +1\right)}{{\tau }^{\alpha }{s}^{\alpha +1}}};t\right],$

i.e.,

(78)$G_{\alpha }\left(t\right)\approx \mu \left[1-{\left({\displaystyle \frac{t}{\tau }}\right)}^{\alpha }\right],t\to 0^{+}.$

Again, taking more terms in the expansion of ${\tilde{G}}_{\alpha }\left(s\right)$ as $s\to +\infty$, we can calculate more terms of $G_{\alpha }\left(t\right)$ as $t\to 0^{+}$ by using (57). Thereby, we obtain

(79)$G_{\alpha }\left(t\right)\approx \mu \left[1-{\left({\displaystyle \frac{t}{\tau }}\right)}^{\alpha }+{\displaystyle \frac{{\mathrm{\Gamma }}^2\left(1+\alpha \right)}{\mathrm{\Gamma }\left(2\alpha -1\right)}}{\left({\displaystyle \frac{t}{\tau }}\right)}^{2\alpha }+\dots \right],t\to 0^{+}.$

Note that the particular case $\alpha ={\displaystyle \frac{1}{3}}$ in (78) yields

(80)$G_{1/3}\left(t\right)\approx \mu \left[1-{\left({\displaystyle \frac{t}{\tau }}\right)}^{1/3}\right],t\to 0^{+}.$

As a consistency test, we can obtain (80) from the the expression given in (50) for $G_{m/n}\left(t\right)$ with $m=1$ and $n=3$, and the definition of the Mittag–Leffler Function (45). Furthermore, the asymptotic formula given in (78) can be obtained from the Volterra integral Equation (4). Indeed, taking into account (2) and (3) and performing the change of variables $x=t-t^{\prime }$, this integral equation reads as

(81)$G_{\alpha }\left(t\right)=\mu -{\displaystyle \frac{1}{\tau }}{\int }_0^t\left[1+\alpha {\left({\displaystyle \frac{t-x}{\tau }}\right)}^{\alpha -1}\right]G_{\alpha }\left(x\right)dx,$

thus,

(82)$\underset{t\to 0^{+}}{lim}{G}_{\alpha }\left(t\right)=\mu .$

According to (82), we can take the approximation $G_{\alpha }\left(x\right)\approx \mu$ as $t\to 0^{+}$ in (81), thus

(83)$\begin{array}{ccc}\hfill G_{\alpha }\left(t\right)& \approx & \mu -{\displaystyle \frac{\mu }{\tau }}{\int }_0^t\left[1+\alpha {\left({\displaystyle \frac{t-x}{\tau }}\right)}^{\alpha -1}\right]dx\hfill \\ & =& \mu \left[1-{\displaystyle \frac{t}{\tau }}-{\left({\displaystyle \frac{t}{\tau }}\right)}^{\alpha }\right]\hfill \end{array}$

Recalling that $\alpha \in \left(0,1\right)$, we recover (78), i.e.,

(84)$G_{\alpha }\left(t\right)\approx \mu \left[1-{\left({\displaystyle \frac{t}{\tau }}\right)}^{\alpha }\right],t\to 0^{+}.$

Figure 1 presents the graph of $G_{1/3}\left(t\right)$ for $\mu =1$ and different values of $\tau$. Figure 2 shows the asymptotic behaviours given in (80) and (69) for $G_{1/3}\left(t\right)$ with $\mu =1$ and $\tau ={\displaystyle \frac{1}{2}}$.

5. Numerical Results

5.1. Volterra Integral Equation

The quadrature formulas to numerically solve the Volterra integral equation [18] that satisfies $G_{\alpha }\left(t\right)$, i.e., Equation (81), fail because from (80), recalling that $\alpha \in \left(0,1\right)$ and $\mu ,\tau >0$, we have that

(85)$\underset{t\to 0^{+}}{lim}{\displaystyle \frac{d{G}_{\alpha }\left(t\right)}{dt}}=-\infty .$

However, we can apply a successive approximation method in order to numerically compute $G_{\alpha }\left(t\right)$ ([19] [Sect. 2.1]). This method states that if we have the Volterra integral equation of the second kind

(86)$u\left(t\right)=f\left(t\right)+{\int }_0^{t}K\left(x,t\right)u\left(x\right)dx,$

we take as zeroth approximation

(87)$u^{\left(0\right)}\left(t\right)=f\left(t\right),$

and for the successive approximations $j=1,2,\dots$

(88)$u^{\left(j\right)}\left(t\right)=f\left(t\right)+{\int }_0^{t}K\left(x,t\right)u^{\left(j-1\right)}\left(x\right)dx.$

Figure 3 shows the application of this successive approximation method to the solution of (81) (i.e., taking as kernel $K\left(x,t\right)=-{\displaystyle \frac{1}{\tau }}\left[1+\alpha {\left({\displaystyle \frac{t-x}{\tau }}\right)}^{\alpha -1}\right]$, and $f\left(t\right)=\mu$ in (86)), for $\mu =\tau =1$ and $\alpha =1/3$. It is apparent that as the order of approximation increases, we get a better approximation to the analytical solution $G_{1/3}\left(t\right)$ obtained in (37) or (54). Similar graphs are obtained for other rational values of $\alpha \in \left(0,1\right)$ compared to the analytical solution $G_{\alpha }\left(t\right)$ obtained in (50).

Note that this successive approximation method has been successfully applied in Section 4 in order to derive the first order asymptotic formula (78).

5.2. Inverse Laplace Transform

According to our numerical experiments, the relative error between the analytical formulas of $G_{\alpha }\left(t\right)$ and the numerical Laplace inversion of ${\tilde{G}}_{\alpha }\left(s\right)$, never exceeds the value of ${10}^{-9}$ in the time interval $t\in \left[0,5\right]$. Below we present some of these numerical experiments.

Figure 4 shows the relative error ${\mathrm{\Delta }}_{\mathrm{MR}}\left(t\right)$ between $G_{1/3}^{\mathrm{MR}}\left(t\right)$ and $G_{1/3}^{\mathrm{num}}\left(t\right)$, i.e., the numerical Laplace inversion of ${\tilde{G}}_{1/3}\left(s\right)$ using Talbot’s method [20],

(89)${\mathrm{\Delta }}_{\mathrm{MR}}\left(t\right)=\left|1-{\displaystyle \frac{G_{1/3}^{\mathrm{MR}}\left(t\right)}{G_{1/3}^{\mathrm{num}}\left(t\right)}}\right|.$

Figure 5 shows the relative error ${\mathrm{\Delta }}_{\mathrm{R}}\left(t\right)$ between $G_{1/3}^{\mathrm{R}}\left(t\right)$ and $G_{1/3}^{\mathrm{num}}\left(t\right)$,

(90)${\mathrm{\Delta }}_{\mathrm{R}}\left(t\right)=\left|1-{\displaystyle \frac{G_{1/3}^{\mathrm{R}}\left(t\right)}{G_{1/3}^{\mathrm{num}}\left(t\right)}}\right|.$

6. Conclusions

Considering the Andrade model in linear viscoelasticity, we have derived for the first time an analytical expression for the relaxation modulus in the time domain $G_{\alpha }\left(t\right)$ considering a rational parameter $\alpha \in \left(0,1\right)$ in terms of Mittag–Leffler functions (or equivalently, as a linear combination of Rabotnov functions). For the original parameter $\alpha =1/3$ of the Andrade model, we have derived a particular expression for $G_{1/3}\left(t\right)$ in terms of Miller-Ross functions. It turns out that this last expression is numerically more efficient (approximately twice faster) than the equivalent one in terms of Rabotnov functions.

Furthermore, we have obtained the asymptotic behaviour of $G_{\alpha }\left(t\right)$ for $t\to 0^{+}$ and $t\to +\infty$ using the Tauberian theorem. We have derived the same expression for the asymptotic behaviour as $t\to 0^{+}$ by using the Volterra integral equation of the second kind that $G_{\alpha }\left(t\right)$ satisfies.

Finally, numerical computations for particular values of the parameters have been performed in order to verify the analytical solutions obtained. For this purpose, we have used Talbot’s method for the numerical computation of the inverse Laplace transform, and the method of successive approximations for the numerical evaluation of the Volterra integral equation of the second kind.

Author Contributions

Conceptualization, G.S. and F.M.; methodology, J.L.G.-S. and A.A.; software, J.L.G.-S.; validation, J.L.G.-S. and A.A.; formal analysis, J.L.G.-S.; investigation, G.S.; writing—original draft preparation, J.L.G.-S.; writing—review and editing, G.S.; supervision, G.S. and A.A.; project administration, G.S. and F.M. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Graph of [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.].

Figure 2. Asymptotic behaviour of [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 3. Successive approximations [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.].

Figure 4. Graph of [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.].

Figure 5. Graph of [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.].

References

1. Christensen, R.M. Theory of Viscoelasticity; Courier Corporation: North Chelmsford, MA, USA, 2003.

2. Andrade, E.N.D.C. On the viscous flow in metals, and allied phenomena. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character; Generic: Homebush West, NSW, Australia, 1910; Volume 84, pp. 1-12.

3. Andrade, E.N.D.C. The validity of the t^1/3 law of flow of metals. Philos. Mag.; 1962; 7, pp. 2003-2014. [DOI: https://dx.doi.org/10.1080/14786436208214469]

4. Walterová, M.; Plesa, A.C.; Wagner, F.W.; Breuer, D. Andrade Rheology in Planetary Science. Authorea Preprints. 2023. Available online: https://essopenarchive.org/doi/full/10.22541/essoar.169008299.99203183 (accessed on 25 July 2024).

5. Rosti, J.; Koivisto, J.; Laurson, L.; Alava, M. Fluctuations and scaling in creep deformation. Phys. Rev. Lett.; 2010; 105, 100601. [DOI: https://dx.doi.org/10.1103/PhysRevLett.105.100601] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20867504]

6. Pandey, V. Hidden jerk in universal creep and aftershocks. Phys. Rev.; 2023; 107, L022602. [DOI: https://dx.doi.org/10.1103/PhysRevE.107.L022602] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36932618]

7. Bagheri, A.; Khan, A.; Al-Attar, D.; Crawford, O.; Giardini, D. Tidal response of Mars constrained from laboratory-based viscoelastic dissipation models and geophysical data. J. Geophys. Res. Planets; 2019; 124, pp. 2703-2727. [DOI: https://dx.doi.org/10.1029/2019JE006015]

8. Castillo-Rogez, J.C.; Efroimsky, M.; Lainey, V. The tidal history of Iapetus: Spin dynamics in the light of a refined dissipation model. J. Geophys. Res. Planets; 2011; 116, E09008. [DOI: https://dx.doi.org/10.1029/2010JE003664]

9. Mainardi, F. Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models; World Scientific: Singapore, 2022.

10. Mainardi, F.; Spada, G. On the viscoelastic characterization of the Jeffreys–Lomnitz law of creep. Rheol. Acta; 2012; 51, pp. 783-791. [DOI: https://dx.doi.org/10.1007/s00397-012-0634-x]

11. Cosorzi, A.; Mellini, D.; González-Santander, J.; Spada, G. On the Love numbers of an Andrade planet. Earth Space Sci.; 2024.

12. Apelblat, A. Laplace Transforms and Their Applications; Mathematics Research Developments Series; Nova Science Publishers: Hauppauge, NY, USA, 2012.

13. Prudnikov, A.P.; Brychkov, Y.A.; Marichev, O.I. Integrals and Series: Inverse Laplace Transforms; CRC Press: Boca Raton, FL, USA, 1986; Volume 5.

14. Olver, F.W.J.; Lozier, D.W.; Boisvert, R.F.; Clark, C.W. NIST Handbook of Mathematical Functions; Cambridge University Press: Cambridge, UK, 2010.

15. Oldham, K.B.; Myland, J.; Spanier, J. An Atlas of Functions: With Equator, the Atlas Function Calculator; Springer: Berlin/Heidelberg, Germany, 2009.

16. Sandev, T.; Tomovski, Z. Fractional equations and models. Theory and Applications; Springer Nature Switzerland AG: Cham, Switzerland, 2019.

17. Erdélyi, A.; Magnus, W.; Oberhettinger, F.; Tricomi, F.G. Higher Transcendental Functions; McGraw-Hill: New York, NY, USA, 1955; Volume 3.

18. Rakotch, E. Numerical solution of Volterra integral equations. Numer. Math.; 1972; 20, pp. 271-279. [DOI: https://dx.doi.org/10.1007/BF01407369]

19. Tricomi, F.G. Integral Equations; Courier Corporation: North Chelmsford, MA, USA, 1985; Volume 5.

20. Talbot, A. The accurate numerical inversion of Laplace transforms. IMA J. Appl. Math.; 1979; 23, pp. 97-120. [DOI: https://dx.doi.org/10.1093/imamat/23.1.97]