Nonlinear waves in solids with slow dynamics: an internal-variable model

TL;DR

This paper develops a thermodynamically consistent 3D internal-variable model for nonlinear wave softening in solids with slow dynamics, capturing hysteresis and strain-dependent sound speed in heterogeneous materials.

Contribution

A new 3D thermodynamically admissible model incorporating internal variables for slow dynamics in solids is introduced, extending previous 1D phenomenological models.

Findings

Model reproduces main experimental features of slow dynamics.

The model is thermodynamically consistent and dissipative.

Analytical results confirm qualitative agreement with experiments.

Abstract

In heterogeneous solids such as rocks and concrete, the speed of sound diminishes with the strain amplitude of a dynamic loading (softening). This decrease known as "slow dynamics" occurs at time scales larger than the period of the forcing. Also, hysteresis is observed in the steady-state response. The phenomenological model by Vakhnenko et al. is based on a variable that describes the softening of the material [Phys. Rev. E 70-1, 2004]. However, this model is 1D and it is not thermodynamically admissible. In the present article, a 3D model is derived in the framework of the finite strain theory. An internal variable that describes the softening of the material is introduced, as well as an expression of the specific internal energy. A mechanical constitu-tive law is deduced from the Clausius-Duhem inequality. Moreover, a family of evolution equations for the internal variable is proposed. Here, an evolution equation with one relaxation time is chosen. By construction, this new model of continuum is thermodynamically admissible and dissipative (inelas-tic). In the case of small uniaxial deformations, it is shown analytically that the model reproduces qualitatively the main features of real experiments.

Figures

Tables

Table 1. Table 1: Nonzero Fourier coefficients ( 39 ) in the case of Landau’s law ( 35 ).

$n$	0	1	2	3	4
$a_n$	$\frac{1}{2}-\frac{3}{16}\delta V^2$		$-\frac{1}{4}+\frac{1}{8}\delta V^2$		$-\frac{1}{32}\delta V^2$
$b_n$		$-\frac{1}{4}\beta V$		$\frac{1}{12}\beta V$

Table 2. Table 2: Physical parameters.

${\rho }_0$ (kg.m-3)	$E$ (GPa)	$\gamma$ (J.m-3)	$\tau$ (s)
$2600$	$10$	$20$	$7\times {10}^{-3}$

Equations

$$\mathbf{\sigma }=\left(1-g\right)\bar{\mathbf{\sigma }}(\mathbf{\chi }),$$

$$\dot{g}=\mathcal{S}(\mathbf{\chi },g)$$

$$\mathbf{F}=\mathbf{grad}{\mathbf{x}}_t=\mathbf{G}+\mathbf{grad}\mathbf{u},$$

$$\dot{\mathbf{\psi }}=\frac{\partial \mathbf{\psi }}{\partial t}.$$

$$\dot{\mathbf{F}}=\mathbf{grad}\mathbf{v},$$

$$\frac{{\rho }_0}{\rho }=\det (\mathbf{F}),$$

$${\rho }_0\dot{\mathbf{v}}=\mathbf{div}\left(\det (\mathbf{F})\mathbf{\sigma }\cdot {\mathbf{F}}^{-\top }\right),$$

$$\rho \dot{e}=\mathbf{\sigma }:\mathbf{D},$$

$$\rho \dot{s}⩾0,$$

$$\dot{e}=T\dot{s}+{\frac{\partial e}{\partial \mathbf{\chi }}|}_{s,g}:\dot{\mathbf{\chi }}+{\frac{\partial e}{\partial g}|}_{s,\mathbf{\chi }}\dot{g},$$

$$\mathcal{D}=\mathbf{\sigma }:\mathbf{D}-\rho {\frac{\partial e}{\partial \mathbf{\chi }}|}_{s,g}:\dot{\mathbf{\chi }}-\rho {\frac{\partial e}{\partial g}|}_{s,\mathbf{\chi }}\dot{g}⩾0,$$

$${\rho }_{0}e={\varphi }_1(g)W(\mathbf{\chi })+{\varphi }_2(g),$$

$$\rho {\frac{\partial e}{\partial \mathbf{\chi }}|}_{s,g}:\dot{\mathbf{\chi }}={\varphi }_1\frac{\rho }{{\rho }_0}\frac{\partial W}{\partial \mathbf{\chi }}:\dot{\mathbf{\chi }},$$

$$\rho {\frac{\partial e}{\partial g}|}_{s,\mathbf{\chi }}\dot{g}=\frac{\rho }{{\rho }_0}\left({\varphi }_1^{\prime }W+{\varphi }_2^{\prime }\right)\dot{g},$$

$$\mathbf{T}:({\mathbf{F}}^{\top }\cdot \mathbf{D}\cdot \mathbf{F})$$

$$=\mathrm{tr}(\mathbf{D}\cdot \mathbf{F}\cdot (\mathbf{F}\cdot \mathbf{T}{)}^{\top })$$

$$=\mathrm{tr}(\mathbf{D}\cdot (\mathbf{F}\cdot \mathbf{T}\cdot {\mathbf{F}}^{\top }{)}^{\top })$$

$$=(\mathbf{F}\cdot \mathbf{T}\cdot {\mathbf{F}}^{\top }):\mathbf{D}.$$

