Models of fractional viscous stresses for incompressible materials

Abstract

We present and review several models of fractional viscous stresses from the literature, which generalise classical viscosity theories to fractional orders by replacing total strain derivatives in time with fractional time derivatives. We also briefly introduce Prony-type approximations of these theories. Here, we investigate the issues of material frame-indifference and thermodynamic consistency for these models and find that on these bases, some are physically unacceptable. Next, we study elementary shearing and tensile motions, observing that some models are more convenient to use than others for the analysis of creep and relaxation. Finally, we compute the incremental stresses due to small-amplitude wave propagation in a deformed material, with a view to establish acoustoelastic formulas for prospective experimental calibrations.

1. Introduction

Differential operators based on the Riemann–Liouville integral are commonly used to generalise differentiation of integer order to fractional orders, thus laying the foundations of fractional calculus [1]. This branch has found various applications over time [2, 3], most notably in electronics and in material rheology. In the latter case [4, 5], fractional calculus provides accurate predictions of the time-dependent mechanical behaviour with a limited number of parameters. Figure 1 lists a few materials whose stress relaxation response can be described by such theories, e.g., xanthan gum, bread dough, and nylon. One property of these materials is the power-law time evolution of the stress in response to a sudden deformation (see the review by Bonfanti et al. [6] for details). A possible explanation of the micromechanical origins of fractional behaviour is provided by Brenner [7].

Figure 1.

Stress relaxation data sets show that the response of many materials follows a power-law behaviour. Source: Reproduced from Bonfanti et al. [6] under the Creative Commons Attribution 3.0 Unported Licence.

Classical rheological models are described by an assembly of one-dimensional elements such as springs and dashpots. Within this framework, fractional viscous elements are represented by springpots which are intermediates between springs (accounting for elasticity) and dashpots (accounting for viscous dissipation), see Mainardi [4] for an overview of the linear theory valid at small strains. The simplest such model is based on a single springpot element, which can be found under the name of the “Kjartansson constant-Q model” in geophysics [8].

The mechanical analogue of Newtonian or Kelvin–Voigt viscoelastic materials is a dashpot connected to a spring in parallel. In this case, elastic stresses are proportional to the scalar strain $\epsilon$ , whereas viscous stresses are proportional to the velocity gradient or strain rate $\stackrel{·}{\epsilon }$ (in what follows, $\stackrel{·}{\epsilon }$ denotes the material time derivative of $\epsilon$ ). The strain rate can be generalised to fractional differential orders $\alpha \in ]0,1[$ based on the Caputo time-derivative ${\mathfrak{D}}^{\alpha }$ defined by:

$${\mathfrak{D}}^{\alpha }\epsilon (t)=\frac{1}{\mathrm{\Gamma }(1-\alpha )}{\int }_0^t\frac{\stackrel{·}{\epsilon }(s)}{{(t-s)}^{\alpha }}\mathrm{d}s,$$

(1)

where $\epsilon$ is causal and $\mathrm{\Gamma }(z)={\int }_0^{\infty }{t}^{z-1}{\mathrm{e}}^{-t}\mathrm{d}t$ defines Euler’s Gamma function. This step amounts to replacing the dashpot with a springpot in the diagram (see Figure 2).

Figure 2.

Schematic representation of the fractional Kelvin–Voigt model.

With the above definition of the fractional derivative, we have the limits ${\mathfrak{D}}^0\epsilon =\epsilon$ and, by differentiation, ${\mathfrak{D}}^1\epsilon =\stackrel{·}{\epsilon }$ . Thus, $\alpha =0$ corresponds to an elastic model, whereas $\alpha =1$ recovers the classical Kelvin–Voigt theory. Generalisation of equation (1) to arbitrary orders $\alpha \ge 0$ can be carried out [9], as well as generalisation to non-causal functions (e.g., periodic ones [10]). In this latter case, the Caputo derivative may be replaced by the Weyl derivative, which amounts to evaluating the integral (1) from $-\infty$ to $t$ .

Fractional calculus has also been used to model the time-dependent mechanical response of highly deformable soft materials for which the small strain assumption is no longer valid. Examples of such materials include amorphous polymers, organic glasses [11], rubber, elastomers [12, 13], tendons [14], liver tissue [15], and other similar materials. While the generalisation of the spring to three-dimensional finite strain is rather straightforward (see the textbook by Holzapfel [16] for a description of hyperelasticity), the crux of the matter is a proper definition of the fractional viscous stress. In this study, we consider the finite motion of incompressible viscous materials. More specifically, we describe nonlinear three-dimensional generalisations of the fractional Kelvin–Voigt rheology depicted in Figure 2.

To describe their motion, we introduce $\mathit{F}(t)$ , the deformation gradient tensor at time t, which is defined as the gradient $\partial \mathit{x}/\partial \mathit{X}$ of the current position $\mathit{x}=\mathit{x}(t)$ of a particle with respect to its position $\mathit{X}=\mathit{x}\left(0\right)$ in the reference configuration. Incompressible materials do not allow for volume change. Hence, isochoricity is enforced at all times $(det\mathit{F}\equiv 1)$ , and the mass density $\rho$ is constant.

In a similar fashion to the linear scalar case, we assume that the second Piola–Kirchhoff stress tensor $\mathit{S}$ and the Cauchy stress tensor $\mathit{\sigma }=\mathit{FS}{\mathit{F}}^{\mathrm{T}}$ can be split additively into elastic and viscous parts, i.e.,

$$\mathit{S}=-p{\mathit{C}}^{-1}+{\mathit{S}}^{\mathrm{e}}+{\mathit{S}}^{\mathrm{v}},\mathit{\sigma }=-p\mathit{I}+{\mathit{\sigma }}^{\mathrm{e}}+{\mathit{\sigma }}^{\mathrm{v}},$$

(2)

where $\mathit{C}={\mathit{F}}^{\mathrm{T}}\mathit{F}$ is the right Cauchy–Green deformation tensor, $\mathit{I}$ is the identity tensor, and the scalar p is a Lagrange multiplier accounting for the incompressibility constraint. The stresses with exponents e and v are elastic and viscous contributions, respectively. Elastic contributions vanish in the fluid limit, which can therefore be viewed as a special case of the present theories.

Up to a suitable redefinition of the pressure coefficient p, we note that the partial stress tensors ${\mathit{S}}^{•}$ , ${\mathit{\sigma }}^{•}$ in equation (2) can be replaced with their deviatoric counterparts ${\mathit{S}}_{\mathrm{D}}^{•}$ , ${\mathit{\sigma }}_{\mathrm{d}}^{•}$ . Here, the deviators are denoted by the subscripts “D” and “d,” such that:

$$(•{)}_{\mathrm{D}}=(•)-\frac{1}{3}\mathrm{tr}[(•)\mathit{C}]{\mathit{C}}^{-1},(•{)}_{\mathrm{d}}=(•)-\frac{1}{3}\mathrm{tr}(•)\mathit{I},$$

(3)

and tr is the trace. With these definitions, we note that the tensor $(•{)}_{\mathrm{d}}$ is trace-free.

In classical Newtonian viscosity theories, the viscous Cauchy stress is typically expressed as ${\mathit{\sigma }}^{\mathrm{v}}=2\eta \mathit{D}$ , where $\eta >0$ is the dynamic viscosity. The tensor $\mathit{D}=\frac{1}{2}(\mathit{L}+{\mathit{L}}^{\mathrm{T}})$ with $\mathit{L}=\stackrel{·}{\mathit{F}}{\mathit{F}}^{-1}$ is the Eulerian strain rate tensor, which is trace-free by virtue of incompressibility—in other words, ${\mathit{\sigma }}^{\mathrm{v}}={\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{v}}$ . In terms of the second Piola–Kirchhoff stress, the Newtonian viscosity is then expressed by the relationship:

$${\mathit{S}}^{\mathrm{v}}=2\eta \stackrel{·}{\mathbf{\Pi }}=2\eta {\mathit{C}}^{-1}\stackrel{·}{\mathit{E}}{\mathit{C}}^{-1},$$

(4)

where

$$\mathbf{\Pi }=\frac{1}{2}(\mathit{I}-{\mathit{C}}^{-1}),\mathit{E}=\frac{1}{2}(\mathit{C}-\mathit{I}),$$

(5)

denote a Piola-type deformation tensor and the Green–Lagrange strain tensor $\mathit{E}=\mathit{C}\mathbf{\Pi }$ , respectively. The equality $\stackrel{·}{\mathbf{\Pi }}={\mathit{C}}^{-1}\stackrel{·}{\mathit{E}}{\mathit{C}}^{-1}$ used in equation (4) follows from differential rules.

As an alternative to Newtonian viscosity (4), one might define the viscous stress as:

$${\mathit{S}}^{\mathrm{v}}=2\eta \stackrel{·}{\mathit{E}},$$

(6)

which is proportional to the rate of Green–Lagrange strain. In Physical Acoustics, this viscous term is sometimes preferred to the Newtonian viscous term, although care must be taken that it is added to the second, and not the first, Piola–Kirchhoff stress tensor [17].

In terms of the Cauchy stress tensor, the viscous stress (6) reads ${\mathit{\sigma }}^{\mathrm{v}}=2\eta \mathit{BDB}$ , where we have used the identity $\stackrel{·}{\mathit{E}}={\mathit{F}}^{\mathrm{T}}\mathit{D}\mathit{F}$ and the definition $\mathit{B}=\mathit{F}{\mathit{F}}^{\mathrm{T}}$ of the left Cauchy–Green deformation tensor. Here, we note that the viscous Cauchy stress is not necessarily trace-free (the Cauchy–Schwarz inequality does not apply). Nevertheless, the latter can still be replaced by its deviatoric part ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{v}}$ up to a suitable redefinition of the arbitrary Lagrange multiplier in equation (2).

