Nicolas Favrie
email: nicolas.favrie@univ-amu.fr

Bruno Lombard
email: lombard@lma.cnrs-mrs.fr

Cédric Payan
email: cedric.payan@univ-amu.fr

Fast and slow dynamics in a nonlinear elastic bar excited by longitudinal vibrations

Keywords: Nonlinear acoustics, time-dependent materials, viscoelasticity, acoustic conditioning, numerical methods, hyperbolic system

The dynamics of heterogeneous materials, like rocks and concrete, is complex. It includes such features as nonlinear elasticity, hysteresis, and long-time relaxation. This dynamics is very sensitive to microstructural changes and damage. The goal of this paper is to propose a physical model describing the longitudinal vibrations in heterogeneous material, and to develop a numerical strategy to solve the evolution equations. The theory relies on the coupling of two processes with radically different time scales: a fast process at the frequency of the excitation, governed by nonlinear elasticity and viscoelasticity, and a slow process, governed by the evolution of defects. The evolution equations are written as a nonlinear hyperbolic system with relaxation. A timedomain numerical scheme is developed, based on a splitting strategy. The features observed by numerical simulations show qualitative agreement with the features observed experimentally by Dynamic Acousto-Elastic Testing.

Introduction

Understanding the mechanisms of acoustic nonlinearity in heterogeneous materials is an object of intensive studies [START_REF] Guyer | Nonlinear mesoscopic elasticity: Evidence for a new class of materials[END_REF][START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF][START_REF] Guyer | Nonlinear Mesoscopic Elasticity: The Complex Behaviour of Rocks, Solis, Concrete[END_REF][START_REF] Lebedev | A unified model of hysteresis and long-time relaxation in heterogeneous materials[END_REF]. Experimental evidence has shown that media such as rocks and concrete possess an anomalously strong acoustic nonlinearity, which is of great importance for the description of ultrasonic phenomena including damage diagnostics. Besides the widely-studied nonlinear and hysteretic stress-strain relation [START_REF] Kadish | Frequency spectra of nonlinear elastic pulse-mode waves[END_REF], a long-time relaxation is also reported by most of the authors [START_REF] Ten Cate | Slow dynamics in the nonlinear elastic response of Berea sandstone[END_REF][START_REF] Tencate | Universal slow dynamics in granular solids[END_REF]. This slow dynamics is typically observed in experiments of softening / hardening [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF][START_REF] Riviere | Pump and probe waves in dynamic acousto-elasticity: Comprehensive description and comparison with nonlinear elastic theories[END_REF], where a bar is forced by a monochromatic excitation on a time interval, before the source is switched-off. During the experiment, the elastic modulus is measured by Dynamic Acousto-Elastic Testing methods. It can be observed that the elastic modulus decreases gradually (softening), and then it recovers progressively its initial value after the extinction of the source (hardening). The time scales of each stage is much longer than the time scale of the forcing, which justifies the term "slow dynamics".

The modelling of this slow dynamic effect has been investigated by many authors. An essentially phenomenological model is widely used for this purpose: the Preisach-Mayergoyz model (P-M model) based on the integral action of hysteretic elements connecting stress and strain [START_REF] Scalerandi | Preisach-Mayergoyz approach to fatigue-induced irreversibility[END_REF][START_REF] Scalerandi | Nonequilibrium and hysteresis in solids: disentangling conditioning from nonlinear elasticity[END_REF][START_REF] Lebedev | A unified model of hysteresis and long-time relaxation in heterogeneous materials[END_REF]. This model initially arose from the theory of magnetism, where the "hysteron" has a clear physical significance. In elasticity, such a physical interpretation is not available. To overcome this limitation and to develop a rigorous theory, various authors have proposed alternative models based on clear mechanical concepts. To our knowledge, the first physical model of slow dynamics was described in [START_REF] Tencate | Universal slow dynamics in granular solids[END_REF], where the relaxation was related to the recovery of microscopic contact impeded by a smooth spectrum of energy barriers. This theory was extended in [START_REF] V. Aleshin | Microcontact-based theory for acoustics in microdamaged materials[END_REF][START_REF] V. Aleshin | Friction in unconforming grain contacts as a mechanism for tensorial stressstrain hysteresis[END_REF], and recently improved based on the analysis of inter-grain contacts and the resulting surface force potential with a barrier [START_REF] Lebedev | A unified model of hysteresis and long-time relaxation in heterogeneous materials[END_REF]. Another approach was followed in [START_REF] Pecorari | Adhesion and nonlinear scattering by rough surfaces in contact: beyond the phenomenology of the Preisach-Mayergoyz framework[END_REF], where the author shows that two rough surfaces interacting via adhesion forces yield dynamics similar to that of the fictitious elements of the Preisach-Mayergoyz space [START_REF] Pecorari | Adhesion and nonlinear scattering by rough surfaces in contact: beyond the phenomenology of the Preisach-Mayergoyz framework[END_REF].

Here, we present an alternative mechanical description of slow dynamics based on the works of Vakhnenko and coauthors [START_REF] Vakhnenko | Strain-induced kinetics of intergrain defects as the mechanism of slow dynamics in the nonlinear resonant response of humid sandstone bars[END_REF][START_REF] Vakhnenko | Soft-ratchet modeling of end-point memory in the nonlinear resonant response of sedimentary rocks[END_REF], where the following scenario is proposed:

• the Young's modulus E varies with time. One can write E(g), where g is a time-dependent concentration of defects. It is closely related to the notion of damage in solids mechanics. But contrary to what happens in this irreversible case, where g strictly increases with time, the evolution of g is reversible. Waiting a sufficiently long time, the initial material properties are recovered;

• at equilibrium, stress σ yields a concentration of defects g σ . The dependence of g σ with respect to σ is monotonic;

• out of equilibrium, relaxation times are required for g to reach g σ . Whether g < g σ (increase in the number of defects) or g > g σ (decrease in the number of defects), Vakhnenko et al state that the time scales differ. The argument is given in section III of [START_REF] Vakhnenko | Soft-ratchet modeling of end-point memory in the nonlinear resonant response of sedimentary rocks[END_REF]: "there are various ways for an already existing crack in equilibrium to be further expanded when surplus tensile load is applied. However, under compressive load a crack, once formed, has only one spatial way to be annhilated or contracted". In both cases, these relaxation times are much longer than the time scale of the excitation, which explains the slow dynamics.

Comparisons with experimental data are given in section V of [START_REF] Vakhnenko | Soft-ratchet modeling of end-point memory in the nonlinear resonant response of sedimentary rocks[END_REF], where the authors reproduced experiments done on Berea sandstone [START_REF] Ten Cate | Slow dynamics in the nonlinear elastic response of Berea sandstone[END_REF]. One current weakness is that no micro-mechanical description of the involved defects has been proposed so far. A possible analogy may be found with populations of open / closed cracks filled with air, equivalent to a population of bubbles that relax towards an equilibrium state, depending on the applied stress [START_REF] Emelianov | Nonlinear dynamics of a gas bubble in an incompressible elastic medium[END_REF][START_REF] Gavrilyuk | Generalized vorticity for bubbly liquid and dispersive shallow water[END_REF]. In counterpart, one attractive feature of Vakhnenko's model is that it combines hyperbolic equations and relaxation terms, which constitutes a sound basis of physical phenomena [START_REF] Godunov | Elements of Continuum Mechanics and Conservation Laws[END_REF].

