Ngoc-Son Nguyen
email: ngoc-son.nguyen@inrialpes.fr

Bernard Brogliato
email: bernard.brogliato@inrialpes.fr

N.-S Nguyen

Shock dynamics in granular chains: numerical simulations and comparison with experimental tests

Keywords: Multiple impacts, LZB model, Tapered chain, Stepped chain, Solitary wave interaction

The aim of this paper is to simulate the nonlinear wave propagation in granular chains of beads using a recently introduced multiple impact model and to compare numerical results to experimental ones. Different kinds of granular chains are investigated: monodisperse chains, tapered chains and stepped chains. Particular attention is paid to the dispersion effect, and the wave propagation in tapered chains, the interaction of two solitary waves in monodisperse chains, and the formation of solitary wave trains in stepped chains. We show that the main features of the wave propagation observed experimentally in these granular chains are very well reproduced. This proves that the considered multiple impact model and numerical scheme are able to encapsulate the main physical effects that occur in such multibody systems.

Introduction

The dynamics of chains of beads subjected to collisions has been the object of many studies for a long time. The most well-known example is the so-called Newton's cradle (see e.g. [3, §6.5.6] for references). Chains of beads are of interest for researchers in Solid Mechanics because it is a nice example of a system with multiple impacts (the system is subjected to several simultaneous impacts), and for research-ers in Physics because it is the simplest example of a granular material.

Multiple impacts are known to be a challenging issue in the impact modeling, even for frictionless impacts (like in chains of aligned beads) [START_REF] Glocker | Concepts for modeling impacts without friction[END_REF][START_REF] Payr | An experimental and theoretical study of perfect multiple contact collisions in linear chains of bodies[END_REF]. Despite their apparent simplicity, chains of beads possess a quite complex collision dynamics that involves two main physical effects: dissipation and dispersion of energy. The dissipation is mainly due to the local irreversible deformation at the contact/impact areas and may have many sources: viscosity, plasticity, sound, hysteretic effects, etc. The dispersion rather quantifies the way the initial energy of the chain (before the impact) is spread within the chain after the impact has occurred. The wave effects are responsible for the dispersion, which is due to the discrete nature of multi-body systems. For this reason wave phenomena have been extensively analyzed [START_REF] Nesterenko | Dynamics of Heterogeneous Materials[END_REF][START_REF] Sen | Solitary waves in the granular chain[END_REF]. The dissipation may have, on the other hand, an influence on the dispersion.

Obtaining good models of multiple impacts with reliable numerical methods is an important issue in the field of granular matter in general, and for chains of beads as a particular case. Models based on kinematic restitution coefficients (Newton's like) or on kinetic restitution coefficients (Poisson's like) are known not to be rich enough to correctly represent the dispersion effects. Moreover, they present serious deficiencies like non-uniqueness of the restitution coefficients for a given energetic behavior, or the necessity to estimate the restitution matrix for each chain and each initial data. Models based on Routh's impact dynamics [3, §4.2.13] and an impulse correlation ratio (ICR) have been proposed [START_REF] Ceanga | A new look at an old problem: Newton's cradle[END_REF][START_REF] Gharib | Shock absorption using linear particle chains with multiple impacts[END_REF]. However, it seems that the assumption that the ICR is constant during the impact process and can be estimated from experiments between triplets of beads may fail [43, §6.4]. Let us mention that kinematic laws and binary collisions (the contact gaps are assumed to be all open) are sometimes used [START_REF] Harbola | Pulse propagation in tapered granular chains: an analytic study[END_REF][START_REF] Harbola | Pulse propagation in decorated granular chains: An analytical approach[END_REF][START_REF] Harbola | Pulse propagation in randomly decorated chains[END_REF][START_REF] Pöschel | Extremal collision sequences of particles on a line: optimal transmission of kinetic energy[END_REF][START_REF] Wu | Conservation principles in solitary impulse propagation through granular chains[END_REF]. This is also prone to some fundamental issues like the fact that it is not guaranteed, in general, that letting the gaps tend to zero (hence recovering the real system where the beads touch each other) yields a unique limit, because the trajectories may be discontinuous with respect to the initial data [START_REF] Paoli | Continuous dependence on data for vibro-impact problems[END_REF].

Models based on the second order equation of motion with a linear or nonlinear spring modeling the interaction at each contact point have been extensively used [START_REF] Doney | Impulse absorption by tapered horizontal alignments of elastic spheres[END_REF][START_REF] Falcon | Collision of a 1-D column of beads with a wall[END_REF][START_REF] Manciu | Impulse propagation in dissipative and disordered chains with power-law repulsive potentials[END_REF][START_REF] Mishra | On the determination of contact parameters for realistic DEM simulations of ball mills[END_REF][START_REF] Sen | Thermalizing an impulse[END_REF][START_REF] Zhang | The calculation of contact forces between particles using spring and damping models[END_REF]. It is well known that this kind of models gives satisfactory results for elastic systems, in particular for chains of elastic balls. Therefore, most of studies on the wave propagation in granular chains carried out by physicists are based on this kind of models. However, the main difficulty of these models lies in the modeling of the dissipative behavior of mechanical systems. In fact, the loss of energy during the collision process results from many phenomena which can be ranged into two distinct categories: local and global phenomena. The local phenomena concern inelastic behaviors at the contact points between colliding bodies such as plasticity, viscosity, friction etc., while the global phenomena concern the conversion of energy into noise, light, vibrational waves inside the bodies, etc. The vibrational waves may be damped by the viscosity, plasticity inside the bodies afterward. In the case of massive bodies like balls, the global phenomena are very small compared to the local ones [START_REF] Johnson | Contact Mechanics[END_REF], thus one needs to focus on the modeling of the dissipation of energy caused by the local phenomena. Let us assume that the local dissipation is mainly due to the viscosity and plasticity at the contact points, i.e. the friction is neglected, and discuss how the local dissipation is modeled.

Many works introduced the viscous dissipation with a linear or nonlinear dashpot coupled to a linear or Hertz's spring [START_REF] Carretero-González | Dissipative solitary waves in granular crystals[END_REF][START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Job | How hertzian solitary waves interact with boundaries in a 1D granular medium[END_REF][START_REF] Mishra | On the determination of contact parameters for realistic DEM simulations of ball mills[END_REF][START_REF] Rosas | Observation of two-wave structures in strongly nonlinear dissipative granular chains[END_REF][START_REF] Zhang | The calculation of contact forces between particles using spring and damping models[END_REF]. This kind of models is usually called "spring-dashpot models". In this respect, there exist a variety of spring-dashpot models: different assemblages of linear springs and linear dashpots (Kelvin-Voigt model [START_REF] Yigit | Spring-dashpot models for the dynamics of a radially rotating beam with impact[END_REF], Maxwell model [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF], Zener model [START_REF] Zener | The intrinsic inelasticity of large plates[END_REF]), different assemblages of Hertz springs and nonlinear dashpots [START_REF] Brilliantov | Model for collisions in granular gases[END_REF][START_REF] Carretero-González | Dissipative solitary waves in granular crystals[END_REF][START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF][START_REF] Kuwabara | Restitution coefficient in a collision between two spheres[END_REF]. The assemblages of linear springs and linear dashpots may produce some physical inconsistencies such as the jump of the contact force at the beginning of collision or the traction force at the end of collision. Except the model proposed by Kuwabara and Kono [START_REF] Kuwabara | Restitution coefficient in a collision between two spheres[END_REF] that is based on the theory of viscous-elastic continuum materials, the choice of the assemblages of nonlinear springs and nonlinear dashpots is, in general, arbitrary. Moreover, in some models the parameters do not possess clear physical meanings, and in this case a calibration procedure is needed to identify the parameters, i.e. the values of the parameters are chosen so that the numerical results fit the best to the experimental ones, for example, the models described in [START_REF] Carretero-González | Dissipative solitary waves in granular crystals[END_REF][START_REF] Vergara | Model for dissipative highly nonlinear waves in dry granular systems[END_REF]. It is worth noting that viscous dissipation models may be valid for viscous-elastic systems such as chains of polytetraflouroethylene (PTFE) balls, how-ever they may produce some physical inconsistencies for the systems exhibiting mainly the plastic-elastic behavior such as chains of metal balls.

To our best knowledge, there are few models treating the plastic dissipation during the impact process. In [START_REF] Doney | Impulse absorption by tapered horizontal alignments of elastic spheres[END_REF][START_REF] Manciu | Impulse propagation in dissipative and disordered chains with power-law repulsive potentials[END_REF], the local dissipation at contact points is modeled by a restitution coefficient defined as the ratio of the unloading force (during the expansion phase) to the loading force (during the compression phase). This model does not represent even qualitatively the plastic-elastic behavior observed experimentally: it includes a jump of the contact force at the end of the compression phase and the deformation is entirely recovered at the end of the expansion phase despite the fact that energy is dissipated during the collision process. Walton and Braun [START_REF] Walton | Viscosity, granular-temperature, and stress calculations for shearing assemblies of inelastic, frictional disks[END_REF] introduced a bistiffness contact model to take the plastic dissipation caused by the normal interaction between two disks into account. Afterward, this bistiffness model has been used in many works, for example in the discrete element method for granular materials [START_REF] Luding | Introduction to discrete element methods: basis of contact force models and how to perform the micro-macro transition to continuum theory[END_REF][START_REF] Sadd | Contact law effects on wave propagation in particulate materials using distinct element modeling[END_REF]. This kind of normal contact model seems to be appropriate to model the plastic-elastic behavior. However, the parameters in these works are generally identified by a calibration procedure as in the case of spring-dashpot models.

