Nonlinear waves in solids with slow dynamics: an internalvariable model
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 steadystate response. The phenomenological model by Vakhnenko et al. is based on a variable that describes the softening of the material [Phys. Rev. E 701, 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 ClausiusDuhem 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.
keywords:
dynamic acoustoelasticity; softening; hysteresis; NDE1 Introduction
Rocks and concrete are known to have a strong nonlinear behaviour. Quasistatic compression or traction tests show a nonlinear stressstrain 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).
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 stressstrain relationship (6). Another approach is related to the PreisachMayergoyz model, which is based on a discrete representation of hysteresis (7); (4); (8). The softratchet 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 softratchet 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 softratchet 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 ClausiusDuhem 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 finitestrain 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 softratchet 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, selfgravitation 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 secondorder tensor defined by (see e.g. (14); (13); (15))
(1) 
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
(2) 
In particular, the material derivative of the deformation gradient satisfies
(3) 
where is the velocity field. The conservation of mass implies
(4) 
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
(5) 
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:
(6) 
where is the strainrate tensor. The second principle of thermodynamics reads
(7) 
where is the specific entropy.
2.2 The model
Preliminaries.
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
(8) 
where is the absolute temperature. Multiplying (8) by , the local equations of thermodynamics (6)(7) yield the ClausiusDuhem inequality
(9) 
for all state and all evolution . The lefthand 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
(10) 
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 OgdenRoxburgh 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):
(11)  
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 CauchyGreen 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 CauchyGreen tensor . The material derivative of the strain tensor is . For any secondorder tensors , and , we recall that
(12)  
Therefore, the ClausiusDuhem inequality (9) with the substitutions (11) reduces to
(13) 
where the hyperelastic stress
(14) 
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 ClausiusDuhem inequality (13) for all yields the constitutive law
(15) 
where the hyperelastic stress is defined in (14).
Now, the ClausiusDuhem inequality (13) reduces to , for all state and all . Therefore, is either dependent on or equal to zero. We choose the simplest nontrivial dependence:
(16) 
where J/m^{3} 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
(17) 
Otherwise (), the internal variable satisfies , i.e. where
(18) 
In this case, the internal variable is instantaneously modified when the strain varies: no slow dynamics occurs.
The previous choice ensures that the ClausiusDuhem inequality is satisfied, independently of the sign of . Indeed, with the assumption (16), the dissipation per unit volume in the material (13) is
(19) 
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
(20) 
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
(21)  
(22) 
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
(23) 
where and are specified by (20) and (21)(22), respectively. The expression of the hyperelastic stress is specified by (14) if the right CauchyGreen 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
(24)  
In particular, the conservation of mass (4) rewrites as . The hyperelastic stress satisfies (14), where
(25) 
with the tensor derivatives (13)
(26)  
Thus, the following substitution
(27) 
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 MooneyRivlin model (18)
(28) 
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 GreenLagrange strain tensor (see e.g. (15)). An example of strain energy density in terms of the invariants of is the Murnaghan’s law (19)
(29) 
where are the Lamé parameters and are the Murnaghan coefficients. The latter are thirdorder 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 :
(30)  
Infinitesimal strain.
The GreenLagrange strain tensor is linearised with respect to the displacement:
(31) 
where is the infinitesimal strain tensor. Murnaghan’s law is used and the expression of the first PiolaKirchhoff stress tensor is linearised with respect to the coordinates of :
(32) 
The equations of motion (23) reduce to
(33) 
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:
(34) 
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):
(35) 
where is the Young’s modulus and are higherorder elastic constants. When and are zero, Hooke’s law of linear elasticity
(36) 
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,
(37) 
where . By identification with Landau’s law (35), the parameters can be expressed in terms of the Lamé and Murnaghan parameters:
(38) 
A similar calculus can be performed with the MooneyRivlin 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
(39)  
where is a time constant and are the Fourier coefficients of the normalized strain energy .
The solution of the ordinary differential equation (39) is
(40)  
The first term in (40) decreases exponentially in time with constant . The second term is the steadystate term, which oscillates at the frequency around its average value
(41) 
where .
In the case of Landau’s law (35), the nonzero Fourier coefficients are given in table 1. At small strain amplitudes, and , the highorder 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:
(42) 
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.
In the 1D case (34), the speed of sound is
(43) 
If the material is linearelastic without slow dynamics, the speed of sound reduces to . It is easier for the analysis to introduce the elastic modulus and its variation
(44) 
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
(45) 
which reduces to if and equal zero. The average of over a period of forcing is deduced from (40) and (45):
(46) 
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) 

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 steadystate 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 stressstrain relationship is outlined.
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.
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
(47) 
which vanishes at high , low and high frequency , and low and high . The maximum value reached by the steadystate solution is
(48) 
The strain amplitude for which the material is destroyed satisfies :
(49) 
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
(50) 
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
(51) 
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.
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.
Dissipation.
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 ClausiusDuhem inequality. Also, one can observe that no dissipation occurs if , which corresponds to the curve .
4 Conclusion
A new model for the dynamic behaviour of solids is proposed. The following features are common with the softratchet 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 nonclassical effects is required.
In comparison with the softratchet model, several differences can be outlined:

the new model satisfies the second principle of thermodynamics;

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

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 longtime 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 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 stressstrain relationship.
Pseudoelastic 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
(52) 
is proposed in (16). From (18), one deduces the expression of the internal variable
(53) 
This expression satisfies if . In particular, along the primary loading path. In the case of the endpoint memory phenomenon which is observed in rocks (1), the pseudoelastic 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 partiallyimmersed 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 quasistatically, 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)
(54) 
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
(55)  
Due to the geometry of the meniscus and the law of hydrostatics, one has
(56) 
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 softratchet model
Thermodynamical analysis.
The softratchet 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 softratchet model introduces a concentration of activated defects , which modifies the stress according to
(57) 
This constitutive law is the same as (15) with the softening function (20). In one space dimension, the strain rate satisfies , where . The ClausiusDuhem inequality (9) rewrites as
(58) 
for all state and all evolution. Due to the constitutive law (57), the specific internal energy must satisfy
(59) 
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 ClausiusDuhem inequality (58) implies
(60) 
for all and all .
In the softratchet model, the evolution equation for has the form
(61) 
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
(62) 
where is a stress and is the value of at zero stress. This expression is modified in (10) to ensure :
(63) 
Injecting (61) in (60) yields the condition
(64) 
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 softratchet 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
(65) 
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
(66) 
It rewrites as a cubic equation:
(67) 
which may have multiple solutions.
When using Cardano’s method, the discriminant
(68) 
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)):
(69) 
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 softratchet 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.
References
 Guyer RA, Johnson PA. 1999 Nonlinear mesoscopic elasticity: Evidence for a new class of materials. Phys. Today 52, 30–36.
 TenCate JA. 2011 Slow dynamics of earth materials: An experimental overview. Pure Appl. Geophys.