Mathematical Modelling of Viscoelastic Media Without Bulk Relaxation Via Fractional Calculus Approach

1. Introduction

It is known [1] that pure elastic bodies do not exist in nature. All media possess viscoelastic features to one degree or another, and their main physical-and-mechanical properties are time-dependent. Due to wide application of the theory of elasticity in studies of advanced and traditional materials, much attention is given to modelling and investigative techniques of viscoelastic media and bodies subjected to various types of loading [1,2,3,4,5].

During the last three decades Fractional Calculus has gained wide acceptance in modelling such viscoelastic bodies as beams, plates, and shells [6,7]. Their damping features are described most often by defining the Young’s operator by the simplest fractional derivative models, namely: Kelvin-Voight model, Maxwell model, and standard linear solid model [8,9,10]. As this takes place, the Poisson’s ratio of a viscoelastic material is frequently assumed to be a constant [11,12]. However, it has been emphasized in [13] that the fractional derivative Kelvin-Voigt model with a time-independent Poisson’s ratio is only acceptable for the description of the dynamic behaviour of elastic bodies in a viscoelastic medium [6,14,15] or on a viscoelastic foundation[16].

As experimental data have shown [17], the Poisson’s ratio is always a time-dependent operator [18,19,20], and only the bulk extension-compression operator could be considered as a constant value, since for the most viscoelastic materials it varies weakly during deformation [1,2].

The detailed reviews of `traditional’ fractional calculus models in viscoelasticity (`traditional’ in the sense that such models consider time-independent Poisson’s ratios) are given in [6,10,13,21]. In the present paper, the fractional derivative models involving the time-dependent Poisson’s operators will be studied, which allows one to reveal rather interesting properties of advanced viscoelastic materials, among them auxetic materials possessing negative Poisson’s ratios [22,23,24,25].

For solving different dynamic problems of mechanics of viscoelastic solids and structures, it is essential to know the form of viscoelastic operators entering in governing equations. For example, Poisson’s ratio and Young’s modulus are involved in the cylindrical rigidity of plates and shells, and therefore in the Hertz’s law for the solution of the problem of viscoelastic contact interaction during impact [26,27,28]. Therefore, their pinpointing is of paramount importance in the impact response analysis.

It is well known [29,30] that each isotropic elastic material possesses only two independent material constants, and all others are expressed in terms of two constants that should be preassigned or determined experimentally. Possible combinations are presented in Table 1, where two prescribed constants are shown in the first column, while the others are determined in terms of two given constants according to formulas, which are located at the intersections of the corresponding lines and columns.

From Table 1 it is seen that for materials, the bulk relaxation of which could be neglected, i.e., to consider K as a constant value, two independent material moduli could be assigned via four ways.

Similarly, in the case of isotropic viscoelastic media, material properties are time-dependent and are described by operators, which should be expressed in terms of two preassigned ones (or determined from experimental data) utilizing the correspondence principle and relationships given in Table 1.

For solving one-dimensional dynamic problems of viscoelasticity, it is a need to know the Young’s and shear operators defining phase velocities and coefficients of attenuation of viscoelastic longitudinal and shear waves. Thus, modelling of these time-dependent operators without volume relaxation (K = const) has been considered using the following viscoelastic fractional derivative models:

• Young’s modulus is preassigned by the fractional derivative Kevin-Voigt model [13,31,32,33];
• Young’s modulus is relaxed according to the fractional derivative Maxwell model [33];
• Young’s modulus is relaxed according to the fractional derivative standard linear solid model [13,31,33,34];
• Shear modulus is defined by Scott Blair model [13];
• Shear modulus is preassigned by the fractional derivative Kelvin-Voigt model [13,31,32,33,35];
• Shear modulus is relaxed according the Maxwell model [13,33];
• Shear modulus is relaxed according to the standard linear solid model [13,33,36,37].

Interestingly to note that only the first and second Lamé constants,λ and µ, or the bulk and shear moduli, K and µ, appear in Hooke’s law for three-dimensional media, but not Young’s modulus E, or Poisson’s ratio $\nu$. This indicates that K, $\lambda$ and µ are the most intrinsic operators to express stress in terms of strain when studying wave propagation in 3D viscoelastic media [38,39,40].

That is why in the present paper, the emphasis will be make on the comprehensive analysis of time-dependent operators for Lamé parameters. Using the procedure for the study of viscoelastic operators suggested in Rossikhin and Shitikova [13,34], below for the first time the models based on the application of fractional derivatives will be studied for the cases when the first Lame parameter or the P-wave modulus is given a priori, without considering the bulk relaxation.

2. Models of Viscoelasticity Involving Fractional Order Operators With Time-Dependent Poisson’s Ratio

2.1. Preliminary Remarks

The rheological equations of the simplest fractional derivative models of viscoelasticity widely used in mechanics are the following [10,13]:

1.

Scott Blair model (Figure 1A)

$$\sigma =E{\tau }^{\gamma }{D}^{\gamma }\epsilon ,$$

(1)