Another multiple impact model, called "LZB multiple impact model", has been introduced recently in [START_REF] Liu | Energy dissipation and dispersion effects in a granular media[END_REF][START_REF] Liu | Frictionless multiple impacts in multibody systems: Part I. Theoretical framework[END_REF][START_REF] Liu | Frictionless multiple impacts in multibody systems: Part II. Numerical algorithm and simulation results[END_REF]. Similarly to the model introduced by Walton and Braun [START_REF] Walton | Viscosity, granular-temperature, and stress calculations for shearing assemblies of inelastic, frictional disks[END_REF], this model takes the compliance effect and the dissipation effect of elasto-plastic colliding bodies into account by a bistiffness compliant contact model. However, the bistiffness contact model is modified in such a way that it represents the nonlinear nature of the interaction between colliding bodies, for example Hertz's interaction. LZB model differs from other compliant models in the following points. Firstly, it is based on Darboux-Keller's approach [START_REF] Darboux | Etude géométrique sur les percussions et les chocs des corps[END_REF][START_REF] Keller | Impact with friction[END_REF] to replace the time scale by the normal impulse scale. In this approach, the evolution of a system during impact is described at the impulsevelocity level rather than at the force-acceleration level as in other compliant contact models. As a consequence, the unimportant information for an impulsive process such as time, displacements and forces is completely filtered out from the model. It is worth noting that the latter information is needed in other compliant contact models (for example, spring-dashpot models), thus some numerical difficulties may be raised due to the stiff force-displacement relation at each contact point. From the point of view of impact modeling, Darboux-Keller's approach is consistent with the fact that, for a rigid body system, the impact duration can be considered as infinitesimal, during which the configuration of the system can be considered as constant. As a consequence, the change of velocities and impulses are important information to be determined. Secondly, the interaction at each contact point is described in terms of the evolution of energy rather than the force-displacement relation: during the compression phase, the kinetic energy of the colliding bodies is transferred to the potential energy, whose part is released during the expansion phase, respecting Stronge's energetic constraint [START_REF] Stronge | Impact Mechanics[END_REF]. Consequently, the dissipation of energy during the collision at each contact point is taken into account by Stronge's restitution coefficient, also called "energetic restitution coefficient". This restitution coefficient possesses a clear physical meaning, and can be properly identified from independent tests, for example from a central collision between two balls. Thirdly, the dispersion of energy throughout the whole system during the impact is described by a rule distributing normal impulse increments among the contact network, depending on the stiffness ratio and the potential energy ratio between various contacts. Consequently, the absolute stiffness at each contact point is not a factor affecting the integration process for which the choice of the step size is independent of the absolute stiffness. This is not the case for other compliant models for which a high stiffness at the contact points requires an extremely small step size. These are the main reasons why we use LZB model to simulate impact problems within chain of balls in this paper. It is worth noting that this impact model is suitable for elasto-plastic multi-body systems. Many comparisons between numerical results obtained with LZB multiple impact model and experimental data have been presented in [START_REF] Liu | Energy dissipation and dispersion effects in a granular media[END_REF] for column of beads (experiments of [START_REF] Falcon | Collision of a 1-D column of beads with a wall[END_REF]), in [START_REF] Zhao | Planar dynamics of a rigid body system with frictional impacts. II. Qualitative analysis and numerical simulations[END_REF] for the bouncing dimer (experiments of [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF]), Newton's cradle in [START_REF] Liu | Frictionless multiple impacts in multibody systems: Part II. Numerical algorithm and simulation results[END_REF] (experiments of [START_REF] Ceanga | A new look at an old problem: Newton's cradle[END_REF]), and more recently some results on the rocking block system used in the Earthquake Engineering literature [START_REF] Zhang | The planar rocking block: analysis of kinematic restitution laws, and a new rigid-body impact model with friction[END_REF].

Scientists have paid much attention in investigating the energy propagation in granular chains in order to find out shock protection devices that are able to attenuate efficiently the energy induced in a shock process. Various kinds of granular chains have been investigated: (1) monodisperse chains, i.e. chains of identical beads [START_REF] Falcon | Collision of a 1-D column of beads with a wall[END_REF][START_REF] Job | How hertzian solitary waves interact with boundaries in a 1D granular medium[END_REF][START_REF] Nesterenko | Dynamics of Heterogeneous Materials[END_REF][START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF][START_REF] Sen | Solitary waves in the granular chain[END_REF];

(2) monodisperse chains with defects [START_REF] Hascoë | Shock in non-loaded bead chains with impurities[END_REF][START_REF] Job | Wave localization in strongly nonlinear hertzian chains with mass defect[END_REF]; (3) tapered chains, i.e. chains with decreasing size of the beads [START_REF] Doney | Impulse absorption by tapered horizontal alignments of elastic spheres[END_REF][START_REF] Harbola | Pulse propagation in tapered granular chains: an analytic study[END_REF][START_REF] Pfannnes | Energy absorption and recovery in tapered granular chains: small chains and low tapering[END_REF][START_REF] Sen | Thermalizing an impulse[END_REF][START_REF] Sokolow | Absorption of short duration pulses by small, scalable, tapered granular chains[END_REF][START_REF] Wu | Conservation principles in solitary impulse propagation through granular chains[END_REF]; (4) stepped chains, i.e. chains composed of a large monodisperse section followed by a small monodisperse section [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF][START_REF] Nesterenko | Solitary waves in discrete media with anomalous compressibility and similar to "sonic vacuum[END_REF]; (5) decorated chains with small masses placed regularly or randomly among larger masses [START_REF] Gharib | Shock absorption using linear particle chains with multiple impacts[END_REF][START_REF] Harbola | Pulse propagation in decorated granular chains: An analytical approach[END_REF][START_REF] Harbola | Pulse propagation in randomly decorated chains[END_REF]; (6) composite chains, i.e. chains composed of beads made of different materials periodically or randomly distributed [START_REF] Daraio | Energy trapping and shock disintegration in a composite granular medium[END_REF][START_REF] Nesterenko | Anomalous wave reflection at the interface of two strongly nonlinear granular media[END_REF][START_REF] Ponson | Nonlinear waves in disordered diatomic granular chains[END_REF][START_REF] Porter | Highly nonlinear solitary waves in heterogeneous periodic granular media[END_REF][START_REF] Rosas | Observation of two-wave structures in strongly nonlinear dissipative granular chains[END_REF]; [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] disordered chains with beads of masses randomly distributed [START_REF] Manciu | Impulse propagation in dissipative and disordered chains with power-law repulsive potentials[END_REF][START_REF] Pöschel | Extremal collision sequences of particles on a line: optimal transmission of kinetic energy[END_REF]. Monodisperse chains, tapered chains, stepped chains are the object of this paper. More specifically, LZB multiple impact model is used to simulate the dispersion effect and the wave propagation in tapered chains, the formation of solitary wave trains in stepped chains and the interaction of two solitary waves in monodisperse chains. The numerical results obtained will be compared to the experimental data available in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF][START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF][START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF][START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF]. This paper is organized as follows: in Sect. 2, LZB multiple impact model will be briefly described. Section 3 will be dedicated to investigate the dispersion effect and the wave propagation in tapered chains and compare the numerical results to the experimental results presented in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF][START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF]. The formation of solitary wave trains in stepped chains will be the subject of Sect. 4, and the numerical results will be compared to the experimental ones shown in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF]. In Sect. 5, the numerical results concerning the interaction of two solitary waves in monodisperse chains will be presented and compared to the experimental data available in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF]. Finally, some conclusions about the numerical simulations with LZB multiple impact model will be drawn.

Brief presentation of LZB multiple impact model

LZB multiple impact model has been introduced in [START_REF] Liu | Energy dissipation and dispersion effects in a granular media[END_REF][START_REF] Liu | Frictionless multiple impacts in multibody systems: Part I. Theoretical framework[END_REF][START_REF] Liu | Frictionless multiple impacts in multibody systems: Part II. Numerical algorithm and simulation results[END_REF] to solve impact dynamics in a multi-rigid body system where various contacts are established simultaneously. Granular chains are typical examples of such mechanical systems. In this section, we present briefly how to carry out numerical simulations of impact dynamics in a granular chain with LZB multiple impact model. For other systems, for example a rocking block or a bouncing dimer, see [START_REF] Zhang | The planar rocking block: analysis of kinematic restitution laws, and a new rigid-body impact model with friction[END_REF][START_REF] Zhao | Planar dynamics of a rigid body system with frictional impacts. II. Qualitative analysis and numerical simulations[END_REF]. As illustrated in Fig. 1, a granular chain is composed of beads whose center is aligned on a horizontal axis. This alignment aims at ensuring colinear collisions between beads. The last bead in the chain may move freely or may be constrained by a rigid wall placed at the end of the chain. This rigid wall may represent the wall sensor used in experimental setups in order to measure the force felt at the end of granular chains. The presence of a rigid wall at the end of the chain is considered here. The beads in the chain are numbered as increasing integers from 1 for the first bead to N for the last bead (N is the number of beads in the chain). For a bead i, the position of its center is defined by x i , and its radius is R i . The beads in the chain are initially stationary and barely touch each other. The chain is then impacted by a bead, called the striker, numbered as 0 with an impact velocity V s . The considered system (striker + chain) has N +1 degrees of freedom (when none of contacts is closed) and N + 1 unilateral contacts. The contact between beads i and i + 1 is numbered as i and the contact between the last bead N and the wall is numbered as N . Let us assume that there is no friction at the contacts between beads and rails that are used in some experimental setups to align the beads in the chain. Because of colinear collisions between beads, there is no bead rotation during the impact process. Therefore, the state of this system can be described by the following generalized coordinate q(t) = [x 0 (t), x 1 (t), . . . , x N (t)] T . Dynamics of this system is described by the equation of motion and a set of complementarity conditions as follows:

M(q) q = F ext (t) + W (q)λ(t) 0 ≤ λ(t) ⊥ δ(q) ≥ 0, (1)
where:

-M(q) ∈ IR (N +1)×(N +1) is the inertia matrix:

M = ⎡ ⎢ ⎢ ⎢ ⎣ m 0 0 • • • 0 0 m 1 • • • 0 . . . . . . . . . . . . 0 0 • • • m N ⎤ ⎥ ⎥ ⎥ ⎦ ; (2)
δ(q) = [δ 0 (q) δ 1 (q)...δ N (q)] T collects the gaps at the contact points in the system defined as:

δ i = x i+1 -x i -(R i + R i+1 ), ∀i = N δ N = x wall -x N -R N , (3)
and λ(t) = [λ 0 (t) λ 1 (t)...λ N (t)] T collects the normal forces at the contact points in the system; -F ext (t) is the external loading, for instance the force of gravity; -W (q) ∈ IR (N +1)×(N +1) is the Jacobian matrix defined as:

W = ∂δ(q) ∂q T = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -1 0 0 • • • 0 0 1 -1 0 • • • 0 0 0 1 -1 . . . . . . 0 . . . . . . . . . . . . 0 . . . 0 0 0 1 -1 0 0 0 0 0 1 -1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ; (4)
-the symbol ⊥ means that λ i = 0 if δ i > 0, and

λ i > 0 if δ i = 0, for all 0 ≤ i ≤ N .
The link between the relative velocity at the contact points δ and the generalized velocity q is given by:

δ = W T q. ( 5
)
When the striker touches the chain, the impact process occurs. To solve the impact process, the Darboux-Keller's dynamical equation [START_REF] Darboux | Etude géométrique sur les percussions et les chocs des corps[END_REF][START_REF] Keller | Impact with friction[END_REF] is considered to change the time scale to the impulse scale with a basic assumption that the configuration q of the system is constant during the impact process, i.e. M(q) and W (q) are constant. It is worth noting that M(q) and W (q) for a granular chain shown in (2) and (4) are constant so the Darboux-Keller's basic assumption mentioned above is automatically satisfied. The Darboux-Keller's dynamical equation is as follows: Md q = W d P, [START_REF] Coste | On the validity of Hertz contact law for granular material acoustics[END_REF] where P = [P 0 P 1 . . . P N ] T collects the normal impulses at the contact points in the system. The interaction at each contact point is assumed to respect Stronge's energetic constraint [START_REF] Stronge | Impact Mechanics[END_REF]:

e 2 s, j = - W e, j W c, j = - P f j P c j δ j d P j P c j 0 δ j d P j ( 7
)
where W c, j and W e, j are the works done by the force at contact j during the compression and expansion phases, respectively, and P c j and P f j are the impulses at the end of the compression phase ( δ j (P c j ) = 0) and of the expansion phase, respectively. The constant e s is called the "energetic restitution coefficient" which is restricted between 0 and 1: e s = 1 (resp. e s = 0) corresponds to a purely elastic (resp. plastic) contact. It is worth noting that for a single impact between two bodies, Eqs. ( 6) and ( 7) are enough to solve the impact problem with the energetic restitution coefficient e s . This means that one does not need to model what exactly takes place during the impact process by using a contact model. However, this is not the case for the multiple impact problem for which the coupling between various contacts requires an additional information about the distribution of the impulse change d P i among the contacts in the system. A distributing rule is aimed at giving the information required, and it can be formulated by using a compliant model at each contact point in the system. In the following, we will present how to formulate the distributing rule from the chosen compliant contact model.

The compliance effect is added into LZB model by using a compliance contact model such as the mono-stiffness or bi-stiffness compliance models. We limit ourselves to the bi-stiffness compliance model. The reader can refer to [START_REF] Liu | Energy dissipation and dispersion effects in a granular media[END_REF][START_REF] Liu | Frictionless multiple impacts in multibody systems: Part I. Theoretical framework[END_REF] for the details about the mono-stiffness model. As illustrated in Fig. 2, for the bi-stiffness compliance model the force-indentation relationship at a given contact point j for the compression phase expressed as:

λ c, j = K j (δ j ) η , ( 8
)
is different from that for the expansion phase expressed as:

λ e, j = λ M, j δ j -δ r, j δ M, j -δ r, j η , ( 9
)
Fig. 2 Bi-stiffness contact compliant model where δ r, j is the plastic deformation, and λ M, j and δ M, j are respectively the maximum values of the normal contact force and of the normal deformation at the contact point j at the end of the compression phase (when δ j = 0). The elasticity coefficient η is equal to 1 for the linear elasticity and to 3/2 for the Hertz's contact. In this paper, the Hertz's elasticity coefficient is considered. It is worth noting that the dissipated energy during the inelastic collision between two particles is taken into account by the area enclosed by the compression and expansion curves as seen in Fig. 2. Applying the energetic constraint [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] to the bi-stiffness compliance model described in ( 8) and ( 9), the relation of the plastic parameter δ r, j to the energetic restitution coefficient can be found as δ r, j = δ M, j (1e 2 s, j ). The stiffness K j is computed according to the Hertz's theory as follows:

K j = 4 3 R j E j , (10)
where R j and E j are defined in [START_REF] Darboux | Etude géométrique sur les percussions et les chocs des corps[END_REF] for a contact point between two beads and in [START_REF] Doney | Impulse absorption by tapered horizontal alignments of elastic spheres[END_REF] for the contact between the last bead and the wall.

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1 R j = 1 R j + 1 R j+1 1 E j = 1 -ν 2 j E j + 1 -ν 2 j+1 E j+1 , ( 11
) ⎧ ⎨ ⎩ R N = R N 1 E N = 1 -ν 2 N E N + 1 -ν 2 w E w . ( 12
)
Here E j (E w ) and ν j (ν w ) are the Young modulus and the Poisson coefficient of bead j (of the wall, respectively).

According to this compliance contact model, the potential energy at a contact point is computed in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF] for the com-pressions phase and in [START_REF] Falcon | Collision of a 1-D column of beads with a wall[END_REF] for the expansion phase:

E j (P j ) = P j 0 δ j (P j )d P j , ( 13
)
E j (P j ) = E M, j + 1 e 2 s, j P j P c j δ j (P j )d P j , ( 14
)
where E M, j is the potential energy accumulated at the end of the compression phase (at the point M). The force at a contact point is computed from the potential energy as follows:

λ j = (1 + η) η η+1 K 1 1+η j (E j (P j (t))) η η+1 . ( 15
)
From ( 15) and with the relation d P j = λ j dt, we can obtain the distributing law:

d P j d P * = Γ j, * = K j K * 1 1+η E j (P j ) E * (P * ) η η+1 . ( 16
)
The symbol (*) in ( 16) depicts the primary contact point that is defined as the contact point at which the potential energy is maximal among all the contacts in the system. As shown in ( 16), the distributing law expresses that the change of the impulse at a given contact d P j during the multiple impact process is related to the change of the principal impulse d P * (defined at the primary contact), to the relative potential energy E j /E * and to the relative stiffness K j /K * between the considered contact and the primary contact. Combining ( 16) with ( 6), the Darboux-Keller's dynamical equation can be rewritten as:

M d q d P * = W Γ , ( 17
)
where Γ = [Γ 0, * Γ 1, * ...Γ N , * ] T is computed with [START_REF] Glocker | Concepts for modeling impacts without friction[END_REF]. It can be noted that ( 17) is an Ordinary Differential Equation (ODE) whose the dumb integration variable is the principal impulse P * defined at the primary contact (*). It is worth noting that the distributing law ( 16) results naturally from the compliant contact model chosen, however it is needed to solve the multiple impact process with Darboux-Keller dynamics ( 6) at the velocity-impulse level instead of the acceleration-force level considered in the second order dynamics. On another point, the time scale t is related to the principal impulse scale P * by:

dt = d P * F * , ( 18
)
where F * is the force computed at the primary contact point using [START_REF] Gharib | Shock absorption using linear particle chains with multiple impacts[END_REF]. In some situations, a contact point in the system may take some initial potential energy due to the precompression or may experience multiple compression and expansion phases. Therefore, one needs to extend the above formulations to Fig. 3 Illustration of a contact experiencing multiple compression and expansion phases these complex phenomena. A full description of this extension can be found in [START_REF] Liu | Frictionless multiple impacts in multibody systems: Part I. Theoretical framework[END_REF]. We present briefly here how to treat these phenomena. For such an extension, the relationship between the contact force and the indentation expressed in [START_REF] Daraio | Strongly nonlinear wave dynamics in a chain of polymer coated beads[END_REF] and in ( 9) is assumed to be unchanged for interactions with an initial potential energy or with multiple compression and expansion phases. Consequently, the distributing law ( 16) is still valid for these cases. In the case of the precompression, the potential energy during the compression phase is computed as:

E j (P j ) = E o, j + P j 0 δ j (P j )d P j , ( 19
)
where E o, j is the initial potential energy at the contact under consideration, and the release of the potential energy during the expansion phase must respect the energetic constraint extended as:

e 2 s, j = - W e, j W c, j + E o, j , (20)
where W c, j is the work done by the contact force during the compression phase, excluding the precompression phase.

With the extended energetic constraint [START_REF] Hascoë | Shock in non-loaded bead chains with impurities[END_REF], it is easy to prove that the potential energy during the expansion phase still follows relation [START_REF] Falcon | Collision of a 1-D column of beads with a wall[END_REF]. When a contact experiences multiple cycles of compression and expansion phases, secondary cycles can be treated as the case of single cycles with an initial potential energy. The initial potential energy of each secondary cycle is the one computed for the preceding cycle at the point where the secondary cycle under consideration starts. To illustrate this, we consider a contact experiencing two compression and expansion phases shown in Fig. 3 as an example. The secondary cycle starts at the point R, follows the compression curve R M 2 and then the expansion curve M 2 B, and finally stops at the point B where the potential energy is entirely released or dissipated. This secondary cycle is treated as a single cycle whose initial potential energy E o, j is the potential energy at the point R computed for the cycle O M 1 R.