Furthermore, it is worth pointing out that the viscosity theories (4)–(6) have fundamentally different mathematical properties. In fact, contrary to the Newtonian viscosity case (4), viscoelastic shearing motions might be limited to a finite time when using equation (6), which is potentially problematic for the purpose of experimental characterisation involving long-time relaxation processes [18].

In the next section, we present straightforward generalisations of the viscous stresses (4)–(6) to fractional orders (1) and establish connections with the literature (section 2). We briefly discuss approximations of the fractional derivative to be used in computational applications (section 3). Then, we discuss the physical properties of the models at hand, more specifically in relation to the objectivity requirement and thermodynamic consistency (section 4). It appears that some theories from the literature are physically unacceptable in those respects. Furthermore, we study elementary shearing and tensile motions (section 5). Finally, we compute incremental stresses via a “small-on-large” linearisation and obtain acoustoelastic formulas (section 6). The results of this study might be used in experimental setups or in other applications.

2. Constitutive models

In this section, we introduce straightforward generalisations of the viscous stresses (4)–(6) to fractional differential orders (1). A summary of these constitutive laws is given in Table 1.

Table 1.

Expression of the viscous Cauchy stress ${\mathit{\sigma }}^{\mathrm{v}}=\mathit{F}{\mathit{S}}^{\mathrm{v}}{\mathit{F}}^{\mathrm{T}}$ for the Models A, B, and C, including the elastic limit $\alpha \to 0$ and the viscous limit $\alpha \to 1$ .

Model		Fractional viscous stress ${\mathit{\sigma }}^{\mathrm{v}}$	Elastic limit	Viscous limit
A	(8)	$2{\eta }_{\alpha }\mathit{F}\left({\mathfrak{D}}^{\alpha }\mathbf{\Pi }\right){\mathit{F}}^{\mathrm{T}}$	$\mu (\mathit{B}-\mathit{I})$	$2\eta \mathit{D}$
B	(11)	$2{\eta }_{\alpha }{\mathit{F}}^{-\mathrm{T}}\left({\mathfrak{D}}^{\alpha }\mathit{E}\right){\mathit{F}}^{-1}$	$\mu (\mathit{I}-{\mathit{B}}^{-1})$	$2\eta \mathit{D}$
C	(12)	$2{\eta }_{\alpha }\mathit{F}\left({\mathfrak{D}}^{\alpha }\mathit{E}\right){\mathit{F}}^{\mathrm{T}}$	$\mu ({\mathit{B}}^2-\mathit{B})$	$2\eta \mathit{BDB}$

The fractional time derivative ${\mathfrak{D}}^{\alpha }$ is defined in equation (1) and the fractional viscosity ${\eta }_{\alpha }$ is expressed in $\mathrm{Pa}.{\mathrm{s}}^{\alpha }$ .

2.1. Model A

We introduce a fractional time derivative (1) of the Piola strain $\mathbf{\Pi }$ in the definition of the Newtonian viscous Piola–Kirchhoff stress (4), as follows: ${\mathit{S}}^{\mathrm{v}}=2{\eta }_{\alpha }{\mathfrak{D}}^{\alpha }\mathbf{\Pi }$ , where ${\eta }_{\alpha }>0$ is a fractional dynamic viscosity (in $\mathrm{Pa}.{\mathrm{s}}^{\alpha }$ ). This expression involves only two independent parameters—the coefficient ${\eta }_{\alpha }$ and the differential order $\alpha$ . For convenience, we introduce the redundant parameter $\mu ={\eta }_{\alpha }/{\tau }^{\alpha }$ , where $\tau >0$ is a characteristic time. In other words, we have ${\mathit{S}}^{\mathrm{v}}=2\mu {\tau }^{\alpha }{\mathfrak{D}}^{\alpha }\mathbf{\Pi }$ . Without loss of generality, the parameter $\mu$ can be chosen in such a way that it equals the initial shear modulus in solid materials (in Pa).

According to the definition of the fractional derivative (1), this expression might be rewritten as a time-domain convolution product with kernel $\kappa$ :

$${\mathit{S}}^{\mathrm{v}}=2\mu \kappa *\stackrel{·}{\mathbf{\Pi }},\mathrm{where}\kappa (t)=\frac{{(t/\tau )}^{-\alpha }}{\mathrm{\Gamma }(1-\alpha )}\mathrm{H}(t),$$

(7)

and H represents the Heaviside step function. In agreement with previous definitions, we thus have the relationship $\kappa *\stackrel{·}{\mathbf{\Pi }}={\tau }^{\alpha }{\mathfrak{D}}^{\alpha }\mathbf{\Pi }$ . In terms of the Cauchy stress tensor, we find:

$${\mathit{\sigma }}^{\mathrm{v}}=2\mu \mathit{F}(\kappa *\stackrel{·}{\mathbf{\Pi }}){\mathit{F}}^{\mathrm{T}}=2\mu \kappa *\left({\mathit{F}}_{t|s}\mathit{D}{\mathit{F}}_{t|s}^{\mathrm{T}}\right),$$

(8)

due to the identity $\stackrel{·}{\mathbf{\Pi }}={\mathit{F}}^{-1}\mathit{D}{\mathit{F}}^{-\mathrm{T}}$ . Here, the tensor

$${\mathit{F}}_{t|s}=\mathit{F}(t){\mathit{F}}^{-1}(s),$$

(9)

is the relative deformation gradient from the configuration at the integration time $s\in [0,t]$ to the configuration at the current time t. We emphasise that other expressions for the convolution kernel are possible (see the review by Freed and Diethelm [19] for alternatives).

With the above expressions of the viscous stress, we recover the Newtonian viscous stress ${\mathit{S}}^{\mathrm{v}}\to 2\eta \stackrel{·}{\mathbf{\Pi }}$ , i.e., ${\mathit{\sigma }}^{\mathrm{v}}\to 2\eta \mathit{D}$ , in the limit of integer differentiation $\alpha \to 1$ , where $\eta =\mu \tau$ . In the limit of no differentiation $\alpha \to 0$ , we obtain the extra elastic contribution ${\mathit{S}}^{\mathrm{v}}\to 2\mu \mathbf{\Pi }$ , i.e., ${\mathit{\sigma }}^{\mathrm{v}}\to \mu (\mathit{B}-\mathit{I})$ , which is of neo-Hookean type—this property can be inferred from a redefinition of the Lagrange multiplier in equation (2)b, see also the expression (23) of the Mooney–Rivlin stress with the second Mooney parameter ${\mathfrak{C}}_{-}$ equal to zero. Therefore, this model manages to cover and connect an elastic theory and a viscous theory.

The present theory is strongly related to other models found in the literature. Indeed, equation (8) matches equation (2.43) of Shen [20], where it is linked to the theory by Drozdov [12] (cf. next paragraph). It is also of the general form proposed in Capilnasiu et al. [15], equation (5) therein, although it does not match later propositions from that study. Furthermore, the above expression is found in Palade et al. [11] as a special case of equation (16) therein. It is also aligned with “Model A” of Haupt and Lion [13], see equation (5.12) therein. In particular, if the kernel $\kappa$ is chosen exponential instead of the power-law expression (7), then we recover the “upper-convected” Maxwell model which involves the Oldroyd rate of Cauchy stress (also equivalent to the Truesdell rate in the incompressible case).

