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

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


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 constitutive 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 (inelastic). In the case of small uniaxial deformations, it is shown analytically that the model reproduces qualitatively the main features of real experiments.

dynamic acoustoelasticity; softening; hysteresis; NDE

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

Figure 1: Dynamic acoustoelasticity measurement. (a) Evolution of the axial strain at the location of the probe over time. (b) Relative variation of the sound speed with respect to its initial value over time. (c) Hysteresis loop: versus in steady state. Reproduced from (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 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 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 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

where is the Cauchy stress and 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

is obtained, where denotes the material derivative of . Here, both and are possible. If , the sound speed proportional to decreases (softening). Inversely, increases the sound speed (recovery). In the choice of the evolution equation , particular care is taken to ensure that the Clausius-Duhem inequality is satisfied whatever the sign of . This is a major difference with damage modelling, where the internal variable describes an irreversible process, so that only 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 of the reference configuration moves to a position of the current configuration. The deformation gradient is a second-order tensor defined by (see e.g. (14); (13); (15))


where denotes the displacement field and is the gradient with respect to the material coordinates (Lagrangian gradient). In the reference configuration, the deformation gradient (1) is equal to the metric tensor . If the Euclidean space is described by an orthonormal basis and a Cartesian coordinate system, the matrix of the coordinates of 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 of any field is


In particular, the material derivative of the deformation gradient satisfies


where is the velocity field. The conservation of mass implies


where denotes the mass density in the deformed configuration, and denotes the mass density in the reference configuration. The motion is also driven by the conservation of momentum


where denotes the divergence with respect to the material coordinates. The expression of the Cauchy stress tensor 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 . The conservation of total energy writes:


where is the strain-rate tensor. The second principle of thermodynamics reads


where is the specific entropy.

2.2 The model


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


where is the absolute temperature. Multiplying (8) by , the local equations of thermodynamics (6)-(7) yield the Clausius-Duhem inequality


for all state and all evolution . The left-hand term in (9) is the dissipation 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


where is the strain energy density function, expressed in terms of the strain tensor . The function has dimensionless values, and is a storage energy. If and for all , then the classical case of hyperelasticity is recovered, where . 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):


where and denote the derivatives of and , respectively. The final constitutive laws are obtained for a given choice of strain tensor . In the next paragraph, the right Cauchy-Green tensor is used. For many other strain tensors, the constitutive laws can be deduced from , and similar derivations can be done.

Constitutive laws.

We choose the right Cauchy-Green tensor . The material derivative of the strain tensor is . For any second-order tensors , and , we recall that


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


where the hyperelastic stress


depends on .

The stress is a state function: it does not dependent on , which is not a variable of state. Thus, the term in the dissipation (13) is a scalar product between and a tensor which does not depend on . Moreover, the term does not depend on . Therefore, the Clausius-Duhem inequality (13) for all yields the constitutive law


where the hyperelastic stress is defined in (14).

Now, the Clausius-Duhem inequality (13) reduces to , for all state and all . Therefore, is either dependent on or equal to zero. We choose the simplest nontrivial dependence:


where  J/m3 and is a relaxation time. The parameter may be variable, e.g. dependent on the sign of , temperature, or any desired parameter. If , equation (16) gives the evolution equation


Otherwise (), the internal variable satisfies , i.e. where


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 . Indeed, with the assumption (16), the dissipation per unit volume in the material (13) is


If or , then no dissipation occurs: the thermodynamic process is reversible. If , the thermodynamic process is irreversible, which is the origin of hysteresis curves under a dynamic loading.

The effect of on the stress (15) is specified through . If for all , then no stress softening occurs. Indeed, classical hyperelasticity is recovered. If for all , 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 . Moreover, we assume that entails no stress softening: . A natural choice satisfying these requirements is


where .

We require that is an equilibrium point (18) if no strain is applied. Hence, one must have . If the softening function (20) is chosen, the convexity of ensures that the equilibrium point (18) is unique. Simple choices for are


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

To summarize, the equations of motion in Lagrangian coordinates are


where and are specified by (20) and (21)-(22), respectively. The expression of the hyperelastic stress is specified by (14) if the right Cauchy-Green tensor 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 of the internal energy can be replaced by a dependence to the invariants


In particular, the conservation of mass (4) rewrites as . The hyperelastic stress satisfies (14), where


with the tensor derivatives (13)


Thus, the following substitution


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


where 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 (see e.g. (15)). An example of strain energy density in terms of the invariants of is the Murnaghan’s law (19)