The present paper is a contribution to the theoretical analysis of this model and to its practical implementation to describe wave motion in damaged media. First, we point out that no mechanisms prevents the concentration of defects from exceeding 1, which is physically unrealistic.

We fix this problem by proposing another expression for the equilibrium concentration. Second, the Stokes model describing viscoelasticity behaviour in [START_REF] Vakhnenko | Soft-ratchet modeling of end-point memory in the nonlinear resonant response of sedimentary rocks[END_REF] poorly describes the attenuation in real media, and it is badly suited to time-domain simulations of wave propagation. Instead, we propose a new nonlinear version of the Zener model. This viscoelastic model degenerates correctly towards a pure nonlinear elasticity model when attenuation effects vanish. Moreover, the usual Zener model in the linear regime is recovered [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF]. In practice, this model only requires one physical parameter under the assumption of constant quality factor. Third, hyperbolicity is analyzed. Depending on the chosen model of nonlinear elasticity, a real sound speed may obtained only on a finite interval of strains; this is true in particular with the widely-used Landau's model.

The main effort of Vakhnenko et al was devoted to the construction of a model of slow dynamics. The resolution of the involved equations was quite rudimentary and not satisfying. Indeed, the equilibrium concentration of defects g σ was assumed to be known and was imposed (eq [START_REF] Landau | Theory of Elasticity[END_REF] in [START_REF] Vakhnenko | Soft-ratchet modeling of end-point memory in the nonlinear resonant response of sedimentary rocks[END_REF]), while it depends on σ. But treating the full coupled nonlinear equations is out of reach of a semi-analytical approach, which explains the strategy of these authors. On the contrary, we propose here a numerical method to integrate the full system of equations, involving the nonlinear elasticity, the hysteretic terms of viscoelasticity, and the slow dynamics. Due to the existence of different time scales, a splitting strategy is followed, ensuring the optimal time step for integration. The full system is split into a propagative hyperbolic part (resolved by a standard scheme for conservation laws) and into a relaxed part (resolved exactly).

Our numerical model is very modular. The various bricks (nonlinear elasticity, viscoelasticity, slow dynamics) can be incorporated easily. Numerical tests validate each part separately. When all the whole bricks are put together, typical features of wave motion in damaged media are observed. The softening / hardening experiments are qualitatively reproduced.

Physical modeling

In this section, we write the basic components describing the wave motion in a 1D material with damage. The fundations rely on linear elastodynamics, whose equations are recalled in section 2.1. Then, the soft-ratchet model of Vakhnenko and coauthors is introduced and enhanced in section 2.2. The fast dynamics is described in section 2.3, where various known models of nonlinear elasticity are presented, and a nonlinear model of viscoelasticity is proposed. This latter degenerates correctly in the limit cases of linear elasticity or null attenuation.

Linear elastodynamics

In the case of small deformations, the propagation of 1D elastic waves can be described by the following system [START_REF] Achenbach | Wave Propagation in Elastic Solids[END_REF]:

                 ∂v ∂t - 1 ρ ∂σ ∂x = γ, (1a)
∂ε ∂t - ∂v ∂x = 0, ( 1b
)
where t is the time, x is the spatial coordinate, γ is a forcing term, u is the displacement, v = ∂u ∂t is the velocity, ε = ∂u ∂x is the strain, and σ is the stress. The latter is a function of strain: σ = σ(ε).

In the linear case, Hooke's law writes σ = E ε, where E is the Young's modulus, which is assumed to be constant over time. In the particular case where γ is a Dirac source at x s with time evolution G(t), then the exact solution of (1) is straightforward

ε = - sgn(x -x s ) 2 c 2 G t - |x -x s | c , ( 2
)
where sgn is the sign distribution, and c = 1 ρ ∂σ ∂ε ≡ E/ρ is the speed of sound.

The goal of the forthcoming sections is to extend the model (1) in three ways:

• time variations of E due to the stress;

• nonlinear Hooke's law;

• hereditary effects (viscoelasticity).

The time scales for the first effect (variation of E) are much greater than for the second and third effect. This is consequently referred to as slow dynamics.

Slow dynamics: soft-ratchet model

Here we follow the approach taken from [START_REF] Vakhnenko | Strain-induced kinetics of intergrain defects as the mechanism of slow dynamics in the nonlinear resonant response of humid sandstone bars[END_REF][START_REF] Vakhnenko | Soft-ratchet modeling of end-point memory in the nonlinear resonant response of sedimentary rocks[END_REF] with some modifications. The slow dynamics of the medium is assumed to rely on the concentration of activated defects g, which varies with σ. In the lowest approximation, the Young's modulus is written:

E = 1 - g g cr E + , (3)
where g cr and E + are the critical concentration of defects and the maximum possible value of Young's modulus, respectively (figure 1-(a)). The following constraints hold: The concentration g is assumed to evolve to its stress-dependent equilibrium value g σ at a rate f r if g > g σ (restoration), or f d if g < g σ (destruction). This mechanism can be modeled by the ordinary differential equation

0 ≤ g ≤ g cr ≤ 1. (4)
dg dt = -( f r H(g -g σ ) + f d H(g σ -g)) (g -g σ ), ( 5
)
where H is the Heaviside step distribution. The frequencies f r and f d differ substantially:

f r ≪ f d ≪ f c , (6)
where f c is a typical frequency of the excitation. It remains to define the evolution of g σ with σ. In [START_REF] Vakhnenko | Strain-induced kinetics of intergrain defects as the mechanism of slow dynamics in the nonlinear resonant response of humid sandstone bars[END_REF][START_REF] Vakhnenko | Soft-ratchet modeling of end-point memory in the nonlinear resonant response of sedimentary rocks[END_REF], the authors propose the expression

g σ = g 0 exp(σ/σ), σ = kT υ , (7)
where g 0 is the unstrained equilibrium concentration of defects, k is the Boltzmann constant, T is the temperature, and υ is a typical volume accounting for a single defect. If σ > σ ln g cr /g 0 , then g σ > g cr ; in this case, the concentration may evolve to g > g cr due to equation [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF], which contradicts the second assumption in (4). To remove this drawback and to build a physically realistic expression of g σ , we enforce (4) together with the following requirements:

                               0 ≤ g σ < g cr , (8a)
g σ (0) = g 0 , (8b) lim σ→-∞ g σ = 0, (8c)
lim σ→+∞ g σ = g cr , (8d)
∂g σ ∂σ > 0. ( 8e
)
The simplest smooth function satisfying [START_REF] Emelianov | Nonlinear dynamics of a gas bubble in an incompressible elastic medium[END_REF] is

g σ = g cr 2 1 + tanh σ -σ c σ , (9)
where the central stress is

σ c = σ tanh -1 1 -2 g 0 g cr . ( 10
)
Figure 2-(a) illustrates the two expressions of the stress-dependent equilibrium value g σ : the "exponential model" [START_REF] Davis | Simplified second-order Godunov-type methods[END_REF], and the "tanh model" ( 9)- [START_REF] Godunov | Elements of Continuum Mechanics and Conservation Laws[END_REF]. The numerical values are g 0 = 0.3 and σ = 10 5 Pa. The two expressions are the same at null stress. But for tractions greater than 230 kPa, the value of g σ deduced from (7) exceeds 1, leading to non-physical negative Young's modulus. Figure 2-(b) illustrates the influence of σ in [START_REF] Gavrilyuk | Generalized vorticity for bubbly liquid and dispersive shallow water[END_REF]. As σ decreases, g σ may evolve more easily towards the extreme values 0 and g cr , and hence the damage may increase thanks to [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF].