2.

Maxwell model (Figure 1B)

$$\sigma +{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }\sigma =E_{\infty }{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }\epsilon ,$$

(2)

or

$$\begin{array}{c}\sigma \left(t\right)=E_{\infty }{\displaystyle \frac{{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}\epsilon \left(t\right)=E_{\infty }\left(1-{\displaystyle \frac{1}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}\right)\epsilon \left(t\right)=E_{\infty }\left(1-{϶}_{\gamma }^{*}\left({\tau }_i^{\gamma }\right)\right)\epsilon \left(t\right)\hfill \\ =E_{\infty }\left[\epsilon \left(t\right)-{\int }_0^t{϶}_{\gamma }(t^{\prime }/{\tau }_{\epsilon })\epsilon (t-t^{\prime })d{t}^{\prime }\right],\hfill \end{array}$$

(3)

3.

Kelvin-Voigt model (Figure 1C)

$$\sigma =E_0\left(\epsilon +{\tau }_{\sigma }^{\gamma }{D}^{\gamma }\epsilon \right)=\tilde{E}\epsilon ,\tilde{E}=E_0\left(1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }\right),$$

(4)

4.

standard linear solid model (Figure 1D and Figure 1E)

$$\sigma +{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }\sigma =E_0\left(\epsilon +{\tau }_{\sigma }^{\gamma }{D}^{\gamma }\epsilon \right),$$

(5)

or

$$\sigma \left(t\right)=E_0{\displaystyle \frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}\epsilon \left(t\right)=E_{\infty }\left(1-{\nu }_{\epsilon }{\displaystyle \frac{1}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}\right)\epsilon \left(t\right)=\tilde{E}\left(t\right)\epsilon \left(t\right),\tilde{E}=E_{\infty }[1-{\nu }_{\epsilon }{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)],$$

(6)

where σandε are the stress and the strain, respectively,${\tau }_{\epsilon }$ and ${\tau }_{\sigma }$ are the relaxation and retardation (or creep) times, respectively, $E_{\infty }$ and $E_0$ are the nonrelaxed (instantaneous) and relaxed (prolonged) moduli of elasticity, ${\nu }_{\epsilon }=\Delta E{E}_{\infty }^{-1}$, $\Delta E=E_{\infty }-E_0$ is the defect of the modulus, i.e. the value characterizing the decrease in the elastic modulus from its nonrelaxed value to its relaxed magnitude, $D^{\gamma }$ is the Riemann-Liouville fractional derivative [41] of the order $0<\gamma \le 1$

$$D^{\gamma }\sigma ={\displaystyle \frac{d}{dt}}{\int }_0^t{\displaystyle \frac{\sigma \left(t^{\prime }\right)}{\Gamma (1-\gamma ){(t-t^{\prime })}^{\gamma }}}d{t}^{\prime },$$

(7)

$\Gamma$ is the Gamma-function, ${϶}_{\gamma }^{*}\left({\tau }_i^{\gamma }\right)={\displaystyle \frac{1}{1+{\tau }_i^{\gamma }{D}^{\gamma }}}$ ($i=\sigma ,\epsilon$) is the dimensionless Rabotnov’s operator [9,34], and ${϶}_{\gamma }(t/{\tau }_i)={\displaystyle \frac{t^{\gamma -1}}{{\tau }_i^{\gamma }}}{\sum }_{n=0}^{\infty }{\displaystyle \frac{{(-1)}^n{(t/{\tau }_i)}^{\gamma n}}{\Gamma \left[\gamma \right(n+1\left)\right]}}$ is the fractional exponential function [1,2], which at $\gamma =1$ goes over into a conventional exponential function.

It is known [13] that both versions of the standard linear solid model shown in Figure 1D and Figure 1E are described by one and the same equation (4), wherein model parameters are interconnected by the following relationship:

$${\left({\displaystyle \frac{{\tau }_{\epsilon }}{{\tau }_{\sigma }}}\right)}^{\gamma }={\displaystyle \frac{E_0}{E_{\infty }}}.$$

(8)

Equations (2)-(6) govern stress-strain relationships for one-dimensional viscoelastic media. The constitutive equation connecting the strain and stress in a linear viscoelastic isotropic medium has the form

$${\sigma }_{ij}=\tilde{\lambda }\left(t\right){\epsilon }_{kk}{\delta }_{ij}+2\tilde{\mu }\left(t\right){\epsilon }_{ij},$$

(9)

or

$${\sigma }_{ij}=\tilde{K}\left(t\right){\epsilon }_{kk}{\delta }_{ij}+2\tilde{\mu }\left(t\right)\left({\epsilon }_{ij}-{\displaystyle \frac{1}{3}}{\epsilon }_{kk}{\delta }_{ij}\right),$$

(10)

where ${\sigma }_{ij}$ and ${\epsilon }_{ij}$ are stress and strain tensors components, $\tilde{\lambda }\left(t\right)$, $\tilde{\mu }\left(t\right)$ and $\tilde{K}\left(t\right)$ are, respectively, time-dependent Lamé and bulk operators, and ${\delta }_{ij}$ is the Kronecker delta.

Thus, for solving three-dimensional dynamic problems of viscoelasticity, it is a need to know the form of two time-dependent operators: $\tilde{\lambda }\left(t\right)$ and $\tilde{\mu }\left(t\right)$, or $\tilde{K}\left(t\right)$ and $\tilde{\mu }\left(t\right)$.

Below we will consider models involving time-dependent Poisson’s ratio but without volume relaxation, i.e. when bulk modulus is time-independent (this assumption is due to the fact that for many viscoelastic materials volumetric relaxation is much smaller than the shear relaxation)

$$\tilde{K}=K\tilde{I},$$

(11)

where $K=const$ is a certain constant which could take on the value of nonrelaxed bulk modulus $K_{\infty }$ or relaxed bulk modulus $K_0$, and $\tilde{I}$ is the identity operator.

Therefore, it is necessary to assign the second operator, i.e., one of the Lamé operators $\tilde{\lambda }\left(t\right)$ or $\tilde{\mu }\left(t\right)$, or P-wave operator $\tilde{M}=\tilde{\lambda }+2\tilde{\mu }$, what could be done using the fractional derivative Maxwell, Scott Blair, Kelvin-Voigt, or standard linear solid models (2)-(6). Knowing the form of two viscoelastic operators, it is possible to define the form of the time-dependent Poisson’s operator and all other operators.

2.2. Modelling of the Shear Operator $\tilde{\mu }$ Using the Fractional Derivative Kelvin-Voigt Model

The shear operator $\tilde{\mu }$ is most frequently preassigned via the fractional derivative Kelvin-Voigt model (4) as

$$\tilde{\mu }={\mu }_0\left[1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }\right],$$