where are the Lamé parameters and 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 and :


Infinitesimal strain.

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


where is the infinitesimal strain tensor. Murnaghan’s law is used and the expression of the first Piola-Kirchhoff stress tensor is linearised with respect to the coordinates of :


The equations of motion (23) reduce to


which is nonlinear due to the slow dynamics. Classical elastodynamics are recovered if in (33).

Uniaxial strain.

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


where is the space derivative, is the strain and is the particle velocity. The hyperelastic stress satisfies , where is the derivative of the strain energy density function (22).

The functions and are specified by (20) and (21)-(22), respectively. An example of strain energy density function is given by Landau’s law (4); (5); (10):


where is the Young’s modulus and are higher-order elastic constants. When and are zero, Hooke’s law of linear elasticity


is recovered.

The relationship between Murnaghan’s law (29) and Landau’s law (35) is the following. If the uniaxial approximation is made, then and in (29). A polynomial expression of the strain energy density function with respect to is obtained,


where . By identification with Landau’s law (35), the parameters can be expressed in terms of the Lamé and Murnaghan parameters:


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 is driven by (17), where the strain energy is a constant. With the quadratic expression (21) of , the internal variable evolves exponentially in time towards , which is defined in (18). The corresponding relaxation time is .

Now, the case of uniaxial strain is considered. A sinusoidal strain with frequency  kHz and amplitude is applied locally. With the quadratic expression (21) of , the evolution equation (17) writes


where is a time constant and are the Fourier coefficients of the normalized strain energy .

The solution of the ordinary differential equation (39) is


The first term in (40) decreases exponentially in time with constant . The second term is the steady-state term, which oscillates at the frequency around its average value


where .

In the case of Landau’s law (35), the nonzero Fourier coefficients are given in table 1. At small strain amplitudes, and , the high-order terms in table 1 can be neglected. Thus, the case of Hooke’s law (36) is recovered, where and are the only nonzero Fourier coefficients. In particular, the value of the average of (41) is very close to the value obtained in the case of Hooke’s law:


From a practical point of view, if the Young’s modulus is known and the constants and are deduced from measurements at small sinusoidal loadings, then the parameters and of the model can be estimated.

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

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


If the material is linear-elastic without slow dynamics, the speed of sound reduces to . It is easier for the analysis to introduce the elastic modulus and its variation


On figure 1, the experimental variation in speed of sound is represented instead.

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


which reduces to if and equal zero. The average of over a period of forcing is deduced from (40) and (45):


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

(kg.m) (GPa)  (J.m-3) (s)
Table 2: Physical parameters.

On figure 2-(a), is represented up to  ms in the case of Hooke’s law (36) with the parameters from table 2. In this softening phase, decreases and reaches the steady state. At  ms, the excitation is stopped. Thus, goes to infinity in (40). During the recovery, increases exponentially in time towards zero with time constant  ms.

Figures 2-(b) and 2-(c) show the steady-state solution. On figure 2-(b), is represented with respect to the strain for several forcing amplitudes, according to equation (45) with . 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.




Figure 2: Analytical computation in the case of Hooke’s law. (a) Evolution of the relative variation in elastic modulus with respect to its initial value, when a sinusoidal strain is applied until  ms (40). (b) Hysteresis curves versus in steady state ( ms); (c) effect of hysteresis on the stress-strain relationship, where the stress is represented with respect to .

On figure 3, the behaviour of our model with Landau’s law (35) and is compared to the previous case of Hooke’s law (36). At strain amplitudes , the contribution of and 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 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 with the strain, when and are nonzero.



Figure 3: Comparison of the analytical computations in the cases of Hooke’s law and Landau’s law. (a) Evolution of the variation in elastic modulus when a sinusoidal strain with amplitude is applied until  ms (40)-(45). (b) Hysteresis curves versus in steady state.

Supplementary analytical results can be obtained in the case of Hooke’s law (36). In this case, the variation in elastic modulus (45) is , and the only nonzero Fourier coefficients in table 1 are and . The surface area of the hysteresis loops in figure 2-(b) is


which vanishes at high , low and high frequency , and low and high . The maximum value reached by the steady-state solution is


The strain amplitude for which the material is destroyed satisfies :


In the present configuration, . Thus, if the quadratic expression (21) of the storage energy is chosen, the model is only valid for small strains. In the case of the logarithmic expression (22) of , 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 , the other term corresponds to the storage energy . When , the internal energy is only elastic. As increases at constant strain, the strain energy decreases and the storage energy increases. Therefore, the internal energy is transferred from the strain to when increases, and inversely.

