Phenomenological isotropic visco-hyperelasticity: a differential model based on fractional derivatives

Abstract

The rheological equation of a standard linear solid, i.e., the Zener model, is thermodynamically consistent. Thus, it was often used as a starting point for the development of nonlinear viscoelastic models, especially for elastomers. The basic idea of this paper is a generalization of the one-dimensional fractional constitutive equation of the Zener model to large strains. To reduce the number of material parameters of differential models based on the concept of the internal variables and to avoid integral constitutive equations, we develop a differential model based on the concept of dual stress and strain tensors and their derivatives. To this end, we select two couples of dual stress and strain tensors that have been used in finite elasticity. Then we obtain two constitutive models of incompressible isotropic materials called M1 and M2. We show that the M1 model is not suitable for describing the viscoelastic behavior of elastomers. To improve the predictions of the M2 model, we assume that the material is thixotropic. Therefore, the ratio of the relaxation and creep time depends on deformation. Experimental results show that this ratio may be represented as a function of the first invariant of the Cauchy–Green strain tensor. This yields a new constitutive equation whose material parameters were identified using experimental data on relaxation loadings in the literature. Next, we show that the model is able to predict the experimental data for combined loads of tension–torsion. Consequently, the model seems to be efficient at predicting the multiaxial visco-hyperelastic behavior of elastomers. The main advantage of the current model is that it has a differential form with relatively few parameters and is mathematically convenient.

Appendices

Appendix 1: Resolution of fractional differential equation

The analytical solution of linear fractional differential equations (FDEs) with constant coefficients can be obtained using the Laplace transform technique.

In the case of relaxation, the Laplace transform formula for the Caputo fractional derivatives involves the value of the function f(t) at the lower terminal $t=0$:

$$\begin{aligned} \mathcal {L}[{D}^{\alpha }f(t)]=p^{\alpha }\mathcal {L}[f(t)]-f(0)p^{\alpha -1}. \end{aligned}$$

(68)

Let us consider the following FDE:

$$\begin{aligned} {\mathbf {S}}(t)+\tau _\mathrm{R}^\alpha {D}^\alpha {\mathbf {S}}(t)=\frac{\partial {W}}{\partial {\mathbf {E}(t)}}-2\mu _0\tau _\mathrm{R}^\alpha \xi ^{-1}(t) {D}^\alpha {\mathbf {e}}(t), \end{aligned}$$

(69)

with the initial condition ${\mathbf {S}}(0)=0$.

Applying the Laplace transform to both sides of this equation, and using the linearity of the Laplace transform, we obtain the following result:

$$\begin{aligned} \mathcal {L}\left[ {\mathbf {S}}(t)\right] +\tau _\mathrm{R}^\alpha \mathcal {L}\left[ {D}^\alpha {{\mathbf {S}}}(t)\right] = \mathcal {L}\left[ \frac{\partial {W}}{\partial {\mathbf {E}(t)}}\right] -2\mu _0\tau _\mathrm{R}^\alpha \mathcal {L}\left[ \xi ^{-1}(t){D}^\alpha {{\mathbf {e}}}(t)\right] . \end{aligned}$$

(70)

Then

$$\begin{aligned} {\hat{\mathbf {S}}}(p) + \tau _\mathrm{R}^\alpha p^\alpha \hat{\mathbf {S}}(p)={\hat{\mathbf {F}}}_1(p)-2\mu _0\tau _\mathrm{R}^ \alpha {\hat{\mathbf {F}}}_2(p), \end{aligned}$$

(71)

where ${\hat{\mathbf {F}}}_1(p)$ and ${\hat{\mathbf {F}}}_2(p)$ denote the Laplace transforms of $\partial {W}/\partial {\mathbf {E}(t)}$ and $\xi ^{-1}(t){D}^\alpha {{\mathbf {e}}}(t)$.

The equation’s solution is

$$\begin{aligned} {\hat{\mathbf {S}}}(p) =\left( \frac{\tau _\mathrm{R}^{-\alpha }}{p^\alpha +\tau _\mathrm{R}^{-\alpha }}\right) {\hat{\mathbf {F}}}_1(p)-2\mu _0\left( \frac{1}{p^\alpha +\tau _\mathrm{R}^{-\alpha }}\right) {\hat{\mathbf {F}}}_2(p). \end{aligned}$$

(72)

The solution in the time–space domain is obtained using the inverse Laplace transform, and we obtain