Equations (27)–(32) of Drozdov [12] introduce a viscous stress based on a relative deformation gradient tensor and the strain-rate tensor $\mathit{D}$ . However, the conventions therein differ from the present ones. To reconcile the two, we take the transpose of the deformation gradients in Drozdov [12], see also equation (4) and later sections therein. Thus, equation (3) of Drozdov [12] becomes ${\mathit{F}}_{t|s}^{\mathrm{T}}={\mathit{F}}^{-\mathrm{T}}(s){\mathit{F}}^{\mathrm{T}}(t)$ if our notations (9) are used, and the viscous stress (32) proposed therein takes the form of equation (8)b.

2.2. Model B

We introduce a fractional time derivative (1) of the Green–Lagrange strain $\mathit{E}$ in the definition of the viscous stress (4) as follows: ${\mathit{S}}^{\mathrm{v}}=2{\eta }_{\alpha }{\mathit{C}}^{-1}\left({\mathfrak{D}}^{\alpha }\mathit{E}\right){\mathit{C}}^{-1}$ , i.e.,

$${\mathit{S}}^{\mathrm{v}}=2\mu {\mathit{C}}^{-1}(\kappa *\stackrel{·}{\mathit{E}}){\mathit{C}}^{-1},$$

(10)

$${\mathit{\sigma }}^{\mathrm{v}}=2\mu {\mathit{F}}^{-\mathrm{T}}(\kappa *\stackrel{·}{\mathit{E}}){\mathit{F}}^{-1}=2\mu \kappa *\left({\mathit{F}}_{t|s}^{-\mathrm{T}}\mathit{D}{\mathit{F}}_{t|s}^{-1}\right),$$

(11)

where we have used the same notations as for Model A, as well as the identity $\stackrel{·}{\mathit{E}}={\mathit{F}}^{\mathrm{T}}\mathit{D}\mathit{F}$ . Thus, the Cauchy stress has a similar expression as in equation (8), up to the fact that the relative deformation gradient tensor ${\mathit{F}}_{t|s}$ has been replaced with its inverse transpose.

With the above expression, we recover the same viscous limit as for Model A when $\alpha \to 1$ . In the limit of no differentiation $\alpha \to 0$ , equation (10) yields the elastic stress ${\mathit{S}}^{\mathrm{v}}\to 2\mu \mathbf{\Pi }{\mathit{C}}^{-1}$ , i.e., ${\mathit{\sigma }}^{\mathrm{v}}\to \mu (\mathit{I}-{\mathit{B}}^{-1})$ . Thus, up to a redefinition of the Lagrange multiplier in equation (2)b, we observe that the elastic limit of this theory corresponds to Mooney–Rivlin elasticity (23) with the first Mooney parameter ${\mathfrak{C}}_{+}$ equal to zero. The viscous stress (10) corresponds to Model B of the study by Haupt and Lion [13], see equation (5.22) therein. In particular, if the kernel $\kappa$ is chosen exponential instead of the power-law expression (7), then we recover the “lower-convected” Maxwell model involving the Cotter–Rivlin rate of Cauchy stress.

2.3. Model C

We introduce a fractional time derivative (1) of the Green–Lagrange strain $\mathit{E}$ in the definition of the second form of viscous stress (6) as follows: ${\mathit{S}}^{\mathrm{v}}=2{\eta }_{\alpha }{\mathfrak{D}}^{\alpha }\mathit{E}$ , i.e.,

$${\mathit{S}}^{\mathrm{v}}=2\mu \kappa *\stackrel{·}{\mathit{E}},{\mathit{\sigma }}^{\mathrm{v}}=2\mu \mathit{F}(\kappa *\stackrel{·}{\mathit{E}}){\mathit{F}}^{\mathrm{T}},$$

(12)

where we have used the same notations as for Model A. The only difference with respect to Model B is the multiplication of the second Piola–Kirchhoff stress tensor by ${\mathit{C}}^{-1}$ on the left and on the right in the latter case.

With the above expression, we recover ${\mathit{S}}^{\mathrm{v}}\to 2\eta \stackrel{·}{\mathit{E}}$ , i.e., ${\mathit{\sigma }}^{\mathrm{v}}\to 2\eta \mathit{BDB}$ , in the limit of integer differentiation $\alpha \to 1$ . In the limit of no differentiation $\alpha \to 0$ , equation (12) yields the elastic stress ${\mathit{S}}^{\mathrm{v}}\to 2\mu \mathit{E}$ , i.e., ${\mathit{\sigma }}^{\mathrm{v}}\to \mu ({\mathit{B}}^2-\mathit{B})$ . We note that the proposed viscous stress is of the general form found in Capilnasiu et al. [15]. Using the identity $\stackrel{·}{\mathit{E}}={\mathit{F}}^{\mathrm{T}}\mathit{D}\mathit{F}$ , we observe that the viscous stress (12) is included in equation (18) of Palade et al. [11] as a special case. In particular, if the kernel $\kappa$ is chosen exponential instead of the power-law expression (7), then we recover a Maxwell-type model involving the material rate of second Piola–Kirchhoff stress. A scalar compressible model of this form was also proposed by Sugimoto [21], equation (3.7) therein.

2.4. Other models

To facilitate comparisons, the viscous stresses for the Models A, B, and C are summarised in Table 1. Alternative modelling approaches can lead to different definitions of the fractional viscous stress. In the following, we list other propositions found in the literature:

• Freed and Diethelm [19] propose to exploit the K-BKZ hypothesis [22,23] (after Kaye, Bernstein, Kearsley, Zapas), which consists of a viscoelastic extension of elastic constitutive models based initially on exponential relaxation kernels. Variations of this theory are provided by Coleman and Noll [24] (see equation (5.20) therein). Based on more general expressions of the stored viscous energy than for Model A, see Rao and Rajagopal [25], these models can lead to rather complex expressions of the fractional viscous stress.

• Another theory by Freed and Diethelm [9] includes formally similar viscous stresses to equation (8)b up to the substitution of the relative deformation gradient ${\mathit{F}}_{t|s}$ through the relative rotation ${\mathit{R}}_{t|s}=\mathit{R}(t){\mathit{R}}^{\mathrm{T}}(s)$ , where $\mathit{R}$ is the proper orthogonal tensor in the polar decomposition of $\mathit{F}$ (see, e.g., Holzapfel [16]).

• A given elastic stress ${\mathit{S}}^{\mathrm{r}}$ can be used to express the fractional viscous stress as ${\mathit{S}}^{\mathrm{v}}=\kappa *{\stackrel{·}{\mathit{S}}}_{\mathrm{D}}^{\mathrm{r}}$ , which involves the material time derivative of the deviator ${\mathit{S}}_{\mathrm{D}}^{\mathrm{r}}$ , see [15, 26]. This way, we have the limit ${\mathit{S}}^{\mathrm{v}}\to {\mathit{S}}_{\mathrm{D}}^{\mathrm{r}}$ as $\alpha \to 0$ . Here, the elastic stress ${\mathit{S}}^{\mathrm{r}}$ is not necessarily identical to the elastic response ${\mathit{S}}^{\mathrm{e}}$ . For instance, one might set ${\mathit{S}}^{\mathrm{r}}={\mu }^{\mathrm{r}}\mathit{I}$ if this elastic stress is chosen neo-Hookean, where ${\mu }^{\mathrm{r}}$ is a given shear modulus. The present approach is further investigated in a recent preprint [27].

• Based on the same definition of the relative deformation gradient tensor ${\mathit{F}}_{t|s}$ as in equation (9), Zhao et al. [28] switch the position of the transposition symbols in equation (8)b to mimic Drozdov [12] (see equations (2)–(10) of Zhao et al. [28] and also equations (6)–(8) of Gao et al. [29]).

• Delory et al. [30] introduce pseudo-Newtonian stresses based on a fractional rate of deformation gradient $\stackrel{·}{\mathit{F}}$ , invoking Zhang et al. [31] where fractional viscosity is incorporated in a linear fashion.

3. Approximation of the fractional derivative

For practical use, we represent the convolution kernel of equation (7) by a continuous sum of exponentials:

$$\kappa (t)=\mathrm{H}(t){\int }_0^{\infty }\varphi \left(\zeta \right){\mathrm{e}}^{-t/\zeta }\mathrm{d}\zeta ,\varphi \left(\zeta \right)=\frac{{\tau }^{\alpha }\sin \left(\alpha \pi \right)}{\pi {\zeta }^{\alpha +1}},$$