System of Eqs. ( 5), ( 13), ( 14), ( 15), ( 16), ( 17) and ( 18) is integrated with respect to the primary impulse P * . The primary impulse P * is first discretized into small steps of size ΔP and then the numerical algorithm described in [START_REF] Liu | Frictionless multiple impacts in multibody systems: Part II. Numerical algorithm and simulation results[END_REF][START_REF] Nguyen | Shock dynamics in granular chains: numerical simulations and comparison with experimental tests[END_REF] is applied for the numerical integration. LZB multiple impact model is coupled to an event-driven scheme implemented in the open source Siconos available at http://siconos.gforge. inria.fr. All the numerical simulations shown in this paper are performed with a step size ΔP = 10 -6 N s which gives satisfactory numerical solutions.

Dry friction can be incorporated in LZB impact model. This has been done in [START_REF] Zhang | The planar rocking block: analysis of kinematic restitution laws, and a new rigid-body impact model with friction[END_REF][START_REF] Zhao | Planar dynamics of a rigid body system with frictional impacts. II. Qualitative analysis and numerical simulations[END_REF]. For chains of balls, the friction is present at contacts between ball and ball and between ball and rail used to align balls. These two frictions can be either static or dynamic and they are a priori different. Consequently, taking into account the friction requires some more parameters. However, these friction parameters are not available in the experimental data that we extract from some papers. This is the reason why we neglect the friction in the numerical simulations presented in this paper, which may lead to some discrepancies between numerical and experimental results as we will see later.

Numerical tests on tapered chains

A tapered chain is composed of beads whose diameter is progressively decreased such that Φ i+1 = (1q)Φ i (Φ i is the diameter of bead i and q is the tapering factor). This kind of granular chains has been extensively investigated analytically, numerically and experimentally [START_REF] Doney | Impulse absorption by tapered horizontal alignments of elastic spheres[END_REF][START_REF] Harbola | Pulse propagation in tapered granular chains: an analytic study[END_REF][START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF][START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF][START_REF] Pfannnes | Energy absorption and recovery in tapered granular chains: small chains and low tapering[END_REF][START_REF] Sen | Thermalizing an impulse[END_REF][START_REF] Sokolow | Absorption of short duration pulses by small, scalable, tapered granular chains[END_REF][START_REF] Wu | Conservation principles in solitary impulse propagation through granular chains[END_REF]. Tapered chains show interesting mechanical features such as the dispersion effect and the capability of shock absorption. Nakagawa et al. [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] and Melo et al. [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] have conducted experimental tests on some tapered chains to put in evidence the above features. The authors in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] investigate the outcome of the impact (post-impact velocity, kinetic energy of the beads, etc.), whereas the authors in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] investigate the force pulse propagation in such media. In this section we show the numerical results obtained from the numerical simulations carried out with LZB multiple impact model on the tapered chains considered in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF][START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] and compare with the corresponding experimental results.

In order to obtain the experimental data shown in a figure of a paper, we use the following data extracting technique. First the considered figure is extracted from an electronic version of the paper. Then we set the resolution of the

Description of the experimental test

The configuration of the experimental tests performed by Nakagawa et al. [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] is illustrated in Fig. 4. The considered tapered chain is composed of 19 beads (numbered from 1 for the largest bead to 19 for the smallest bead) with the diameter of the first bead Φ 1 = 9.5 mm and the tapering factor q = 5 %. Initially, the beads in the tapered chain barely touch one another. Then a bead (numbered 0) of diameter Φ 0 = 10 mm strikes the chain at the left end. The beads are made of chrome steel with the following properties: mass density ρ = 7833 kg/m 3 , Young modulus E = 203 GPa, Poisson ratio ν = 0.3. In order to estimate the value of the restitution coefficient for the beads, the authors have performed binary collision experiments with various pairs of identical beads. Three values are found for the restitution coefficient: 0.947, 0.965 and 0.955 corresponding to three values of the bead velocity averaged over 5, 15 and 20 ms, respectively. The velocity of the beads at the end of the impact (defined at the moment when the smallest bead leaves the chain) is measured by a high-speed digital image analysis.

Numerical simulations

The parameters used in the numerical simulations presented here are the same as those presented above. In fact, the restitution coefficient estimated in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] is the kinematic restitution coefficient (Newton's restitution coefficient). In the case of colinear collisions the kinematic restitution coefficient is Fig. 5 Velocity of different beads, normalized by the impact velocity, versus time during the impact process resulting from the numerical simulation of the tapered chain considered in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] with e s = 0.965. The number above each curve corresponds to the bead number equal to the energetic restitution coefficient. Therefore, we can use the restitution coefficient estimated in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] as the energetic restitution coefficient e s . Since there is not a unique value presented in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF], we choose to use all the three values presented above: e s = 0.947, 0.955, 0.965.

The evolution in time of the velocity of various beads in the considered tapered chain, which is obtained from the numerical simulation with e s = 0.965, is plotted in Fig. 5. We can see that the effect of tapering the bead diameter of the chain leads to an increase in the peak velocity of the beads as the wave propagates. Moreover, a tail is left after the wave has gone through each bead. This means that the velocity of the beads in the tapered chain does not go back to zero after the wave has left them. This demonstrates the dispersion effect of the tapered chain. These results are similar to those shown in [START_REF] Sen | Thermalizing an impulse[END_REF]. It can be noted that the multiple impact process in the tapered chain is over at instant of about 0.18 ms. According to LZB model, the multiple impact process is considered to be over when the relative velocity is negative at all the contacts, i.e. the contacting bodies move away from each other, and the potential energy at all the contacts is entirely released or dissipated. The velocity of the beads at the end of the multiple impact process (at instant of about 0.18 ms) is called the "post-impact velocity of the beads".

Figure 6 shows the post-impact velocity of the beads in the considered tapered chain obtained from the numerical simulations with the three values of the restitution coefficient (0.947, 0.965, 0.955), compared to the experimental results presented in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF]. In general, the numerical results are close to the experimental ones. One can note that the numerical simulations slightly underestimate and overestimate the postimpact velocity of first beads (beads 1-5) and of last beads (beads [START_REF] Gharib | Shock absorption using linear particle chains with multiple impacts[END_REF][START_REF] Glocker | Concepts for modeling impacts without friction[END_REF][START_REF] Harbola | Pulse propagation in tapered granular chains: an analytic study[END_REF][START_REF] Harbola | Pulse propagation in decorated granular chains: An analytical approach[END_REF][START_REF] Harbola | Pulse propagation in randomly decorated chains[END_REF][START_REF] Hascoë | Shock in non-loaded bead chains with impurities[END_REF], respectively (we will see these discrepancies more clearly in Fig. 7). It can be observed that the numerical results obtained with e s = 0.965 match the best the experimental ones. In [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF], the authors have compared the post-Fig. 6 Post-impact velocity of the beads, normalized by the impact velocity, obtained in the experimental tests and in the numerical simulations with different values for the restitution coefficient e s . The experimental data are extracted from figure 6.a in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] Fig. 7 Post-impact velocity (V ) and kinetic energy (KE) of the beads (normalized by the respective initial values of the impactor) plotted in the semilog scale, obtained in the experimental tests and in the numerical simulations with the restitution coefficient e s = 0.965. The experimental data are extracted from figure 6.b in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] impact velocity of the beads in the chain obtained from the experimental tests and from an independent-collision model (sequence of binary collisions introduced in [START_REF] Wu | Conservation principles in solitary impulse propagation through granular chains[END_REF]). The comparison has shown that the independent-collision model significantly underestimates the dispersion effect in the chain (see figure 6.a in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF]). On the other hand, the numerical results obtained with LZB multiple impact model can reproduce much better this dispersion effect.

The numerical and experimental results can be more closely compared when plotting the post-impact velocity and kinetic energy of the beads in the semilog scale as shown in Fig. 7. In this figure, only the numerical results obtained with e s = 0.965, which yields the best results in Fig. 6, are shown. With this plot, we can see more clearly that the numerical simulation overestimates the post-impact velocity of first beads, and thus their post-impact kinetic energy. Nevertheless, the numerical simulation gives, on the whole, satisfactory results.

In conclusion, the numerical simulations with LZB multiple impact model reproduce well the dispersion effect Fig. 8 Configuration of the granular chains considered by Melo et al. [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] observed in the experimental tests in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF]. It can be noted that the parameters used for these numerical simulations are all given in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF]. We do not proceed any parameter fitting. The numerical results in Figs. 6 and7 clearly show that LZB multiple impact model gives results much better than the binary collision model employed in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF][START_REF] Wu | Conservation principles in solitary impulse propagation through granular chains[END_REF] (see figure 6 in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF]). The post-impact velocities are not well predicted for the first 5 beads. However for all the other beads the prediction is very good (for the last two beads the experimental and numerical results are too close to be distinguishable).

Tapered chains considered in [34]

Description of the experimental tests

The granular chains considered in the experimental tests carried out by Melo et al. [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] are illustrated in Fig. 8. A monodisperse chain composed of 16 beads of diameter Φ o = 26 mm is placed in front of a tapered chain. Two tapered chains are considered: the first one is composed of 14 beads with the tapering factor q 1 = 5.6 % and the second one is composed of 12 beads with q 2 = 8.27 %. Initially, the beads in the chains barely touch one another. Then a smaller bead of diameter Φ s = 8 mm strikes at the left extremity of the monodisperse chain in order to generalize a solitary wave propagating in the chain. The beads are made of high carbon hardened steel whose properties are as follows: mass density ρ = 7780 kg/m 3 , Young modulus E = 203 GPa, Poisson ratio ν = 0.3.