$$\begin{aligned} {\mathbf {S}}(t)= & {} \tau _\mathrm{R}^{-\alpha }\int \limits _{0}^{t}\frac{\partial W}{\partial \mathbf {E}(s)}(t-s)^{\alpha -1}E_\alpha ,_\alpha \Bigl (-\Bigl (\frac{t-s}{\tau _\mathrm{R}}\Bigr )^\alpha \Bigr ) {\mathrm {d}}s\nonumber \\&-2\mu _0\int \limits _{0}^{t}\xi ^{-1}(s){D}^\alpha {\mathbf {e}}(s)(t-s)^{\alpha -1}E_\alpha ,_\alpha \Bigl (-\Bigl (\frac{t-s}{\tau _\mathrm{R}}\Bigr )^\alpha \Bigr ) {\mathrm {d}}s. \end{aligned}$$

(73)

Application of the Laplace transform to the FDE, in the case of $\xi$ equal to a constant, gives

$$\begin{aligned} {\hat{\mathbf {S}}}(p) =\Bigl (\frac{\tau _\mathrm{R}^{-\alpha }}{p^\alpha +\tau _\mathrm{R}^{-\alpha }}\Bigr ){\hat{\mathbf {F}}}_1(p)-2\mu _0\xi ^{-1}\Bigl (\frac{p^\alpha }{p^\alpha +\tau _\mathrm{R}^{-\alpha }}\Bigr )\hat{ {\mathbf {e}}}(p). \end{aligned}$$

(74)

We have

$$\begin{aligned} \mathcal {L}^{-1}\left[ \left( \frac{p^\alpha }{p^\alpha +\tau _\mathrm{R}^{-\alpha }}\right) {\hat{\mathbf {e}}}(p)\right] =\mathcal {L}^{-1}\left[ \left( \frac{p^{\alpha -1}}{p^\alpha +\tau _\mathrm{R}^{-\alpha }}\right) \Bigr (p {\hat{\mathbf {e}}}(p)\Bigr )\right] . \end{aligned}$$

(75)

We obtain the solution as follows:

$$\begin{aligned} {\mathbf {S}}(t)=\tau _\mathrm{R}^{-\alpha }\int \limits _{0}^{t}\frac{\partial W}{\partial \mathbf {E}(s)}(t-s)^{\alpha -1}E_\alpha ,_\alpha \left( -\Bigl (\frac{t-s}{\tau _\mathrm{R}}\right) ^\alpha \Bigr ) {\mathrm {d}}s -2\mu _0\xi ^{-1}\int \limits _{0}^{t}{\dot{\mathbf {e}}}(s)E_\alpha ,_1\left( -\Bigl (\frac{t-s}{\tau _\mathrm{R}}\right) ^\alpha \Bigr ) {\mathrm {d}}s. \end{aligned}$$

(76)

Appendix 2: Mittag–Leffler function

The function $E_\alpha$, with $\alpha >0$ defined by $E_\alpha (z)=\sum _{k=0}^{\infty }z^k/\varGamma (\alpha k+1)$, is called the Mittag–Leffler function of order $\alpha$. This function provides a simple generalization of the exponential function, i.e., $E_1(z)=\mathrm {exp}(z)$.

When $0<\alpha \le 1$, the Mittag–Leffler function on the negative axis is completely monotone and is capable of representing phenomena of relaxation. The graph of the function $E_\alpha (-x^\alpha )$, with $x=(t/\tau _\mathrm{R})\ge 0$, is illustrated in Fig. 11 for different values of $\alpha$.

Fig. 11

Mittag–Leffler function for different values of $\alpha$

Appendix 3: Lianis material functions for SBR circular bar at $0\,^{\circ }$C [73]

The Lianis constitutive equation for an isotropic incompressible viscoelastic material is given by