(13)

with a suitable expression of the spectrum $\varphi$ (see Lion [32] and Euler’s reflection formula). This particular form of $\kappa$ is a diffusive representation of the fractional derivative (see section 7.4.1 of Matignon [1] where the same expression is proposed up to a change of variable in the integral).

In practice, the continuous spectrum of relaxation (13) might be approximated by a discrete one, which leads to a Prony-type theory with exponential relaxation [26]. A straightforward discretisation of this kind is obtained based on Gauss–Laguerre quadrature, which consists in evaluating the integrand of equation (13) at the roots of the Laguerre polynomial of degree N with suitable weights. Faster convergence can be obtained based on a change of variable followed by Gauss–Jacobi quadrature (see Diethelm [33] as well as Birk and Song [34] for more elaborated techniques).

For the purpose of illustration, we consider a straightforward Prony series approximation here, and we discuss its suitability. Thus, we write:

$$\kappa (t)=\frac{\sin \left(\pi \alpha \right)}{\pi \alpha }\mathrm{H}(t){\int }_0^1\frac{{\mathrm{e}}^{-t/\zeta }}{\theta }\mathrm{d}\theta \simeq {\kappa }_N(t),$$

(14)

where

$${\kappa }_N(t)=\mathrm{H}(t)\sum _{n=1}^N{w}_n{\mathrm{e}}^{-t/{\zeta }_n},w_n=\frac{\sin \left(\pi \alpha \right)}{\pi \alpha N{\theta }_n}.$$

(15)

Here, we have used the change of variables $\zeta =\tau {(-\ln \theta )}^{-1/\alpha }$ in the integral (13) over the time coordinate $\zeta$ . The resulting integral over $\theta =\exp (-{(\zeta /\tau )}^{-\alpha })$ was then approximated as a discrete sum based on the extended midpoint rule with N points ${\theta }_n=\frac{2n-1}{2N}$ for $1\le n\le N$ , which correspond to the relaxation times ${\zeta }_n=\zeta {|}_{\theta ={\theta }_n}$ .

Such an approximation of the kernel $\kappa$ as a finite sum of exponentials is rather accurate over a broad range of times and frequencies provided a sufficient number N of relaxation mechanisms is included (see the example in Figure 3 where we have used $N=3$ and $N=6$ ). Therein, we display also the modulus of the Fourier transform $\hat{\kappa }(\omega )=\int \kappa (t){\mathrm{e}}^{-\mathrm{i}\omega t}\mathrm{d}t$ of $\kappa$ and ${\kappa }_N$ , which satisfies:

$$\hat{\kappa }(\omega )/\tau ={\left(\mathrm{i}\omega \tau \right)}^{\alpha -1},{\hat{\kappa }}_N(\omega )=\sum _{n=1}^N\frac{w_n{\zeta }_n}{1+\mathrm{i}\omega {\zeta }_n}.$$

(16)

Figure 3.

Exact convolution kernel (full line) and its diffusive Prony series approximation for $\alpha =0.3$ with $N=3$ (dotted line) and $N=6$ terms (dashed lines). (a) Impulse response. (b) Fourier spectrum.

In particular, the figure illustrates the error introduced by the bounded Prony series approximation (15) at short times where $\kappa$ is singular. Nevertheless, the figure shows that the Prony series approximation obtained for $N=6$ relaxation mechanisms is much more accurate than that obtained for $N=3$ relaxation mechanisms at long times, despite the singularity of $\hat{\kappa }$ at low frequency.

Based on a representation of the form (13)–(14), the following identity holds for any causal tensor field $\mathit{T}$ :

$$\kappa *\stackrel{·}{\mathit{T}}={\int }_0^{\infty }\varphi \left(\zeta \right){\mathit{T}}_{\zeta }^{\mathrm{v}}\mathrm{d}\zeta \simeq \sum _{n=1}^N{w}_n{\mathit{T}}_{{\zeta }_n}^{\mathrm{v}}.$$

(17)

The tensors:

$${\mathit{T}}_{\zeta }^{\mathrm{v}}={\int }_0^t{\mathrm{e}}^{-(t-s)/\zeta }\stackrel{·}{\mathit{T}}(s)\mathrm{d}s,$$

(18)

are memory variables governed by a differential equation of the form:

$${\stackrel{·}{\mathit{T}}}_{\zeta }^{\mathrm{v}}=\stackrel{·}{\mathit{T}}-{\mathit{T}}_{\zeta }^{\mathrm{v}}/\zeta .$$

(19)

This result can be obtained by application of the Leibniz integral rule to equation (18)b. Under this form, evaluation of the fractional derivative amounts to solving a linear system of differential equations.

Such approximations of the fractional derivative are therefore convenient from a computational point of view as they avoid the storage of the history of $\mathit{T}$ to evaluate the current viscous stresses. Nevertheless, it is worth pointing out that special care should be taken as the system (19) might become stiff. In fact, if we consider the approximation (15) at large N, the smallest dimensionless relaxation time ${\zeta }_1/\tau =\ln {\left(2N\right)}^{-1/\alpha }$ approaches zero, whereas the largest relaxation time ${\zeta }_N/\tau =\ln {\left(\frac{2N}{2N-1}\right)}^{-1/\alpha }$ approaches infinity. Further improvements of the approximation (15) can be obtained through optimisation of the $2N$ coefficients ${\zeta }_n$ , $w_n$ (see, for instance, Blanc et al. [35]).

4. Properties

4.1. Material frame-indifference

Here, we examine the acceptability of the above-mentioned models in terms of the material frame-indifference principle, which stipulates that the mechanical response of a material should not be affected by a change of observer (Holzapfel [16], sections 5.2–5.4). Doing so, we show that some of the presented models do not comply with the material frame-indifference principle.

We consider the superimposed rigid-body motion defined by ${\mathit{x}}^{+}=\mathit{c}+\mathit{Qx}$ , where $\mathit{c}(t)$ is a vector and $\mathit{Q}(t)$ is a proper orthogonal tensor. Then, the deformation gradient tensor transforms according to ${\mathit{F}}^{+}=\mathit{QF}$ . The kinematic variables $\mathbf{\Pi }$ , $\mathit{E}$ and their material time-derivatives are unaffected by the rigid-body motion: $\mathbf{\Pi }={\mathbf{\Pi }}^{+}$ , $\mathit{E}={\mathit{E}}^{+}$ , as are second Piola–Kirchhoff stresses: ${\mathit{S}}^{+}=\mathit{S}$ . Using the product rule, we derive the change of observer formula: ${\stackrel{·}{\mathit{F}}}^{+}=\stackrel{·}{\mathit{QF}}+\mathit{Q}\stackrel{·}{\mathit{F}}$ , leading to: ${\mathit{L}}^{+}=\mathbf{\Omega }+\mathit{Q}\mathit{L}{\mathit{Q}}^{\mathrm{T}}$ , where $\mathbf{\Omega }=\stackrel{·}{\mathit{Q}}{\mathit{Q}}^{\mathrm{T}}$ is skew-symmetric, showing that $\mathit{L}$ is not an objective tensor. The Eulerian strain rate tensor and the Cauchy stress tensor are objective, because they satisfy the change of observer formulas ${\mathit{D}}^{+}=\mathit{QD}{\mathit{Q}}^{\mathrm{T}}$ and ${\mathit{\sigma }}^{+}=\mathit{Q}\mathit{\sigma }{\mathit{Q}}^{\mathrm{T}}$ .

Now introduce the polar decomposition of the deformation gradient tensor: $\mathit{F}=\mathit{RU}=\mathit{VR}$ , where the stretch tensors $\mathit{U}$ and $\mathit{V}$ are positive definite and symmetric. The tensor $\mathit{R}$ is a proper orthogonal tensor which satisfies the transformation rule ${\mathit{R}}^{+}=\mathit{Q}\mathit{R}$ . Based on the definition of the relative deformation gradient tensor ${\mathit{F}}_{t|s}$ in equation (9) and of the relative rotation tensor ${\mathit{R}}_{t|s}$ , we derive the following transformation rules for these quantities: ${\mathit{F}}_{t|s}^{+}=\mathit{Q}(t){\mathit{F}}_{t|s}{\mathit{Q}}^{\mathrm{T}}(s)$ and ${\mathit{R}}_{t|s}^{+}=\mathit{Q}(t){\mathit{R}}_{t|s}{\mathit{Q}}^{\mathrm{T}}(s)$ .