(12)

where ${\mu }_0$ is the relaxed shear modules, ${\tau }_{\sigma }$ is the retardation time during shear deformations, while the bulk operator is assumed to be constant according to (11).

In order to evaluate dynamic response of viscoelastic bodies, for example, impact response, it is necessary to calculate the Young’s operator. For this purpose using the Volterra principle and formula from the third line in Table 1:

$$\tilde{E}={\displaystyle \frac{9K_0\tilde{\mu }}{3K_0+\tilde{\mu }}}$$

(13)

First we could write the operator

$$3K_0+\tilde{\mu }=(3K_0+{\mu }_0)\left(1+t_{\sigma }^{\gamma }{D}^{\gamma }\right)$$

(14)

where $t_{\sigma }^{\gamma }={\mu }_0{\tau }_{\sigma }^{\gamma }{(3K_0+{\mu }_0)}^{-1}$.

Then we find the operator reverse to (14), i.e.,

$${\left(3K_0+\tilde{\mu }\right)}^{-1}={\displaystyle \frac{1}{3K_0+{\mu }_0}}{\displaystyle \frac{1}{1+t_{\sigma }^{\gamma }{D}^{\gamma }}}={\left(3K_0+{\mu }_0\right)}^{-1}{϶}_{\gamma }^{*}\left(t_{\sigma }^{\gamma }\right).$$

(15)

Substituting (12) and (15) in (13) and considering the formula [31,34]

$$t_{\sigma }^{\gamma }{D}^{\gamma }·{϶}_{\gamma }^{*}\left(t_{\sigma }^{\gamma }\right)={\displaystyle \frac{t_{\sigma }^{\gamma }{D}^{\gamma }}{1+t_{\sigma }^{\gamma }{D}^{\gamma }}}=1-{϶}_{\gamma }^{*}\left(t_{\sigma }^{\gamma }\right)$$

(16)

yield

$$\tilde{E}=9K_0\left[1-{\displaystyle \frac{E_0}{3{\mu }_0}}{϶}_{\gamma }^{*}\left(t_{\sigma }^{\gamma }\right)\right].$$

(17)

Now the Poisson’s operator could be calculated via the following formula:

$${\displaystyle \frac{\tilde{E}}{1-2\tilde{\nu }}}=3K_0.$$

(18)

Substituting (17) in (18), we find

$$\begin{array}{c}\tilde{\nu }=-1+{\displaystyle \frac{E_0}{2{\mu }_0}}{϶}_{\gamma }^{*}\left(t_{\sigma }^{\gamma }\right)=-1+\left(1+{\nu }_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{{\mu }_0}{3K_0+{\mu }_0}}{\tau }_{\sigma }^{\gamma }\right)\hfill \\ \hfill =-1+\left(1+{\nu }_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{1-2{\nu }_0}{3}}{\tau }_{\sigma }^{\gamma }\right),\end{array}$$

(19)

or