During the experimental tests, the force at a given contact in the monodisperse chain is directly measured by a sensor that is inserted inside one of two adjacent beads at the considered contact. This measurement technique was developed by Job et al. [START_REF] Job | How hertzian solitary waves interact with boundaries in a 1D granular medium[END_REF][START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF][START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] in order to achieve the actual force felt exactly at the interface between two beads. The force at a given contact in a tapered chain, for example the right contact of kth bead in the tapered chain, is indirectly measured as follows. First, the same test as the previous ones is performed on the monodisperse chain followed by a tapered sub-chain composed of 1st bead up to kth bead in the tapered chain. The force at the end of the tapered sub-chain is recorded by a wall sensor. This wall sensor is made of the same material as the beads. Then the force pulse at the right contact of the kth bead in the tapered chain (denoted as

F wall m,k = 1.94F chain m,k , τ wall k = 1.09τ chain k . ( 21
)
Scaling relation ( 21) is presented in [START_REF] Job | How hertzian solitary waves interact with boundaries in a 1D granular medium[END_REF], that is obtained from numerical simulations of the interaction between monodisperse chains of perfectly elastic beads and an elastic wall made of the same material as the beads. We have also performed the same numerical simulations as those in [START_REF] Job | How hertzian solitary waves interact with boundaries in a 1D granular medium[END_REF] with LZB multiple impact model for two cases without and with local dissipation. Scaling relation ( 21) is verified in our simulations. However, such a scaling relation might be not suitable for tapered chains for which impact dynamics is much more complex than for monodisperse chains. This point will be analyzed later.

Numerical simulations

The granular chains considered in the numerical simulations presented in this section are the same as those in the experimental tests in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]. In order to analyze the wave propagation in tapered chains and compare with the experimental results shown in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF], we perform two types of numerical tests:

-Type 1: granular chains for this type of tests are composed of the monodisperse chain followed by a tapered sub-chain whose number of beads varies from n 1 = 0 to n 1 = 14 for q 1 = 5.6 %, and from n 2 = 0 to n 2 = 12 for q 2 = 8.27 %. A rigid wall made of the same material as the beads is placed at the end of the tapered sub-chain. This type of tests is illustrated in Fig. 9. In the following, we refer to the numerical tests of this type as the numerical tests T 1 ; -Type 2: granular chains for this type of tests are composed of the monodisperse chain followed by a whole tapered chain (chain of 14 beads with q 1 = 5.6 % or chain of 12 beads with q 2 = 8.27 %). No rigid wall is placed at the end of the tapered chain. Figure 8 illustrates this type of tests. In the following, we refer to the numerical tests of this type as the numerical tests T 2 .

The force pulse at a contact in the monodisperse chain is directly computed for both types of numerical tests, whereas the force pulse at a contact in the tapered chains is computed by the two following methods: (i) the force pulse at the considered contact is extrapolated from the force pulse computed at the end of the corresponding tapered sub-chain in the tests T 1 with the method described in Sect. 3.2.1 using the scaling relation ( 21), (ii) the force pulse at the considered contact is directly computed in the tests T 2 . In the following, the first method is referred to as the "extrapolation method", whereas the second one is referred to as the "direct computation method". It is worth mentioning that the numerical tests T 1 reproduce exactly what has been performed in the experimental process so they allow us a direct comparison between the numerical and the experimental results. The numerical tests T 2 are complementary to the tests T 1 , and they allow us to investigate properly the real force pulse propagating in tapered chains without being disturbed by the presence of a rigid wall. Moreover, the tests T 2 allows us to evaluate the accuracy of the technique of measuring the force pulse at the contacts in a tapered chain carried out in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF].

The bead and wall properties used in the numerical simulations are the same as those given in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]. Concerning the impact velocity and the restitution coefficient, they are not given in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] so it is necessary to determine these parameters by a fitting procedure. The experimental data used for the fitting procedure is the force pulse amplitude at the contacts in the monodisperse chain, which is extracted from figure 5 in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]. It is worth noting that this experimental data are credible because they are obtained by a direct measurement at the contacts in the monodisperse chain without any extrapolation. We carry out the numerical tests T 2 on the monodisperse chain followed by the first tapered chain (with 14 beads and q 1 = 5.6 %) with different values of the restitution coefficient e s and of the impact velocity V s . With e s = 0.965 and V s = 0.62 m/s the numerical result fits the best to the experimental data. It is interesting to note that the restitution coefficient e s = 0.965 determined by the above fitting procedure is one of the three values of the restitution coefficient estimated in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF]. The bead material presented in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] is the same as that presented in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] (stainless steel with the same properties). Moreover, this value of the restitution coefficient gives the best numerical results compared to the experimental data in [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] as analyzed in Sect. 3.1.2. Hereafter, we carry out all the numerical simulations with the restitution coefficient e s = 0.965 and the impact velocity V s = 0.62 m/s.

Force pulses felt at the wall

In this section, we compare the force pulses felt at the contacts between tapered sub-chains of increasing length (number of beads n 1 for q 1 = 5.6 % and n 2 for q 2 = 8.27 %) and the rigid wall obtained from the numerical tests T 1 and from the experimental tests. Figures 10 and11 show the force pulses obtained for q 1 = 5.6 % and for q 2 = 8.27 %, respectively. The experimental data are extracted from figure 1 in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] It can be noted that these force pulses are obtained directly without any extrapolation in the numerical and experimental tests. One can observe firstly that the decreasing tendency of the force pulse amplitude with the increase in the length of the tapered sub-chains is quantitatively well predicted. The numerical amplitudes are quite close to the experimental [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] ones. This will be compared more closely in Sect. 3.2.4. The decrease in amplitude of the force pulses is due to two main effects: the dissipation effect when the collisions between beads are not purely elastic and the dispersion effect introduced by tapering the diameter of the beads. Both effects become more important when the length of the tapered subchains increases. The dispersion effect increases also with the tapering factor so the force pulse amplitudes obtained for q 2 = 8.27 % decrease more quickly than for q 1 = 5.6 %. Besides, the dispersion effect leads to force pulse tails that are formed behind force pulse fronts, at which the contact force oscillates significantly as seen in Figs. 10 and11. The reasons why these force pulse tails are oscillating will be discussed in Sect. 3.2.8. The oscillation at the force pulse tails increases with the length of the tapered sub-chains and with the tapering factor q. On the whole, these oscillating force pulse tails are well reproduced by the numerical simulations, particularly for q 2 = 8.27 % and n 2 ≥ 8. Nevertheless, the force pulse tails are sharper (their base is narrow) in the numerical simulations than in the experimental tests. This can be clearly observed in Fig. 12 where all the force pulses in Figs. 10 and 11 are overlapped by normalizing them by their amplitude F m and their duration τ . One can also observe that the force pulse fronts remain almost unchanged with the length of the tapered sub-chains and with the tapering factor. This suggests that the the force pulses observed in Figs. 10 and 11 might consist of a solitary wave (a soliton observed in monodisperse chains) with decreasing amplitude and duration followed by an oscillating tail.

Force pulse amplitudes

As presented in Sect. 3.2.1, the force pulse amplitude at a given contact in a tapered chain (F chain m,k ) in the experimental ) in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] tests is obtained by the extrapolation method using scaling relation [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF] from that measured at the end of the corresponding tapered sub-chain by the wall sensor (F wall m,k ). In the numerical simulations, the force pulse amplitude at the considered contact can be obtained either by the same extrapolation method as in the experimental tests (with the numerical tests T 1 ) or by the direct computation method (the force pulse amplitude is directly computed at the considered contact in the numerical tests T 2 ). The numerical results are compared to the experimental ones in Fig. 13 for the two considered tapered chains with q 1 = 5.6 % (Fig. 13a) and q 2 = 8.27 % (Fig. 13b). In this figure and the next Figs. 15, 16 and 17, the contacts in the monodisperse chain (resp. tapered chains) correspond to positions < 0.416 m (resp. ≥ 0.416 m).

As shown in Fig. 13, the force pulse amplitude at the contacts in the considered tapered chains obtained from the numerical simulations are quite close to the experimental results for both extrapolation and direct computation methods. However, it can be observed that the extrapolation method leads to a slight overestimation of the force pulse Fig. 14 Scaling factor for the force amplitude versus contact position obtained from the numerical simulation for the two tapered chains with q 1 = 5.6 % and q 2 = 8.27 %, compared to the one used in [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF] amplitude computed directly at the contacts in the tapered chains. This discrepancy can be attributed to the reduction factor of 1.94 used in scaling relation [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF]. In fact, as analyzed in Sect. 3.2.1, this value is determined through numerical simulations performed on non-dissipative monodisperse chains. Considering the same value for tapered chains might be not relevant because dynamics of tapered chains is very different from dynamics of monodisperse chains. To demonstrate this, we compute the scaling factor of force amplitude (F wall m /F chain m ) for the contacts in the two considered tapered chains from the numerical simulations. As shown in Fig. 14, the scaling factor obtained is significantly higher than that used in [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF], in particular for q 2 = 8.27 %. Moreover, the obtained scaling factor increases with the position of the contacts in the tapered chains. Consequently, a constant scaling factor of force amplitude for any tapering factor q and for any contact position may not be relevant to determine the force pulse in a tapered chain by using the extrapolation method presented in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF].

Force pulse duration

In [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] the duration τ of a force pulse is determined by fitting the Nesterenko solution [START_REF] Nesterenko | Dynamics of Heterogeneous Materials[END_REF]

F(t) = F m cos 6 t -T τ , ( 22
)
to the experimental data for the front of the considered force pulse. The force pulse duration at a given contact in a tapered chain (τ chain k

) is experimentally extrapolated using the scaling relation [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF] from that measured at the end of the corresponding tapered sub-chain by the wall sensor (τ wall k ). In the numerical simulations the force pulse duration at a given contact in a tapered chain can be obtained either by the same extrapolation method as in the experimental tests (with the numerical tests T 1 ), or by the direct computation method [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] (with the numerical tests T 2 ). The numerical results corresponding to both methods are compared to the experimental results in Fig. 15 for the two considered tapered chains with q 1 = 5.6 % (Fig. 15a) and q 2 = 8.27 % (Fig. 15b).

