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

Abstract

In the present paper, several viscoelastic models are studied for cases when time-dependent viscoelastic operators of Lamé’s parameters are represented in terms of the fractional derivative Kelvin–Voigt, Scott-Blair, Maxwell, and standard linear solid models. This is practically important since precisely these parameters define the velocities of longitudinal and transverse waves propagating in three-dimensional media. Using the algebra of dimensionless Rabotnov’s fractional exponential functions, time-dependent operators for Poisson’s ratios have been obtained and analysed. It is shown that materials described by some of such models are viscoelastic auxetics because the Poisson’s ratios of such materials are time-dependent operators which could take on both positive and negative magnitudes.

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 the wide application of the theory of elasticity in studies of advanced and traditional materials, much attention is paid to the 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 (whose theory has been presented in numerous monographs, for example, in [6,7,8,9,10,11]) has gained wide acceptance in the modelling of viscoelastic bodies such as beams, plates, and shells [12,13,14,15,16]. Their damping features are described most often by defining the Young’s operator by the simplest fractional derivative models, namely: the Kelvin–Voigt model, Maxwell model, and standard linear solid model [17,18,19]. As this happens, the Poisson ratio of a viscoelastic material is often assumed to be constant [20,21]. However, it has been emphasized in [22] 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 [12,23,24] or on a viscoelastic foundation [25].

As experimental data have shown [26], the Poisson’s ratio is always a time-dependent operator [27,28,29], 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 [12,19,22,30]. 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 [31,32,33,34,35].

For solving different dynamic problems of the mechanics of viscoelastic solids and structures, it is essential to know the form of viscoelastic operators entering the governing equations. For example, Poisson’s ratio $\nu$ and Young’s modulus E are involved in the cylindrical rigidity of plates and shells $D={\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}}$, where h is the thickness of a plate or shell. Poisson’s ratios and Young’s moduli of colliding bodies define the Hertz law for the solution of the problem of viscoelastic contact interaction during impact [36,37,38]. Therefore, their pinpointing is of paramount importance in the impact response analysis.

It is well known [39,40] 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. The interrelationship of elastic material constants: bulk modulus K, Young’s modulus E, Lamé’s parameters $\lambda$ and $\mu$, Poisson’s ratio $\nu$, and P-wave modulus $M=\lambda +2\mu$.

Material Constants	Young’s Modulus E	1st Lamé’s Parameter $\mathit{\lambda }$	Shear Modulus $\mathit{\mu }$	Poisson’s Ratio $\mathit{\nu }$	P-Wave Modulus M
$(K,E)$	-	${\displaystyle \frac{3K(3K-E)}{9K-E}}$	${\displaystyle \frac{3KE}{9K-E}}$	${\displaystyle \frac{3K-E}{6K}}$	${\displaystyle \frac{3K(3K+E)}{9K-E}}$
$(K,\lambda )$	${\displaystyle \frac{9K(K-\lambda )}{3K-\lambda }}$	-	${\displaystyle \frac{3(K-\lambda )}{2}}$	${\displaystyle \frac{\lambda }{3K-\lambda }}$	$3K-2\lambda$
$(K,\mu )$	${\displaystyle \frac{9K\mu }{3K+\mu }}$	$K-{\displaystyle \frac{2\mu }{3}}$	-	${\displaystyle \frac{3K-2\mu }{2(3K+\mu )}}$	$K+{\displaystyle \frac{4\mu }{3}}$
$(K,\nu )$	$3K(1-2\nu )$	${\displaystyle \frac{3K\nu }{1+\nu }}$	${\displaystyle \frac{3K(1-2\nu )}{2(1+\nu )}}$	-	${\displaystyle \frac{3K(1-\nu )}{1+\nu }}$
$(K,M)$	${\displaystyle \frac{9K(M-K)}{3K+M}}$	${\displaystyle \frac{3K-M}{2}}$	${\displaystyle \frac{3(M-K)}{4}}$	${\displaystyle \frac{3K-M}{3K+M}}$	-

, 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.

In [39,40], it was shown that all models of isotropic elastic materials should satisfy the following thermodynamic constraints:

$$K>0,\mu >0,-1\le \nu \le {\displaystyle \frac{1}{2}}.$$

(1)

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 necessary to know the Young’s and shear operators defining phase velocities and coefficients of attenuation of viscoelastic longitudinal ($c_l=\sqrt{E/\rho }$, where $\rho$ is the material density) and shear ($c_t=\sqrt{\mu /\rho }$) waves. The longitudinal wave propagates in plates with velocity $c_l=\sqrt{E/\rho (1-{\nu }^2)}$; thus, in this case, it is also needed to define the time-dependent Poisson’s operator. Therefore, 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 Kelvin–Voigt model [22,41,42,43];
• Young’s modulus is relaxed according to the fractional derivative Maxwell model [43];
• Young’s modulus is relaxed according to the fractional derivative standard linear solid model [22,41,43,44];
• Shear modulus is defined by the Scott-Blair model [22];
• Shear modulus is preassigned by the fractional derivative Kelvin–Voigt model [22,41,42,43,45];
• Shear modulus is relaxed according to the Maxwell model [22,43];
• Shear modulus is relaxed according to the standard linear solid model [22,43,46,47].

It is interesting to note that in contrast with 1D and 2D dynamic problems for rods and plates, 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$ [40]:

$${\sigma }_{ij}=\lambda {\delta }_{ij}{\epsilon }_{ll}+2\mu {\epsilon }_{ij},$$

(2)

or

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

(3)

where ${\sigma }_{ij}$ and ${\epsilon }_{ij}$ are the components of the stress and strain tensors, respectively, and ${\delta }_{ij}$ is the Kronecker tensor.

This indicates that K, $\lambda$, and $\mu$, or K and $M=\lambda +2\mu$, are the most intrinsic moduli to express stress in terms of strain when studying wave propagation in 3D elastic media. Therefore, it is extremely important to define their time-dependent analogues for viscoelastic media [48,49,50].

The parameter $M=\lambda +2\mu$ is called the modulus of the P-wave (primary wave, or compressional wave, or dilatational wave), since it governs the velocity of its propagation as $c_l=\sqrt{(\lambda +2\mu )/\rho }$. The velocity of propagation of the S-wave (secondary wave, rotational, equivoluminal, distortional, or shear wave) through the elastic medium is given by $c_t=\sqrt{\mu /\rho }$ [51].

That is why in the present article, the emphasis will be on a complete analysis of time-dependent operators for Lamé parameters, which would allow us to study the propagation of dilatational and distortional waves through viscoelastic media, the properties of which are described by various fractional derivative models. Using the procedure for the mathematical treatment of viscoelastic operators suggested in Rossikhin and Shitikova [22,44] below for the first time, the models based on the application of fractional derivatives will be studied for the cases where the first Lamé parameter or the P-wave modulus is given a priori, without considering 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 [19,22]:

• Scott-Blair model (Figure 1A)

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

  (4)
• Maxwell model (Figure 1B)

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

  (5)

  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}$$

  (6)
• 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),$$

  (7)
• Standard linear solid model (Figure 1D,E)

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

(8)

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)],$$

(9)

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 [6] 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 },$$

(10)

$\Gamma \left(\gamma \right)$ is the Gamma-function,

$${϶}_{\gamma }^{*}\left({\tau }_i^{\gamma }\right)={\displaystyle \frac{1}{1+{\tau }_i^{\gamma }{D}^{\gamma }}}(i=\epsilon \mathrm{or}\sigma )$$

(11)

is the dimensionless Rabotnov’s operator [18,44], 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$ becomes a conventional exponential function.

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

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

(12)

Equations (4)–(9) govern stress–strain relationships for one-dimensional viscoelastic media. The constitutive equations connecting the strain and the stress in a three-dimensional linear viscoelastic isotropic medium could be written using the corresponding equations for the elastic medium (2) and (3) and the Volterra correspondence principle in the following form:

$${\sigma }_{ij}=\tilde{\lambda }{\delta }_{ij}{\epsilon }_{ll}+2\tilde{\mu }{\epsilon }_{ij},$$