$$\begin{aligned} {\varvec{\sigma }}(t)= & {} \left[ a+\Bigl (I_1({\mathbf {B}})-3\Bigl )b+cI_1({\mathbf {B}})\right] {\mathbf {B}}-c{\mathbf {B}}^2+2\int \limits _{-\infty }^{t}\varphi _0(t-\tau ){\dot{\mathbf {C}}}_t(\tau ){\mathrm {d}}\tau \nonumber \\&+\int \limits _{-\infty }^{t}\varphi _1(t-\tau )[{\mathbf {B}}{\dot{\mathbf {C}}}_t(\tau )+{\dot{\mathbf {C}}}_t(\tau ){\mathbf {B}}]{\mathrm {d}}\tau + \int \limits _{-\infty }^{t}\varphi _2(t-\tau ) [{\mathbf {B}}^2{\dot{\mathbf {C}}}_t(\tau )+ {\dot{\mathbf {C}}}_t(\tau ){\mathbf {B}}^2]{\mathrm {d}}\tau \nonumber \\&+{\mathbf {B}}\int \limits _{-\infty }^{t}\varphi _3(t-\tau ) I_1({\mathbf {B}}{\dot{\mathbf {C}}}_t(\tau )){\mathrm {d}}\tau - P(t)\mathbf {I}, \end{aligned}$$

(77)

where a, b, c are constants corresponding to the equilibrium stress, and the $\varphi _k(u)$ are relaxation functions that approach zero as $u\rightarrow \infty$.

For a SBR at $0\,^{\circ }$C, the constants a, b, c were found to be $a=29.5$ psi, $b=0$, $c=51.07$ psi, where 1 psi = $6.8948 \times 10^3$ Pa, and the Lianis material functions were related by

$$\begin{aligned} \varphi _1(t)=\varphi _3(t)=-2\varphi _2(t). \end{aligned}$$

(78)

The relaxation functions $\varphi _k(t)$ are tabulated in Table 3.

Table 3 Lianis material functions for SBR at $0\,^{\circ }$C

Appendix 4: Tension–torsion of a circular cylinder

In cylindrical coordinates and for the extension–torsion of a solid circular cylinder made of an incompressible isotropic material, we have

$$\begin{aligned} {\mathbf {F}}= & {} \left( \begin{array}{c@{\quad }c@{\quad }c} \frac{1}{\sqrt{\lambda }} &{} 0 &{} 0 \\ 0 &{} \frac{1}{\sqrt{\lambda }} &{} \sqrt{\lambda }\psi R \\ 0 &{} 0 &{} \lambda \end{array} \right) , \end{aligned}$$

(79)

$$\begin{aligned} {\mathbf {B}}= & {} \left( \begin{array}{c@{\quad }c@{\quad }c} \lambda ^{-1} &{} 0 &{} 0 \\ 0 &{} \lambda ^{-1}+\lambda (\psi R)^2 &{} \lambda ^{3/2}\psi R \\ 0 &{} \lambda ^{3/2}\psi R &{} \lambda ^2 \end{array}\right) , \end{aligned}$$

(80)

$$\begin{aligned} {\mathbf {C}}= & {} \left( \begin{array}{c@{\quad }c@{\quad }c} \lambda ^{-1}&{} 0 &{} 0 \\ 0 &{} \lambda ^{-1} &{} \psi R \\ 0 &{} \psi R &{} \lambda ^2+\lambda (\psi R)^2 \end{array}\right) , \end{aligned}$$

(81)

$$\begin{aligned} {\mathbf {e}}= & {} \frac{1}{2} \left( \begin{array}{c@{\quad }c@{\quad }c} \lambda -1&{} 0 &{} 0 \\ 0 &{} \lambda +(\psi R)^2-1 &{} -\lambda ^{-1}\psi R \\ 0 &{} -\lambda ^{-1}\psi R &{} \lambda ^{-2}-1 \end{array}\right) . \end{aligned}$$

(82)

The equation of equilibrium is given by

$$\begin{aligned} \mathrm{div}\,{\varvec{\sigma }}=0. \end{aligned}$$

(83)

Equations (79)–(82), together with the constitutive equation of the proposed model, imply that

$$\begin{aligned} {\varvec{\sigma }}= \left( \begin{array}{c@{\quad }c@{\quad }c} \sigma _{rr} &{} 0 &{} 0 \\ 0 &{} \sigma _{\theta \theta } &{} \sigma _{\theta z} \\ 0 &{} \sigma _{\theta z} &{} \sigma _{zz} \end{array}\right) , \end{aligned}$$

(84)

where the nonzero components of the Cauchy stress tensor are given by