The viscous Piola–Kirchhoff stress ${\mathit{S}}^{\mathrm{v}}$ defined in equations (7)–(12) must remain invariant. Because $\stackrel{·}{\mathbf{\Pi }}$ and $\stackrel{·}{\mathit{E}}$ are unaffected by the superimposed rigid-body motion, the present constitutive laws are frame-indifferent. An alternative proof for Models A and B is to use the transformation rule for ${\mathit{F}}_{t|s}$ to show that equations (8)b and (11)b are frame-indifferent.

Remark 1. Using the transformation rule for the relative deformation gradient tensor, one shows that the K-BKZ viscous stress is frame-indifferent. The transformation rule for ${\mathit{R}}_{t|s}$ yields the acceptability of the theory by Freed and Diethelm [9] from the point of view of material frame-indifference. Given that the models found in Capilnasiu et al. [15] and Zhang et al. [26] involve only invariant quantities (e.g., second Piola–Kirchhoff stresses and their material rates), the material frame-indifference property is straightforwardly satisfied for these theories.

Remark 2. Zhao et al. [28] propose ${\mathit{\sigma }}^{\mathrm{v}}=2\mu \kappa *\left({\mathit{F}}_{t|s}^{\mathrm{T}}\mathit{D}{\mathit{F}}_{t|s}\right)$ . Then, superimposition of the rigid-body motion ${\mathit{x}}^{+}=\mathit{c}+\mathit{Qx}$ provides:

$${\mathit{\sigma }}^{\mathrm{v}+}=2\mu \mathit{Q}(\kappa *\left({\mathbf{\Theta }}^{\mathrm{T}}{F}_{t|s}^{\mathrm{T}}{\mathbf{\Theta }}^{\mathrm{T}}\mathit{D}\mathbf{\Theta }{\mathit{F}}_{t|s}\mathbf{\Theta }\right)){\mathit{Q}}^{\mathrm{T}},\mathbf{\Theta }={\mathit{Q}}^{\mathrm{T}}(s)\mathit{Q}(t).$$

(20)

Delory et al. [28] propose ${\mathit{\sigma }}^{\mathrm{v}}=\mu (\mathit{l}+{\mathit{l}}^{\mathrm{T}})$ with $\mathit{l}=(\kappa *\stackrel{·}{\mathit{F}}){\mathit{F}}^{-1}$ . This way:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{\sigma }}^{\mathrm{v}+}=\mu \mathit{Q}(\kappa *({\mathbf{\Theta }}^{\mathrm{T}}{\mathit{Q}}^{\mathrm{T}}{\stackrel{·}{\mathit{QF}}}_{t|s}^{-1}-{\mathit{F}}_{t|s}^{-\mathrm{T}}{\mathit{Q}}^{\mathrm{T}}\stackrel{·}{\mathit{Q}}\mathbf{\Theta })){\mathit{Q}}^{\mathrm{T}}}\\ \multicolumn{1}{c}{+\mu \mathit{Q}((\kappa *{\mathbf{\Theta }}^{\mathrm{T}}\stackrel{·}{\mathit{F}}){\mathit{F}}^{-1}+{\mathit{F}}^{-\mathrm{T}}(\kappa *{\stackrel{·}{\mathit{F}}}^{\mathrm{T}}\mathbf{\Theta })){\mathit{Q}}^{\mathrm{T}},}\end{array}$$

(21)

where the skew-symmetry of $\mathbf{\Omega }$ was used. As shown in the above computations, we note that ${\mathit{\sigma }}^{\mathrm{v}+}\ne \mathit{Q}{\mathit{\sigma }}^{\mathrm{v}}{\mathit{Q}}^{\mathrm{T}}$ in general, for both theories. Therefore, these constitutive laws are not frame-indifferent.

4.2. Thermodynamic consistency

Thermodynamic consistency of Models A and B was proved by Haupt and Lion [13]. For Model C, thermodynamic consistency can be obtained in a similar way to the linear case [32]. To do so, we use the diffusive representation (13) to rewrite the viscous stress (12)a as ${\mathit{S}}^{\mathrm{v}}=2\mu \int \varphi (\zeta ){\mathit{E}}_{\zeta }^{\mathrm{v}}\mathrm{d}\zeta$ by reversing the order of integration, where ${\mathit{E}}_{\zeta }^{\mathrm{v}}$ is a memory variable defined in equation (18) with $\mathit{T}=\mathit{E}$ .

In an isothermal modelling framework, the free energy per unit of reference volume is then defined as:

$$\mathrm{\Psi }={\mathrm{\Psi }}^{\mathrm{e}}+{\mathrm{\Psi }}^{\mathrm{v}},$$

(22)

where ${\mathrm{\Psi }}^{\mathrm{e}}$ is the strain energy function associated with the elastic stress contribution, ${\mathrm{\Psi }}^{\mathrm{v}}=\mu {\int }_0^{\infty }\varphi (\zeta )\parallel {\mathit{E}}_{\zeta }^{\mathrm{v}}{\parallel }^2\mathrm{d}\zeta$ corresponds to the viscous part, and $\parallel ·\parallel$ is the Frobenius norm. According to the first and second principles of thermodynamics, the dissipation $\mathcal{D}=\mathit{S}:\stackrel{·}{\mathit{E}}-\stackrel{·}{\mathrm{\Psi }}$ per unit of reference volume must remain non-negative at all times, where the colon symbol denotes the double contraction, aka the Frobenius inner product. Using the above expression of $\mathrm{\Psi }$ , the inequality $\mathcal{D}\ge 0$ is obtained straightforwardly (see section 6.3 of Holzapfel [16] for the inclusion of incompressibility).

The dissipative behaviour of the discretised version (15) of this theory follows immediately. In fact, it suffices to define the thermodynamic potential ${\mathrm{\Psi }}^{\mathrm{v}}=\mu \sum _n{w}_n\parallel {\mathit{E}}_{{\zeta }_n}^{\mathrm{v}}{\parallel }^2$ for the viscous part along with the stress ${\mathit{S}}^{\mathrm{v}}=2\mu \sum _n{w}_n{\mathit{E}}_{{\zeta }_n}^{\mathrm{v}}$ . Under this form, similarities with Prony series theories from the literature can be identified [36].

Remark 3. The thermodynamic admissibility of the K-BKZ theory was studied by Bernstein et al. [37] (see also Rao and Rajagopal [38]). Freed and Diethelm [9] leave thermodynamic consistency to the reader’s curiosity, and Capilnasiu et al. [15] and Zhang et al. [26] do not provide the thermodynamic potentials related to their model either.

5. Elementary motions

Various illustrations are provided in the following sections, including simple shear and uniaxial tensile motions. The Cauchy stress tensor $\mathit{\sigma }$ is deduced from equation (2)b with suitable constitutive assumptions for the elastic and viscous parts. Here, we assume that the elastic response is of Mooney–Rivlin type, i.e.,

$${\mathit{\sigma }}^{\mathrm{e}}=2{\mathfrak{C}}_{+}\mathit{B}-2{\mathfrak{C}}_{-}{\mathit{B}}^{-1},2{\mathfrak{C}}_{\pm }=\frac{1}{2}\mu (1\pm \beta ),$$

(23)

where $\mathit{B}=\mathit{F}{\mathit{F}}^{\mathrm{T}}$ is the right Cauchy–Green strain tensor, $\mu >0$ is the shear modulus, and the parameters ${\mathfrak{C}}_{\pm }$ are the Mooney coefficients. The parameter $-1\le \beta \le 1$ is introduced in such a way that $\beta =1$ entails neo-Hookean material behaviour. For the sake of clarity, we now reduce the discussion to Models A, B, and C. Thus, the viscous stress ${\mathit{\sigma }}^{\mathrm{v}}$ is deduced from equations (8)–(12), respectively (see also the expressions in Table 1).

5.1. Simple shear

We consider general simple shear motions described by the deformation gradient tensor $\mathit{F}=\mathit{I}+\gamma ({\mathit{e}}_x\otimes {\mathit{e}}_z)$ , where $\gamma (z,t)$ is the shear strain, and the vectors ${\mathit{e}}_x,{\mathit{e}}_y,{\mathit{e}}_z$ form an orthonormal basis. With the present assumptions, we obtain the corresponding stress–strain relationships:

$$\mathrm{\Sigma }=ϵ+{\tau }^{\alpha }{\mathfrak{D}}^{\alpha }ϵ,\mathrm{\Sigma }=ϵ+{\tau }^{\alpha }({\mathfrak{D}}^{\alpha }ϵ+K^2ϵ{\mathfrak{D}}^{\alpha }{ϵ}^2),$$