As can be seen in Fig. 15, the numerical results obtained with both extrapolation and direct computation methods match generally well the experimental results. One can observe that the numerical simulations slightly underestimate the pulse duration obtained in the experimental tests. On another point, the direct computation and extrapolation methods give close results in terms of force pulse duration. In fact, the scaling factor obtained from the numerical simulations for the force pulse duration (τ wall /τ chain ) is about 1.02, which is slightly lower than the value 1.09 used in [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF], and is approximately constant for any tapering factor q and for any contact position. Consequently, the extrapolation method mentioned above might be relevant to determine the force pulse duration in a tapered chain.

Wave speed

The speed of the wave when it passes through a given bead i is computed as:

v i = Φ i T i -T i-1 , ( 23
)
where T i-1 and T i are respectively the times of flight at the left and right contacts of bead i whose diameter is Φ i . Experimentally, the time of flight at a contact in the monodisperse chain is directly measured by the sensor inserted inside the left bead at the considered contact, and the time of flight at a contact in a tapered chain is measured by the wall sensor in the test with the corresponding tapered sub-chain. Numerically, the time of flight of the force pulse at a given contact in a tapered chain can be either computed at the rigid wall in the tests T 1 or computed directly at the considered contact in the tests T 2 . The numerical results are compared to the experimental results in Fig. 16 for the two considered tapered chains with q 1 = 5.6 % (Fig. 16a) and q 2 = 8.27 % (Fig. 16b).

We can see in Fig. 16 that the wave speed obtained numerically with the time of flight computed at the rigid wall is very close to that obtained with the time of flight computed directly at the contacts in the tapered chains. This is not surprising since the time of flight for a pulse to reach a given contact in a tapered chain is very close to the time of flight for the same pulse to reach the rigid wall that is placed at the end of the corresponding tapered sub-chain. Indeed, to reach the considered contact, the waves in both cases have to travel on the same monodisperse chain and then on the same corresponding tapered sub-chain. It can be observed that the numerical simulations can predict qualitatively well the acceleration of the wave when traveling in the tapered chains. Nevertheless, the discrepancy between the wave speed obtained from the numerical simulations and from the experimental tests is significant particularly for the tapered chain with q 2 = 8.27 %. It is worth noting that the experimental data on the wave speed presented in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] are not suitable for a quantitative comparison, due to very large relative errors when estimating the wave speed from the time of flight difference (see figure 6 in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]). In order to capture a central tendency from the experimental data, the authors in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] use a smoothing procedure. The smoothing procedure consists in fitting the relation v i = Qv i-1 (Q is a constant to be determined) to the experimental data. This relation results from the so-called quasisolitary wave approximation combined with binary collisions [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]. The accuracy of such a relation is not guaranteed. The experimental data that we show in Fig. 16 are obtained with this smoothing procedure so the results should be considered at a qualitative level rather than a quantitative level concerning the wave speed.

Front and tail impulses

Impulses for the front (P F ) and for the tail (P T ) of a force pulse are defined as follows:

P F = T +τ -∞ F(t)dt, ( 24
)
P T = +∞ T +τ F(t)dt, ( 25
)
where the force pulse F(t) is obtained experimentally by the extrapolation method with scaling relation ( 21) and numerically by either the extrapolation method with the numerical tests T 1 or the direct computation method with the numerical tests T 2 .

Let us first compare the front and tail impulses P F and P T obtained with the extrapolation method from the numerical tests T 1 and from the experimental tests. From a qualitative point of view, the numerical simulations pre-Fig. 17 Front and tail impulses P F and P T versus the contact position obtained from the numerical simulations with the extrapolation and the direct computation method (represented by symbols diamond and square, respectively) compared to the experimental data (represented by symbol circle). a for q 1 = 5.6 % and b for q 2 = 8.27 %. The experimental data are extracted from figure 5 (d1 and d2) in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] dict rather well the transfer of the impulse from the front pulse to the tail pulse, as shown in Fig. 17. However, there is a discrepancy between the numerical and experimental results: the numerical simulations underestimate the front impulse P F and overestimate the tail impulse P T . It is worth mentioning that this discrepancy can be observed even for the contacts in the monodisperse chain for which the force pulse amplitude and duration fit very well to the experimental ones (see Figs. 13 and15). This is likely due to the difference between the numerical pulse shape and the experimental one. This discrepancy becomes more significant at the contacts in the tapered chains, which is partly due to the discrepancies in terms of force pulse amplitude and duration obtained from the numerical simulations. Experimentally, force pulses are recorded using a high frequency measurement with a fixed set of points [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]. This technique allows pulse fronts to be fully recorded (high precision) but not pulse tails (the very end of pulse tails is not recorded). Consequently, it is difficult to compare quantitatively the tail impulses P T for the experimental pulses that are not fully recorded and for the numerical pulses that are fully computed in the numerical tests T 1 .

Now we compare the front and tail impulses P F and P T at the contacts in the tapered chains obtained from the numerical simulations with the extrapolation and direct computation methods. It can be observed in Fig. 17 that the front impulses P F obtained with both methods are quite close, whereas the tail impulses P T are very different. The tail impulse P T obtained with the direct computation method is very small with respect to that obtained with the extrapolation method. This result shows that the force pulse tail at a contact in a tapered chain obtained by extrapolating the force pulse computed at the rigid wall in the numerical tests T 1 is much more pronounced than that obtained by computing directly at the considered contact in the numerical tests T 2 . This point will be analyzed in the next section.

Force pulses obtained by different methods

A question that arises here is whether or not the force pulse at a contact in a tapered chain can be experimentally captured by the extrapolation method with scaling relation [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF] (described in Sect. 3.2.1). In other words, whether or not placing a rigid wall at the end of tapered sub-chains might result in some undesirable effects when measuring force pulses in a tapered chain? In Sects. 3.2.4 and 3.2.5 we have shown that the extrapolation method allows us to obtain approximately the force pulse amplitude and duration in a tapered chain. In relation to the force pulse shape, particularly the force pulse tail, does the extrapolation method work as well? This section is dedicated to analyze these questions.

For this purpose, in Fig. 18 we compare normalized force pulses at some contacts in the considered tapered chains obtained by the extrapolation and direct computation methods. It can be observed that the force pulse fronts obtained with the extrapolation and direct computation methods are similar. However, the force pulse tails obtained by the extrapolation method are much more pronounced than those obtained with the direct computation method. Moreover, the tail of the force pulses computed directly is very small compared to their front. This is the reason why the tail impulse P T obtained with the direct computation method is negligible with respect to that obtained with the extrapolation method as shown in Fig. 17 In fact, the pulse in the numerical tests T 1 will hit the rigid wall placed at the end of the tapered sub-chain and then reverses the propagating direction. The pulse propagates now in the direction of increasing bead diameter (we can consider this case as a pulse propagating in an antitapered chain). In this case, part of the pulse propagates toward the largest end and part is reflected once again backward the rigid wall, and so on. This phenomenon was also observed in [START_REF] Doney | Impulse absorption by tapered horizontal alignments of elastic spheres[END_REF]. Consequently, the force pulse computed at the rigid wall has a well pronounced tail with multiple oscillations (multiple compression and expansion phases). It is worth mentioning that LZB multiple impact model is able to handle the phenomenon of multiple compression-expansion phases occurring at certain contact points (see [START_REF] Liu | Frictionless multiple impacts in multibody systems: Part I. Theoretical framework[END_REF][START_REF] Liu | Frictionless multiple impacts in multibody systems: Part II. Numerical algorithm and simulation results[END_REF] for details). On the other hand, the pulse in the numerical tests T 2 propagates in the chain without any reflection. These are the reasons why the force pulses obtained from the direct computation method are different from those obtained from the extrapolation method, particularly for force pulse tails. It seems, from the above numerical results, that the well pronounced force pulse tails obtained in the experimental tests shown in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF], particularly the oscillation at the force pulse tails, result, on one hand, from the mechanical effects of the tapered chains and, on the other hand, from the presence of a rigid wall placed at the end of the tapered subchains.

Summary

In Sects. 3.2.3, 3.2.4, 3.2.5, 3.2.6, 3.2.7 and 3.2.8, we compared the main features of force pulses propagating in two tapered chains with tapering ratio q 1 = 5.6 % and q 2 = 8.27 % such as the pulse shape, the pulse amplitude, the pulse duration, the propagation speed, the front and tail impulses obtained from the experimental tests presented in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] and from the numerical simulations with LZB multiple impact model. The following conclusions can be drawn from the comparison:

-The numerical simulations can reproduce rather well the main features of force pulses in the tapered chains observed in the experimental tests. More precisely, the force pulse amplitude and duration are quantitatively well predicted. -The wave speed measured experimentally is qualitatively well predicted by the numerical simulations. However, it is difficult to compare quantitatively because the experimental technique of measuring the wave speed gives large errors. -The impulse transfer from the front pulse to the tail pulse is qualitatively well predicted by the numerical simulations. Quantitatively, the numerical simulations underestimate the front impulse and overestimate the tail impulse obtained from the experimental tests. The overestimation of the tail impulse by the numerical simulations might result from the fact that force pulse tails are not fully recorded in the experimental tests, whereas they are fully computed in the numerical simulations. -The force pulse tails obtained numerically with the direct computation method are very small compared to those obtained with the extrapolation method. The well pronounced force pulse tails obtained with the extrapolation method in the numerical simulations and in the experimental tests as well, might be due to partly the presence of a rigid wall placed at the end of tapered sub-chains. Therefore, it is likely that the measurement technique described in [START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF] is not suitable to capture the force pulses that actually propagate through a free tapered chain.

Numerical tests on stepped chains