$$\begin{aligned} \sigma _{rr}(t)= & {} 2\tau _\mathrm{R}^{-\alpha }B_{rr}(t)\left[ \int \limits _{0}^{t}\left( c_{10}+c_{01}I_1({\mathbf {C}})\right) K(t-s) {\mathrm {d}}s \right] -2\tau _\mathrm{R}^{-\alpha }c_{01}F_{rr}^2(t)\left[ \int \limits _{0}^{t}C_{rr}(s) K(t-s) {\mathrm {d}}s \right] \nonumber \\&-2\mu _0F_{rr}^2(t)\left[ \int \limits _{0}^{t}\xi ^{-1}(s) {D}^\alpha e_{rr}(s) K(t-s) {\mathrm {d}}s \right] -P(t), \end{aligned}$$

(85)

$$\begin{aligned} \sigma _{\theta \theta }(t)= & {} 2\tau _\mathrm{R}^{-\alpha }\left\{ B_{\theta \theta }(t)\left[ \int \limits _{0}^{t}\left( c_{10}+c_{01}I_1({\mathbf {C}})\right) K(t-s) {\mathrm {d}}s \right] -2\tau _\mathrm{R}^{-\alpha }c_{01}\left\{ F_{\theta \theta }^2(t) \left[ \int \limits _{0}^{t}C_{\theta \theta }(s) K(t-s) {\mathrm {d}}s \right] \right. \right. \nonumber \\&\left. +2F_{\theta \theta }(t)F_{\theta z}(t)\left[ \int \limits _{0}^{t}C_{\theta z}(s) K(t-s) {\mathrm {d}}s \right] +F_{\theta z}^2(t)\left[ \int \limits _{0}^{t}C_{zz}(s) K(t-s) {\mathrm {d}}s \right] \right\} \nonumber \\&-2\mu _0\left\{ F_{\theta \theta }^2(t)\left[ \int \limits _{0}^{t} \xi ^{-1}(s) {D}^\alpha e_{\theta \theta }(s)K(t-s){\mathrm {d}}s \right] +2F_{\theta \theta }(t)F_{\theta z}(t)\left[ \int \limits _{0}^{t} \xi ^{-1}(s) {D}^\alpha e_{\theta z}(s) K(t-s) {\mathrm {d}}s \right] \right. \nonumber \\&\left. +F_{\theta z}^2(t)\left[ \int \limits _{0}^{t}\xi ^{-1}(s) {D}^\alpha e_{zz}(s) K(t-s) {\mathrm {d}}s \right] \right\} -P(t), \end{aligned}$$

(86)

$$\begin{aligned} \sigma _{zz}(t)= & {} 2\tau _\mathrm{R}^{-\alpha }B_{zz}(t)\left[ \int \limits _{0}^{t} \Bigl (c_{10}+c_{01}I_1({\mathbf {C}})\Bigr )K(t-s) {\mathrm {d}}s \right] -2\tau _\mathrm{R}^{-\alpha }c_{01}F_{zz}^2(t)\left[ \int \limits _{0}^{t} C_{zz}(s) K(t-s) {\mathrm {d}}s \right] \nonumber \\&-2\mu _0F_{zz}^2(t)\left[ \int \limits _{0}^{t}\xi ^{-1}(s) {D}^\alpha e_{zz}(s) K(t-s) {\mathrm {d}}s \right] -P(t), \end{aligned}$$

(87)

$$\begin{aligned} \sigma _{\theta z}(t)= & {} 2\tau _\mathrm{R}^{-\alpha }B_{\theta z}(t) \left[ \int \limits _{0}^{t}\Bigl (c_{10}+c_{01}I_1({\mathbf {C}})\Bigr ) K(t-s) {\mathrm {d}}s \right] -2\tau _\mathrm{R}^{-\alpha }c_{01}\left\{ F_{\theta \theta }(t)F_{zz}(t) \left[ \int \limits _{0}^{t}C_{\theta z}(s) K(t-s) {\mathrm {d}}s \right] \right. \nonumber \\&\left. +F_{\theta z}(t)F_{zz}(t)\left[ \int \limits _{0}^{t}C_{zz}(s) K(t-s) {\mathrm {d}}s \right] \right\} -2\mu _0\left\{ F_{\theta \theta }(t)F_{zz}(t)\left[ \int \limits _{0}^{t} \xi ^{-1}(s) {D}^\alpha e_{\theta z}(s) K(t-s) {\mathrm {d}}s \right] \right. \nonumber \\&\left. +F_{\theta z}(t)F_{zz}(t)\left[ \int \limits _{0}^{t}\xi ^{-1}(s) { D}^\alpha e_{zz}(s) K(t-s) {\mathrm {d}}s \right] \right\} , \end{aligned}$$

(88)

where $I_1({\mathbf {C}})=\lambda ^2+2\lambda ^{-1}+\lambda (\psi R)^2$.

Substituting (84) into (83), one obtains

$$\begin{aligned}&\frac{\partial \sigma _{rr}}{\partial r}+\frac{\partial \sigma _{rz}}{\partial z}+\frac{\partial \sigma _{\theta r}}{r\partial \theta }+\frac{\sigma _{rr}-\sigma _{\theta \theta }}{r}= \frac{\partial \sigma _{rr}}{\partial r}+\frac{\sigma _{rr}-\sigma _{\theta \theta }}{r}=0, \end{aligned}$$

(89)

$$\begin{aligned}&\frac{\partial \sigma _{r\theta }}{\partial r}+\frac{\partial \sigma _{z\theta }}{\partial z}+\frac{\partial \sigma _{\theta \theta }}{r\partial \theta }+2\frac{\sigma _{r\theta }}{r}= \frac{\partial P}{\partial \theta }=0, \end{aligned}$$

(90)

$$\begin{aligned}&\frac{\partial \sigma _{rz}}{\partial r}+\frac{\partial \sigma _{z\theta }}{r\partial \theta }+\frac{\partial \sigma _{zz}}{\partial z}+\frac{\sigma _{rz}}{r}=\frac{\partial P}{\partial z}=0. \end{aligned}$$

(91)

We see that Eq. (90) implies that ${\partial } P/{\partial } \theta =0$ and Eq. (91) implies that $\partial P/\partial z=0$. This means that P is independent of $\theta$ and z and is a function of t and r only. Then (89) becomes

$$\begin{aligned} \frac{\partial \sigma _{rr}}{\partial r}+\frac{1}{r}(\sigma _{rr}- \sigma _{\theta \theta })=0. \end{aligned}$$

(92)

The integration of (92) leads to

$$\begin{aligned} \int \limits _{r}^{r_0}{\mathrm {d}}\sigma _{rr}=\sigma _{rr}(r_0,t)- \sigma _{rr}(r,t)=-\int \limits _{r}^{r_0}\frac{1}{r}(\sigma _{rr}- \sigma _{\theta \theta }){\mathrm {d}}r. \end{aligned}$$

(93)

We have

$$\begin{aligned} \sigma _{rr}(r=r_0,t)=0. \end{aligned}$$

(94)

One can deduce from Eq. (93) that

$$\begin{aligned} \sigma _{rr}(r,t)=\int \limits _{r}^{r_0}\frac{1}{r}(\sigma _{rr}- \sigma _{\theta \theta }){\mathrm {d}}r. \end{aligned}$$

(95)

According to Eqs. (85) and (86), one can write the following relations:

$$\begin{aligned} \sigma _{rr}(r,t)= & {} T_1(r,t)-P(r,t), \end{aligned}$$

(96)

$$\begin{aligned} \sigma _{\theta \theta }(r,t)= & {} T_2(r,t)-P(r,t), \end{aligned}$$

(97)

where

$$\begin{aligned} T_1(r,t)= & {} 2\tau _\mathrm{R}^{-\alpha }B_{rr}(t)\left[ \int \limits _{0}^{t} \Bigl (c_{10}+c_{01}I_1({\mathbf {C}})\Bigr )K(t-s) {\mathrm {d}}s\right] -2\tau _\mathrm{R}^{-\alpha }c_{01}F_{rr}^2(t)\left[ \int \limits _{0}^{t}C_{rr} (s) K(t-s) {\mathrm {d}}s \right] \nonumber \\&-2\mu _0F_{rr}^2(t)\left[ \int \limits _{0}^{t}\xi ^{-1}(s) {D}^\alpha e_{rr}(s) K(t-s) {\mathrm {d}}s \right] , \end{aligned}$$

(98)

$$\begin{aligned} T_2(r,t)= & {} 2\tau _\mathrm{R}^{-\alpha }B_{\theta \theta }(t)\left[ \int \limits _{0}^{t}\Bigl (c_{10}+ c_{01}I_1({\mathbf {C}})\Bigr )K(t-s) {\mathrm {d}}s\right] -2\tau _\mathrm{R}^{-\alpha }c_{01}\left\{ F_{\theta \theta }^2(t) \left[ \int \limits _{0}^{t}C_{\theta \theta }(s) K(t-s) {\mathrm {d}}s \right] \right. \nonumber \\&\left. +2F_{\theta \theta }(t)F_{\theta z}(t)\left[ \int \limits _{0}^{t}C_{\theta z}(s) K(t-s) {\mathrm {d}}s \right] +F_{\theta z}^2(t)\left[ \int \limits _{0}^{t}C_{zz}(s) K(t-s) { \mathrm {d}}s \right] \right\} \nonumber \\&-2\mu _0\left\{ F_{\theta \theta }^2(t)\left[ \int \limits _{0}^{t} \xi ^{-1}(s) {D}^\alpha e_{\theta \theta }(s)K(t-s){\mathrm {d}}s \right] +2F_{\theta \theta }(t)F_{\theta z}(t)\left[ \int \limits _{0}^{t}\xi ^{-1}(s) {D}^ \alpha e_{\theta z}(s)K(t-s){\mathrm {d}}s \right] \right. \nonumber \\&\left. +F_{\theta z}^2(t)\left[ \int \limits _{0}^{t}\xi ^{-1}(s) {D}^\alpha e_{zz}(s) K(t-s) {\mathrm {d}}s \right] \right\} . \end{aligned}$$

(99)

Thus, Eq. (95) gives

$$\begin{aligned} \sigma _{rr}(r,t)=\int \limits _{r}^{r_0}\frac{1}{r}(T_1-T_2){\mathrm {d}}r. \end{aligned}$$

(100)

Substituting Eqs. (100) into (96), we can easily find the Lagrange multiplier, i.e., the pressure P(r, t), in terms of the material parameters $\lambda (t)$ and $\psi (t)$:

$$\begin{aligned} P(r,t)=T_1(r,t)-\int \limits _{r}^{r_0}\frac{1}{r} [T_1(r,t)-T_2(r,t)]{\mathrm {d}}r. \end{aligned}$$

(101)

After calculations we obtain

$$\begin{aligned} P(r,t)= & {} 2c_{10}\lambda ^{-1}(t)(1-g(t)) +2c_{01}\tau _\mathrm{R}^{-\alpha }\lambda ^{-1}(t)\left[ \int \limits _{0}^{t}\Bigl (\lambda ^2(s)+\lambda ^{-1}(s)+[r\lambda (s)\psi (s)]^2\Bigr )K(t-s){\mathrm {d}}s \right] \nonumber \\&-\,2\mu _0\lambda ^{-1}(t)\left[ \int \limits _{0}^{t}\xi ^{-1}(s) {D}^\alpha e_{rr}(s)K(t-s){\mathrm {d}}s \right] +\int \limits _{r}^{r_0}\Biggl \{2c_{10}[\lambda (t)\psi (t)]^2 \Bigl (1-g(t)\Bigr )r\nonumber \\&+\,4c_{01}\tau _\mathrm{R}^{-\alpha }[\lambda (t)\psi (t)]^2r\int \limits _{0}^{t} \lambda ^{-1}(s)K(t-s){\mathrm {d}}s -4c_{01}\lambda ^{1/2}(t)\psi (t)\tau _\mathrm{R}^{-\alpha }\int \limits _{0}^{t}[\lambda ^{1/2}(s)\psi (s)r]K(t-s){\mathrm {d}}s \Biggr \}{\mathrm {d}}r\nonumber \\&+\,\int \limits _{r}^{r_0}\left\{ 2\mu _0r^{-1}\lambda ^{-1}(t) \left[ \int \limits _{0}^{t}\xi ^{-1}(s) {D}^\alpha e_{rr}(s)K(t-s){\mathrm {d}}s -\int \limits _{0}^{t}\xi ^{-1}(s) {D}^\alpha e_{\theta \theta }(s)K(t-s){\mathrm {d}}s\right] \right. \nonumber \\&-\,4\mu _0\lambda ^{1/2}(t)\psi (t)\left[ \int \limits _{0}^{t}\xi ^{-1}(s) {D}^\alpha e_{\theta z}(s)K(t-s){\mathrm {d}}s\right] \nonumber \\&\left. -\,2\mu _0r[\lambda (t)\psi (t)]^2\left[ \int \limits _{0}^{t}\xi ^{-1}(s) {D}^ \alpha e_{zz}(s) K(t-s) {\mathrm {d}}s\right] \right\} {\mathrm {d}}r. \end{aligned}$$

(102)