(24)

for Models A, B and for Model C, which correspond to equations (24)a and (24)b, respectively. Here, $\mathrm{\Sigma }={\sigma }_{13}/(\mu K)$ is a non-dimensional measure of the shear stress, $ϵ=\gamma /K$ is a rescaled shear strain, and $K>0$ is a given strain magnitude. We observe that Models A and B lead to a linear fractional Kelvin–Voigt rheology. Note in passing that Models A, B, and C produce the same expression of the shear stress in the limit of infinitesimal shear strains $K\to 0$ , namely, equation (24)a.

5.1.1. Shear creep

In a standard fashion [8], we now assume that the material is initially at rest and that it is suddenly subjected to a step shear stress $\mathrm{\Sigma }=\mathrm{H}(t)$ which entails a simple shear deformation. The evolution of the strain is governed by the fractional differential equations resulting from the above constitutive relationships.

Models A and B. In this case, the relationship (24)a with $\mathrm{\Sigma }=\mathrm{H}$ leads to the linear fractional differential equation $ϵ+{\tau }^{\alpha }{\mathfrak{D}}^{\alpha }ϵ=\mathrm{H}$ . Solutions are given by the creep function [39]:

$$ϵ(t)=1-{\mathrm{E}}_{\alpha }(-{(t/\tau )}^{\alpha }),{\mathrm{E}}_{\alpha }(z)=\sum _{n=0}^{\infty }\frac{z^n}{\mathrm{\Gamma }(\alpha n+1)},$$

(25)

where ${\mathrm{E}}_{\alpha }$ is the one-parameter Mittag–Leffler function (see also Podlubny [2]). Illustrations are provided in Figure 4(a) for several values of $\alpha$ . In the elastic limit $\alpha \to 0$ , the creep response is a step shear strain (dotted line), whereas in the viscous limit $\alpha \to 1$ , we recover an exponential creep response [8] (dashed line).

Figure 4.

Time-dependent shearing motion for Models A and B. (a) Creep function, i.e., evolution of the strain for a constant applied stress. (b) Relaxation function, i.e., evolution of the stress for a constant applied strain.

Model C. Here, we obtain a nonlinear fractional differential equation deduced from equation (24)b with $\mathrm{\Sigma }=\mathrm{H}$ . In the present case, exact resolution of the creep problem seems hardly tractable. Approximate resolution might be performed, for instance, based on the representation (13)–(14) of the fractional derivative or on a perturbation approach involving the small parameter K (see, for instance, [40]).

5.1.2. Shear stress relaxation

Initially at rest, the material is suddenly subjected to a step shear strain $ϵ=\mathrm{H}(t)$ . The evolution of the stress is deduced from equation (24).

Models A and B. The stress–strain relationship (24)a produces $\mathrm{\Sigma }=ϵ+\kappa *\stackrel{·}{ϵ}$ , where we have rewritten the fractional derivative (1) in the form of a convolution product. At positive times, we therefore find $\mathrm{\Sigma }=1+\kappa$ , where the kernel $\kappa (t)$ in equation (7) follows a power-law evolution. Note in passing that the first (unit) term vanishes in fluid materials. Illustrations are provided in Figures 3 and 4(b), which can be compared to experimental results from the literature (see Figure 1). In the elastic limit $\alpha \to 0$ , the relaxation response is a step shear stress (dotted line), whereas in the viscous limit $\alpha \to 1$ , the relaxation response is singular (dashed line). Unlike their fractional counterpart, classical Kelvin–Voigt models cannot account for stress relaxation [8].

Model C. Here, direct evaluation of the stress is not straightforward (the nonlinear term of equation (24)b with $ϵ=\mathrm{H}$ seems not well-defined). Nevertheless, noting that ${\mathrm{H}}^2=\mathrm{H}$ in the weak sense, one would obtain $\mathrm{\Sigma }=1+(1+K^2)\kappa$ at positive times, provided that every step of this computation is mathematically meaningful. This way, the curves displayed in Figure 4(b) undergo a vertical dilation as K is increased.

5.2. Uniaxial tension-compression

We consider a state of uniaxial tension described by the diagonal deformation gradient tensor $\mathit{F}=\mathrm{diag}[\lambda ,{\lambda }^{-1/2},{\lambda }^{-1/2}]$ at all times. Consequently, the Cauchy stress tensor $\mathit{\sigma }$ is diagonal and its lateral components are equal, ${\sigma }_{22}={\sigma }_{33}$ . By making these lateral stresses vanish, the constitutive law (2) yields an additive decomposition of the tensile Cauchy stress:

$${\sigma }_{11}={\sigma }_{11}^{\mathrm{e}}+{\sigma }_{11}^{\mathrm{v}},$$

(26)

where the elastic part ${\sigma }_{11}^{\mathrm{e}}=2(\lambda {\mathfrak{C}}_{+}+{\mathfrak{C}}_{-})(\lambda -{\lambda }^{-2})$ follows from the Mooney–Rivlin model (23). The tensile component of the viscous stress deduced from Models A, B, and C satisfies:

$$\begin{array}{c}\multicolumn{1}{c}{{\sigma }_{11}^{\mathrm{v}}/{\eta }_{\alpha }={\lambda }^{-1}{\mathfrak{D}}^{\alpha }\lambda -{\lambda }^2{\mathfrak{D}}^{\alpha }{\lambda }^{-2},}\\ \multicolumn{1}{c}{{\sigma }_{11}^{\mathrm{v}}/{\eta }_{\alpha }={\lambda }^{-2}{\mathfrak{D}}^{\alpha }{\lambda }^2-\lambda {\mathfrak{D}}^{\alpha }{\lambda }^{-1},}\\ \multicolumn{1}{c}{{\sigma }_{11}^{\mathrm{v}}/{\eta }_{\alpha }={\lambda }^2{\mathfrak{D}}^{\alpha }{\lambda }^2-{\lambda }^{-1}{\mathfrak{D}}^{\alpha }{\lambda }^{-1},}\end{array}$$

(27)

respectively.

In the limit of small tensile strains $ϵ=3(\lambda -1)/K\simeq 0$ , we obtain the linearised expression (24)a of the dimensionless tensile stress $\mathit{\Sigma }={\sigma }_{11}/\left(\mu K\right)$ for all models. Therefore, the creep and relaxation behaviour in shear and tension-compression are formally equivalent at small strain amplitudes (see illustrations in Figure 4). Furthermore, we observe that ${\sigma }_{11}^{\mathrm{v}}/{\eta }^{\alpha }$ has identical limits for Models A and B as $\alpha \to 0$ or $\alpha \to 1$ (namely, 0 and $3\stackrel{·}{\lambda }/\lambda$ ), but that Model C has different limits (namely, ${\lambda }^4-{\lambda }^{-2}$ and $(2{\lambda }^3+{\lambda }^{-3})\stackrel{·}{\lambda }$ ). The above variability of the tensile viscous stress for large stretches provides a potential means of selecting practically relevant constitutive theories based on tensile stress relaxation experiments.

Remark 4. Upon division of ${\sigma }_{11}^{\mathrm{v}}$ in equation (27)a by the stretch $\lambda$ , we recover the expression of the “11”-component of the first Piola–Kirchhoff stress ${\mathit{P}}^{\mathrm{v}}={\mathit{\sigma }}^{\mathrm{v}}{\mathit{F}}^{-\mathrm{T}}$ provided in equation (52) of Zhao et al. [28] (see also equation (26) of Gao et al. [29]). This equivalence is caused by the symmetry of the relative deformation gradient tensor ${\mathit{F}}_{t|s}$ in the uniaxial tensile case (see definition in equation (9)). Thus, this remark applies also to the pure shear and equibiaxial tensile motions for which the deformation gradient tensor can be chosen diagonal. We conclude that the experimental results obtained in Zhao et al. [28] and Gao et al. [29] are consistent with Model A, which is a frame-indifferent version of the theory proposed therein.

6. Incremental stress and acoustoelasticity

In this section, the material is subjected to an infinitesimal perturbation $\tilde{\mathit{u}}=\mathit{x}-\overline{\mathit{x}}$ of a motionless equilibrium $\overline{\mathit{x}}$ , whose stress $\overline{\mathit{\sigma }}=-\overline{p}\mathit{I}+{\overline{\mathit{\sigma }}}^{\mathrm{e}}$ is assumed homogeneous. Hence, the equilibrium equation $\overline{\nabla }·\overline{\mathit{\sigma }}=0$ for the pre-deformation $\overline{\mathit{x}}=\overline{\mathit{FX}}$ is naturally satisfied. Here, quantities with an overbar are related to the statically pre-deformed configuration, whereas tildes mark infinitesimal increments (cf. Figure 5).