In this section, we turn our attention to impact dynamics within stepped chains which are composed of a monodisperse section of large beads followed by a monodisperse section of small beads. A typical feature of these granular chains is that a solitary wave (SW) propagating in the large monodisperse section turns into a solitary wave train (SWT) composed of many single solitary waves with decreasing amplitudes when it passes through the small monodisperse section. Such a SWT has been first observed in the experimental tests presented in [START_REF] Lazaridi | Observation of a new type of solitary waves in a one-dimensional granular medium[END_REF][START_REF] Nesterenko | Propagation of nonlinear compression pulses in granular media[END_REF][START_REF] Nesterenko | Solitary waves in discrete media with anomalous compressibility and similar to "sonic vacuum[END_REF] when a large enough striker impacts a small monodisperse chain. The formation of the SWT has been numerically confirmed in [START_REF] Lazaridi | Observation of a new type of solitary waves in a one-dimensional granular medium[END_REF][START_REF] Nesterenko | Propagation of nonlinear compression pulses in granular media[END_REF][START_REF] Sokolow | Solitary wave train formation in hertzian chains[END_REF]. Job et al. [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF] have carried a series of experimental tests on stepped chains and have put in evidence the formation of SWTs within such granular chains. In the following, we present the numerical results obtained from the simulations of the tests in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF] performed with LZB multiple impact model and compare with the experimental results available in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF].

4.1 Description of the experimental tests in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF] The experimental tests carried out in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF] are illustrated in Fig. 19. The first stepped chain is composed of two mono- In fact, the restitution coefficient e s and the impact velocity V s are not given in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF] so we apply a fitting procedure to the experimental data to determine these parameters. It consists in changing both the restitution coefficient and the impact velocity such that the amplitude of the incident SW and the first peak of each SWT shown Fig. 20b are the best reproduced from the numerical simulations. We obtain e s = 0.99 and V s = 0.34m/s. The value of the restitution coefficient obtained here is reasonable for steel beads. However, it is higher than the value that we obtained from the fitting procedure presented in Sect. 3.2 for beads made of the same material (e s = 0.965). Figure 20 shows the numerical results obtained for a monodisperse chain of 7 beads of radius R 1 = 13 mm and the two stepped chains mentioned above, compared to the experimental data in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF]. It can be observed that the SWTs observed in the experimental tests in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF] are well reproduced by the numerical simulations. The numerical times of flight of single SWs in each SWT are in good accordance with the experimental data. Moreover, the duration of single SWs is close to the experimental duration. However, the number of significant SWs in each SWT obtained from the numerical simulations is higher than that found in the experi-Fig. 22 The incident force and the interfacial force obtained from a the numerical simulations and b from the experimental tests performed in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF]. The experimental data are extracted from figure 3.a in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF] mental tests (about 7 in the numerical simulations and 5 in the experimental tests). This might be due to the friction effect between the beads and the rail (used in the experimental setup to align the beads) that is not taken into account in the numerical simulations. As can be seen in Fig. 20b, the SWT that propagates in the stepped chain with 25 small beads is attenuated when it propagates in the longer stepped chain with 50 small beads. This attenuation is certainly due to inelastic collisions between beads (dissipation effect). With LZB multiple impact model, the dissipation effect is consistently taken into account so the numerical simulations reproduce well the attenuation of the SWT propagating in the stepped chains as seen in Fig. 20a. However, a discrepancy between the numerical and experimental results can be observed: although the amplitude of the incident SW is overestimated (numerical value of about 119 N compared to the experimental value of about 110 N), the first peak of the first SWT shown in figure 20b is underestimated by the numerical simulations (numerical value of about 63 N compared to the experimental value of about 68 N). One can see in Fig. 21 that the numerical simulations predict well the decrease in amplitude of single SWs in each SWT observed experimentally.

It is interesting to compare the incident and interfacial force pulses. It can be observed in Fig. 22 that the numerical simulations predict well the interfacial force pulse: the front and tail of the interfacial force pulse are very well reproduced, and the numerical force pulse amplitude is close to the experimental one (the numerical value of about 19 N compared to the experimental value of about 20 N). Nevertheless, the incident force pulse amplitude is significantly overestimated as Fig. [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF] Illustration of the experimental tests performed in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF] mentioned previously. Therefore, the ratio of the interfacial amplitude to the incident amplitude obtained from the numerical simulations (about 0.31) is lower than the experimental value (about 0.36).

In conclusion, the solitary wave train that is generated when a solitary wave propagates in a stepped chain is well reproduced by the numerical simulations. The time of flight and the decrease in amplitude of single solitary waves in solitary wave trains obtained from the numerical simulations are in good accordance with the experimental data. The force pulse at the interface is well predicted. However, the number of single solitary waves in the numerical simulations is higher than in the experimental tests. In addition, the incident solitary wave is significantly overestimated by the numerical simulations.

Interaction of two solitary waves in monodisperse chains

It is well known that waves propagating in granular chains are nonlinear waves, due to the nonlinear nature of the interaction between spherical particles [START_REF] Nesterenko | Dynamics of Heterogeneous Materials[END_REF]. One of the main features of the nonlinear waves is that the linear superposition principle does not hold when two nonlinear waves interact each other, i.e. the amplitude of the wave resulting from the interaction is larger than the sum of two individual wave amplitudes. The interaction of two solitary waves (SWs) in monodisperse chains has been numerically investigated in [START_REF] Ávalos | How solitary waves collide in discrete granular alignements[END_REF][START_REF] Liu | Energy dissipation and dispersion effects in a granular media[END_REF][START_REF] Nesterenko | Propagation of nonlinear compression pulses in granular media[END_REF][START_REF] Nesterenko | Dynamics of Heterogeneous Materials[END_REF], and experimentally investigated in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF]. Some interesting phenomena have been discovered: (1) the SWs are recovered after interaction and propagate faster than the respective SWs without interaction, although the latter ones have a larger amplitude than the former ones, (2) two secondary solitary waves (SSWs) are formed after interaction, and (3) the resulting SSWs for a chain of even number of beads (even chain) differ from the case for a chain of odd number of beads (odd chain). In this section, we present numerical simulations of the experimental tests carried out in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF] and compare the numerical results to the experimental data available in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF].

Description of the experimental tests in [50]

The experimental tests in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF] are illustrated in Fig. S 1 is inserted inside bead 9 for the odd chain and inside bead 10 for the even chain, S 2 is inserted inside bead 17 for both chains. In order to compare the SWs after interaction to the SWs without interaction, two reference SWs are initiated in the considered monodisperse chains by striking only at one end with the same striker and the same impact velocity (the other end is free).

Numerical simulations

In order to compare the numerical results to the experimental ones, it is necessary to estimate the values of the restitution coefficient e s and of the impact velocity V s of each striker, that are not given in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF]. For this purpose, we apply a fitting procedure to the experimental data available for the reference SWs shown in figure 2 in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF] for both even and odd chains. The fitting procedure consists in varying the restitution coefficient e s and the impact velocity V s in such a way that the experimental reference SWs are reproduced. From the fitting procedure, we obtained V s = 0.46 m/s for the odd chain, V s = 0.49 m/s for the even chain (little higher than for the odd chain) and e s = 0.95 for both chains. It is worth mentioning that the value of the restitution coefficient obtained here is still reasonable for steel beads but it is a little lower than that obtained in Sect. 3.2 (e s = 0.965). This is likely due to the friction effect between beads and the rail in the experimental setup that is not taken into account in LZB multiple impact model for the frictionless case. In fact, in Sect. 3.2 the fitting procedure is applied over a monodisperse section of about 8 beads, whereas the monodisperse section considered here is about 17 beads. Consequently, the frictional dissipation might be more significant for the case considered here than for the case considered in Sect. 3.2 due to a longer traveling distance of the SW in the first case than in the last case. The numerical simulations presented in this section are performed with the bead properties given in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF] and the restitution coefficient and the impact velocity obtained by the fitting procedure.

Figure 24 shows the SWs before and after interaction obtained from the numerical simulations for the odd and even chains, compared to the experimental data presented in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF]. The SWs that are initiated at the left and right ends of the chains are called LSTW and RSTW (the left-hand and righthand side traveling waves), respectively. It can be observed that the numerical simulations reproduce fairly well what occurs when two SWs collide in a monodisperse chain. For the odd chain, the LSTW arrives to bead 17 earlier than the reference SW with a time difference Δ ≈ 29 μs (Fig. 24aB). This numerical time difference is close to the experimental value Δ ≈ 30 μs (Fig. 24bB). The LSTW and RSTW amplitudes after interaction are both about 12 N (Figs. 24aA, aB), compared to the experimental values of about 12 N (Figs. 24bA, bB). For the even chain, the time difference of the LSTW is Δ ≈ 28 μs (Fig. 24aD), close to the experimental value Δ ≈ 29 μs (Fig. 24bD). The RSTW and LSTW amplitudes after interaction are both about 13 N (Figs. 24aC, aD), compared to the experimental values of about 13 N and 11 N (Figs. 24bC, bD), respectively (so higher difference is observed for the LSTW). The reason why the time difference Δ for the odd chain is higher than for the even chain is likely that the impact velocity V s for the odd chain (V s = 0.46 m/s) is lower than for the even chain (V s = 0.49 m/s) so the wave propagation in the even chain is faster than in the odd chain, which reduces the time difference Δ in the even chain. On another point, in the numerical simulations two secondary solitary waves (SSWs) with small amplitude, that are formed when two incident SWs collide each other, are observed for the even chain (Figs. 24aC, aD) but not for the odd chain, whereas SSWs are also observed experimentally for the odd chain (Figs. 24bA, bB). The experimental observation shown in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF] (Fig. 24b) and the numerical observation in [START_REF] Ávalos | How solitary waves collide in discrete granular alignements[END_REF] have pointed out that the SSWs for the even chain are much stronger than for the odd chain. On the other hand, the amplitude of the SSWs observed in the numerical simulations for the even chain (about 0.4 N) is much smaller than the experimental value (about 1.6 N). A question that arises here is whether or not the dissipation affects the SSWs. To answer this question, we perform numerical simulations on the even chain with different values of the restitution coefficient e s varying from 1.0 to 0.9. As can be seen in Fig. 25, the restitution coefficient e s affects significantly the formation of SSWs in the even chain. When e s = 1.0 (no dissipation), there is a disturbance at the wave tail, due to the collision of two SWs. However, no SSW is observed in this case. When the dissipation increases (e s decreases), the SSW is more clearly formed and the amplitude increases. We perform the same numerical simulations on the odd chain and find out that despite the fact that the dissipation is changed, the SSW is not observed for such a chain. These results show that using the energetic restitution coefficient to take into account the dissipation effect allows us to reproduce the formation of SSWs when two SWs collide in an even monodisperse chain. One may expect to observe more clearly SSWs in monodisperse chains when the friction is incorporated in the numerical simulation.