$$\mathcal{D}=\underset{{\mathcal{D}}_{el}}{\underbrace{\left(\mathbf{\sigma }-{\varphi }_1\bar{\mathbf{\sigma }}\right):\mathbf{D}}}\underset{{\mathcal{D}}_{inel}}{\underbrace{-\rho /{\rho }_0\left({\varphi }_1^{\prime }W+{\varphi }_2^{\prime }\right)\dot{g}}}⩾0,$$

$$\bar{\mathbf{\sigma }}=\frac{\rho }{{\rho }_0}\mathbf{F}\cdot 2\frac{\partial W}{\partial \mathbf{C}}\cdot {\mathbf{F}}^{\top }=\frac{1}{\det (\mathbf{F})}\mathbf{F}\cdot 2\frac{\partial W}{\partial \mathbf{C}}\cdot {\mathbf{F}}^{\top }$$

$$\mathbf{\sigma }={\varphi }_1(g)\bar{\mathbf{\sigma }}(\mathbf{\chi }),$$

$${\varphi }_1^{\prime }W+{\varphi }_2^{\prime }=-{\tau }_1\dot{g},$$

$$\dot{g}=-\frac{1}{{\tau }_1}\left({\varphi }_1^{\prime }(g)W(\mathbf{\chi })+{\varphi }_2^{\prime }(g)\right).$$

$$g_{eq}(\mathbf{\chi })={\left({\varphi }_2^{\prime }/{\varphi }_1^{\prime }\right)}^{-1}\left(-W(\mathbf{\chi })\right).$$