Fast dynamics: nonlinear viscoelasticity

Nonlinear elasticity.

The stress-strain relation is given by a smooth function

s ≡ s(ǫ, K, p), ( 11
)
where s is the stress, ǫ is the strain, K is a stiffness, and p is a set of parameters governing the nonlinearity. No pre-stress is considered; K is the slope of s at the origin; lastly, s is homogeneous of degree 1 in K. In other words, s satisfies the following properties:

s(0, K, p) = 0, ∂s ∂ǫ (0, K, p) = K, s(ǫ, αK, p) = α s(ǫ, K, p). (12)
Three models of nonlinear elasticity (11) satisfying ( 12) are now given and illustrated in figure 3.

Model 1. This model is from [START_REF] Vakhnenko | Soft-ratchet modeling of end-point memory in the nonlinear resonant response of sedimentary rocks[END_REF] and mimics the Lennard-Jones potential describing the interaction between a pair of neutral atoms:

s(ǫ, K, p) = K d r -a               1 1 + ǫ d a+1 - 1 1 + ǫ d r+1               , p = (r, a, d) T . ( 13
)
The nonlinear parameters are the repulsion and attraction coefficients r and a (0 < a < r). The strain is bounded below by the maximal allowable closure d. The function [START_REF] Hamilton | Nonlinear Acoustics[END_REF] has an extremal point ǫ c > 0, and then it decreases asymptotically towards 0 when ǫ > ǫ c (figure 3-(a)).

Model 2. A third-order Taylor expansion of the model 1 (13) yields

s(ǫ, K, p) = K ǫ 1 - 1 2 (r + a + 3) ǫ d + 1 6 r 2 + ra + a 2 + 6r + 6a + 11 ǫ d 2 , p = (r, a, d) T . (14)
The nonlinear parameters are the same than in model 1. But contrary to what happened in model 1, the function ( 14) is a strictly monotonically increasing function without extremal point (figure 3-(a)). Moreover, the strain is not bounded below. Model 3. The most widely used law in ultrasonic NonDestructive Testing is the so-called Landau's model [START_REF] Landau | Theory of Elasticity[END_REF]

s(ǫ, K, p) = K ǫ 1 -β ǫ -δ ǫ 2 , p = (β, δ) T . ( 15
)
The Viscoelasticity.

To incorporate attenuation, the following criteria are used as a guideline:

C 1 : when the viscous effects are null, the nonlinear elasticity must be recovered (11);

C 2 : when a linear stress-strain relation holds, it is necessary to recover the standard linear solid model (or generalized Zener model), which accurately represents the behavior of usual solids [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF]. For this purpose, a system with N Zener elements connected in parallel is considered (figure 4).

The total stress acting on the system is

σ = N ℓ=1 σ 1ℓ = N ℓ=1 (σ 2ℓ + σ 3ℓ ), (16)
where the index 1 refers to the springs in series, and indices 2-3 refer to the springs and dashpots in parallel. The strain ε is

ε = ε 1ℓ + ε 2ℓ , ℓ = 1, • • • , N. ( 17
)
The index 1 springs satisfy nonlinear stress-strain relations [START_REF] Guyer | Nonlinear mesoscopic elasticity: Evidence for a new class of materials[END_REF] with stiffnesses K 1ℓ . The parameters p governing the nonlinearity (for instance β and δ in model 3 [START_REF] Kadish | Frequency spectra of nonlinear elastic pulse-mode waves[END_REF]) are assumed to be constant and identical for each element. The index 2 springs satisfy linear stress-strain relations with stiffnesses K 2ℓ . Lastly, the dashpots satisfy linear Maxwell laws with coefficients of viscosity η ℓ . These laws are summed up as follows:

                 σ 1ℓ (ε 1ℓ ) = s(ε 1ℓ , K 1ℓ , p), (18a)
σ 2ℓ (ε 2ℓ ) = s(ε 2ℓ , K 2ℓ , 0), (18b)
σ 3ℓ (ε 2ℓ ) = η ℓ ∂ε 2ℓ ∂t . (18c)
To determine the parameters K 1ℓ , K 2ℓ and η ℓ , the relaxation times τ σℓ , τ εℓ and the relaxed modulus E R are introduced:

τ σℓ = η ℓ K 1ℓ + K 2ℓ , τ εℓ = η ℓ K 2ℓ , E R N = K 1ℓ K 2ℓ K 1ℓ + K 2ℓ . ( 19
)
On the one hand, a procedure is given in Appendix A to compute the relaxation times in terms of the quality factor Q. On the other hand, E R is related to the unrelaxed Young's modulus E [START_REF] V. Aleshin | Friction in unconforming grain contacts as a mechanism for tensorial stressstrain hysteresis[END_REF] and to the relaxation times previsously determined (see [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF]):

E R = N N ℓ=1 τ εℓ τ σℓ E. ( 20
)
Once τ σℓ , τ εℓ and E R are determined, inverting [START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF] provides the values of the viscoelastic model in terms of relaxed modulus and relaxation times (ℓ = 1, • • • , N):

K 1ℓ = τ εℓ τ σℓ E R N , K 2ℓ = τ εℓ τ εℓ -τ σℓ E R N , η ℓ = τ 2 εℓ τ εℓ -τ σℓ E R N . (21)
From ( 20) and ( 21), it follows that the viscoelastic parameters depend indirectly on the Young's modulus E, and thus depend on g. In other words, the proposed model of viscoelasticity evolves with the concentration of defects and thus with the applied stress.

In the inviscid case, the stress-strain relation deduced from ( 16)-( 18) makes it possible to recover the nonlinear elasticity [START_REF] Guyer | Nonlinear mesoscopic elasticity: Evidence for a new class of materials[END_REF], whatever the number N of relaxation mechanisms:

σ = s(ε, E, p). ( 22
)
This property is proven in Appendix B.

Mathematical modeling

In this section, the basic components describing wave motion in damaged media are put together and analysed. Section 3.1 collects the various mechanisms (nonlinear elastodynamics, slow dynamics, hysteresis) into a single system of first-order equations. Two important properties of this system are addressed in section 3.2: hyperbolicity (finite sound velocity) and decrease in energy.

First-order system

The conservation of momentum (1a) writes

∂v ∂t = 1 ρ ∂σ ∂x + γ, ( 23
)
where γ is a forcing term, and σ is given by ( 16). The hypothesis of small deformations (1b)

gives

∂ε ∂t = ∂v ∂x . (24)
Lastly, manipulations on ( 16), ( 17) and (18c) yield

∂ε 1ℓ ∂t = ∂v ∂x + σ 2ℓ (ε -ε 1ℓ ) -σ 1ℓ (ε 1ℓ ) η ℓ , ℓ = 1, • • • , N. (25)
In [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF], ε 1ℓ takes the place of the memory variables proposed in [START_REF] Lombard | Numerical modeling of transient two-dimensional viscoelastic waves[END_REF] and is better suited to nonlinear elasticity. Putting together ( 23)-( 25) and the relaxation equation ( 5) leads to the firstorder system of N + 3 evolution equations

                                         ∂v ∂t - 1 ρ ∂σ ∂x = γ, (26a)