In [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF] the authors have also performed numerical simulations based on the second order equation of motion and the Hertz's contact model. The local dissipation at the contact points is taken into account by using a nonlinear viscoelastic solid model [START_REF] Kuwabara | Restitution coefficient in a collision between two spheres[END_REF]. Comparing the numerical results shown in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF] (see figure 3 in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF]) to those presented in this section with LZB multiple impact model which takes into account the local dissipation by using the energetic restitution coefficient (see Fig. 24a), we can see that the last ones match better the experimental data than the first ones.

As can be seen in Fig. 24b, the amplitude of the LSTW before interaction is significantly lower than the amplitude of the RSTW, although the LSTW and RSTW are measured at two sensors that are placed in symmetrical positions with respect to the middle of the chains. Therefore, the experimental loading on the chains is not perfectly symmetrical. This means that the impact velocity of the left striker may be lower than that of the right striker. In the numerical simulations, the impact velocity is the same for both strikers. This difference between the numerical and experimental tests is certainly a source of the discrepancy between the numerical and experimental results. Consequently, we should consider the above comparison at the qualitative level rather than at the quantitative level.

In conclusion, the numerical simulations performed with LZB multiple impact model are able to well reproduce the phenomena when two SWs (solitary waves) collide in monodisperse chains. The time shift between the SWs after interaction and the respective SWs without interaction, and the amplitude of the SWs after interaction obtained from the numerical simulations match fairly well the experimental data shown in [START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF]. The secondary solitary waves (SSWs) that are generated after interaction are not observed in the numerical simulations for the odd monodisperse chain but they are observed for the even monodisperse chain. However, the SSWs are significantly smaller than those observed in the experimental tests.

Conclusions

In this paper, we have presented the numerical simulations of impact dynamics for granular chains performed with LZB multiple impact model. The main advantage of using such an impact model is that the dissipation effect resulting from many different sources can be consistently taken into account by using the Stronge's energetic restitution coefficient at each contact point. Different kinds of granular chains have been considered: monodisperse chains, tapered chains, stepped chains. We have carefully compared the numerical results to the experimental results presented in [START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF][START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF][START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF][START_REF] Santibanez | Experimental evidence of solitary wave interaction in hertzian chains[END_REF]. In some cases, some parameters necessary to perform numerical simulations (impact velocity, restitution coefficient) are not given in the experimental tests. In such cases, some fitting procedures have been applied to determine these parameters. The comparison between the numerical and experimental results has shown that the numerical simulations with LZB multiple impact model are able to predict qualitatively and quantitatively well shock dynamics of granular chains. More precisely, the dispersion effect and the force pulse propagation in tapered chains with the main features such as the decreasing amplitude, the decreasing duration, the formation of oscillating pulse tails, the increasing propagation speed, etc. are well predicted. Solitary wave trains in stepped chains are well predicted with the time of flight and the amplitude of single solitary waves in good agreement with the experimental observation. Secondary solitary waves are clearly observed in the numerical simulations when two solitary waves collide in monodisperse chains. Moreover, the time delay and the amplitude of the solitary waves after collision are in good agreement with the experimental data. However, some significant discrepancies between the numerical and experimental results are still observed. These discrepancies might result, on one hand, from the limitation of measurement techniques used in the experimental tests and, on the other hand, from the limitation of the numerical simulations. In fact, the friction effect between beads and the rail used in some experimental setups to align the beads is not taken into account in the numerical simulations, which might be a source of discrepancy. Therefore, the friction effect should be taken into account in the numerical simulation.

The granular chains considered in the present paper are all made of steel beads whose elastic behavior is predominant (dissipation effect is weak). In fact, the dissipation effect in granular chains was investigated in some experimental works, in which beads are made of strongly dissipative materials such as brass, lead, glass, nylon, teflon, etc. [START_REF] Coste | On the validity of Hertz contact law for granular material acoustics[END_REF][START_REF] Daraio | Strongly nonlinear wave dynamics in a chain of polymer coated beads[END_REF][START_REF] Daraio | Strongly nonlinear waves in a chain of teflon beads[END_REF][START_REF] Nesterenko | Dynamics of Heterogeneous Materials[END_REF]. Therefore, it would be interesting to simulate the wave propagation in this kind of granular chains and compare with experimental results in order to know whether or nor the energetic restitution coefficient introduced in LZB model is suitable to take into account the dissipation effect. Moreover, analyzing the dispersion effect in dissipative chains and the influence of the dissipation effect on the dispersion effect will also be the subject of our future works.

Fig. 1

 1 Fig. 1 Configuration of a granular chain with possibly a rigid wall at the end, impacted by a striker (bead 0)

Fig. 4

 4 Fig.[START_REF] Carretero-González | Dissipative solitary waves in granular crystals[END_REF] Configuration of the experimental tests carried out by Nakagawa et al.[START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF]

Fig. 9

 9 Fig. 9 Illustration of the first type of numerical tests

Fig. 10

 10 Fig. 10 Force pulses at the rigid wall obtained from the numerical tests T 1 for different tapered sub-chains with n 1 = 0, 2, 4, 6, 8, 10, 12 (left column) and from the experimental tests (right column) for q 1 = 5.6 %.The experimental data are extracted from figure1in[START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]

Fig. 11

 11 Fig. 11 Force pulses at the rigid wall obtained from the numerical tests T 1 for different tapered sub-chains with n 2 = 0, 2, 4, 6, 8, 10, 12 (left column) and from the experimental tests (right column) for q 2 = 8.27 %. The experimental data are extracted from figure 2 in[START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]

Fig. 12

 12 Fig. 12 Collapse of the pulses presented in Figs. 10 and 11 when the force is normalized to the amplitude F m , and the time scale is (t -T )/τ where T and τ are the measured time of flight and pulse duration, respectively. a for the numerical data and b for the experimental data extracted from figure 4 in[START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]

Fig. 13

 13 Fig.[START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF] Force pulse amplitude versus the contact position obtained from the numerical simulations by both extrapolation and direct computation methods (represented by symbols diamond and open square, respectively), compared to the experimental results (represented by symbol filled circle). a for q 1 = 5.6 % and b for q 2 = 8.27 %. The experimental data are extracted from figure 5 (b1 and b2) in[START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]

Fig. 15

 15 Fig.[START_REF] Gharib | Shock absorption using linear particle chains with multiple impacts[END_REF] Force pulse duration versus the contact position obtained from the numerical simulations by the extrapolation and direct computation methods (represented by symbols diamond and open square, respectively), compared to the experimental results (represented by symbol filled circle). a for q 1 = 5.6 % and b for q 2 = 8.27 %. The experimental data are extracted from figure 5 (c1 and c2) in[START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]

Fig. 16

 16 Fig.[START_REF] Glocker | Concepts for modeling impacts without friction[END_REF] Wave speed versus the contact position obtained from the numerical simulations with the extrapolation and direct computation methods (represented by symbols diamond and open square, respectively), compared to the experimental results (represented by symbol filled circle). a for q 1 = 5.6 % and b for q 2 = 8.27 %. The experimental data are extracted from figure6(c1 and c2) in[START_REF] Melo | Experimental evidence of shock mitigation in a hertzian tapered chain[END_REF]

  .

Fig. 18

 18 Fig.[START_REF] Harbola | Pulse propagation in decorated granular chains: An analytical approach[END_REF] Overlap of normalized force pulses at the right contact of beads 16, 20, 24 in the considered tapered chains obtained by the extrapolation and direct computation methods. Force pulses are normalized by their amplitude and the time scale is (t -T )/τ . a for q 1 = 5.6 % and b for q 2 = 8.27 %

Fig. 19 Fig. 20

 1920 Fig.[START_REF] Harbola | Pulse propagation in randomly decorated chains[END_REF] Impact tests on stepped chains in[START_REF] Job | Solitary wave trains in granular chains: experiments, theory and simulations[END_REF]

Fig. 21

 21 Fig. 21 Amplitude of single solitary waves in the solitary wave trains obtained from the numerical simulations, compared to the experimental data extracted from figure 2 in [23] for the solitary wave trains at the end of the stepped chains with a 25 and b 50 small beads

  23. Two monodisperse chains composed of 25 and of 26 steel beads of radius R = 13 mm are considered. Bead properties are as follows: density ρ = 7780 kg/m 3 , Young's modulus E = 203 GPa and Poisson ratio ν = 0.3. The beads in the chains are initially stationary and barely touch each other. Then the chains are impacted simultaneously at the two ends by two identical beads of radius R s = 4 mm. Force pulses induced by the shocks are measured by two force sensors S 1 and S 2 :

Fig. 24

 24 Fig. 24 The numerical results (a) compared to the experimental data extracted from figure 2 in [50] (b). Left and right columns correspond to the chains of 25 and 26 beads, respectively. In each column, the top panel shows the force felt at the right contact of bead 9 for the odd chain and bead 10 for the even chain, and the bottom panel shows the force felt at the left contact of bead 17 for both chains. LSTW and RSTW are respectively the left-hand and right-hand side traveling waves (initiated on the left and right of the chains), and REF is the reference solitary waves

Fig. 25

 25 Fig. 25 Secondary solitary waves (SSWs) obtained numerically for the even chain with different values of the restitution coefficient e s . The force F is computed at the left contact of bead 17 and is normalized by the the maximal value. The time origin corresponds to the maximal force