Let us assume that . The internal variable satisfies (18). With the quadratic expression (21) of , the internal variable is equal to


The value , which corresponds to a destructed material, is reached for strain energies . In the case of Hooke’s law (36) with the parameters from table 2, the maximum admissible strain is . This value is recovered by setting in equation (49). The logarithmic expression (22) of yields


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 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 increases. If , the strain energy does not depend on the strain anymore, which illustrates the destruction of the material.

Figure 4: Sketch of the strain energy per unit volume with respect to the strain , for several values of the internal variable .

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



Figure 5: Sketch of the internal energy per unit volume with respect to (10). It is represented (a) for several values of when  J.m-3; (b) for several values of (J.m-3) when .


In one space dimension and small strain, depends on and . The dissipation per unit volume (19) is a surface in - coordinates (figure 6). The expression of is deduced from the softening function (20), the quadratic storage energy (21), Hooke’s law (36) and the conservation of mass . This figure illustrates that the dissipation is positive, in agreement with the Clausius-Duhem inequality. Also, one can observe that no dissipation occurs if , which corresponds to the curve .

Figure 6: View of the dissipation in - coordinates (19). The black line marks the curve , i.e. the locus where no dissipation occurs (18).

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 describes the softening of the material;

  • an evolution equation for with a relaxation time is given;

  • 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. the new model satisfies the second principle of thermodynamics;

  2. the new model does not require an expression for the equilibrium value of , but an expression of the storage energy ;

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


.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 in (17). The internal variable is then deduced from the strain through (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 . For example, is used in (16) to describe the Mullins effect, which is observed in cyclic loading of filled rubber. An expression of the form


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


This expression satisfies if . In particular, 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).

Figure 7: Sketch of two sticking fibers of length with initial spacing , when withdrawn from a wetting liquid (grey). The height of fluid between the fibers is .

System of wet fibers.

A formal analogy with a system of two partially-immersed fibers of length can be made (figure 7). Initially, their spacing is . Then, the fibers are immersed in a fluid with surface tension . 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)


where is the wet length of the fibers. In the case of a system of cylindrical elastic fibers with radius and Young’s modulus , the expressions in (54) are


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


where is the altitude in the fluid, is the mass density of the fluid and 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 , which modifies the stress according to


This constitutive law is the same as (15) with the softening function (20). In one space dimension, the strain rate satisfies , where . The Clausius-Duhem inequality (9) rewrites as


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


When integrating (59) with respect to the strain , an integration constant appears, which we denote by . Thus, the internal energy per unit volume (10) is recovered, where . The Clausius-Duhem inequality (58) implies


for all and all .

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


where is a variable relaxation time and is the value of at equilibrium for a given stress. Various expressions of are proposed in the literature. In (9), reads


where is a stress and is the value of at zero stress. This expression is modified in (10) to ensure :


Injecting (61) in (60) yields the condition


for all in and all in .

In particular, (64) must hold for all when . In this case, the condition (64) reduces to for all such that . We deduce that must be negative or zero, i.e. . The expressions (62)-(63) of imply that 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 .

Modified model.

Expressions of must be chosen carefully. The condition (64) imposes that and have the same sign. Both functions of and 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 and equal zero for the same values of and . Combining both equalities, the condition


is deduced from the constitutive law (57). An expression of which satisfies (65) is not necessarily thermodynamically admissible. Moreover, one can note that such an expression depends on the strain energy density and on the storage energy .

Now, we examine the existence of a thermodynamically admissible expression of 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 . The necessary condition (65) writes


It rewrites as a cubic equation:


which may have multiple solutions.



Figure 8: Graph of the equilibrium value in the soft-ratchet model. (a) Roots of (67). The solid line corresponds to a thermodynamically admissible expression of . (b) Classical expressions “exp” (62) and “tanh” (63) of .

When using Cardano’s method, the discriminant


of the cubic function in (67) is positive if . 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 are displayed on figure 8-(b), where and  GPa. Among the three real roots of (67), only one satisfies (solid line on figure 8-(a)):


This thermodynamically admissible expression of is only defined when the discriminant (68) is positive, i.e. for strains smaller than . 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 (65). The expression of 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.


  1. Guyer RA, Johnson PA. 1999 Nonlinear mesoscopic elasticity: Evidence for a new class of materials. Phys. Today 52, 30–36.
  2. TenCate JA. 2011 Slow dynamics of earth materials: An experimental overview. Pure Appl. Geophys.