Figure 5.

Acoustoelasticity. Combination of a large static deformation and a small incremental perturbation.

The total deformation gradient reads $\mathit{F}=\overline{\mathit{F}}+\mathit{H}\overline{\mathit{F}}$ , where $\mathit{H}=\partial \tilde{\mathit{u}}/\partial \overline{\mathit{x}}$ is the incremental displacement gradient, and the total particle velocity reduces to the incremental part: $\stackrel{·}{\mathit{x}}=\stackrel{·}{\stackrel{~}{\mathit{u}}}=\tilde{\mathit{v}}$ . Based on the decomposition $\mathit{\sigma }=\overline{\mathit{\sigma }}+\tilde{\mathit{\sigma }}$ of the stress, linearisation of Cauchy’s equation of motion with respect to $\tilde{\mathit{u}}$ yields the incremental equations of motion [36]:

$$\rho {\partial }_t\tilde{\mathit{v}}=\nabla ·\tilde{\mathit{\sigma }},$$

(28)

along with the linearised incompressibility constraint $\nabla ·\tilde{\mathit{v}}=0$ .

By linearisation of equation (2)b, we arrive at the expression of the incremental stress $\tilde{\mathit{\sigma }}=-\tilde{p}\mathit{I}+{\tilde{\mathit{\sigma }}}^{\mathrm{e}}+{\tilde{\mathit{\sigma }}}^{\mathrm{v}}$ , whose elastic part satisfies:

$${\tilde{\mathit{\sigma }}}^{\mathrm{e}}=2{\mathfrak{C}}_{+}(\mathit{H}\overline{\mathit{B}}+\overline{\mathit{B}}{\mathit{H}}^{\mathrm{T}})+2{\mathfrak{C}}_{-}({\mathit{H}}^{\mathrm{T}}{\overline{\mathit{B}}}^{-1}+{\overline{\mathit{B}}}^{-1}\mathit{H}),$$

(29)

see equation (23). Based on the relationship ${\tilde{\mathit{\sigma }}}^{\mathrm{v}}=\overline{\mathit{F}}{\tilde{\mathit{S}}}^{\mathrm{v}}{\overline{\mathit{F}}}^{\mathrm{T}}$ , Models A, B, and C produce:

$${\tilde{\mathit{\sigma }}}^{\mathrm{v}}=2\mu \kappa *\tilde{\mathit{D}},{\tilde{\mathit{\sigma }}}^{\mathrm{v}}=2\mu \overline{\mathit{B}}(\kappa *\tilde{\mathit{D}})\overline{\mathit{B}},$$

(30)

where $\tilde{\mathit{D}}=\mathrm{sym}\nabla \tilde{\mathit{v}}$ is the symmetric part of the incremental velocity gradient. Note that Models A and B yield the same viscous stress increment, equation (30)a.

Remark 5. We note that the non-objective fractional velocity gradient used by Delory et al. [30] produces the same incremental stress as in equation (30)a. Therefore, the Models A and B can be used as a frame-indifferent substitute for the velocity gradient in the discussed study.

Now consider body waves propagating in a material subject to a homogeneous tri-axial stretch $\overline{\mathit{F}}=\mathrm{diag}[{\overline{\lambda }}_1,{\overline{\lambda }}_2,{\overline{\lambda }}_3]$ , with ${\overline{\lambda }}_1{\overline{\lambda }}_2{\overline{\lambda }}_3=1$ by incompressibility. The principal Cauchy stresses ${\overline{\sigma }}_i$ required to effect the pre-deformation are such that:

$${\overline{\sigma }}_i-{\overline{\sigma }}_j=2{\mathfrak{C}}_{+}({\overline{\lambda }}_i^2-{\overline{\lambda }}_j^2)-2{\mathfrak{C}}_{-}({\overline{\lambda }}_i^{-2}-{\overline{\lambda }}_j^{-2}),$$

(31)

where the coefficients ${\mathfrak{C}}_{\pm }$ are defined in equation (23).

We study harmonic principal body waves of the form $\tilde{\mathit{u}}=\hat{\mathit{u}}{\mathrm{e}}^{\mathrm{i}(\omega t-\mathrm{kx})}$ , where $\hat{\mathit{u}}$ is the constant amplitude vector, $\omega$ is the angular frequency, k is the wavenumber, and x is the direction of propagation. We see from $\nabla ·\tilde{\mathit{v}}=0$ that the polarisation of the wave must be transverse to accommodate incremental incompressibility: $\hat{\mathit{u}}·{\mathit{e}}_x=0$ .

Consider the wave polarised along $y:\hat{\mathit{u}}=\hat{u}{\mathit{e}}_y$ . The harmonic amplitudes of ${\partial }_t\tilde{\mathit{v}}$ and $\kappa *\tilde{\mathit{D}}$ are $-{\omega }^2\hat{u}{\mathit{e}}_y$ and $-{\left(\mathrm{i}\omega \tau \right)}^{\alpha }\mathrm{i}k\hat{u}\mathrm{sym}({\mathit{e}}_y\otimes {\mathit{e}}_x)$ , respectively. Using the incremental equation of motion (28), we obtain the following dimensionless dynamic moduli $M^{*}=\rho {\omega }^2/(\mu k^2)$ for Models A–B and for Model C:

$$\begin{array}{c}{M}^{*}=\frac{1+\beta }{2}{\overline{\lambda }}_1^2+\frac{1-\beta }{2}{\overline{\lambda }}_2^{-2}+{\left(\mathrm{i}\omega \tau \right)}^{\alpha },\hfill \\ M^{*}=\frac{1+\beta }{2}{\overline{\lambda }}_1^2+\frac{1-\beta }{2}{\overline{\lambda }}_2^{-2}+{\overline{\lambda }}_1^2{\overline{\lambda }}_2^2{\left(\mathrm{i}\omega \tau \right)}^{\alpha },\hfill \end{array}$$

(32)

respectively. The real and imaginary parts of $M^{*}$ represent the storage modulus $\Re \mathfrak{e}{M}^{*}$ and the loss modulus $\Im \mathfrak{m}{M}^{*}$ , respectively. We note that Models A and B cannot capture a change of loss modulus with the pre-stretch.

We can deduce other dispersion properties from equation (32) (see section 2.3 of Carcione [8]). For instance, the dissipation factors $d^{*}=\Im \mathfrak{m}{M}^{*}/\Re \mathfrak{e}{M}^{*}$ for Models A–B and for Model C are given by:

$$\begin{array}{c}\multicolumn{1}{c}{d^{*}=\frac{{\left(\omega \tau \right)}^{\alpha }\sin \left(\alpha \frac{\pi }{2}\right)}{\frac{1+\beta }{2}{\overline{\lambda }}_1^2+\frac{1-\beta }{2}{\overline{\lambda }}_2^{-2}+{\left(\omega \tau \right)}^{\alpha }\cos \left(\alpha \frac{\pi }{2}\right)},}\\ \multicolumn{1}{c}{d^{*}=\frac{{\overline{\lambda }}_1^2{\overline{\lambda }}_2^2{(\omega \tau )}^{\alpha }\sin \left(\alpha \frac{\pi }{2}\right)}{\frac{1+\beta }{2}{\overline{\lambda }}_1^2+\frac{1-\beta }{2}{\overline{\lambda }}_2^{-2}+{\overline{\lambda }}_1^2{\overline{\lambda }}_2^2{\left(\omega \tau \right)}^{\alpha }\cos \left(\alpha \frac{\pi }{2}\right)},}\end{array}$$

(33)

respectively. In the special case of neo-Hookean behaviour $\beta =1$ , we see that the dissipation factor for Models A and B is sensitive to the stretch along the propagation direction, while for Model C, it depends on the stretch along the polarisation direction. Upon nondimensionlisation, the phase velocity $c=\omega /\Re \mathfrak{e}k$ and the attenuation factor $a=-\Im \mathfrak{m}k$ satisfy:

$$c^2=\frac{2(1+d^{*2})}{1+\sqrt{1+d^{*2}}}\Re \mathfrak{e}{M}^{*},a^2=\frac{{(\omega \tau )}^2}{\Re \mathfrak{e}{M}^{*}}\frac{\sqrt{1+d^{*2}}-1}{2(1+d^{*2})},$$