(13)

or

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

(14)

where $\tilde{\lambda }$, $\tilde{\mu }$, and $\tilde{K}$ are, respectively, hereditary Lamé and bulk operators.

Thus, for solving three-dimensional dynamic viscoelasticity problems, the form of two hereditary operators is needed: $\tilde{\lambda }$ and $\tilde{\mu }$, or $\tilde{K}$ and $\tilde{\mu }$.

Below we will consider models involving a time-dependent Poisson’s ratio but without volume relaxation, i.e., when the 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},$$

(15)

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, that is, one of the Lamé operators $\tilde{\lambda }$ or $\tilde{\mu }$, or P-wave operator $\tilde{M}=\tilde{\lambda }+2\tilde{\mu }$, which could be done using the fractional derivative Maxwell, Scott-Blair, Kelvin–Voigt, or standard linear solid models (4)–(9). Knowing the form of two viscoelastic operators, it is possible to define the form of the hereditary 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 (7) as

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

(16)

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

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

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

(17)

Now, substituting the shear operator $\tilde{\mu }$ (16) in (17), we obtain

$$\begin{array}{c}\tilde{E}={\displaystyle \frac{9K_0{\mu }_0\left(1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }\right)}{3K_0+{\mu }_0\left(1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }\right)}}=9K_0{\displaystyle \frac{1+\frac{{\mu }_0}{3K_0+{\mu }_0}{\tau }_{\sigma }^{\gamma }{D}^{\gamma }-\frac{3K_0}{3K_0+{\mu }_0}}{1+\frac{{\mu }_0}{3K_0+{\mu }_0}{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =9K_0\left[1-{\displaystyle \frac{3K_0}{\left(3K_0+{\mu }_0\right)\left(1+\frac{{\mu }_0}{3K_0+{\mu }_0}{\tau }_{\sigma }^{\gamma }{D}^{\gamma }\right)}}\right].\end{array}$$

(18)

Using the definition of the dimensionless Rabotnov’s operator (11), the Young’s operator (18) could be reduced to the following form:

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

(19)

where $E_0={\displaystyle \frac{9K_0{\mu }_0}{3K_0+{\mu }_0}}$ is the relaxed Young’s modulus.

Now, the Poisson’s operator could be calculated either via the formula [22]:

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

(20)

or by the formula from the first line of Table 1:

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

(21)

Substituting (19) in (21), we find

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

(22)

where ${\nu }_0={\displaystyle \frac{\frac{3}{2}{K}_0-{\mu }_0}{3K_0+{\mu }_0}}={\displaystyle \frac{E_0}{2{\mu }_0}}-1$ is the relaxed Poisson’s ratio.

In order to investigate the time dependence of the Poisson’s ratio during the relaxation process, it is necessary to apply the Poisson’s operator $\tilde{\nu }$ (22) to the unit Heaviside function $H\left(t\right)$

$$\begin{array}{c}\hfill \nu \left(t\right)=\tilde{\nu }H\left(t\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],\end{array}$$

(23)

where $e_{\gamma }\left(t\right)\equiv E_{\gamma }\left(-t^{\gamma }\right)$[50], 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 [7,10].

From the relationship (23), the limiting magnitudes of the Poisson’s ratio 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{\frac{3}{2}{K}_0-{\mu }_0}{3K_0+{\mu }_0}}={\nu }_0.$$

(24)

From (24), it follows that the model (16) with $K_0=const$ could describe the behaviour of isotropic viscoelastic auxetics (materials with negative Poisson ratios), and so the Poisson’s ratio varies from −1 to its relaxed magnitude ${\nu }_0$, which does not violate the laws of thermodynamics [39,40], as follows from the relationship (1).

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

If the fractional derivative Maxwell model (6) is applied to describe 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],$$

(25)

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

Using the procedure described above for the Kelvin–Voigt model in Section 2.2, we could similarly obtain for the Maxwell model

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

(26)

$$\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],$$

(27)

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{3/\phantom{32}2}{{\nu }_{\infty }+1}}{\tau }_{\epsilon }^{\gamma }\right),$$