∂ε ∂t - ∂v ∂x = 0, ( 26b
)
∂ε 1ℓ ∂t - ∂v ∂x = σ 2ℓ (ε -ε 1ℓ ) -σ 1ℓ (ε 1ℓ ) η ℓ , ℓ = 1, • • • , N, ( 26c
)
dg dt = -( f r H(g -g σ ) + f d H(g σ -g)) (g -g σ ). ( 26d
)
To close the system (26), the following equations are recalled:

• The total stress σ in (26a) depends on ε 1ℓ via ( 16), (18a), and a nonlinear law [START_REF] Guyer | Nonlinear mesoscopic elasticity: Evidence for a new class of materials[END_REF]:

σ = N ℓ=1 s(ε 1ℓ , K 1ℓ , p). (27)
• The stress components σ 1ℓ and σ 2ℓ in (26c) depend on the stifnesses K 1ℓ and K 2ℓ (18a) and (18b). The latter, as well as the viscosity coefficients η ℓ , depend on the Young modulus E via ( 20)-( 21), and thus on g:

E = 1 - g g cr E + . (28)
• The equilibrium value of the defect concentration g σ in (26d) satisfies ( 9) and [START_REF] Godunov | Elements of Continuum Mechanics and Conservation Laws[END_REF]:

g σ = g cr 2 1 + tanh σ -σ c σ . ( 29
)
The system [START_REF] Riviere | Pump and probe waves in dynamic acousto-elasticity: Comprehensive description and comparison with nonlinear elastic theories[END_REF], together with equations ( 27)-( 29), generalizes the standard equations of linear elastodynamics [START_REF] Achenbach | Wave Propagation in Elastic Solids[END_REF]. It accounts for softening / recovering of Young's modulus, nonlinearity and viscoelasticity.

For the sake of clarity, the vector of N + 3 variables is introduced

U = (v, ε, ε 11 , • • • , ε 1N , g) T . (30)
Then the system (26) can be put in the form

∂ ∂t U + ∂ ∂x F(U) = R(U) + Γ. (31)
The flux function F, the relaxation term R, and the forcing Γ are

F(U) = - σ ρ , -v, -v, • • • , -v, 0 T , R(U) = (0, 0, ∆ 1 , • • • , ∆ N , -( f r H(g -g σ ) + f d H(g σ -g)) (g -g σ )) T , Γ = (γ, 0, • • • , 0, 0) T , (32)
where

∆ ℓ = σ 2ℓ (ε -ε 1ℓ ) -σ 1ℓ (ε 1ℓ ) η ℓ . ( 33
)
To conclude, let us consider the limit-case where the viscoelastic attenuation is neglected. In this case, equation [START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF] states that the stress-strain relations degenerate rigorously towards pure nonlinear elasticity, whatever N.

Properties

Hyperbolicity is a crucial issue in wave problems -physically, mathematically, and numerically. It amounts to saying that there exists a real and finite sound velocity c. This property was analysed in [START_REF] Ndanou | Criterion of hyperbolicity in hyperelasticity in the case of the stored energy in separable form[END_REF] for a particular nonlinear stress-strain relation in 3D. In 1D, it reduces to a simpler case detailed as follows. Let us define the sound speed c by

c 2 = N ℓ=1 c 2 ℓ = 1 ρ N ℓ=1 ∂σ 1ℓ ∂ε 1ℓ . ( 34
)
The system ( 31) is hyperbolic if and only if c 2 > 0 in [START_REF] Whitham | Linear and Nonlinear Waves[END_REF]. The proof, as well as sufficient conditions on the strain to ensure hyperbolicity, is given in Appendix B. From [START_REF] Whitham | Linear and Nonlinear Waves[END_REF], the local elastic modulus M can be deduced:

M = ρ c 2 = N ℓ=1 ∂σ 1ℓ ∂ε 1ℓ . (35)
Note that the Stokes viscoelastic model used in [START_REF] Vakhnenko | Soft-ratchet modeling of end-point memory in the nonlinear resonant response of sedimentary rocks[END_REF] introduces a term ∂ 2 v ∂x 2 in the right-hand side of (26c). This Laplacian term destroys the hyperbolic character of the system [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction[END_REF]. The viscoelastic model used here has therefore better mathematical properties. Now let us examine the spectrum of the relaxation function in [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction[END_REF]. Let us consider linear stress-strain relations. The parameters K 1ℓ , K 2ℓ and η ℓ are "freezed" in ( 20)-( 21), so that they do not depend on g via E (3). Then, the eigenvalues of the Jacobian matrix J = ∂R ∂U are

Sp(J) = 0 2 , -f ξ , - K 1ℓ + K 2ℓ η ℓ = 0 2 , -f ξ , - 1 τ σℓ , ℓ = 1, • • • , N, (36)
(see [START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF]), with

f ξ = f r if g > g σ , f ξ = f d if g < g σ , f ξ = 0 else. The proof is detailed in Appendix C
. Two observations can be made:

• J is definite-negative if the relaxation frequencies τ σℓ are positive. The latter parameters are deduced from an optimization process based on the quality factor (Appendix A). To ensure the energy decrease, it is therefore crucial to perform nonlinear optimization with constraint of positivity.

• The optimization procedure detailed in Appendix A is performed on the frequency range [ f min , f max ] surrounding the excitation frequency f c . These frequencies satisfy

f min ≈ 1 max τ σℓ < f c < f max ≈ 1 min τ σℓ . (37)
In (37), ≈ are replaced by equalities if a linear optimisation is used [START_REF] Lombard | Numerical modeling of transient two-dimensional viscoelastic waves[END_REF]. From [START_REF] Ben Jazia | Wave propagation in a fractional viscoelastic Andrade medium: diffusive approximation and numerical modeling[END_REF], it follows the spectral radius of J

̺(J) = 1 min τ σℓ ≫ f ξ , (38)
so that the system (31) is stiff.

Numerical modeling

In this section, a numerical strategy is proposed to integrate the first-order equations [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction[END_REF].

For the sake of efficiency, a splitting approach is followed in section 4.1. The original equations are splitted into two parts, solved successively: a propagative part (section 4.2) and a relaxation part (section 4.3).

Splitting

To integrate [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction[END_REF], a uniform spatial mesh ∆x and a variable time step ∆t (n) ≡ ∆t are introduced. An approximation U n i of the exact solution U(

x i = i ∆x, t n = t n-1 + ∆t) is sought. A first
strategy is to discretize explicitly the non-homogeneous system [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction[END_REF]. But numerical stability implies a bound of the form

∆t ≤ min ∆x c max , 2 ̺(J) , (39)
where c max = max c n i is the maximal sound velocity at time t n , and ̺(J) is the spectral radius of the Jacobian of the relaxation term. As deduced from (38), the second bound in (39) is penalizing compared with the standard CFL condition ∆t ≤ ∆x/c max .

Here we follow another strategy: equation ( 31) is split into a hyperbolic step

∂ ∂t U + ∂ ∂x F(U) = 0 (40)
and a relaxation step

∂ ∂t U = R(U) + Γ. (41)
The discrete operators associated with the discretization of (40) and (41) are denoted H h and H r , respectively. The second-order Strang splitting is used, solving successively (40) and (41) with adequate time increments:

                             U (1) i = H r ∆t 2 U n i , (42a)
U (2) i = H h (∆t) U (1) i , (42b)
U n+1 i = H r ∆t 2 U (2) i . (42c)
Provided that H h and H r are second-order accurate and stable operators, the time-marching (42)

gives a second-order accurate approximation of the original equation (31) [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF].

Hyperbolic step

The homogeneous equation ( 40) is solved by a conservative scheme for hyperbolic systems [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]]

U n+1 i = U n i - ∆t ∆x F i+1/2 -F i-1/2 . ( 43
)
Many sophisticated schemes can be used for this purpose [START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF]. For the sake of simplicity and robustness, the Godunov scheme is used here. The numerical flux function F i+1/2 is computed using the Rusanov method [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction[END_REF]

F i+1/2 = 1 2 F(U n i+1 ) + F(U n i ) -λ n i+1/2 (U n i+1 -U n i ) , ( 44
)
where F is the flux function [START_REF] Vakhnenko | Strain-induced kinetics of intergrain defects as the mechanism of slow dynamics in the nonlinear resonant response of humid sandstone bars[END_REF], and the diffusion parameter λ n i+1/2 is given by the Davis approximation [7]

λ n i+1/2 = max c n i , c n i+1 . (45)
The Godunov scheme is first-order accurate and stable under the usual Courant-Friedrichs-Lewy

(CFL) condition ∆t = α ∆x c max , with α ≤ 1. ( 46
)

Relaxation step

Let us denote U = (ε, ε 11 , • • • , ε 1N ) and R the restriction of R(U) to the strain components (32)- [START_REF] Vakhnenko | Soft-ratchet modeling of end-point memory in the nonlinear resonant response of sedimentary rocks[END_REF]. The ordinary differential equation (41) can then be written

                           ∂v ∂t = γ, (47a)
∂ ∂t U = R(U), (47b)
dg dt = -( f r H(g -g σ ) + f d H(g σ -g)) (g -g σ ), (47c)
The viscoelastic parameters in the relaxation function R depend implicitly on g (see section 2.3), which complicates the resolution of (47a). However, one can take advantage of the scaling (6).

Indeed, ε and ε 1ℓ evolve much faster than g, so that the viscoelastic parameters K 1ℓ , K 2ℓ , η ℓ are almost constant on a time step. Consequently, they are "freezed" and the three equations in (47) can be solved separately.

The half-time step in the relaxation steps (42a)-(42c) is denoted by τ = ∆t 2 . One details the time-stepping from t n to the first intermediate step (42a); adaptation to the third intermediate step

(42c) is straightforward.
The first equation (47a) is integrated using the Euler method:

v n+1 i = v (1) i + ∆t γ(i, t n ). ( 48
)
To integrate the second equation (47b), a first-order Taylor expansion of R(U) is performed

∂ ∂t U ≈ R(0) + ∂R ∂U (0) U = J U, ( 49
)
where J is the Jacobian matrix (C.2); the nullity of stress at zero strain has been used (18a). Then (49) is solved exactly, leading to the relaxation operator

U (1) i = e J τ U n i ( 50
)
with the matrix exponential

e J τ =                              1 0 • • • 0 E 21 E 11 + E 21 1 -e - E 11 +E 21 η 1 τ e - E 11 +E 21 η 1 τ . . . . . . E 2N E 1N + E 2N 1 -e - E 1N +E 2N η N τ e - E 1N +E 2N η N τ                              . ( 51
)
Lastly, the third equation ( 47c) is solved exactly. The grid value g σi is evaluated thanks to [START_REF] Gavrilyuk | Generalized vorticity for bubbly liquid and dispersive shallow water[END_REF].

Setting

f ξ =        f r if g n i ≥ g n σi , f d if g n i < g n σi , (52)
leads to

g (1) i = g n σi + g n i -g n σi e -f ξ τ . (53)
The integrations (50), ( 48) and ( 53) are unconditionally stable. As a consequence, the splitting (42) is stable under the CFL condition (46).

Summary of the algorithm

The numerical method can be divided in two parts:

1. initialisation

• bulk modulus ρ, Young's modulus E = E 0 = ρ c 2 ∞ ;
• soft-ratchet coefficients g cr = 1, g = g 0 , f r , f d , σ;

• maximum Young's modulus E + (3)

• nonlinear coefficients (e.g. β and δ in (15);

• quality factor Q, frequency range of optimization [ f min , f max ], number of relaxation mechanisms N;

• optimization of the viscoelastic coefficients (Appendix A);

2. time-marching t n → t n+1 , x i = i ∆x (n = 0, • • • , N t , i = 1, • • • , N x )
• physical and numerical parameters -Young's modulus E (3), viscoelastic parameters E R [START_REF] Lombard | Numerical modeling of transient two-dimensional viscoelastic waves[END_REF], K 1ℓ , K 2ℓ and η ℓ (21);

partial stresses σ 1ℓ (18a) and total stress σ (16);

sound velocity c (34) and (B.6), maximal velocity c max ;

time step ∆t (46);

• relaxation step H r (42a)
strains (50) and (51);

velocity v (48);

concentration of defects at equilibrium g σ (9) and out of equilibrium g (53);

• hyperbolic step H h (42b)
coefficient λ i+1/2 of Davis (45);

computation of the flux F (32), e.g., by the Rusanov flux F i+1/2 (44);

time-marching of the conservative scheme (43);

• relaxation step H r (42c). The physical parameters are detailed in table 1. Depending on the test, some of these parameters are modified. In the limit-case of linear elasticity, the sound velocity is c = E/ρ = 3280 m/s. The maximal CFL number is α = 0.95 in (46). The mesh size is ∆x = 4 10 -3 m. Depending on the test, two lengths of domain are considered. For each test, a receiver put at x r = 0.2 m stores the numerical solution at each time step.

Numerical experiments

The wave fields are excited by a punctual source at x s = 10 -2 m, with a central frequency f c = 10 kHz. Depending on the expression of the forcing γ in (26c), it is possible to deduce the magnitude of the maximal strain ε max emitted by the source in the limit-case of linear elasticity

(2):

ε max = 1 2 c 2 max G(t). ( 54
)
The Landau model for nonlinear elasticity is used [START_REF] Kadish | Frequency spectra of nonlinear elastic pulse-mode waves[END_REF]. The coefficient β is much smaller than δ. The critical value of strain that ensures hyperbolicity (B.5) is ε c = 3.08 10 -4 . The viscoelastic effects are described by N = 4 relaxation mechanisms. The relaxation times τ σℓ and τ εℓ [START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF] are computed by optimization on the frequency range [ f min = f c /10, f max = f c × 10] (see Appendix A); they are given in table 2.

ℓ = 1 ℓ = 2 ℓ = 3 ℓ = 4 τ σℓ (s) 1.

Test 1: nonlinear elastodynamics

In the first test, the viscoelasticity is neglected, and the activation / restoration of defects is annihilated: f r = f d = 0 Hz. This test corresponds to the example 12 of [START_REF] Whitham | Linear and Nonlinear Waves[END_REF]. Our goal is to show typical features of wave propagation in purely nonlinear elastic media. The source is a monochromatic excitation:

ε max = 10 -5 ε max = 2.
G(t) = A sin(ω c t) H(t), ( 55
)
where A is the magnitude of the forcing, and ω c = 2 π f c . From ( 54) and (55), it is possible to estimate the maximal strain ε max emitted by the source in the linear elastic case. The domain of propagation is L x = 2 m long and is discretized onto 400 grid nodes.

Figure 5 displays the spatial evolution of ε after 400 time steps. For ε max = 10 -5 , almost no distorsion of the wave is seen. On the contrary, ε max = 2.0 10 -4 yields a high distorsion as the wave propagates. Shocks, as well as the attenuation due to the intersection of characteristic curves [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF], are observed.

Figure 6 displays the time evolution of the strain recorded at the receiver (vertical dotted line in Figure 5) for ε max = 2.0 10 -4 . The normalized amplitudes of the Fourier series decomposition show a typical feature of cubic nonlinear elasticity: the spectrum involves mainly odd harmonics [START_REF] Hamilton | Nonlinear Acoustics[END_REF].

Test 2: linear viscoelasticity

The goal of the second test is to validate the numerical modeling of attenuation. For this purpose, a linear stress-strain relation is chosen (β = δ = 0), and the activation / restoration of defects is still annihilated ( f r = f d = 0 Hz). Consequently, the system (26) simplifies into The domain of propagation is L x = 2 m long and is discretized onto 400 grid nodes. The time evolution of the source is a truncated combination of sinusoids with C 6 smoothness: The attenuation is slightly overestimated by the scheme, due to the numerical diffusion of the Godunov scheme. This numerical artifact can be fixed by choosing a higher-order scheme [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction[END_REF].

                             ∂v ∂t - 1 ρ ∂σ ∂x = γ, ( 56a
) ∂ε ∂t - ∂v ∂x = 0, ( 56b
)
∂ε 1ℓ ∂t - ∂v ∂x = K 2ℓ η ℓ (ε -ε 1ℓ ) - K 1ℓ η ℓ ε 1ℓ . ( 56c
G(t) =            4 m=1 a m sin (b m ω c t) if 0 ≤ t ≤ 1 f c , 0 otherwise, (57)

Test 3: softening / recovering

The goal of the third test is to illustrate the softening / recovering of the elastic modulus, and to validate the numerical modeling of this phenomenon. For this purpose, linear elasticity is assumed and the viscoelasticity is neglected (β = δ = 0, Q = +∞). Even if a linear stress-strain relation is used, the evolution problem ( 26) is nonlinear by virtue of (26d), ( 28) and [START_REF] Ten Cate | Slow dynamics in the nonlinear elastic response of Berea sandstone[END_REF]. Like in test 1, the source is monochromatic; but is is switched off after a time t * : As long as the source is switched on (0 < t < t * ), the equilibrium concentration of defects increases from the initial value g 0 up to g * = g(t * ). At the same time, the Young's modulus

G(t) = A sin(ω c t) (H(t) -H(t * )) . (58) (a)
decreases from E 0 to E * via (3).

For t > t * , the waves go out of the domain, and the elastodynamic fields vanish. From [START_REF] Ten Cate | Slow dynamics in the nonlinear elastic response of Berea sandstone[END_REF] and [START_REF] Godunov | Elements of Continuum Mechanics and Conservation Laws[END_REF], σ = 0 implies that the equilibrium concentration of defects becomes g σ = g 0 . As a consequence, the ordinary differential equation (ODE) (26d) describing the evolution of defects simplifies into

         dg dt = -f r (g -g 0 ), g(t * ) = g * . (59)
The solution of (59) is

g(t) = g 0 + (g * -g 0 ) e -f r (t-t * ) . (60)
Equation ( 60) is injected into (3), which gives the time evolution of the Young's modulus during the recovering process (t ≥ t * ):

E(t) = E 0 - 1 g cr (g * -g 0 ) e -f r (t-t * ) E + . (61)
The domain of propagation is L x = 0.4 m long and is discretized onto 100 grid nodes. The maximal strain is ε max = 10 -5 . Time integration is performed up to t = 460 ms. Pa or 10 7 Pa (the other parameters are those of table 1). According to the Vakhnenko's expression [START_REF] Davis | Simplified second-order Godunov-type methods[END_REF], these values correspond to spherical defects of radius 2.13 10 -10 m and 4.59 10 -10 m, respectively. In both cases, equilibrium has been reached at t * . The lower value of σ yields a greater variation of the elastic modulus. This property follows from (9): as σ decreases, the curve g → g σ stiffens and tend towards a Heaviside step function. Consequently, greater values of g σ are obtained when σ is smaller. This implies a greater evolution of g [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF], and hence of E

(3). 1). The lowest value of f r yields a greater variation of the elastic modulus. This is a consequence of the competition between restoration (with frequency f r ) and destruction (with frequency f d ). When f r is too low compared with f d , restoration has almost no time to occur during one period T = 1/ f c , and destruction plays a preponderant role.

Test 4: full model

The fourth and last test incorporates all the physical mechanisms of the model: nonlinear stress-strain law, viscoelasticity, activation / restoration of defects. The domain is L x = 0.4 m long and is discretized onto 100 grid nodes. The source is a monochromatic excitation (55).

Time integration is performed during 5 10 4 time steps. The fields are recorded at x r . First, nonlinear curves are obtained, which is a signature of the nonlinear stress-strain relation.

Second, ∆M increases with ε max : softening increases monotonically with the forcing. Third and last, loops are obtained if and only if viscoelasticity is incorporated (c-d). These three features are qualitatively similar to those obtained experimentally [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF][START_REF] Riviere | Pump and probe waves in dynamic acousto-elasticity: Comprehensive description and comparison with nonlinear elastic theories[END_REF].

Conclusion

We have proposed a one-dimensional model that captures the behavior of real media under longitudinal bar excitation, including the following features: nonlinear elasticity; softening / recovering of the elastic modulus; hysteretic evolution of the elastic modulus with the strain. The proposed model is very modular. It involves three different bricks which can be used also independently: see for instance the numerical experiments in section 5, in which are considered various combinations of elasticity, attenuation and slow dynamics. Experimentally, the parameters corresponding to each mechanisms can be identified separately:

• the measure of nonlinear elastic parameters is described in many books [START_REF] Guyer | Nonlinear mesoscopic elasticity: Evidence for a new class of materials[END_REF][START_REF] Hamilton | Nonlinear Acoustics[END_REF];

• the measure of the quality factor must be performed in the linear regime. See the reference book [START_REF] Bourbié | Acoustics of Porous Media[END_REF] for a description of an experimental protocol;

• lastly, measuring the parameters of the slow dynamics is detailed in many papers cited in the bibliography. The current challenge is to link the physical observations to the parameters of Vakhnenko's model. Our ambition, with the present paper, is to provide experi- menters with a tool for testing various sets of parameters, and hence testing the validity of Vakhnenko's model.

(a) (b) T¤¥T @ ¦ § ¦ w ¨ © ¦ w @ 7 ¨0H H 0.12 ε ¥T¤¥T ¥E¤¥T 0 E¤¥T T¤¥T ε ! "#$"5 "%$"5 0 %$"R #$"R 0 1 Q 2 (c) (d) & ' ( ' ) 0 1
A major interest of the numerical approach is the possibility to tackle with variable coefficients in space, which is representative of localized defects [START_REF] Pecorari | Forced nonlinear vibrations of a one-dimensional bar with arbitrary distributions of hysteretic damage[END_REF]. In particular, a random initial distribution of defects g 0 (x) can be considered straightforwardly.

Many improvements can be investigated, to mention but a few. More sophisticated models can be built quite naturally, considering for instance relaxation of the nonlinear coefficients p in (18a), or a nonlinear law in (18b). Concerning the numerical simulations, higher-order schemes (such as WENO schemes [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]) can easily be adapted to the proposed formulation. Lastly, theoretical analyses should be done to prove rigorously the well-posedness of the model and its thermodynamic properties.

Work is currently proceeding along two directions. First, numerical simulations are being done to recover quantitatively the experimental results of the litterature [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF][START_REF] Riviere | Pump and probe waves in dynamic acousto-elasticity: Comprehensive description and comparison with nonlinear elastic theories[END_REF]. Second, the extension of this model to 2D and 3D geometries is under progress.

Appendix A. Parameters of the viscoelastic model

Standard calculations on ( 16), ( 18) and ( 19) yield the reciprocal of the quality factor Q

[5] Q -1 (ω) =        N ℓ=1 ω (τ εℓ -τ σℓ ) 1 + ω 2 τ 2 σℓ        /        N ℓ=1 1 + ω 2 τ εℓ τ σℓ 1 + ω 2 τ 2 σℓ        . (A.1)
Optimizing Q -1 towards a given law (for instance a constant quality factor on a frequency range of interest [ f min , f max ]) provides a means to determine τ σℓ and τ εℓ [START_REF] Lombard | Numerical modeling of transient two-dimensional viscoelastic waves[END_REF]. Here an optimization with constraint is applied to ensure positive values of τ σℓ and τ εℓ , as required by the decrease in energy (see section 3.2). See [START_REF] Ben Jazia | Wave propagation in a fractional viscoelastic Andrade medium: diffusive approximation and numerical modeling[END_REF] for details about such an optimization.

= E R /ρ to c ∞ = E/ρ.
The reader is referred to [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF] for details about these quantities.

Lastly, the consistancy relation ( 22) is proven here. Null attenuation amounts to an infinite quality factor. Equation (A.1) implies that Q = +∞ is obtained if τ εℓ = τ σℓ . In this case, the viscoelastic coefficients [START_REF] Lombard | Numerical modeling of transient two-dimensional viscoelastic waves[END_REF] and ( 21) are

E R = E, K 1ℓ = E N , K 2ℓ = +∞, η ℓ = +∞. (A.2)
To get a bounded stress, (18c) implies ε 2ℓ = 0, and hence ε 1ℓ = ε for ℓ = 1, • • • , N [START_REF] Landau | Theory of Elasticity[END_REF]. Putting together the total stress ( 16), the nonlinear elasticity [START_REF] Guyer | Nonlinear mesoscopic elasticity: Evidence for a new class of materials[END_REF] and the homogeneity property in [START_REF] Guyer | Nonlinear Mesoscopic Elasticity: The Complex Behaviour of Rocks, Solis, Concrete[END_REF],

Appendix B. Analysis of hyperbolicity

The Jacobian A of f (32) is

A(U) =                             0 0 Φ 1 • • • Φ N 0 -1 0 0 • • • 0 0 -1 0 0 • • • 0 0 . . . . . . . . . . . . . . . -1 0 0 • • • 0 0 0 0 0 0 0 0                             , (B.1)
where

Φ ℓ = - 1 ρ ∂σ 1ℓ ∂ε 1ℓ . (B.
2)

The determinant of A writes

P A (λ) = -λ -λ 0 Φ 1 • • • Φ N -1 -λ 0 • • • 0 -1 0 -λ 0 . . . . . . . . . -1 0 -λ (B.
3)

The columns and lines are denoted by C j and L j , respectively. The following algebraic manipulations are performed successively:

(i) C 1 ← λ C 1 , (ii) C 1 ← C 1 -C j , with j = 2, • • • , N + 1,
which yields

λ P A (λ) = -λ -λ 2 - N ℓ=1 Φ ℓ 0 Φ 1 • • • Φ N 0 -λ . . . . . . 0 -λ , = (-1) N+1 λ N+2        λ 2 + N ℓ=1 Φ ℓ        . (B.4)
It follows that the eigenvalues are 0 (with multiplicity N + 1) and ±c, with the sound velocity [START_REF] Whitham | Linear and Nonlinear Waves[END_REF]. From (B.2), real eigenvalues are obtained if and only if c 2 > 0 in [START_REF] Whitham | Linear and Nonlinear Waves[END_REF].

Necessary and sufficient conditions are easily deduced from [START_REF] Whitham | Linear and Nonlinear Waves[END_REF] for the models ( 13)- [START_REF] Kadish | Frequency spectra of nonlinear elastic pulse-mode waves[END_REF] when N = 1: hyperbolicity is satisfied if |ε| < ε c , where Given the nonlinear elastic models ( 13)-( 15), the speed of sound c satisfies:

ε c =                                        d           r + 1 a + 1 1 r -a -1           (model 1), +∞ (model 2), 1 2β if δ = 0, β 3δ         1 + 3δ β 2 -1         otherwise ( model
c 2 =                                                  N ℓ=1 K 1ℓ ρ 1 r -a              r + 1 1 + ε 1ℓ d r - a + 1 1 + ε 1ℓ d a              ( model

Appendix C. Analysis of the relaxation terms

For linear stress-strain relations (18), the relaxation coefficients (33) yield

               ∂∆ ℓ ∂ε (0) = 1 η ℓ σ ′ 2ℓ (0) = K 2ℓ η ℓ , ∂∆ ℓ ∂ε 1ℓ (0) = - 1 η ℓ σ ′ 1ℓ (0) + σ ′ 2ℓ (0) = - 1 η ℓ (K 1ℓ + K 2ℓ ) . (C.1)
The Jacobian matrix of the relaxation function (32) can be obtained

J =                                0 0 • • • 0 0 E 21 η 1 - E 11 + E 21 η 1 0 . . . . . . E 2N η 1 - E 1N + E 2N η N 0
with f ξ = f r if g > g σ , f ξ = f d if g < g σ , f ξ = 0 else. It follows that the eigenvalues are 0,

-K 1ℓ +K 2ℓ η ℓ
, andf ξ .

Appendix D. Semi-analytical solution

The semi-analytical solution of the viscodynamic equations is computed as follows. Fourier transforms in space and time are applied to the system (56). Applying an inverse Fourier transform in space yields v(x, ω) = iω ρ

N ℓ=1 K 1ℓ iω + 1/τ ε ℓ iω + 1/τ σ ℓ Ĝ(ω) 2 π +∞ -∞ 1 k 2 -k 2 0 e -ikx 0 dk, (D.1)
where the hat refers to the Fourier transform, G is the time evolution of the source, the relaxation times τ ε ℓ and τ σ ℓ are defined in [START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF], and k is the wavenumber. The poles ±k 0 satisfy

k 2 0 = ρ ω 2 N ℓ=1 K 1ℓ iω + 1/τ ε ℓ iω + 1/τ σ ℓ (D.2)
with ℑm(k 0 ) < 0. Applying the residue theorem gives the time-domain velocity v(x, t) = ρ

Figure 1 :

 1 Figure 1: parameters of the slow dynamics. (a): Young's modulus E in terms of the concentration of defects g (3), for E + = 14.28 GPa; the vertical dotted line denotes the initial concentration of defects g 0 = 0.3 and the corresponding Young's modulus E 0 = E(g 0 ) = 10 GPa. (b): time evolution of the concentration of defects g given an equilibrium stress σ and two initial values g 0 ; the horizontal dotted line denotes g σ .

Figure 1 -

 1 (b) represents the time evolution of g, given a constant equilibrium concentration g σ = 0.3 denoted by a horizontal dotted line. The restoration and rupturation frequencies are f r = 25 Hz and f d = 250 Hz, respectively. Two initial value of the concentration of defects are considered: g 0 = 0.2 and g 0 = 0.4. In both cases, g tends towards g σ with different rates: destruction is much faster than restoration.

Figure 2 :

 2 Figure 2: equilibrium concentration of defects g σ in terms of the applied stress σ. (a): "exponential model" (7) and "tanh model" (9). (b): "tanh model" (9) with various values of σ. The horizontal dotted line denotes the critical concentration of defects g cr ; the vertical dotted line denotes the central stress σ c .

Figure 3 :

 3 Figure3: Stress-strain relations for the three models[START_REF] Guyer | Nonlinear mesoscopic elasticity: Evidence for a new class of materials[END_REF]. In (a), the dotted lines denote the coordinates of the inflexion point for model 1. The physical parameters are: E = 10 GPa, d = 4.3 10 -4 m, a = 2, r = 4 (models 1 and 2), β = 100, δ = 10 8 (model 3).

Figure 4 :

 4 Figure 4: Rheological model of a generalized Zener material.

5 Table 2 :

 52 [START_REF] Lebedev | A unified model of hysteresis and long-time relaxation in heterogeneous materials[END_REF] 10 -3 2.05 10 -4 4.49 10 -5 7.75 10 -6 τ εℓ (s) 1.53 10 -3 2.49 10 -4 5.50 10 -5 1.06 10 -Relaxation times for a quality factor Q = 20. Optimization with N = 4 relaxation mechanisms on the frequency range [1 kHz, 100 kHz].

Figure 5 :

 5 Figure 5: test 1. Snapshot of the strain after 400 time steps, for two amplitudes of the excitation. The vertical dotted line denotes the location x r of the receiver.

Figure 6 :

 6 Figure 6: test 1. Time history of the strain at the receiver at x r (a), normalized Fourier coefficients (b). The amplitude of the excitation is ε max = 2.0 10 -4 .

with parameters b m = 2 m- 1 , a 1 = 1 , a 2 =

 21112 -21/32, a 3 = 63/768 and a 4 = -1/512. Five receivers are put at abscissae x r = 0.5 + 0.3 ( j -1), with j = 1, • • • 5.

Figure 7 -

 7 Figure 7-(a) shows a seismogram of the velocity recorded at the receivers. Attenuation and dispersion of the waves is clearly observed. Figure 7-(b) compares the numerical solution with the semi-analytical solution after 400 time steps. The computation of the semi-analytical solution is described in Appendix D; it is numerically evaluated with N f = 512 Fourier modes, with a frequency step ∆ f = 200 Hz. Good agreement is observed between numerical and exact values.

Figure 7 :

 7 Figure 7: test 2. Wave propagation in a viscoelastic medium. (a): time evolution of v at a set of receivers; (b): snapshot of v at t = 0.46 ms, and comparison between the numerical and the semi-analytical solution.

  Figure 8 shows the time evolution of the elastic modulus M ≡ E (35); this equality occurs only because a linear stress-strain relation is assumed. The numerical values of M are shown from the beginning of the simulation, whereas the exact values of E (61) are shown from t * . For the sake of clarity, the values are shown only each 5000 time steps. Logically, the elastic modulus decreases as long as the source is switched on (softening), and then increases up to its initial value (recovering).

Figure 8 :

 8 Figure 8: test 3. Time evolution of the elastic modulus M (35) at x r . (a): influence of the central stress σ = 10 8 Pa and 10 7 Pa. (b): influence of the frequency of restoration f r = 2.5 Hz and 100 Hz. The vertical dotted line denotes the time t * where the source is switched off.

Figure 8 -

 8 Figure 8-(a) illustrates the influence of the central stress on the evolution of M: σ = 10 8

Figure 8 -

 8 Figure 8-(b) illustrates the influence of the frequency of restoration on the evolution of M: f r = 2.5 Hz or 100 Hz (the other parameters are those of table1). The lowest value of f r yields

Figure 9 -Figure 9 :

 99 Figure 9-(a) illustrates the influence of viscoelasticity on the stress-strain law. When viscous effects are neglected (Q = +∞, where Q is the quality factor), the behavior induced by the Landau law (15) is observed. Moreover, the scaling (6) induces that the evolution of defects on one cycle is insufficient to provide a measurable hysteretic effect. On the contrary, hysteresis is obtained when viscoelasticity is accounted for (Q = 20). Figure9-(b) mimics the simulation of test 3, where the source a switched-on and off. But contrary to test 3, a nonlinear stress-strain

Figure 10

 10 Figure 10 displays the relative variation of the elastic modulus ∆M = (M -M 0 )/M 0 in terms of the strain, for various amplitudes of the forcing. Three observations can be made.

5 dFigure 10 :

 510 Figure 10: test 4. Relative variations in the elastic modulus M for various amplitudes of forcing ε max , from 10 -6 to 7.5 10 -5 . Top (a-b): without viscoelasticity; bottom (c-d): with viscoelasticity.

Figure A. 11 :

 11 Figure A.11: Properties of the viscoelastic model in the linear regime. (a): reciprocal of the quality factor Q = 20 (A.1). The constant exact value is denoted by a horizontal line; the values obtained after optimization with N = 2 and N = 4 relaxation mechanisms are denoted in blue and red, respectively; the range of optimization [ f min , f max ] is denoted by vertical dotted lines. (b): frequency evolution of the phase velocity; the horizontal dotted lines denote the phase velocity at zero and infinite frequency.

Figure A. 11

 11 Figure A.11 illustrates the properties of the viscoelastic model. Figure A.11-(a) compares the reciprocal of the constant quality factor Q = 20 with the value deduced from (A.1), for N = 2 and N = 4 relaxation mechanisms. Nonlinear optimization is performed from f min = 1 kHz to f max = 100 kHz. Large oscillations are obtained for N = 2; excellent agreement is observed for N = 4. Figure A.11-(b) shows the increase in phase velocity from c 0 = E R /ρ to c ∞ = E/ρ.

  , E, p) = s(ε, E, p), (A.3)which concludes the proof.

  3). (B.5) Model 2 is always hyperbolic. On the contrary, the widely-used Landau model (model 3) is conditionally hyperbolic. When N > 1, the hyperbolicity condition |ε 1ℓ | < ε c is sufficient.

2 +1ℓ d 2 (

 22 ra + a 2 + 6r + 6a + 11 ε

Table 1 :

 1 Physical parameters.

	5.1. Configuration
	ρ (kg/m 3 ) E 0 (GPa)	g 0	f r (Hz) f d (Hz) σ (GPa)	β	δ	Q
	2054	2.21	0.1	25	250	0.1	40 3.5 10 6 20