(34)

for all models.

For illustrative purposes, consider a neo-Hookean viscous body for which $\beta =1$ , which is subjected to uni-axial traction of stretch $\overline{\lambda }$ (see, for instance, the experimental configuration in Delory et al. [41]). Three types of principal waves may propagate in such a deformed body [42]. First, the wave travelling in a direction perpendicular to the uni-axial tension and polarised along that direction (e.g., a wave propagating along the x-axis and polarised along the y-axis, which is the direction of uni-axial tension). Here, ${\overline{\lambda }}_1={\overline{\lambda }}_3={\overline{\lambda }}^{-1/2}$ , ${\overline{\lambda }}_2=\overline{\lambda }$ . Wave dispersion properties are then deduced from equations (32)–(34).

In the first row of Figure 6, we show the variations of the corresponding storage moduli $\Re \mathfrak{e}{M}^{*}$ and loss moduli $\Im \mathfrak{m}{M}^{*}$ with the dimensionless frequency $\omega \tau$ for some given values of pre-stretch, respectively, found from:

$$M^{*}={\overline{\lambda }}^{-1}+(\mathrm{i}\omega \tau {)}^{\alpha },M^{*}={\overline{\lambda }}^{-1}+\overline{\lambda }(\mathrm{i}\omega \tau {)}^{\alpha }.$$

(35)

Figure 6.

Principal plane shear waves propagating in a neo-Hookean viscous material $(\beta =1)$ subject to uni-axial tension. The material has fractional viscosity of order $\alpha =0.3$ and the pre-stretch equals $\overline{\lambda }=\frac{3}{5},1,\frac{5}{3}$ . Dimensionless storage modulus (a) and loss modulus (b) for waves polarised along the stretching direction (top) and for waves with propagation direction and polarisation transverse to the uni-axial tension (middle). Same for shear waves propagating along the stretching axis (bottom).

In particular, the figure shows the evolution of the loss modulus in terms of the frequency (Figure 6(b)). As predicted earlier, the loss modulus is unaffected by variations of the pre-stretch for Models A and B, whereas the loss modulus increases with increasing values of $\overline{\lambda }$ when Model C is used.

For Models A and B—i.e., for $M^{*}$ satisfying (35)a — the relationships (33)–(34) entail the following asymptotic expansions for $0<\alpha <1$ :

$$\begin{array}{c}{d}^{*}\tilde{0}(\omega \tau {)}^{\alpha }\sin \left(\alpha \frac{\pi }{2}\right)\overline{\lambda },d^{*}\tilde{\infty }(1-{(\omega \tau )}^{-\alpha }\sec \left(\alpha \frac{\pi }{2}\right)/\overline{\lambda })\tan \left(\alpha \frac{\pi }{2}\right),\hfill \\ c\tilde{0}(1+{(\omega \tau )}^{2\alpha }\sin {\left(\alpha \frac{\pi }{2}\right)}^2){\overline{\lambda }}^{-\frac{1}{2}},c\tilde{\infty }(\omega \tau {)}^{\frac{\alpha }{2}}\sec (\alpha \frac{\pi }{4}),\hfill \\ a\tilde{0}\frac{1}{2}(\omega \tau {)}^{1+\alpha }\sin \left(\alpha \frac{\pi }{2}\right){\overline{\lambda }}^{\frac{3}{2}},a\tilde{\infty }(\omega \tau {)}^{1-\frac{\alpha }{2}}\sin \left(\alpha \frac{\pi }{2}\right),\hfill \end{array}$$

(36)

at low- and high-dimensionless frequency $\omega \tau$ . The above quantities (33)a–(34) with $\beta =1$ and the approximations (36) are displayed in Figure 7 for several values of the pre-stretch. Although the loss modulus does not vary with pre-stretch, $\overline{\lambda }$ still influences the dispersion properties through modification of the storage modulus. Of course, similar computations can be carried out for viscous Mooney–Rivlin solids $-1\le \beta <1$ , for Model C, and even with other wave polarisations, given that equation (34) is model-independent.

Figure 7.

Dispersion characteristics (33)–(36) in terms of the dimensionless frequency: (a) dissipation factor, (b) normalised phase velocity, and (c) normalised attenuation coefficient. Models A and B with $\alpha =0.3$ , $\beta =1$ , and stretch $\overline{\lambda }=\frac{3}{5},1,\frac{5}{3}$ , where shear waves are polarised along the stretching direction.

The second principal shear wave propagates and is polarised transversely to the direction of stretching. Then, ${\overline{\lambda }}_1={\overline{\lambda }}_2={\overline{\lambda }}^{-1/2}$ , ${\overline{\lambda }}_3=\overline{\lambda }$ , and:

$$M^{*}={\overline{\lambda }}^{-1}+(\mathrm{i}\omega \tau {)}^{\alpha },M^{*}={\overline{\lambda }}^{-1}+{\overline{\lambda }}^{-2}(\mathrm{i}\omega \tau {)}^{\alpha },$$

(37)

for Models A–B and for Model C, respectively. Finally, the third wave propagates along the direction of stretching and ${\overline{\lambda }}_1=\overline{\lambda }$ , ${\overline{\lambda }}_2={\overline{\lambda }}_3={\overline{\lambda }}^{-1/2}$ . In this case, we have the following dimensionless dynamic moduli for Models A–B and for Model C:

$$M^{*}={\overline{\lambda }}^2+(\mathrm{i}\omega \tau {)}^{\alpha },M^{*}={\overline{\lambda }}^2+\overline{\lambda }(\mathrm{i}\omega \tau {)}^{\alpha },$$

(38)

respectively, see the second and third rows of Figure 6 for the variations of the storage and loss moduli with the frequency when $\alpha =0.3$ .

As shown in the figure, the first two waves yield a decrease of the storage modulus with increasing values of the applied stretch (see the first two rows of Figure 6), whereas the third wave leads to the opposite trend (third row of Figure 6). This feature can be explained by the relationship between the propagation direction and the stretching direction. In fact, the first two waves both propagate orthogonally to the direction of stretching. Hence, by virtue of incompressibility, they are subject to a contraction along their propagation direction if $\overline{\lambda }>1$ . In contrast, the third wave travels along the direction of stretching, which implies that this wave is subject to an extension along its propagation direction if $\overline{\lambda }>1$ .

7. Conclusion

Several fractional viscous stresses were investigated. Among others, Models A, B, and C are both physically satisfactory from the points of view of material frame-indifference and thermodynamic consistency. Nevertheless, one might prefer Model A or Model B due to their ease of use (e.g., for the study of shearing and tensile motions in creep and relaxation) and for the mathematical properties of the viscous limit. Although Models A and B both correspond to the classical Newtonian theory of viscosity for $\alpha \to 1$ , their elastic limits $\alpha \to 0$ differ. In passing, we observe that linear combinations of these models produce Mooney–Rivlin stresses (23) in the elastic limit (see the expressions in Table 1).

Despite this observation, it appears that Models A and B are very similar in many respects. In fact, they lead to the same creep and relaxation behaviour in simple shear. They produce also the same incremental stresses, as shown by the theory of acoustoelasticity. Nevertheless, they can be distinguished according to their creep and relaxation behaviour at large stretches (27), which could be useful for experimental calibration purposes.

Interestingly, we note that the dispersion relationships obtained for body waves in section 6 are reminiscent of experimental and theoretical results obtained for Lamb wave propagation in stretched plates [30], as well as for wave propagation in stretched strips [41]. On this basis, an experimentally calibrated theory that complies with the elementary principles of physics has yet to be established. This study provides relevant models to carry out such a task.

On one hand, if the loss modulus does not vary significantly with applied stretches, then Models A or B might be satisfactory. On the other hand, if the experimental data shows that the loss modulus varies with the pre-stretch, then Model C might be preferred, despite its mathematical complexity for elementary shearing and tensile motions. Alternatively, a more general fractional viscosity theory based on Model A or Model B could be derived by including further tensor invariants [43] (e.g., in the framework of the K-BKZ theory), especially if it provides a better fit with experiments. Such an approach was followed by Delory et al. [30] for the modelling of their experimental data.

Some limitations of this study are the restriction to incompressible motions and to fractional viscosity theories with one differential order only. Related compressible theories could be obtained by removing the incompressibility constraint and by incorporating the missing strain and strain-rate tensor invariants. In the same spirit, fractional Maxwell models with two fractional differential orders could be developed [11].