(28)

where ${\nu }_{\infty }={\displaystyle \frac{\frac{3}{2}{K}_{\infty }-{\mu }_{\infty }}{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 \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{1}{2}}-\left({\displaystyle \frac{1}{2}}-{\nu }_{\infty }\right)e_{\gamma }\left[{\left({\displaystyle \frac{{\nu }_{\infty }+1}{3/\phantom{32}2}}\right)}^{1/\phantom{1\gamma }\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right],\end{array}$$

(29)

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 }}}={\nu }_{\infty },\underset{t\to \infty }{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{1}{2}}.$$

(30)

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

2.4. Scott-Blair Model for Shear Relaxation

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

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

(31)

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).$$

(32)

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),$$

(33)

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),$$

(34)

$$\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}}.$$

(35)

From (35), it is evident that according to this model, the Poisson’s ratio could vary in a very broad range, namely, from −1 to 0.5. Thus, this model is thermodynamically admissible for viscoelastic auxetics since all magnitudes of the Poisson’s ratio fall within the domain bounded by the thermodynamic constraints (1).

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

If the fractional derivative standard linear solid model (9) is applied to describe 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 }}},$$

(36)

or

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

(37)

where ${\nu }_{\epsilon }=({\mu }_{\infty }-{\mu }_0){\mu }_{\infty }^{-1}$, and ${\mu }_{\infty }$ is defined similarly to Relationship (12)

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

(38)

In so doing, in this model, operator $\tilde{K}$ is still defined by (15). Therefore, in this case, the 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)}}\{9{K_0}^2-3K_0\left(2{\mu }_{\infty }-{\mu }_0\right)\hfill \\ -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)\}\\ \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}$$

(39)

For the model under consideration, the Poisson’s ratio varies with time according to the following relationship:

$$\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}\\ \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}$$

(40)

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

$$\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 }_{\infty },\underset{t\to \infty }{lim}\tilde{\nu }H\left(t\right)={\displaystyle \frac{\frac{3}{2}{K}_0-{\mu }_0}{3K_0+{\mu }_0}}={\nu }_0.$$

(41)

Thus, Model (37) allows one to describe both the relaxation and creep processes of viscoelastic materials, and it could be applied for solving different dynamics problems, the governing equations of which involve time-dependent viscoelastic operators, including a viscoelastic Poisson’s ratio, for example, in the problems of viscoelastic contact and/or impact interaction [36,37,38].

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 (9):

$$\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)],$$

(42)

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 (12) 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 (42) and (15), respectively, and 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 }}},$$

(43)

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}$$

(44)

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

$$\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}$$

(45)

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

(46)

Reference to (46) 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 for 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 [49]. In this case, the combination of both Lamé parameters, which is called 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 (9):

$$\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)],$$

(47)

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 (12) 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 (15), respectively, the time-dependent Poisson’s operator according to the last line in Table 1 is expressed as follows:

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

(48)

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}$$

(49)

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

$$\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}$$

(50)

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

(51)

From (51), it is seen that similar to the cases when each of the Lamé parameters is defined separately by the fractional derivative standard linear solid model, modelling the P-wave modulus allows one to describe the behaviour of viscoelastic materials, the Poisson’s ratio of which varies over 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., the 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 a time-dependent Poisson’s ratio, seven of which have been considered above in Section 2, and some models are presented in [22,41,42,43,44].

The limiting values of the time-dependent Poisson’s ratio are summarized in

Table 2. Limiting values of the Poisson’s ratio $\nu \left(t\right)=\tilde{\nu }H\left(t\right)$ for isotropic materials, the viscoelastic features of which are described by different fractional-order operators without bulk relaxation (Part 1).