$$\begin{array}{c}\nu \left(t\right)=\tilde{\nu }H\left(t\right)={\nu }_0-{\displaystyle \frac{E_0}{6K_0}}{\tau }_{\sigma }^{\gamma }{\displaystyle \frac{t^{-\gamma }}{\Gamma \left(1-\gamma \right)}}-1+\left(1+{\nu }_0\right)\left\{1-e_{\gamma }\left[{\left({\displaystyle \frac{3K_0}{{\mu }_0}}+1\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right]\right\}\hfill \\ ={\nu }_0-\left(1+{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{3}{1-2{\nu }_0}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right],\hfill \end{array}$$

(20)

where $e_{\gamma }\left(t\right)\equiv E_{\gamma }\left(-t^{\gamma }\right)$ [40], and $E_{\gamma }\left(t\right)\equiv {\sum }_{n=0}^{\infty }{\displaystyle \frac{t^n}{\Gamma \left(\gamma n+1\right)}}$ is the Mittag-Leffler function.

From relationship (20) the limiting magnitudes of the Poisson’s ration could be calculated

$$\underset{t\to 0}{lim}\tilde{\nu }H\left(t\right)=-1,\underset{t\to \infty }{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{3K_0-2{\mu }_0}{2(3K_0+{\mu }_0)}}={\nu }_0.$$

(21)

From (21) it follows that the model (12) with $K_0=const$ could describe the behavior of isotropic viscoelastic auxetics (materials with negative Poisson’s ratios), in so doing the Poisson’s ratio varies from -1 to its relaxed magnitude ${\nu }_0$, what does not violate the laws of thermodynamics [29,30].

2.3. Modelling the Shear Operator $\tilde{\mu }$ Using the Fractional Derivative Maxwell Model

If the fractional derivative Maxwell model (3) is applied for describing viscoelastic bodies, then the shear operator $\tilde{\mu }$ could be written in the form

$$\tilde{\mu }={\mu }_{\infty }\left[1-{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right],$$

(22)

where ${\mu }_{\infty }$ is the nonrelaxed magnitude of the shear modulus, in so doing the volumetric operator is still considered as a constant $K_{\infty }$.

Using the procedure described above for the Kelvin-Voigt model, we could similarly obtain for the Maxwell model

$$\tilde{E}=E_{\infty }\left[1-{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right],$$

(23)

$$\tilde{\nu }={\displaystyle \frac{1}{2}}-{\displaystyle \frac{E_{\infty }}{6K_{\infty }}}\left[1-{϶}_{\gamma }^{*}\left(\left({\displaystyle \frac{3K_{\infty }+{\mu }_{\infty }}{3K_{\infty }}}\right){\tau }_{\epsilon }^{\gamma }\right)\right],$$

(24)

or

$$\tilde{\nu }={\nu }_{\infty }+\left({\displaystyle \frac{1}{2}}-{\nu }_{\infty }\right){϶}_{\gamma }^{*}\left(\left({\displaystyle \frac{3K_{\infty }+{\mu }_{\infty }}{3K_{\infty }}}\right){\tau }_{\epsilon }^{\gamma }\right)={\nu }_{\infty }+\left({\displaystyle \frac{1}{2}}-{\nu }_{\infty }\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{\frac{3}{2}}{{\nu }_{\infty }+1}}{\tau }_{\epsilon }^{\gamma }\right),$$

(25)

where ${\nu }_{\infty }={\displaystyle \frac{3K_{\infty }-{\mu }_{\infty }}{2(3K_{\infty }+{\mu }_{\infty })}}$ is the nonrelaxed magnitude of the Poisson’s ratio.

Then the Poisson’s operator will take the form

$$\begin{array}{c}\nu \left(t\right)=\tilde{\nu }H\left(t\right)={\nu }_{\infty }+\left({\displaystyle \frac{1}{2}}-{\nu }_{\infty }\right)\left\{1-e_{\gamma }\left[{\left({\displaystyle \frac{3K_{\infty }}{3K_{\infty }+{\mu }_{\infty }}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right]\right\}\hfill \\ \hfill ={\displaystyle \frac{1}{2}}-\left({\displaystyle \frac{1}{2}}-{\nu }_{\infty }\right)e_{\gamma }\left[{\left({\displaystyle \frac{{\nu }_{\infty }+1}{\frac{3}{2}}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right],\end{array}$$

(26)

whence the limiting values of the operator are the following:

$$\underset{t\to 0}{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{\frac{3}{2}{K}_{\infty }-{\mu }_{\infty }}{3K_{\infty }+{\mu }_{\infty }}},\underset{t\to \infty }{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{1}{2}}.$$

(27)

From relationships (27) it is seen that for the fractional derivative Maxwell model, Poisson’s ratio could increase from ${\nu }_{\infty }$ to its limiting value of 0.5, what means that this model is suitable for the analysis of viscoelastic rubber-like materials.

2.4. Scott Blair Model for Shear Relaxation

Some authors prefer to use the simplest fractional derivative model, i.e., the Scott Blair element (1) for modelling the shear operator

$$\tilde{\mu }={\mu \tau }^{\gamma }{D}^{\gamma }$$

(28)

and assume that volumetric relaxation is absent, i.e. $\tilde{K}=K\tilde{I}$, where $K=const$.

In this case

$$3\tilde{K}+\tilde{\mu }=3K+{\mu \tau }^{\gamma }{D}^{\gamma }=3K\left(1+{\displaystyle \frac{\mu }{3K}}{\tau }^{\gamma }{D}^{\gamma }\right).$$

(29)

Then the Poisson’s operator takes the form

$$\tilde{\nu }={\displaystyle \frac{3K-2\tilde{\mu }}{2(3K+\tilde{\mu })}}=\left({\displaystyle \frac{1}{2}}-{\displaystyle \frac{\mu }{3K}}{\tau }^{\gamma }{D}^{\gamma }\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{\mu }{3K}}{\tau }^{\gamma }\right)=-1+{\displaystyle \frac{3}{2}}{϶}_{\gamma }^{*}\left({\displaystyle \frac{\mu }{3K}}{\tau }^{\gamma }\right),$$

(30)

whence it follows that

$$\nu \left(t\right)=\tilde{\nu }H\left(t\right)=-1+{\displaystyle \frac{3}{2}}{϶}_{\gamma }^{*}\left({\displaystyle \frac{\mu }{3K}}{\tau }^{\gamma }\right),$$

(31)

$$\underset{t\to 0}{lim}\tilde{\nu }H\left(t\right)=-1,\underset{t\to \infty }{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{1}{2}}.$$

(32)

From (32) it is evident that according to this model the Poisson’s ratio could vary in a very broad range, namely, form -1 to 0.5. Thus, this model is thermodynamically admissible for viscoelastic auxetics.

2.5. Modelling the Shear Operator $\tilde{\mu }$ Via the Fractional Derivative Standard Linear Solid Model

If the fractional derivative standard linear solid model (6) is applied for describing viscoelastic bodies, then the shear operator $\tilde{\mu }$ has the form

$$\tilde{\mu }={\mu }_0{\displaystyle \frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}},$$

(33)

or

$$\tilde{\mu }={\mu }_0{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)+{\mu }_0{\displaystyle \frac{{\tau }_{\sigma }^{\gamma }}{{\tau }_{\epsilon }^{\gamma }}}{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }·{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right).$$

(34)

Substituting (16) in (34) and introducing the notation (similar to relationship (8))

$${\mu }_{\infty }={\mu }_0{\left({\displaystyle \frac{{\tau }_{\sigma }}{{\tau }_{\epsilon }}}\right)}^{\gamma }$$

(35)

yield

$$\tilde{\mu }={\mu }_{\infty }\left[1-{\nu }_{\mu }^{\epsilon }{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right],$$

(36)

where ${\nu }_{\mu }^{\epsilon }=({\mu }_{\infty }-{\mu }_0){\mu }_{\infty }^{-1}$, in so doing operator $\tilde{K}$ is still defined by (11).

For this model the time-dependent Poisson’s operator will take the form

$$\begin{array}{c}\tilde{\nu }={\displaystyle \frac{1}{2\left(3K_0+{\mu }_0\right)\left(3K_0+{\mu }_{\infty }\right)}}\left\{9{K_0}^2-3K_0\left(2{\mu }_{\infty }-{\mu }_0\right)\right.\hfill \\ \left.-2{\mu }_0{\mu }_{\infty }+9K_0\left({\mu }_{\infty }-{\mu }_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{3K_0+{\mu }_{\infty }}{3K_0+{\mu }_0}}{\tau }_{\epsilon }^{\gamma }\right)\right\}\\ \hfill ={\nu }_{\infty }+\left({\nu }_{\infty }-{\nu }_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{{\nu }_0+1}{{\nu }_{\infty }+1}}{\tau }_{\epsilon }^{\gamma }\right)={\nu }_{\infty }+\left({\nu }_{\infty }-{\nu }_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{1-2{\nu }_0}{1-2{\nu }_{\infty }}}{\tau }_{\sigma }^{\gamma }\right),\end{array}$$

(37)

or

$$\begin{array}{c}\nu \left(t\right)=\tilde{\nu }H\left(t\right)={\displaystyle \frac{1}{2\left(3K_0+{\mu }_0\right)\left(3K_0+{\mu }_{\infty }\right)}}\left\{9{K_0}^2-3K_0\left(2{\mu }_{\infty }-{\mu }_0\right)-2{\mu }_0{\mu }_{\infty }\begin{array}{c}\hspace{1em}\\ \hspace{1em}\end{array}\right.\hfill \\ +\left.9K_0\left({\mu }_{\infty }-{\mu }_0\right)\left\{1-e_{\gamma }\left[{\left({\displaystyle \frac{3K_0+{\mu }_0}{3K_0+{\mu }_{\infty }}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right]\right\}\right\}\hfill \\ ={\nu }_0+\left({\nu }_{\infty }-{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{{\nu }_{\infty }+1}{{\nu }_0+1}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right]={\nu }_0+\left({\nu }_{\infty }-{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{1-2{\nu }_{\infty }}{1-2{\nu }_0}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right].\hfill \end{array}$$

(38)

From relationship (38) it is seen that the limiting values of the Poisson’s ratio, i.e. nonrelaxed and relaxed magnitudes, are the following:

$${\nu }_{\infty }=\underset{t\to 0}{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{\frac{3}{2}{K}_0-{\mu }_{\infty }}{3K_0+{\mu }_{\infty }}},{\nu }_0=\underset{t\to \infty }{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{\frac{3}{2}{K}_0-{\mu }_0}{3K_0+{\mu }_0}}.$$

(39)

Thus, the model (36) allows one to describe both the relaxation and creep processes of viscoelastic materials, and it could be applied for solving different dynamics problems.

2.6. Modelling the Relaxation of the First Lamé parameter Via the Fractional Derivative Standard Linear Solid Model

Now let us consider the case when the first Lamé parameter $\lambda$ is preassigned by the fractional derivative standard linear solid model (6):

$$\tilde{\lambda }={\lambda }_0{\displaystyle \frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}={\lambda }_{\infty }[1-{\nu }_{\epsilon }{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)],$$

(40)

where the nonrelaxed ${\lambda }_{\infty }$ and relaxed ${\lambda }_0$ moduli of the first Lamé parameter are connected with the relaxation ${\tau }_{\epsilon }$ and retardation ${\tau }_{\sigma }$ times similar to (8) by the relationship

$${\left({\displaystyle \frac{{\tau }_{\epsilon }}{{\tau }_{\sigma }}}\right)}^{\gamma }={\displaystyle \frac{{\lambda }_0}{{\lambda }_{\infty }}},$$

and ${\nu }_{\epsilon }=({\lambda }_{\infty }-{\lambda }_0){\lambda }_{\infty }^{-1}$.

For the model under consideration with $\tilde{\lambda }$ and K defined by relationships (40) and (11), respectively, in so doing $K=K_0=K_{\infty }$, the time-dependent Poisson’s operator according to line 2 in Table 1 is expressed as follows:

$$\tilde{\nu }={\displaystyle \frac{\tilde{\lambda }}{3K-\tilde{\lambda }}},$$

(41)

or

$$\begin{array}{c}\tilde{\nu }={\displaystyle \frac{{\lambda }_0\frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}{3K-{\lambda }_0\frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}}={\displaystyle \frac{1}{3K-{\lambda }_{\infty }}}\left\{{\lambda }_{\infty }-{\displaystyle \frac{3K\left({\lambda }_{\infty }-{\lambda }_0\right)}{3K-{\lambda }_0}}{\displaystyle \frac{1}{1+\frac{3K-{\lambda }_{\infty }}{3K-{\lambda }_0}{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}\right\}=\hfill \\ ={\displaystyle \frac{1}{3K-{\lambda }_{\infty }}}\left\{{\lambda }_{\infty }-{\displaystyle \frac{3K\left({\lambda }_{\infty }-{\lambda }_0\right)}{3K-{\lambda }_0}}{϶}_{\gamma }^{*}\left({\displaystyle \frac{3K-{\lambda }_{\infty }}{3K-{\lambda }_0}}{\tau }_{\epsilon }^{\gamma }\right)\right\}=\hfill \\ ={\nu }_{\infty }-\left({\nu }_{\infty }-{\nu }_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{1+{\nu }_0}{1+{\nu }_{\infty }}}{\tau }_{\epsilon }^{\gamma }\right)={\nu }_{\infty }-\left({\nu }_{\infty }-{\nu }_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{{\nu }_0}{{\nu }_{\infty }}}{\tau }_{\sigma }^{\gamma }\right).\hfill \end{array}$$

(42)

Then the time-dependence of the Poisson’s ratio and its nonrelaxed and relaxed magnitudes could be obtained from relationship (42) in the form

$${\nu }_{\infty }=\underset{t\to 0}{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{{\lambda }_{\infty }}{3K-{\lambda }_{\infty }}},{\nu }_0=\underset{t\to \infty }{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{{\lambda }_0}{3K-{\lambda }_0}},$$

(43)

$$\begin{array}{c}\nu \left(t\right)=\tilde{\nu }H\left(t\right)={\displaystyle \frac{1}{3K-{\lambda }_{\infty }}}\left\{{\lambda }_{\infty }-{\displaystyle \frac{3K\left({\lambda }_{\infty }-{\lambda }_0\right)}{3K-{\lambda }_0}}\left\{1-e_{\gamma }\left[{\left({\displaystyle \frac{3K-{\lambda }_0}{3K-{\lambda }_{\infty }}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right]\right\}\right\}\hfill \\ ={\nu }_0+\left({\nu }_{\infty }-{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{1+{\nu }_{\infty }}{1+{\nu }_0}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right]={\nu }_0+\left({\nu }_{\infty }-{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{{\nu }_{\infty }}{{\nu }_0}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right]\hfill \end{array}$$

(44)

Reference to (44) shows that this model describes the behaviour of viscoelastic materials, the Poisson’s ratio of which varies with time from ${\nu }_{\infty }$ to ${\nu }_0$.

2.7. Modelling the P-Wave Modulus $\tilde{M}\equiv \tilde{\lambda }+2\tilde{\mu }$ Via Fractional Derivative Standard Linear Solid Model

One of the most efficient methods of the reconstruction of the material parameters of a linear isotropic viscoelastic structure is from time-dependent measurements of a viscoelastic wave on the surface of a bounded domain of propagation [39]. In this case, the combination of both Lamé parameters, which is called as the modulus of the longitudinal wave, or P-wave, could be preassigned $\tilde{M}=\tilde{\lambda }+2\tilde{\mu }$.

Assume that operator $\tilde{M}$ is given by the fractional derivative standard linear solid model (4):

$$\tilde{M}=M_0{\displaystyle \frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}=M_{\infty }[1-{\nu }_{\epsilon }{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)],$$

(45)

where the nonrelaxed $M_{\infty }$ and relaxed $M_0$ magnitudes of the P-wave modulus are connected with the relaxation ${\tau }_{\epsilon }$ and retardation ${\tau }_{\sigma }$ times similar to (8) by the relationship

$${\left({\displaystyle \frac{{\tau }_{\epsilon }}{{\tau }_{\sigma }}}\right)}^{\gamma }={\displaystyle \frac{M_0}{M_{\infty }}},$$

and ${\nu }_{\epsilon }=(M_{\infty }-M_0)M_{\infty }^{-1}$.

For the model under consideration with $\tilde{M}$ and K defined by relationships (47) and (11), respectively, the time-dependent Poisson’s operator according to the last line in Table 1 is expressed as follows:

$$\nu ={\displaystyle \frac{3K-\tilde{M}}{3K+\tilde{M}}},$$

(46)

or

$$\begin{array}{c}\tilde{\nu }={\displaystyle \frac{3K-M_0\frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}{3K+M_0\frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}}=\hfill \\ ={\displaystyle \frac{1}{\left(3K+M_0\right)\left(3K+M_{\infty }\right)}}\left\{9K^2-3K\left(M_{\infty }-M_0\right)-M_0{M}_{\infty }+6K\left(M_{\infty }-M_0\right){\displaystyle \frac{1}{1+\frac{3K+M_{\infty }}{3K+M_0}{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}\right\}\hfill \\ ={\displaystyle \frac{1}{\left(3K+M_0\right)\left(3K+M_{\infty }\right)}}\left\{9K^2-3K\left(M_{\infty }-M_0\right)-M_0{M}_{\infty }+6K\left(M_{\infty }-M_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{3K+M_{\infty }}{3K+M_0}}{\tau }_{\epsilon }^{\gamma }\right)\right\}\hfill \\ ={\nu }_{\infty }-\left({\nu }_{\infty }-{\nu }_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{{\nu }_0+1}{{\nu }_{\infty }+1}}{\tau }_{\epsilon }^{\gamma }\right)={\nu }_{\infty }-\left({\nu }_{\infty }-{\nu }_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{1-{\nu }_0}{1-{\nu }_{\infty }}}{\tau }_{\sigma }^{\gamma }\right)\hfill \end{array}$$

(47)

Then the time-dependence of the Poisson’s ratio and its nonrelaxed and relaxed magnitudes could be obtained from relationship (47) in the form

$${\nu }_{\infty }=\underset{t\to 0}{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{3K-M_{\infty }}{3K+M_{\infty }}},{\nu }_0=\underset{t\to \infty }{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{3K-M_0}{3K+M_0}},$$

(48)

$$\begin{array}{c}\nu \left(t\right)=\tilde{\nu }H\left(t\right)={\displaystyle \frac{1}{\left(3K+M_0\right)\left(3K+M_{\infty }\right)}}\hfill \\ \times \left\{9K^2-3K\left(M_{\infty }-M_0\right)-M_0{M}_{\infty }+6K\left(M_{\infty }-M_0\right)\left\{1-e_{\gamma }\left[{\left({\displaystyle \frac{3K+M_0}{3K+M_{\infty }}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right]\right\}\right\}\hfill \\ ={\nu }_0+\left({\nu }_{\infty }-{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{{\nu }_{\infty }+1}{{\nu }_0+1}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right]={\nu }_0+\left({\nu }_{\infty }-{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{1-{\nu }_{\infty }}{1-{\nu }_0}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right].\hfill \end{array}$$

(49)

From (49) it is seen that similar to the cases, when each of the Lamé parameters is defined separately by the fractional derivative standard linear solid, modelling the P-wave modulus allows one to describe the behaviour of viscoelastic materials, Poisson’s ratio of which varies with time from its nonrelaxed value to its relaxed value.

3. Analysis of the Fractional Derivative Models of Viscoelasticity Involving Time-Dependent Poisson’s Ratio and Without Volume Relaxation

From Table 1 it is seen that for viscoelastic materials, the volume relaxation of which could be neglected, i.e., bulk modulus could be considered as a constant value $K=\mathrm{const}$, as a second time-dependent operator which should be given together with $K=\mathrm{const}$ one of the four material characteristics could be preassigned: E, $\lambda$, $\mu$, or $M=\lambda +2\mu$. In its turn, each of these operators could be described by four fractional derivative models: Scott Blair, Kelvin-Voigt, Maxwell, or standard linear solid models. Thus, there could be 16 different variants of viscoelastic models involving time-dependent Poisson’s ratio, seven of which have been considered above in Section 2, and some models are presented in [13,31,32,33,34].

The limiting values of the time-dependent Poisson’s ratio are summarized in Table 2Table 3 for all 16 models, what will allow one to classify the models constructed.

Reference to Table 2 and Table 3 shows that all models could be divided into three groups:

1.

models describing the behaviour of ’traditional’ viscoelastic isotropic materials, i.e., materials with positive magnitudes of Poisson’s ratio within the thermodynamically admissible range $0<\nu \le 1/2$ – models No. 3, 4, 7, 8, 11, 12, 16;

2.

models describing the behaviour of isotropic viscoelastic materials with negative Poisson’s ratios within the thermodynamically admissible range $-1\le \nu \le 1/2$ – models No. 5, 6, 9, 10, 14;

3.

physically meaningless models, i.e., models with Poisson’s ratios lying without the thermodynamically admissible domain $-1\le \gamma \le 1/2$ either from the left with $\nu <-1$ – models 1 and 2, or from the right with $\nu >1/2$ – models 13 and 15.

The fractional derivative standard linear solid model, which could describe both the relaxation and creep phenomena occurring during deformation of viscoelastic materials, provides the variation of Poisson’s ratio with time from its nonrelaxed (instantaneous) magnitude ${\nu }_{\infty }$ to its relaxed (prolonged) magnitude ${\nu }_0$ according to the relationship similar in the form to all models:

$$\nu \left(t\right)={\nu }_0-\left({\nu }_0-{\nu }_{\infty }\right)e_{\gamma }\left[{\left({\displaystyle \frac{{\nu }_{1+\infty }}{1+{\nu }_0}}\right)}^{1/\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right],$$

in so doing the limiting values of Poisson’s ratios are calculated for each model individually.

The time-dependence of the Poisson’s ratio for fractional derivative standard linear solid models for Lamé parameters: model 8 for $\mu$ (26), model 12 for $\lambda$ (44), and model 16 for $M=\lambda +2\mu$ (49) are presented in Figure 2 for ${\nu }_{\infty }=0.25$ and ${\nu }_0=0.3$ at different values of the fractional parameter $\gamma =0.4,0.6,0.8$, and 1, whence it is evident that all curves increase monotonically from 0.25 to 0.3, and the curve corresponding to $\gamma =1$ approaches the upper limit more rapidly than the curves for fractional magnitudes of $\gamma$.

The Maxwell models for Young’s operator (model 3) and for shear operator (model 7) behave in a similar way: Poisson’s ratio varies with time from its nonrelaxed magnitude ${\nu }_{\infty }$ to 1/2, as it is shown in Figure 3.

Scott Blair models for the shear operator (model 5) and the first Lamé operator (model 9), Kelvin-Voigt models for the shear operator (model 6), the first Lamé operator (model 10), and the P-wave operator (model 14) describe the behaviour of viscoelastic auxetics with the lower limit of the Poisson’s ratio equal to ${\left.\nu \left(t\right)\right|}_{t\to 0}=-1$. As for the upper limiting magnitude of the Poisson’s ratio, then for the model 5 it is 1/2, for the model 9 it is zero, and for the models 6, 10, and 14, it is a finite value within the range of $0<\nu <1/2$, which depends on the given magnitudes of the preassigned operators. Time-dependence of the Poisson’s operator for the model 14 is shown in Figure 4.

The Maxwell model for the first Lamé operator (model 11) is of particular interest, since in this case $\nu \left(t\right)$ either increases monotonically from a certain negative magnitude up to zero (see Figure 5), or decreases monotonically from a certain positive magnitude to zero (see Figure 6).

As for physically meaningless models, then in models 1 and 2 $\nu \left(t\right)$ monotonically increases from $-\infty$ to 1/2 and a certain positive magnitude (see Figure 7), respectively, or monotonically increases from -1 to 1 in model 13 and from a certain positive magnitude to 1 in model 15 (see Figure 8).

4. Conclusion

In the present paper, several viscoelastic models are studied for the cases when time-dependent viscoelastic operators are represented in terms of the fractional derivative Scott Blair, Kelvin-Voigt, Maxwell, and standard linear solid models. Using the algebra of dimensionless Rabotnov’s fractional exponential functions, time-dependent operators for Poisson’s ratios have been obtained and analyzed. It is shown that materials described by some of such models are viscoelastic auxetics, because Poisson’s ratios of such materials are time-dependent operators which could take on both positive and negative magnitudes.

In the companion paper, it is planned to analyze viscoelastic materials involving time-dependent operators, including the Poisson’s operator with due account for bulk relaxation.