$$0⩽\mathcal{D}=\{\begin{array}{rlrl}& 0& & \text{if }\tau =0,\\ & \frac{\rho }{{\rho }_0}\frac{{\left({\varphi }_1^{\prime }W+{\varphi }_2^{\prime }\right)}^2}{{\tau }_1}& & \text{if }\tau >0.\end{array}$$

$${\varphi }_1(g)=1-g,$$

$${\varphi }_2(g)=\frac{1}{2}\gamma g^2,$$

$${\varphi }_2(g)=-\frac{1}{2}\gamma \ln (1-g^2),$$

$$\{\begin{array}{rl}& \dot{\mathbf{F}}=\mathbf{grad}\mathbf{v},\\ & {\rho }_0\dot{\mathbf{v}}=\mathbf{div}\left({\varphi }_1\det (\mathbf{F})\bar{\mathbf{\sigma }}\cdot {\mathbf{F}}^{-\top }\right),\\ & -{\tau }_1\dot{g}={\varphi }_1^{\prime }W+{\varphi }_2^{\prime },\end{array}$$

$$C_{\mathrm{I}}=\mathrm{tr}(\mathbf{C}),$$

$$C_{\mathrm{II}}=\frac{1}{2}\left(\mathrm{tr}(\mathbf{C}{)}^2-\mathrm{tr}({\mathbf{C}}^2)\right),$$

$$C_{\mathrm{III}}=\det (\mathbf{C}).$$

$$\frac{\partial W}{\partial \mathbf{C}}=\frac{\partial W}{\partial C_{\mathrm{I}}}\frac{\partial C_{\mathrm{I}}}{\partial \mathbf{C}}+\frac{\partial W}{\partial C_{\mathrm{II}}}\frac{\partial C_{\mathrm{II}}}{\partial \mathbf{C}}+\frac{\partial W}{\partial C_{\mathrm{III}}}\frac{\partial C_{\mathrm{III}}}{\partial \mathbf{C}},$$

$$\frac{\partial C_{\mathrm{I}}}{\partial \mathbf{C}}=\mathbf{G},$$

Full text

1 Introduction

Rocks and concrete are known to have a strong nonlinear behaviour. Quasistatic compression or traction tests show a nonlinear stress-strain relationship. A hysteresis loop is observed when the loading is increased and decreased. This phenomenon is interpreted as a memory effect [1]. The longitudinal vibrations of a rod of material also show highly nonlinear features. Indeed, a frequency shift of resonance peaks is observed when the amplitude of the vibration is increased. The frequency shift reveals a global softening of the material with the strain amplitude [1, 2].

In dynamic acoustoelastic testing (DAET), the speed of sound is measured locally over time, when longitudinal vibrations are simultaneously applied to the whole sample. As illustrated on figure 1-(b), a decrease with time of the measured sound speed is observed. This softening occurs over a time scale larger than the period of the dynamic loading, which highlights the phenomenon of slow dynamics. Moreover, the evolution of this speed with respect to the strain presents an hysteresis curve in steady state (figure 1-(c)). When the excitation is stopped ($t\approx 0.08$ s in figure 1-(b)), the sound speed increases, and recovers gradually its initial value (recovery). All these phenomena are accentuated when the strain amplitude is increased [3].

Several dynamic models which reproduce these features can be found in the literature [4, 5]. One approach consists in incorporating a dependency on the strain rate in the stress-strain relationship [6]. Another approach is related to the Preisach-Mayergoyz model, which is based on a discrete representation of hysteresis [7, 4, 8]. The soft-ratchet model of Vakhnenko et al. results from a different approach [9, 10]. A new variable $g$ is introduced so as to describe the softening. Interpreted as a concentration of activated defects, this variable modifies the apparent elastic modulus. Also, an evolution equation for $g$ is provided. In this equation, a relaxation time is incorporated to describe the slow dynamics. Nevertheless, the soft-ratchet model was developed in one space dimension, and does not generalize straightforwardly to higher space dimensions. Moreover, thermodynamical issues are not considered in the construction of this model.

In the present article, a new phenomenological model is proposed in the context of the finite strain theory. Similarly to the soft-ratchet model, a scalar internal variable $g$ is introduced to describe the softening. Our model, derived from beginning in the framework of continuum mechanics with internal variables [11, 12], satisfies by construction the principles of thermodynamics. As shown later in the document, a particular choice of the internal energy yields a separable constitutive law

[TABLE]

where $\bm{\sigma}$ is the Cauchy stress and $\bm{\chi}$ is a strain tensor. Such a constitutive law resembles classical models of irreversible damage. Furthermore, an evolution equation for the internal variable of the form

[TABLE]

is obtained, where $\dot{g}$ denotes the material derivative of $g$. Here, both $\dot{g}\geqslant 0$ and $\dot{g}\leqslant 0$ are possible. If $\dot{g}\geqslant 0$, the sound speed proportional to $\sqrt{1-g}$ decreases (softening). Inversely, $\dot{g}\leqslant 0$ increases the sound speed (recovery). In the choice of the evolution equation $\mathscr{S}(\bm{\chi},g)$, particular care is taken to ensure that the Clausius-Duhem inequality is satisfied whatever the sign of $\dot{g}$. This is a major difference with damage modelling, where the internal variable $g$ describes an irreversible process, so that only $\dot{g}\geqslant 0$ is possible [13].

The article is organized as follows. In section 2, the model is constructed, leading to the constitutive law and the evolution equation. Several examples of finite-strain models are provided for illustration purposes. In particular, the cases of infinitesimal strain and uniaxial strain are addressed. In section 3, the equations are solved analytically in a particular configuration. The three expected phenomena — softening, slow dynamics and hysteresis — are reproduced by the model. In appendix .1, a link is made between the new model and quasistatic models of filled rubber. Moreover, a formal analogy with a system of wet sticking fibers is proposed. Lastly, in appendix .2, we demonstrate that the soft-ratchet model originally proposed by Vakhnenko et al. is not thermodynamically relevant.

2 Construction of the model

2.1 Basic equations

Let us consider an homogeneous continuum on which no external volume force is applied, and no heat transfer occurs. Furthermore, self-gravitation is neglected. A particle initially located at some position $\bm{x}_{0}$ of the reference configuration moves to a position $\bm{x}_{t}$ of the current configuration. The deformation gradient is a second-order tensor defined by (see e.g. [14, 13, 15])

[TABLE]

where $\bm{u}=\bm{x}_{t}-\bm{x}_{0}$ denotes the displacement field and $\mathbf{grad}$ is the gradient with respect to the material coordinates $\bm{x}_{0}$ (Lagrangian gradient). In the reference configuration, the deformation gradient (1) is equal to the metric tensor $\bm{G}$. If the Euclidean space is described by an orthonormal basis and a Cartesian coordinate system, the matrix of the coordinates of $\bm{G}$ is the identity matrix.

The choice of a representation of motion — Eulerian or Lagrangian — does not affect the expressions of the constitutive laws and the evolution equations. However, it affects the expression of the material derivative and the equations of motion. Here, the Lagrangian representation of motion is used. Hence, the material derivative $\dot{\bm{\psi}}$ of any field $\bm{\psi}(\bm{x}_{0},t)$ is

[TABLE]

In particular, the material derivative of the deformation gradient satisfies

[TABLE]

where $\bm{v}(\bm{x}_{0},t)$ is the velocity field. The conservation of mass implies

[TABLE]

where $\rho$ denotes the mass density in the deformed configuration, and $\rho_{0}$ denotes the mass density in the reference configuration. The motion is also driven by the conservation of momentum

[TABLE]

where $\mathbf{div}$ denotes the divergence with respect to the material coordinates. The expression of the Cauchy stress tensor $\bm{\sigma}$ will be specified later on.

As usual in acoustics, the thermodynamic process is assumed to be adiabatic. The first principle of thermodynamics introduces the specific internal energy $e$. The conservation of total energy writes:

[TABLE]

where $\bm{D}=\frac{1}{2}(\dot{\bm{F}}\cdot\bm{F}^{-1}+\bm{F}^{-\top}\!\cdot\dot{\bm{F}}^{\top})$ is the strain-rate tensor. The second principle of thermodynamics reads

[TABLE]

where $s$ is the specific entropy.

2.2 The model

Preliminaries.

We choose the following variables of state: the specific entropy $s$, the strain tensor $\bm{\chi}$, and an additional scalar variable $g$, which is introduced to represent the softening/recovery of the material. Consequently, the Gibbs identity reads

[TABLE]

where $T=\left.{\partial e}/{\partial s}\right|_{\bm{\chi},g}>0$ is the absolute temperature. Multiplying (8) by $\rho$, the local equations of thermodynamics (6)-(7) yield the Clausius-Duhem inequality

[TABLE]

for all state $\{s,\bm{\chi},g\}$ and all evolution $\{\dot{s},\dot{\bm{\chi}},\dot{g}\}$. The left-hand term in (9) is the dissipation $\mathscr{D}$ per unit volume of material (W.m-3).

The main ingredient of the model is an expression of the internal energy per unit volume of the form

[TABLE]

where $W$ is the strain energy density function, expressed in terms of the strain tensor $\bm{\chi}$. The function $\phi_{1}$ has dimensionless values, and $\phi_{2}$ is a storage energy. If $\phi_{1}(g)=1$ and $\phi_{2}(g)=0$ for all $g$, then the classical case of hyperelasticity is recovered, where $\rho_{0}\,e=W(\bm{\chi})$. The expression of the internal energy (10) is analogous to the Ogden-Roxburgh model of filled rubber [16]. It is also formally analogous to a model of wet sticking fibers [17]. These similarities are detailed in appendix .1.

With the assumption (10), the following substitutions are made in the inequality (9):

[TABLE]

where $\phi^{\prime}_{1}$ and $\phi^{\prime}_{2}$ denote the derivatives of $\phi_{1}$ and $\phi_{2}$, respectively. The final constitutive laws are obtained for a given choice of strain tensor $\bm{\chi}$. In the next paragraph, the right Cauchy-Green tensor $\bm{C}=\bm{F}^{\top}\!\cdot\bm{F}$ is used. For many other strain tensors, the constitutive laws can be deduced from $\bm{C}$, and similar derivations can be done.

Constitutive laws.

We choose the right Cauchy-Green tensor $\bm{\chi}=\bm{C}=\bm{F}^{\top}\!\cdot\bm{F}$. The material derivative of the strain tensor is $\dot{\bm{C}}=2\bm{F}^{\top}\!\cdot\bm{D}\cdot\bm{F}$. For any second-order tensors $\bm{T}$, $\bm{F}$ and $\bm{D}$, we recall that

[TABLE]

Therefore, the Clausius-Duhem inequality (9) with the substitutions (11) reduces to

[TABLE]

where the hyperelastic stress

[TABLE]

depends on $\bm{F}$.

The stress $\bm{\sigma}$ is a state function: it does not dependent on $\bm{D}$, which is not a variable of state. Thus, the term $\mathscr{D}_{\mathit{el}}$ in the dissipation (13) is a scalar product between $\bm{D}$ and a tensor which does not depend on $\bm{D}$. Moreover, the term $\mathscr{D}_{\mathit{inel}}$ does not depend on $\bm{D}$. Therefore, the Clausius-Duhem inequality (13) for all $\bm{D}$ yields the constitutive law

[TABLE]

where the hyperelastic stress $\bar{\bm{\sigma}}$ is defined in (14).

Now, the Clausius-Duhem inequality (13) reduces to $\mathscr{D}_{\mathit{inel}}\geqslant 0$, for all state and all $\dot{g}$. Therefore, $\phi_{1}^{\prime}\,W+\phi_{2}^{\prime}$ is either dependent on $\dot{g}$ or equal to zero. We choose the simplest nontrivial dependence:

[TABLE]

where $\tau_{1}=\tau\times 1$ J/m3 and $\tau\geqslant 0$ is a relaxation time. The parameter $\tau$ may be variable, e.g. dependent on the sign of $\dot{g}$, temperature, or any desired parameter. If $\tau\neq 0$, equation (16) gives the evolution equation

[TABLE]

Otherwise ($\tau=0$), the internal variable $g$ satisfies $\phi_{1}^{\prime}\,W+\phi_{2}^{\prime}=0$, i.e. $g=g_{\mathit{eq}}(\bm{\chi})$ where

[TABLE]

In this case, the internal variable is instantaneously modified when the strain varies: no slow dynamics occurs.

The previous choice ensures that the Clausius-Duhem inequality is satisfied, independently of the sign of $\dot{g}$. Indeed, with the assumption (16), the dissipation per unit volume in the material (13) is

[TABLE]

If $\tau=0$ or $\tau\rightarrow{+\infty}$, then no dissipation occurs: the thermodynamic process is reversible. If $0<\tau<{+\infty}$, the thermodynamic process is irreversible, which is the origin of hysteresis curves under a dynamic loading.

The effect of $g$ on the stress (15) is specified through $\phi_{1}$. If $\phi_{1}(g)=1$ for all $g$, then no stress softening occurs. Indeed, classical hyperelasticity is recovered. If $\phi_{1}(g)=0$ for all $g$, then the stress does not depend on the strain any more: the material is destroyed. For the physical relevance of the constitutive law (15), we assume that $\phi_{1}>0$. Moreover, we assume that $g=0$ entails no stress softening: $\phi_{1}(0)=1$. A natural choice satisfying these requirements is

[TABLE]

where $g<1$.

We require that $g=0$ is an equilibrium point (18) if no strain is applied. Hence, one must have $\phi_{2}^{\prime}(0)=0$. If the softening function (20) is chosen, the convexity of $\phi_{2}$ ensures that the equilibrium point (18) is unique. Simple choices for $\phi_{2}$ are

[TABLE]

where $\gamma>0$ is an energy per unit volume. The choice (20)-(22) ensures that $g$ is bounded by 1. In the vicinity of $g=0$, both expressions (21) and (22) are equivalent.

To summarize, the equations of motion in Lagrangian coordinates are

[TABLE]

where $\phi_{1}$ and $\phi_{2}$ are specified by (20) and (21)-(22), respectively. The expression of the hyperelastic stress $\bar{\bm{\sigma}}$ is specified by (14) if the right Cauchy-Green tensor $\bm{\chi}=\bm{C}$ is used. Otherwise, elementary tensor algebra yields the expression of the constitutive law.

In the next section, a few cases are detailed: the isotropic case, the case of infinitesimal strain and the case of uniaxial strain.

2.3 Particular cases

Isotropic case.

The dependence to $\bm{C}$ of the internal energy can be replaced by a dependence to the invariants

[TABLE]

In particular, the conservation of mass (4) rewrites as $\rho_{0}/\rho=\sqrt{C_{\mathrm{III}}}$. The hyperelastic stress $\bar{\bm{\sigma}}$ satisfies (14), where

[TABLE]

with the tensor derivatives [13]

[TABLE]

Thus, the following substitution

[TABLE]

can be made in equation (14). In the literature, several strain energy density functions can be found. In terms of the invariants of $\bm{C}$, a classical example is the compressible Mooney-Rivlin model [18]

[TABLE]

where $(\mathfrak{c}_{1},\mathfrak{c}_{2},\mathfrak{d}_{1})$ are material parameters. This hyperelastic model (28) is classically used in mechanics of elastomers.

Sometimes, the strain energy density function is expressed in terms of the Green-Lagrange strain tensor $\bm{E}=\frac{1}{2}(\bm{C}-\bm{G})$ (see e.g. [15]). An example of strain energy density in terms of the invariants of $\bm{E}$ is the Murnaghan’s law [19]

[TABLE]

where $(\lambda,\mu)$ are the Lamé parameters and $(\mathfrak{l},\mathfrak{m},\mathfrak{n})$ are the Murnaghan coefficients. The latter are third-order elastic constants. The hyperelastic model (29) is widely used in the community of nondestructive testing [20, 21]. For conversions, one has the following relations between the invariants of $\bm{E}$ and $\bm{C}$:

[TABLE]

Infinitesimal strain.

The Green-Lagrange strain tensor is linearised with respect to the displacement:

[TABLE]

where $\bm{\varepsilon}=\frac{1}{2}(\bm{F}+\bm{F}^{\top})-\bm{G}$ is the infinitesimal strain tensor. Murnaghan’s law is used and the expression of the first Piola-Kirchhoff stress tensor $\det(\bm{F})\,\bar{\bm{\sigma}}\cdot\bm{F}^{-\top}\!$ is linearised with respect to the coordinates of $\bm{\varepsilon}$:

[TABLE]

The equations of motion (23) reduce to

[TABLE]

which is nonlinear due to the slow dynamics. Classical elastodynamics are recovered if $\tau_{1}\rightarrow{+\infty}$ in (33).

Uniaxial strain.

In this case, only one component of the displacement field remains. The corresponding coordinate $u$ is assumed to be invariant with respect to the other coordinates. Thus, the equations of motion (23) write now as a $3\times 3$ differential system:

[TABLE]

where $\partial_{x}$ is the space derivative, $\varepsilon=\partial_{x}u$ is the strain and $v=\dot{u}$ is the particle velocity. The hyperelastic stress satisfies $\bar{\sigma}=W^{\prime}(\varepsilon)$, where $W^{\prime}$ is the derivative of the strain energy density function [22].

The functions $\phi_{1}$ and $\phi_{2}$ are specified by (20) and (21)-(22), respectively. An example of strain energy density function is given by Landau’s law [4, 5, 10]:

[TABLE]

where $E$ is the Young’s modulus and $(\beta,\delta)$ are higher-order elastic constants. When $\beta$ and $\delta$ are zero, Hooke’s law of linear elasticity

[TABLE]

is recovered.

The relationship between Murnaghan’s law (29) and Landau’s law (35) is the following. If the uniaxial approximation is made, then $E_{\mathrm{I}}=\varepsilon\left(1+\frac{1}{2}\varepsilon\right)$ and $E_{\mathrm{II}}=E_{\mathrm{III}}=0$ in (29). A polynomial expression of the strain energy density function with respect to $\varepsilon$ is obtained,

[TABLE]

where $\vartheta=(\mathfrak{l}+2\mathfrak{m})/(\lambda+2\mu)$. By identification with Landau’s law (35), the parameters $(E,\beta,\delta)$ can be expressed in terms of the Lamé and Murnaghan parameters:

[TABLE]

A similar calculus can be performed with the Mooney-Rivlin model (28).

3 Analysis of the model

3.1 Analytical results

From now on, the softening function (20) is used. If a strain step is applied locally, then $g$ is driven by (17), where the strain energy $W$ is a constant. With the quadratic expression (21) of $\phi_{2}$, the internal variable $g$ evolves exponentially in time towards $g_{\mathit{eq}}(\bm{\chi})$, which is defined in (18). The corresponding relaxation time is $\tau_{\gamma}=\tau_{1}/\gamma$.

Now, the case of uniaxial strain is considered. A sinusoidal strain with frequency $f_{c}=\omega_{c}/2\pi=10$ kHz and amplitude $V$ is applied locally. With the quadratic expression (21) of $\phi_{2}$, the evolution equation (17) writes

[TABLE]

where $\tau_{EV}=\tau_{1}/(EV^{2})$ is a time constant and $(a_{n},b_{n})$ are the Fourier coefficients of the normalized strain energy $W/(EV^{2})$.

The solution of the ordinary differential equation (39) is

[TABLE]

The first term in (40) decreases exponentially in time with constant $\tau_{\gamma}$. The second term is the steady-state term, which oscillates at the frequency $f_{c}$ around its average value

[TABLE]

where $\tau_{\gamma}/\tau_{EV}=EV^{2}/\gamma$.

In the case of Landau’s law (35), the nonzero Fourier coefficients are given in table 1. At small strain amplitudes, $\beta V\ll 1$ and $\delta V^{2}\ll 1$, the high-order terms in table 1 can be neglected. Thus, the case of Hooke’s law (36) is recovered, where $a_{0}=1/2$ and $a_{2}=-1/4$ are the only nonzero Fourier coefficients. In particular, the value of the average of $g$ (41) is very close to the value obtained in the case of Hooke’s law:

[TABLE]

From a practical point of view, if the Young’s modulus $E$ is known and the constants $\tau_{\gamma}$ and $\left\langle g\right\rangle_{t\gg\tau_{\gamma}}$ are deduced from measurements at small sinusoidal loadings, then the parameters $\tau$ and $\gamma$ of the model can be estimated.

In the 1D case (34), the speed of sound is

[TABLE]

If the material is linear-elastic without slow dynamics, the speed of sound reduces to $c_{0}=\sqrt{E/\rho_{0}}$. It is easier for the analysis to introduce the elastic modulus $M=\rho_{0}\,c^{2}$ and its variation

[TABLE]

On figure 1, the experimental variation in speed of sound $\Delta c/c=(c-c_{0})/c_{0}$ is represented instead.

When Landau’s law (35) is used, the variation in elastic modulus is

[TABLE]

which reduces to $-g$ if $\beta$ and $\delta$ equal zero. The average of $\Delta M/M$ over a period of forcing is deduced from (40) and (45):

[TABLE]

The diminution of the elastic modulus with the square of the strain amplitude is similar to the Payne effect in filled rubber [23].

On figure 2-(a), $\Delta M/M$ is represented up to $t=5$ ms in the case of Hooke’s law (36) with the parameters from table 2. In this softening phase, $\Delta M/M$ decreases and reaches the steady state. At $t=5$ ms, the excitation is stopped. Thus, $\tau_{EV}$ goes to infinity in (40). During the recovery, $\Delta M/M$ increases exponentially in time towards zero with time constant $\tau_{\gamma}=0.35$ ms.

Figures 2-(b) and 2-(c) show the steady-state solution. On figure 2-(b), $\Delta M/M$ is represented with respect to the strain for several forcing amplitudes, according to equation (45) with $\beta=\delta=0$. A hysteretic behaviour caused by the dissipation is observed. Figure 2-(c) is an alternative representation of the phenomenon for several strain amplitudes. Here, the effect of increasing strain levels on the stress-strain relationship is outlined.

On figure 3, the behaviour of our model with Landau’s law (35) and $\{\beta=10^{2},\delta=10^{6}\}$ is compared to the previous case of Hooke’s law (36). At strain amplitudes $V\approx 10^{-5}$, the contribution of $\beta$ and $\delta$ in the Fourier coefficients is not significant (table 1). On figure 3-(a), the softening phases are compared. Figure 3-(b) represents the hysteresis curves. More important variations of $\Delta M/M$ are observed in the case of Landau’s law, as well as a loss of symmetry in the hysteresis curves. These phenomena are due to the dependence (45) of $\Delta M/M$ with the strain, when $\beta$ and $\delta$ are nonzero.

Supplementary analytical results can be obtained in the case of Hooke’s law (36). In this case, the variation in elastic modulus (45) is $\Delta M/M=-g$, and the only nonzero Fourier coefficients in table 1 are $a_{0}$ and $a_{2}$. The surface area of the hysteresis loops in figure 2-(b) is

[TABLE]

which vanishes at high $\gamma$, low and high frequency $f_{c}$, and low and high $\tau$. The maximum value reached by the steady-state solution is

[TABLE]

The strain amplitude $V_{\mathit{max}}$ for which the material is destroyed satisfies $g_{\mathit{max}}=1$:

[TABLE]

In the present configuration, $V_{\mathit{max}}\approx 8.8\times 10^{-5}$. Thus, if the quadratic expression (21) of the storage energy $\phi_{2}$ is chosen, the model is only valid for small strains. In the case of the logarithmic expression (22) of $\phi_{2}$, no strain limit is imposed by the slow dynamics.

3.2 Properties

Internal energy.

According to the equation (10), the internal energy per unit volume is separated into two terms. One term corresponds to the strain energy $\phi_{1}W$, the other term corresponds to the storage energy $\phi_{2}$. When $g=0$, the internal energy is only elastic. As $g$ increases at constant strain, the strain energy decreases and the storage energy increases. Therefore, the internal energy is transferred from the strain to $\phi_{2}$ when $g$ increases, and inversely.

Let us assume that $\tau=0$. The internal variable satisfies $g=g_{\mathit{eq}}(\bm{\chi})$ (18). With the quadratic expression (21) of $\phi_{2}$, the internal variable is equal to

[TABLE]

The value $g=1$, which corresponds to a destructed material, is reached for strain energies $W\geqslant\gamma$. In the case of Hooke’s law (36) with the parameters from table 2, the maximum admissible strain is $\sqrt{2\gamma/E}\approx 6.3\times 10^{-5}$. This value is recovered by setting $\tau_{1}=0$ in equation (49). The logarithmic expression (22) of $\phi_{2}$ yields

[TABLE]

which is always between zero and one. Therefore, there is no strain limit in this case.

Figure 4 represents the strain energy per unit volume $\phi_{1}W$ when the geometry is 1D. The strain energy density function is issued from Hooke’s law (36) and the softening function (20) is used (parameters from table 2). One can observe that the strain energy decreases as $g$ increases. If $g=1$, the strain energy does not depend on the strain anymore, which illustrates the destruction of the material.

On figures 5-(a) and 5-(b), the internal energy is represented with respect to $g$, where the quadratic expression (21) of the storage energy $\phi_{2}$ is used. The values of $g_{\mathit{eq}}$ correspond to the abscissas of the local minima of the curves (50). On figure 5-(a), one can observe an increase in $g_{\mathit{eq}}$ when the strain increases. No asymptote avoids to reach the value $g=1$, which destroys the material. On figure 5-(b), one can observe an increase in $g_{\mathit{eq}}$ when $\gamma$ decreases. Again, no asymptote avoids to reach the value $g=1$, which destroys the material.

Dissipation.

In one space dimension and small strain, $\mathscr{D}$ depends on $\varepsilon$ and $g$. The dissipation per unit volume (19) is a surface in $\varepsilon$-$g$ coordinates (figure 6). The expression of $\mathscr{D}$ is deduced from the softening function (20), the quadratic storage energy (21), Hooke’s law (36) and the conservation of mass $\rho_{0}/\rho=1+\varepsilon$. This figure illustrates that the dissipation is positive, in agreement with the Clausius-Duhem inequality. Also, one can observe that no dissipation occurs if $\tau=0$, which corresponds to the curve $g=g_{\mathit{eq}}(\varepsilon)=E\varepsilon^{2}/(2\gamma)$.

4 Conclusion

A new model for the dynamic behaviour of solids is proposed. The following features are common with the soft-ratchet model of Vakhnenko et al. [9]:

a variable $g$ describes the softening of the material;

• 2.

an evolution equation for $g$ with a relaxation time $\tau$ is given;

• 3.

a low number of extra parameters for the non-classical effects is required.

In comparison with the soft-ratchet model, several differences can be outlined:

1. (i)

the new model satisfies the second principle of thermodynamics; 2. (ii)

the new model does not require an expression for the equilibrium value $g_{\mathit{eq}}(\sigma)$ of $g$, but an expression of the storage energy $\phi_{2}(g)$; 3. (iii)

the new model generalizes naturally to higher space dimensions.

The point (i) is a major difference (see appendix .2), which ensures that our model is thermodynamically relevant. As shown in section 3, the new model reproduces qualitatively the macroscopic behaviour of real media.

Our approach is purely phenomenological. No physical interpretation of $g$ at the microscopic scale is known. To go further, some similarities with other materials are pointed out in appendix .1, in particular with filled rubber. It seems that the dynamic response of rocks is similar to the Payne effect [23], and that the quasistatic response of rocks is similar to the Mullins effect [24, 25]. In mechanics of elastomers, existing quasistatic models have a very similar structure to our dynamic model [16, 26]. By analogy, the coupling of nonlinear viscoelasticity and heat conduction could be a key for future physical modelling (see e.g. [27]). Lastly, from a microscopic point of view, both materials are roughly made of a matrix with particles inside. These similarities may be used for future micromechanical modelling.

Future work will be devoted to 2D and 3D numerical modelling of the nonlinear wave propagation. Since the system of partial differential equations is nonlinear, a mathematical study of the existence and the smoothness of solutions is required. Also, the computation of long-time periodic solutions will be addressed. Lastly, comparisons with real experiments should be done to validate the model.

Appendix

.1 Analogies with other models

Quasistatic loading of filled rubber.

In the case of a quasistatic process, equilibrium is satisfied over the transformation. This is equivalent to have $\dot{g}=0$ in (17). The internal variable is then deduced from the strain through $g=g_{\mathit{eq}}(\bm{\chi})$ (18). Due to the constitutive relation (15), the stress depends explicitly on the strain. Therefore, no hysteresis occurs in the stress-strain relationship.

Pseudo-elastic models are designed to incorporate hysteresis and memory effects. Additional variables which are stored along the loading path can be used in the storage energy $\phi_{2}$. For example, $W_{\mathit{max}}=\max_{t}W(\bm{\chi})$ is used in [16] to describe the Mullins effect, which is observed in cyclic loading of filled rubber. An expression of the form

[TABLE]

is proposed in [16]. From (18), one deduces the expression of the internal variable

[TABLE]

This expression satisfies $g_{\mathit{eq}}=0$ if $W(\bm{\chi})=W_{\mathit{max}}$. In particular, $g_{\mathit{eq}}=0$ along the primary loading path. In the case of the end-point memory phenomenon which is observed in rocks [1], the pseudo-elastic model [16] can be adapted as described in section 4 of [26]. For further reading, a review on existing models of rubber can be found in [24, 25].

System of wet fibers.

A formal analogy with a system of two partially-immersed fibers of length $L$ can be made (figure 7). Initially, their spacing is $d$. Then, the fibers are immersed in a fluid with surface tension $\Upsilon$. When withdrawn quasi-statically, they stick together. The internal energy of this system is the sum of the bending energy in the fibers and the energy due to the surface tension of the fluid. Thus, [17]

[TABLE]

where $L_{\mathit{stick}}$ is the wet length of the fibers. In the case of a system of cylindrical elastic fibers with radius $r$ and Young’s modulus $E$, the expressions in (54) are

[TABLE]

Due to the geometry of the meniscus and the law of hydrostatics, one has

[TABLE]

where $z$ is the altitude in the fluid, $\rho_{\mathrm{f}}$ is the mass density of the fluid and $g_{\mathrm{n}}$ is the standard gravity. A sign mistake has been found in equation (2) of [17]. Equations (55)-(56) are taken from equations (3)-(4) of [17], where the sign is correct. Formally, the energy (54) is similar to the energy (10).

.2 Limitations of the soft-ratchet model

Thermodynamical analysis.

The soft-ratchet model is a particular case of 1D model with internal variable of state [9]. Thus, we carry out the thermodynamical analysis from section 2. The soft-ratchet model introduces a concentration of activated defects $g$, which modifies the stress according to

[TABLE]

This constitutive law is the same as (15) with the softening function (20). In one space dimension, the strain rate satisfies $D=\dot{\varepsilon}/F$, where $F=1+\varepsilon$. The Clausius-Duhem inequality (9) rewrites as

[TABLE]

for all state and all evolution. Due to the constitutive law (57), the specific internal energy must satisfy

[TABLE]

When integrating (59) with respect to the strain $\varepsilon$, an integration constant appears, which we denote by $\phi_{2}(g)$. Thus, the internal energy per unit volume (10) is recovered, where $W^{\prime}(\varepsilon)=\bar{\sigma}(\varepsilon)$. The Clausius-Duhem inequality (58) implies

[TABLE]

for all $\{\varepsilon,g\}$ and all $\dot{g}$.

In the soft-ratchet model, the evolution equation for $g$ has the form

[TABLE]

where $\tau>0$ is a variable relaxation time and $g_{\mathit{eq}}(\sigma)$ is the value of $g$ at equilibrium for a given stress. Various expressions of $g_{\mathit{eq}}$ are proposed in the literature. In [9], $g_{\mathit{eq}}$ reads

[TABLE]

where $\tilde{\sigma}$ is a stress and $g_{0}$ is the value of $g_{\mathit{eq}}$ at zero stress. This expression is modified in [10] to ensure $g_{\mathit{eq}}<1$:

[TABLE]

Injecting (61) in (60) yields the condition

[TABLE]

for all $\varepsilon$ in $\left]{-1},{+\infty}\right[$ and all $g$ in $\left[0,1\right]$.

In particular, (64) must hold for all $g$ when $\varepsilon=0$. In this case, the condition (64) reduces to $g\geqslant g_{0}$ for all $g$ such that $\phi^{\prime}_{2}(g)>0$. We deduce that $g_{0}$ must be negative or zero, i.e. $g_{0}=0$. The expressions (62)-(63) of $g_{\mathit{eq}}$ imply that $g_{\mathit{eq}}$ is always equal to zero, which is not physically relevant. Something must be modified in the soft-ratchet model to satisfy equation (64). Here, we propose to seek thermodynamically admissible expressions of $g_{\mathit{eq}}$.

Modified model.

Expressions of $g_{\mathit{eq}}$ must be chosen carefully. The condition (64) imposes that $\phi_{2}^{\prime}(g)-W(\varepsilon)$ and $g-g_{\mathit{eq}}(\sigma)$ have the same sign. Both functions of $\varepsilon$ and $g$ are smooth. Hence, they equal zero with a change in sign or with a gradient equal to zero. Since the gradient of both functions is nonzero, it implies that $\phi_{2}^{\prime}(g)-W(\varepsilon)$ and $g-g_{\mathit{eq}}(\sigma)$ equal zero for the same values of $\varepsilon$ and $g$. Combining both equalities, the condition

[TABLE]

is deduced from the constitutive law (57). An expression of $g_{\mathit{eq}}$ which satisfies (65) is not necessarily thermodynamically admissible. Moreover, one can note that such an expression depends on the strain energy density $W$ and on the storage energy $\phi_{2}$.

Now, we examine the existence of a thermodynamically admissible expression of $g_{\mathit{eq}}$ in a particular case. To do so, the strain energy density from Hooke’s law (36) is chosen. We select the quadratic expression (21) of the storage energy $\phi_{2}$. The necessary condition (65) writes

[TABLE]

It rewrites as a cubic equation:

[TABLE]

which may have multiple solutions.

When using Cardano’s method, the discriminant

[TABLE]

of the cubic function in (67) is positive if $|\sigma|<\sqrt{8E\gamma/27}$. In this case, the three roots of (67) are real. On figure 8-(a), the three real roots are represented, where the parameters are issued from table 2. For comparison, the classical expressions (62) and (63) of $g_{\mathit{eq}}(\sigma)$ are displayed on figure 8-(b), where $g_{0}=0.1$ and $\tilde{\sigma}=0.1$ GPa. Among the three real roots of (67), only one satisfies $g_{\mathit{eq}}(0)=0$ (solid line on figure 8-(a)):

[TABLE]

This thermodynamically admissible expression of $g_{\mathit{eq}}$ is only defined when the discriminant (68) is positive, i.e. for strains smaller than $\sqrt{8\gamma/(27E)}\approx 2.4\times 10^{-5}$. This bound has the same order of magnitude as (49).

To summarize, we have shown that the soft-ratchet model is not thermodynamically relevant. A modification of this model has been examined, which results in an implicit definition of $g_{\mathit{eq}}$ (65). The expression of $g_{\mathit{eq}}$ is dependent on the choice of a strain energy density function and a storage energy. Furthermore, equation (65) may be hard to solve analytically in some cases. Lastly, the domain of validity of the model may be restricted.