Type of the Model Involving Fractional Derivatives	${\left.\mathit{\nu }\left(\mathit{t}\right)\right|}_{\mathit{t}\to 0}$	${\left.\mathit{\nu }\left(\mathit{t}\right)\right|}_{\mathit{t}\to \mathit{\infty }}$
A. Modelling the Young’s operator $\tilde{E}$
(1) Scott-Blair model
$\tilde{E}=E{\tau }^{\gamma }{D}^{\gamma }$
$\tilde{\nu }={\displaystyle \frac{1}{2}}-{\displaystyle \frac{E}{6K}}{\tau }^{\gamma }{D}^{\gamma }$
$\nu \left(t\right)={\displaystyle \frac{1}{2}}-{\displaystyle \frac{E}{6K}}{\displaystyle \frac{{\left(t/\tau \right)}^{-\gamma }}{\Gamma (1-\gamma )}}$	$-\infty$	${\displaystyle \frac{1}{2}}$
(2) Kelvin–Voigt model
$\tilde{E}=E_0\left(1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }\right)$
$\tilde{\nu }={\nu }_0-{\displaystyle \frac{E_0}{6K_0}}{\tau }_{\sigma }^{\gamma }{D}^{\gamma }$
$\nu \left(t\right)={\nu }_0-{\displaystyle \frac{E_0}{6K_0}}{\displaystyle \frac{{\left(t/{\tau }_{\sigma }\right)}^{-\gamma }}{\Gamma (1-\gamma )}}$	$-\infty$	${\nu }_0={\displaystyle \frac{1}{2}}-{\displaystyle \frac{E_0}{6K_0}}$
(3) Maxwell model
$\tilde{E}=E_{\infty }{\displaystyle \frac{{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}=E_{\infty }\left[1-{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right]$
$\tilde{\nu }={\displaystyle \frac{1}{2}}-{\displaystyle \frac{E_{\infty }}{6K_{\infty }}}\left[1-{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right]$
$\nu \left(t\right)={\displaystyle \frac{1}{2}}-{\displaystyle \frac{E_{\infty }}{6K_{\infty }}}{e}_{\gamma }\left({\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right)$	${\nu }_{\infty }={\displaystyle \frac{1}{2}}-{\displaystyle \frac{E_{\infty }}{6K_{\infty }}}$	${\displaystyle \frac{1}{2}}$
(4) Standard linear solid model
$\tilde{E}=E_0{\displaystyle \frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}=E_0\left[{\displaystyle \frac{{\tau }_{\sigma }^{\gamma }}{{\tau }_{\epsilon }^{\gamma }}}-\left({\displaystyle \frac{{\tau }_{\sigma }^{\gamma }}{{\tau }_{\epsilon }^{\gamma }}}-1\right){϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right]$
$\tilde{\nu }={\nu }_{\infty }-\left({\nu }_{\infty }-{\nu }_0\right){϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)$
$\nu \left(t\right)={\nu }_0+\left({\nu }_{\infty }-{\nu }_0\right)e_{\gamma }\left({\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right)$	${\nu }_{\infty }={\displaystyle \frac{1}{2}}-{\displaystyle \frac{E_{\infty }}{6K}}$	${\nu }_0={\displaystyle \frac{1}{2}}-{\displaystyle \frac{E_0}{6K}}$
B. Modelling the shear operator $\tilde{\mu }$ (second Lamé parameter)
(5) Scott-Blair model
$\tilde{\mu }=\mu {\tau }^{\gamma }{D}^{\gamma }$
$\tilde{\nu }=-1+{\displaystyle \frac{3}{2}}{϶}_{\gamma }^{*}\left({\displaystyle \frac{\mu }{3K}}{\tau }^{\gamma }\right)$
$\nu \left(t\right)=-1+{\displaystyle \frac{3}{2}}\left\{1-e_{\gamma }\left[{\left({\displaystyle \frac{3K}{\mu }}\right)}^{1/\gamma }{\displaystyle \frac{t}{\tau }}\right]\right\}$	$-1$	${\displaystyle \frac{1}{2}}$
(6) Kelvin–Voigt model
$\tilde{\mu }={\mu }_0\left(1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }\right)$
$\tilde{\nu }=-1+\left(1+{\nu }_0\right){϶}_{\gamma }^{*}\left[{\displaystyle \frac{{\mu }_0}{3K_0+{\mu }_0}}{\tau }_{\sigma }^{\gamma }\right]=-1+\left(1+{\nu }_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{1-2{\nu }_0}{3}}{\tau }_{\sigma }^{\gamma }\right)$
$\nu \left(t\right)={\nu }_0-\left(1+{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{3K_0}{{\mu }_0}}+1\right)}^{1/\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right]={\nu }_0-\left(1+{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{3}{1-2{\nu }_0}}\right)}^{1/\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right]$	$-1$	${\nu }_0={\displaystyle \frac{\frac{3}{2}{K}_0-{\mu }_0}{3K_0+{\mu }_0}}$
(7) Maxwell model
$\tilde{\mu }={\mu }_{\infty }{\displaystyle \frac{{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}={\mu }_{\infty }\left[1-{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right]$
$\tilde{\nu }={\nu }_{\infty }+\left({\displaystyle \frac{1}{2}}-{\nu }_{\infty }\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{3K_{\infty }+{\mu }_{\infty }}{3K_{\infty }}}{\tau }_{\epsilon }^{\gamma }\right)={\nu }_{\infty }+\left({\displaystyle \frac{1}{2}}-{\nu }_{\infty }\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{3/2}{{\nu }_{\infty }+1}}{\tau }_{\epsilon }^{\gamma }\right)$
$\nu \left(t\right)={\displaystyle \frac{1}{2}}-\left({\displaystyle \frac{1}{2}}-{\nu }_{\infty }\right)e_{\gamma }\left[{\left({\displaystyle \frac{3K_{\infty }}{3K_{\infty }+{\mu }_{\infty }}}\right)}^{1/\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right]={\displaystyle \frac{1}{2}}-\left({\displaystyle \frac{1}{2}}-{\nu }_{\infty }\right)e_{\gamma }\left[{\left({\displaystyle \frac{{\nu }_{\infty }+1}{3/2}}\right)}^{1/\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right]$	${\nu }_{\infty }={\displaystyle \frac{\frac{3}{2}{K}_{\infty }-{\mu }_{\infty }}{3K_{\infty }+{\mu }_{\infty }}}$	${\displaystyle \frac{1}{2}}$
(8) Standard linear solid model
$\tilde{\mu }={\mu }_0{\displaystyle \frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}={\mu }_0\left[{\displaystyle \frac{{\tau }_{\sigma }^{\gamma }}{{\tau }_{\epsilon }^{\gamma }}}-\left({\displaystyle \frac{{\tau }_{\sigma }^{\gamma }}{{\tau }_{\epsilon }^{\gamma }}}-1\right){϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right]$
$\tilde{\nu }={\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)$
$\nu \left(t\right)={\nu }_0+\left({\nu }_{\infty }-{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{{\nu }_{\infty }+1}{{\nu }_0+1}}\right)}^{1/\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/\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right]$	${\nu }_{\infty }={\displaystyle \frac{\frac{3}{2}K-{\mu }_{\infty }}{3K+{\mu }_{\infty }}}$	${\nu }_0={\displaystyle \frac{\frac{3}{2}K-{\mu }_0}{3K+{\mu }_0}}$

and

Table 3. Limiting values of the Poisson’s ratio $\nu \left(t\right)=\tilde{\nu }H\left(t\right)$ for isotropic materials, the viscoelastic features of which are described by different fractional-order operators without bulk relaxation (Part 2).

Type of the Model Involving Fractional Derivatives	${\left.\mathit{\nu }\left(\mathit{t}\right)\right|}_{\mathit{t}\to 0}$	${\left.\mathit{\nu }\left(\mathit{t}\right)\right|}_{\mathit{t}\to \mathit{\infty }}$
C. Modelling of the first Lamé parameter $\tilde{\lambda }$
(9) Scott-Blair model
$\tilde{\lambda }=\lambda {\tau }^{\gamma }{D}^{\gamma }$
$\tilde{\nu }=-1+{϶}_{\gamma }^{*}\left(-{\displaystyle \frac{\lambda }{3K}}{\tau }^{\gamma }\right)$
$\nu \left(t\right)=-e_{\gamma }\left[{\left(-{\displaystyle \frac{3K}{\lambda }}\right)}^{1/\gamma }{\displaystyle \frac{t}{\tau }}\right]$	$-1$	0
(10) Kelvin–Voigt model
$\tilde{\lambda }={\lambda }_0\left(1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }\right)$
$\tilde{\nu }=-1+\left(1+{\nu }_0\right){϶}_{\gamma }^{*}\left(-{\nu }_0{\tau }_{\sigma }^{\gamma }\right)$
$\nu \left(t\right)={\nu }_0-\left(1+{\nu }_0\right)e_{\gamma }\left[{\left(-{\displaystyle \frac{1}{{\nu }_0}}\right)}^{1/\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right]$	$-1$	${\nu }_0={\displaystyle \frac{{\lambda }_0}{3K_0-{\lambda }_0}}$
(11) Maxwell model
$\tilde{\lambda }={\lambda }_{\infty }\left[1-{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right]$
$\tilde{\nu }={\nu }_{\infty }\left[1-{϶}_{\gamma }^{*}\left({\displaystyle \frac{1}{{\nu }_{\infty }+1}}{\tau }_{\epsilon }^{\gamma }\right)\right]$
$\nu \left(t\right)={\nu }_{\infty }{e}_{\gamma }\left[{\left(1+{\nu }_{\infty }\right)}^{1/\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right]$	${\nu }_{\infty }={\displaystyle \frac{{\lambda }_{\infty }}{3K_{\infty }-{\lambda }_{\infty }}}$	0
(12) Standard linear solid model
$\tilde{\lambda }={\lambda }_0{\displaystyle \frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}={\lambda }_0\left[{\displaystyle \frac{{\tau }_{\sigma }^{\gamma }}{{\tau }_{\epsilon }^{\gamma }}}-\left({\displaystyle \frac{{\tau }_{\sigma }^{\gamma }}{{\tau }_{\epsilon }^{\gamma }}}-1\right){϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right]$
$\tilde{\nu }={\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{{\nu }_0}{{\nu }_{\infty }}}{\tau }_{\sigma }^{\gamma }\right)$
$\nu \left(t\right)={\nu }_0+\left({\nu }_{\infty }-{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{{\nu }_{\infty }+1}{{\nu }_0+1}}\right)}^{1/\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/\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right]$	${\nu }_{\infty }={\displaystyle \frac{{\lambda }_{\infty }}{3K-{\lambda }_{\infty }}}$	${\nu }_0={\displaystyle \frac{{\lambda }_0}{3K-{\lambda }_0}}$
D. Modelling the P-wave modulus operator $\tilde{M}$
(13) Scott-Blair model
$\tilde{M}=M{\tau }^{\gamma }{D}^{\gamma }$
$\tilde{\nu }=-1+2{϶}_{\gamma }^{*}\left(-{\displaystyle \frac{M}{3K}}{\tau }^{\gamma }\right)$
$\nu \left(t\right)=-1+2\left\{1-e_{\gamma }\left[{\left({\displaystyle \frac{3K}{M}}\right)}^{1/\gamma }{\displaystyle \frac{t}{\tau }}\right]\right\}$	$-1$	1
(14) Kelvin–Voigt model
$\tilde{M}=M_0\left(1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }\right)$
$\tilde{\nu }=-1+\left(1+{\nu }_0\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{1-{\nu }_0}{2}}{\tau }_{\sigma }^{\gamma }\right)$
$\nu \left(t\right)={\nu }_0-\left(1+{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{2}{1-{\nu }_0}}\right)}^{1/\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right]$	$-1$	${\nu }_0={\displaystyle \frac{3K_0-M_0}{3K_0+M_0}}$
(15) Maxwell model
$\tilde{M}=M_{\infty }\left[1-{϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right]$
$\tilde{\nu }={\nu }_{\infty }+\left(1-{\nu }_{\infty }\right){϶}_{\gamma }^{*}\left({\displaystyle \frac{2}{{\nu }_{\infty }+1}}{\tau }_{\epsilon }^{\gamma }\right)$
$\nu \left(t\right)=1-\left(1-{\nu }_{\infty }\right)e_{\gamma }\left[{\left({\displaystyle \frac{{\nu }_{\infty }+1}{2}}\right)}^{1/\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right]$	${\nu }_{\infty }={\displaystyle \frac{3K_{\infty }-M_{\infty }}{3K_{\infty }+M_{\infty }}}$	1
(16) Standard linear solid model
$\tilde{M}=M_0{\displaystyle \frac{1+{\tau }_{\sigma }^{\gamma }{D}^{\gamma }}{1+{\tau }_{\epsilon }^{\gamma }{D}^{\gamma }}}=M_0\left[{\displaystyle \frac{{\tau }_{\sigma }^{\gamma }}{{\tau }_{\epsilon }^{\gamma }}}-\left({\displaystyle \frac{{\tau }_{\sigma }^{\gamma }}{{\tau }_{\epsilon }^{\gamma }}}-1\right){϶}_{\gamma }^{*}\left({\tau }_{\epsilon }^{\gamma }\right)\right]$
$\tilde{\nu }={\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)$
$\nu \left(t\right)={\nu }_0+\left({\nu }_{\infty }-{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{1+{\nu }_{\infty }}{1+{\nu }_0}}\right)}^{1/\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/\gamma }{\displaystyle \frac{t}{{\tau }_{\sigma }}}\right]$	${\nu }_{\infty }={\displaystyle \frac{3K-M_{\infty }}{3K+M_{\infty }}}$	${\nu }_0={\displaystyle \frac{3K-M_0}{3K+M_0}}$

for all 16 models, which 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:

• 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;
• 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;
• Physically meaningless models, i.e., models with Poisson’s ratios lying outside of the thermodynamically admissible domain $-1\le \nu \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 in the 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 }_{\infty }-{\nu }_0\right)e_{\gamma }\left[{\left({\displaystyle \frac{1+{\nu }_{\infty }}{1+{\nu }_0}}\right)}^{1/\gamma }{\displaystyle \frac{t}{{\tau }_{\epsilon }}}\right],$$

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

The time dependences of the Poisson ratio for the fractional derivative standard linear solid model for Lamé parameters: model 8 for $\mu$ (40), model 12 for $\lambda$ (45), and model 16 for $M=\lambda +2\mu$ (50), 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. From Figure 2, it is evident that all the 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 at fractional magnitudes of $\gamma$.

The Maxwell models for Young’s operator (model 3) and for the shear operator (model 7) behave in a similar way: the Poisson’s ratio varies with time from its unrelaxed magnitude ${\nu }_{\infty }$ to one-half, as 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 ratio equal to ${\left.\nu \left(t\right)\right|}_{t\to 0}=-1$. As for the upper limiting magnitude of the Poisson ratio, then for model 5 it is one-half, for model 9 it is zero, and for 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. The time dependence of the Poisson operator for 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).

For physically meaningless models, it is evident from the data presented in Table 2 and Table 3 that in models 1 and 2, $\nu \left(t\right)$ monotonically increases from $-\infty$ to 1/2 and to 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). Since the asymptotic magnitudes of the Poisson ratios for these models fall without the domain bounded for isotropic media by the thermodynamic constraints (1), it could be deduced that such models lose the physical meaning and could not be used in dynamic problems of viscoelasticity.

4. Conclusions

In the present paper, a complete analysis of time-dependent operators of Lamé’s parameters are studied for the cases where hereditary operators are represented in terms of the fractional derivative Scott-Blair, Kelvin–Voigt, Maxwell, and standard linear solid models without bulk relaxation. This is practically important since precisely these parameters define the velocities of longitudinal and transverse waves propagating in three-dimensional media.

Using the algebra of dimensionless Rabotnov’s fractional exponential functions, time-dependent operators for Poisson’s ratios have been obtained and analysed. 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.

The suggested models would allow one to solve different dynamic problems dealing with the propagation of dilatational and distortional waves through the viscoelastic media, the properties of which are described by various fractional derivative models.

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.