Multiple Scattering of Elastic Waves in Granular Media: Theory and Experiments

Acoustics is an discipline that deals with many types of fields wave phenomena. Originally the field of Acoustics was consecrated to the sound, that is, the study of small pressure waves in air detected by the human ear. The scope of this field of physics has been extended to higher and lower frequencies and to higher intensity levels. Moreover, structural vibrations are also included in acoustics as a wave phenomena produced by elastic waves. This book is focused on acoustic waves in fluid media and elastic perturbations in media. Many different systems are analyzed in this book like layered media, solitons, piezoelectric substrates, crystalline systems, granular materials, interface waves, phononic crystals, acoustic levitation and soft media. Numerical methods are also presented as a fourth-order Runge-Kutta method and an inverse scattering method.


Introduction
Granular materials consist of a collection of discrete macroscopic solid particles interacting via repulsive contact forces. Classical examples are sand, powders, sugar, salt and gravel, which range from tens of micrometers to the macroscopic scale. Their physical behaviour involves complex nonlinear phenomena, such as non equilibrium configurations, energy dissipation, nonlinear elastic response, and peculiar flow dynamics (Jaeger et al. (1996)). We know from classical continuum dynamics that when a deformable body is under the action of a uniform external load, the force is transmitted to every point inside the body (Landau & Lifshitz (1999)). Conversely, when we deal with a static packing of granular particles, the way by which the forces are transmitted within the packing remains a complex, and still unresolved problem (Luding (2005)). One important aspect lies on the observation that force networks form the skeleton that carries most of the load in a static granular medium (Majmudar & Behringer (2005)). The mechanical response of granular packings to external perturbations plays a major role in numerous scientific endeavours (such as soil mechanics and geophysics), as well as in industry (oil exploration, structural stability, product formulation in pharmacology and domestics, granular composites, heterogeneous materials, etc.) Most laboratory experiments have been carried out in two-dimensional disc packings using the photoelastic (i.e., birefringent under strain) discs, which have allowed the finest visualization of the generation and dynamical evolution of force chains (Majmudar & Behringer (2005)). However, real granular materials are optically opaque, the photoelastic technique becomes difficult to practice. New tools such as pulsed ultrasound propagation through granular beds under stress have been recently developed to probe the elastic response of three-dimensional granular packings (Brunet et al. (2008a;; Jia et al. (1999); Jia (2004); Jia et al (2009) ;Johnson & Jia (2005); Khidas & Jia (2010)). In particular, by studying the 5 www.intechopen.com low-amplitude coherent wave propagation and multiple ultrasound scattering, it is possible to infer many fundamental properties of granular materials such as elastic constants and dissipation mechanisms (Brunet et al. (2008a;; Jia et al. (1999); Jia (2004); Jia et al (2009)). In soil mechanics and geophysics, the effective medium approach (a continuum mechanics, long wavelength description) has been commonly used to describe sound propagation in granular media (Duffy & Mindlin (1957); Goddard (1990)). However, the experimental investigation performed by Liu & Nagel (1992) on the sound propagation in glass bead packings under gravity revealed strong hysteretic behaviour and high sensitivity to the arrangement of the particles in the container. The authors interpreted this as being due to sound propagating within the granular medium predominantly along strong force chains. In recent years, our understanding of wave propagation in granular materials has advanced both experimentally and theoretically, covering topics such as surface elastic waves, booming avalanches (Bonneau et al. (2007;2008)), earthquake triggering (Johnson & Jia (2005)), and coda-like scattered waves (Jia et al. (1999); Jia (2004); Jia et al (2009)). The understanding of wave motion in granular media took a major step forward with the experimental observation of the coexistence of a coherent ballistic pulse travelling through an "effective contact medium" and a multiply scattered signal (Jia et al. (1999)). Such a picture is confirmed by the experiments of sound propagation in two dimensional regular lattice of spheres under isotropic compression (Gilles & Coste (2003)), where the transmitted coherent signal remains almost unchanged for different packing realizations followed by an incoherent tail which depends on the specific packing configuration, in agreement with the completely disordered three-dimensional case (Jia et al. (1999); Jia (2004)). Similar results were obtained with extensive numerical simulations for two and three dimensional confined granular systems by Somfai et al. (2005). This reconciles both points of view in terms of classical wave propagation in a random medium: (i) At low frequencies such that the wavelengths are very long compared with the correlation length of force chains or the spacing between them, the granular medium is effectively homogeneous continuum to the propagating wave. In this case, nonlinear effective medium theories based upon the Hertz-Mindlin theory of grain-grain contacts describe correctly the pressure dependence of sound velocity observed if one includes the increasing number of contacts with the external load (Goddard (1990); Makse et al. (1999)). Most standard measurements of acoustic velocities and attenuation focus on the coherent propagation of effective waves, which provide a means for characterizing the large-scale properties of the granular medium, though micro structural features are not readily resolved; (ii) At high frequencies when the wavelength decreases down to the order of the grain size, scattering effects caused by the spatial fluctuations of force chains become very significant and the effective contact medium is no longer a valid description. One can observe that the continuous wave trains in the tail portion of the transmitted temporal signal through the granular packing have a broadband strongly irregular high frequency spectrum ("coda waves") (Jia et al. (1999)). The energy of a propagating wave spreads in many directions, and strong interference effects occur between scattered wave that have travelled different path through the medium, resulting in a complicated pattern of nodes and antinodes (i.e. acoustic speckles) (Jia et al. (1999)). The speckles are highly sensitive to changes in the granular medium, and configuration specific, i.e. fingerprints of the structure of the force chains. If the multiply scattered waves are excited between or during a temporal change in the granular packing, then one can exploit the sensitivity of the waves to quantify structural variations of the contact network. This opens the possibility to engineer a novel and useful method for investigating the complex response of granular packings under mechanical perturbations. This new technique could be extended to other amorphous materials such as glasses where heterogeneous force chains have been observed (Tanguy et al. (2002)). Before scattered waves are made more quantitative probe, it is essential first to understand the nature of the wave transport in the granular medium. Quite recently, experimental evidence (Jia (2004); Jia et al (2009)) indicates a well-established diffusive behaviour of the elastic waves transport over long distance scales. These experiments were carried out in confined granular packings under uniaxial loading. The system was excited using high-frequency ultrasonic pulses. A qualitative investigation of the statistical characteristics of scattered waves was performed using the intensity evolution in space and time of a wave train injected into the granular medium. A key observation reported in (Jia (2004)) is that the intensity of scattered waves is very similar to the transmitted pulses of classical waves across strongly scattering random media, in accordance with the diffusive field approximation. A further important observation is that under strong static loading, the normal and shear loads of individual grain contacts exhibit a random distribution. Therefore, the topological disorder of the granular medium induces space fluctuations on both density and elastic stiffness. This opens the possibility to interpret the propagation of ultrasonic waves within granular media in terms of random fluctuations of density and elastic stiffness by employing the same framework used to describe the vibrational properties of heterogeneous materials (Frisch (1968); Sheng (2006); Vitelli et al. (2010)). The modelling of multiple scattering and the diffusive wave motion in granular media is by no means a simple task. However, the energy envelope of the waveform can be interpreted within the framework of the Radiative Transport Equation (RTE) developed for modelling wave propagation in random media ; Weaver (1990)) The energy density calculated using RTE for multiple isotropic scattering processes converges to the diffusion solution over long distance scales, and describes adequately the transport of elastic waves dominated by shear waves in granular media (Jia (2004); Jia et al (2009)). This development is rather heuristic and it lacks a rigorous basis on the wave equation. The potential importance of the original calculation (Jia (2004)) stimulated further investigations towards the construction of a self-consistent theory of transport equations for elastic waves in granular media (Trujillo et al. (2010)). The theory is based upon a nonlinear elasticity of granular media developed by Jiang & Liu (2007), which emphasizes the role of intrinsic features of granular dynamics such as volume dilatancy, mechanical yield and anisotropies in the stress distribution. The formalism developed in (Trujillo et al. (2010)) introduces an extension of the Jiang-Liu granular elasticity that includes spatial fluctuations for the elastic moduli and density, providing a characterization at the grain-scale. In this chapter we give an introduction to elastic wave propagation in confined granular systems under external load. Our analysis of elastic wave scattering is developed from both experimental and theoretical viewpoints. The present systematic description and interpretation of multiple scattering of elastic waves in granular media is based on a synthesis between the experiments carried out by Jia and co-workers from 1999, and the theory constructed by Trujillo, Peniche and Sigalotti in 2010. This chapter is structured in four main parts as follows: In section 2 we elaborate a presentation of the principal experimental outcomes of ultrasound propagation through a granular packing. We start with a brief schema of the laboratory setup and experimental protocol. Then, we present the characteristics of the transmitted signals, observing a coexistence of a coherent ballistic and a speckle-like multiply scattered signal. For over long distances scales, the diffusion approximation is shown to describe adequately the transport of elastic waves dominated by shear waves; In section 3 we elaborate a theory for elastic wave propagation in granular media. As granular materials are disordered systems, we have to introduce several concepts that are unknown in the classical theory of linear elastic waves in homogeneous solids. After a presentation of the Jiang-Liu model for nonlinear granular elasticity, we provide a rational path for the choice of the spatial variations of the elastic constants. Furthermore we derive the equation of motion for elastic waves in granular media and present the vector-field mathematical formalism; In section 4 the mathematical formulation of the problem leads to a theoretic formalism analogous to the analytical structure of a quantum field theory. Then, introducing the disorder perturbation as a small fluctuation of the time-evolution operator associated to a Schrödinger-like equation, the RTE and the related diffusion equation are constructed. This result provides the theoretical interpretation, derived from first principles, of the intensity of scattered waves propagating through granular packing; In section 5 we summarize the relevant conclusions.

Experiments of ultrasound propagation in externally stressed granular packings
This section deals with different experimental aspects of the propagation of acoustic waves in granular media. After a short presentation of the experimental protocol, we describe two observations that demonstrate the presence of multiple scattering in granular packings under stress: (i) the coexistence of a coherent ballistic pulse and a multiply scattered signal; (ii) the intensity of scattered waves is described by the diffuse field approximation.

Experimental setup and procedure
A schematic diagram of the experimental setup is shown in Fig. 1. The samples consist of random packings of polydisperse glass beads of diameter d = 0.6-0.8 mm, confined in a duralumin cylinder of inner diameter 30 mm. The container is closed with two pistons and a normal load P is applied to the granular sample across the top and bottom pistons. To minimize the hysteretic behavior and improve the reproducibility of experiments, one cycle of loading and unloading is performed in the granular packs before the ultrasonic measurements. Statistically independent ensembles of the packing configuration are realized by stirring vigorously glass beads after each measurement and repeating carefully the same loading protocol. The volume fraction of our glass beads packs is 0.63 ± 0.01. The height of the granular sample ranges from L = 5 mm to 20 mm. A generating transducer of 30 mm and a small detecting transducer of 2 mm are placed on the axis at the top and the bottom of the cell. Both the source transducer and detector are in direct contact with glass beads. The excitation is realized by using ten-cycle tone burst excitation of 20 µs duration centered at a frequency f (= ω/2π)=500 kHz is applied to the calibrated longitudinal source transducer. This narrow band excitation corresponds to the product of granular skeleton acoustic wave number and bead diameter, kd = ωd/ν ≈ 2.9, with ν ≈ 750 m/s being a typical sound speed in the solid frame. At such a high frequency, one expects to deal with a strongly scattering medium. The transmitted ultrasonic signals are digitized and signal averaged to improve the signal-to-noise ratio and to permit subsequent analysis of data.

Transmitted signals: Coherent propagation and codalike multiple scattering
In Fig. 2 we show the transmitted ultrasonic field through a granular packing of thickness L = 11.4 mm under axial stress of P = 0.75 MPa, a typical pressure at depths of tens of meters in soils due to the weight of the overburden. To ascertain that the ultrasound propagates from one grain to its neighbors only through their mutual contacts and not via air, we have checked that no ultrasonic signal is detected at vanishing external load. A typical record is presented in Fig. 2 (a), where the transmitted pulse exhibits a primary low frequency (LF) coherent component E P . This coherent pulse E P well defined at the leading edge of the transmitted signal corresponds to a self-averaging effective wave propagating ballistically at compressional wave velocity c P ≈ 1000 m/s, and frequency f P ∼ 70 kHz. After the arrival of the coherent pulse, one observes a continuous wave trains in the tail portion of the signal. This wave trains, which are named "coda" (Fehler & Sato (2003)), looks like a random interference pattern having an envelope whose amplitude gradually decreases with increasing time. This coda type high frequency (HF) incoherent signal S, with frequency f S ∼ 500 kHz, is associated with speckle-like scattered waves by the inhomogeneous distribution of force chains (Jia et al. (1999)). Indeed, at such a high frequency, the acoustic wavelength is comparable to the bead diameter, λ/d ∼ 1.5; thus one expects to encounter strong sound scattering in a granular medium. The sensitivity of the coherent and incoherent waves to changes in packing configurations is shown in Fig. 2 (b) over 15 independent granular samples. In contrast to the coherent pulse, which is self-averaged and configuration insensitive, the acoustic speckles are configuration specific, and exhibits a fluctuating behavior due to the random phases of the scattered waves through a given contact force network. Hence, an ensemble average of configurations can cancel scattered wave signals and leave only the coherent wave E P . Moreover, another coherent signal noted as E S survives from this averaging procedure, which propagates ballistically at a shear velocity about c S ≈ 450 m/s. The inset of Fig. 2 (b) shows that the use of a transverse transducer as source can lead to a considerable enhancement of this shear wave excitation without ensemble averaging thanks to the temporal and spatial coherence.

Transmitted intensity and diffusive wave propagation
Now we present one the most relevant result of this experiment which builds a bridge between experiments and the quest towards a theory of elastic waves propagation in granular media. To investigate quantitatively the statistical characteristics of scattered waves we measure the time-resolved transmitted intensity I(t) through the granular sample. For each configuration we subtract the LF coherent pulses E P and E S from the transmitted ultrasonic field by means of a high-pass (HP) filter ( f ≥ 300 kHz) and determine the intensity of the scattered wave by squaring the envelope of the filtered waveform. In the inset of Fig. 3, we present the corresponding average amplitude profile of scattered wave transmission, which rise gradually from an early time value below the noise level to a maximum and decays exponentially at late times. The ensemble-average is performed over fifty independent configurations realized according to the same protocol of sample preparation. In Fig. 3, we show that the average transmitted intensity I(t),a tL = 11 mm, was found to decay exponentially at long times, with the entire time dependence of I(t) being well described by the diffusion model (Jia (2004)). We conclude that multiply scattered ultrasound propagates, for this experiment, in a normal, diffusive way. In what follows, we will provide theoretical foundations to model the propagation of elastic waves in granular packings. In subsection 4.2.1 we derive the analytical expression for the transmitted intensity, Eq. (58). This expression fits the time profile of the average transmitted intensity (solid line in Fig. 3).

Theory of elastic wave propagation in compressed granular media
We now proceed to show that a theory for elastic wave propagation in granular media can be constructed based upon a nonlinear elasticity developed by Jiang & Liu (2007), which emphasizes the role of intrinsic features of granular materials such as volume dilatancy, mechanical yield and anisotropies in the stress distribution.

Granular elasticity
The rigorous passage from a microscopic to a macroscopic (continuum) mechanical description of granular materials is a challenging task due to the intrinsic disorder of these materials (shape, size, density, contact forces and friction) and the apparent lack of well separated scales between the grain-level dynamics and the entire bulk. An important aspect concerns the anisotropic environment at the particle scale where force chains are clearly evidenced, i.e., chains of contact along which the forces are stronger than the mean interparticle force. The presence of these force chains, implying preferred force paths, has served as empirical argument against an isotropic continuum description of granular matter. However, recent experimental findings on the stress distribution response to localized perturbations have shed some light on the validity of using a continuum elastic theory (e.g., Serero et al. (2001)). On the other hand, Goldhirsch & Goldenberg (2002) showed that exact continuum forms of the balance equations (for mass, momentum and energy) can be established as relations between weighted space (and time) averages. In this framework it is possible to derive exact expressions for the elasticity of disordered solids, such as a granular packing. For granular systems under strong static compression, as is the case for experiments presented in section 2, the theory starts with the assumption of small deformations, i.e., for an infinitesimal deformation the displacements field u = r − r ′ and their gradients are small compared to unity (|u| and |∇u| ≪ 1). This assumption is physically reasonable for small amplitude of wave motion. The displacement associated with deforming the grains, stores energy reversible and maintains a static strain. Sliding and rolling lead to irreversible, plastic process that only heat up the system. So, the total strain tensor ǫ ij may be decomposed into elastic and plastic parts ǫ ij = u ij + u p ij . In this work we limit our analysis to granular packings under strong static compression. Therefore, the granular energy is a function of the elastic energy alone and we can neglect the plastic strain contribution. Up to linear order, the strain tensor is The only way to ensure that the elastic strains are indeed reversible is to specify a strain energy potential F (ǫ ij ) (Helmholtz's free energy), from which the stresses are given as functions of the strains: σ ij = ∂F (ǫ ij )/∂ǫ ij . A material which possesses such a strain energy function is said to be hyperelastic. Jiang and Liu have recently proposed the following free energy strain potential for granular materials where K and G are the compressional and shear moduli, respectively, ǫ ll is the trace of the strain tensor ǫ ij , and ǫ 2 s ≡ ǫ 0 ij ǫ 0 ij , with ǫ 0 ij the traceless part, i.e., ǫ 0 ij = ǫ ij − ǫ ll δ ij /3, and δ ij is the Dirac delta function. As always, summation over repeated indices is assumed implicitly. In the presence of external forces the granular packing deforms leading to a change in the density given by δ = 1 − ρ 0 /ρ = −ǫ ll , where ρ 0 is the density in the absence of external forces. In granular elasticity the elastic moduli are modified to take into account the interaction between grains, which may deform as a result of contact with one another. The elastic moduli K and G are assumed to be proportional to the volume compression ǫ ll , i.e., with K, G > 0 for δ ≥ 0 and K, G = 0 for δ < 0 so that the elastic moduli remain finite. The exponents a and b are related to the type of contact between the grains. That is, when a = b = 0 linear elasticity is recovered, whereas a = b = 1/2 implies Hertz contacts (Landau & Lifshitz (1999)). This formulation provides a much better approximation to granular elastic behavior in which we can specify any type of contact by suitably choosing the exponents a and b. The free energy strain potential (1) is stable only in the range of strains values that keeps it convex. Therefore, Eq.(1) naturally accounts for unstable configurations of the system, as the yield, which appears as a phase transition on a potential-strain diagram.
If there are no changes in temperature (or an analogous granular temperature), the stresses are given as derivatives of the Helmholtz free energy with respect to the strains. The constitutive behavior of a granular packing is completely specified by (1), then we get the following stress-strain relation The above equation contains the stress elements of both the linear and Boussinesq elasticity models as we may see from inspection of the first two terms of the right-hand side of Eq.
(3). In this way, the stress is completely defined once we specify a and b together with the stress-dependent elastic moduli, which account for the desired granular behavior. The strain-stress relation (3) also includes the effects of volume dilatancy. These effects are represented by the second term between parentheses on the right-hand side of Eq.(3), where the pure shear stress is proportional to shear strain and the volumetric deformation δ.

Local disorder and randomness
Apart from the validity of continuum elasticity description, granular packings are heterogenous materials due to the intrinsic disorder. In order to capture these heterogeneities we introduce spatially-varying constitutive relations. To do so we need to know how the particle displacements are related to the local disorder in the deformation of the bulk. This can be elucidated with the aid of the stress-strain relation. In particular, we introduce spatial fluctuations for the elastic moduli, i.e., K(r)=λ(r)+ 2 3 µ(r) and G(r)=µ(r), where λ(r) and µ(r) are the Lamé coefficients and assume that these fluctuations can be described by a random process. Moreover, the stress tensor becomes spatial dependentσ →σ(r). For the range of frequency used in the experiments, which is below the acoustic resonances of individual glass beads, the granular network can be modeled as an effective random network (Jia et al. (1999); Jia (2004)). Due to the inhomogeneous distribution of individual bead contact forces, this network exhibits spatial fluctuations of both density and elastic modulus, closely analogous to an amorphous solid. However, it is worth emphasizing the peculiar position of the granular medium amongst randomly scattering media in general. In fact, the topologically disordered granular packing plays a twofold role. It builds the medium in which elastic wave propagates, but being random it is also responsible for the disorder effects. There is no separation of the system into a reference medium and scatterers. Here we formulate a heuristic approach to calculate the local spatial variations. Let us recall that in the Jiang-Liu elastic theory, the elastic moduli K and G are strain-dependent functions given by Eq.
(2), with δ = δ(r)=−Tr [ǫ ij (r)] being the volume compression. Here Tr() denotes the trace. We note that in absence of compression and shear δ = 0 when the grains are in contact. Therefore, we can define the localized compressional fluctuations at position r as where δ 0 = δ(r) is the imposed bulk compression. Angular brackets designate average expectation values with respect to the corresponding probability distribution. Here Δ(r) is assumed to be a delta correlated Gaussian random process with a zero mean and covariance given by where σ is the strength of the delta correlated disorder. We assume that σ ∼ d (Goldenberg et al. (2007)). Our choice of δ(r) is based on the empirical observation that for isotropically compressed systems the mean normal force (at the grain level) is distributed randomly around an average value with short-range correlation (Majmudar & Behringer (2005)). However, when the system is subjected to an external shear, the force correlations are of much larger and longer range and may be characterized by a power law (Majmudar & Behringer (2005)).
Here we restrict only to the case when the system is compressed isotropically and ignore for simplicity the action of external shear. The spatial variations of the local compression can be estimated from the strain field which, up to linear order, can be calculated using the coarse-graining procedure introduced by Goldhirsch & Goldenberg (2002). The fluctuations in the Lamé coefficients can be expressed in terms of the fluctuating local compression by means of the compressional and shear elastic moduli, i.e., For Hertzian interactions a = b = 1/2 and therefore where λ 0 = K − 2 3 G and µ 0 = G. The coupling between the elastic moduli through the compression δ(r)=−Tr [ǫ ij (r)] reduces the number of free parameters needed to characterize the spatial perturbation to three, namely σ, λ 0 , and µ 0 . This avoids other cross couplings between the elastic moduli (Trégourès & van Tiggelen (2002)) and simplifies the analytical calculations. Finally, let us remark that in the present work the word random implies some kind of statistical or ensemble averaging in the theory. Since a granular packing is a nonequilibrium, quenched disordered medium, we do not have access to a true statistical ensemble. Thus, theory and experiment can only be connected through some kind of ergodicity. The equivalence between theoretical and observational averaging is a very difficult task and will not be addressed in this Chapter.

Equation of motion
Now we are ready to formulate the mathematics of elastic wave propagation in granular media, in a way that is suited to apply the methods of multiple scattering of waves. We start with the equation of motion for the elastic displacement field u at time t and position r, where ρ(r) is the local density and f(r, t) is an external force per unit volume. Mathematically, if ∂R is the boundary of a region R occupied by the granular packing, then u is prescribed on ∂R. This assumption is physically reasonable for the experiments reported in this work where the source and detecting transducers are placed at the surface of the granular medium. From the Jiang-Liu elastic model the spatial dependent stress tensor (3) which, by the Hooke's law, is given in terms of the Lamé coefficients by = λ(r)δ ij δ kl + µ(r) δ ik δ jl + δ il δ jk ǫ kl (r), where the second equality applies to an isotropic medium, in which case the fourth-rank stiffness tensor C ijkl (r) can only have two independent contributions, proportional to the Lamé moduli λ(r) and µ(r). Inserting the third equality of Eq.(8) into Eq. (7) and rearranging terms, the equation of motion in index notation is given by where the symbols ∂ t and ∂ i are used to denote the partial derivatives with respect to time and space. Most of the existing theoretical and numerical models on granular materials currently use the Hertzian force law because it simulates the nonlinear elastic contacts between grains with fairly good approximation. In what follows, we will restrict ourselves only to the case of Hertz contacts (i.e., a = b = 1/2). Then, the dynamics of the displacement fields becomes For pedagogical completeness let us remark that setting a = b = 0, as it would be appropriate for linear elasticity, in vector notation, Eq. (9) reduces to this last equality indicates that elastic waves have both dilatational ∇·u(r, t) and rotational deformations ∇×u(r, t). For isotropic homogeneous media, the Lamé coefficients λ(r) and µ(r) are independent of r, and the above equation further simplifies to the well-known wave equation The terms c P = (λ + 2µ)/ρ(r) and c S = µ/ρ(r) are the compressional wavespeed, and shear or transverse wavespeed, respectively. This proves that the Jiang-Liu granular elasticity includes the well-known linear elasticity of isotropic and homogeneous elastic solids.

Total elastic energy
Now we proceed to calculate the total elastic energy for a compressed granular packing. It is well-known that for a deformable granular packing the deformation represented by the strain tensorǫ is caused by the external forces applied to the packing itself. For simplicity, we shall assume that no appreciable changes on temperature occur within the packing owing to its deformation so that the flux of heat across its boundary ∂R can be neglected. Moreover, since we assume that the contact between the grains are governed by a Hertzian law, the assumption is also made that they deform the packing at a sufficiently slow rate. The potential energy U P due to the displacement field u is given by the strain energy (Chou & Pagano (1992)) where in the first equality the superindex T means transposition. In the second equality we have shown the potential energy in terms of the dilatational and rotational deformations. The different terms in the second equality of (13) represent the compressional E P and shear E S energy. Note that this relation (13) is strictly valid only when the integral is independent of the path of deformation (Norris & Johnson (1997)). This provides a fairly good approximation provided that the displacements of the grains in the packing are assumed to be small enough. Keeping in mind that the kinetic energy is just the volume integral of the quantity 1 2 ρ(∂u/∂t) 2 , we may then evaluate the total energy E T of the elastic displacement u as which works for the Jiang-Liu model specialized to Hertzian contact. In subsection 4.2.2, we will discuss that after several scattering, a stationary regime is set and mode conversion equilibrates: the initial energy spreads with equal probability over all modes. The stabilization of the energy partition ratio E S /E P is a strong indication that coda waves are made of multiple scattered waves (Hennino et al. (2001); Jia (2004); Jia et al (2009); Papanicolaou et el. (1996); Weaver (1990)).

Vector-field mathematical formalism
Working with the wave equation (10) for the analysis multiple scattering of elastic waves in granular media is greatly facilitated by introducing an abstract vector space formed by the collection of vector fields where i = √ −1 is the imaginary number. This vector has 13 components for the running indexes j and k. The most convenient way to perform the algebra of this vector field is by using the Dirac bra-ket notation of quantum mechanics, i.e., |Ψ ≡Ψ and, for the complex conjugate field, Ψ|≡Ψ * . A physical interpretation of vector Ψ follows by realizing that the Cartesian scalar inner product Ψ(r, t)|Ψ(r, t) ≡ d 3 rΨ * (r, t) · Ψ(r, t), is exactly to the total energy E T of the elastic displacement u, defined by Eq. (14). This result suggests that Ψ can be viewed as a complex amplitude for the elastic energy of the granular media. After some algebra, we show that Eq.(10) for the elastic displacement field u is indeed equivalent to a Schrödinger-like equation for Ψ, i.e., where |Ψ f (t) is the external force field defined by the 13-component vector and K is the time-evolution operator given by the 13×13 matrix where we have defined the operator P = −i∇ and introduced the third-rank tensor L jkl ≡ 1/2 P j δ kl + P k δ jl . In Eqs. (18) and (19), in order to gain compactness in writing the vector and matrix representations, the notation (0) T 9 is used to signify a column array consisting of nine zeros. Conversely, (0) 9 is employed to denote a row array of nine zeros filling the right top part of the matrix. In addition, (0) 3×3 and (0) 9×9 indicate square arrays of 3 × 3 and 9 × 9 zeros, respectively. A similar time-evolution operator to Eq. (19) was previously obtained by Trégourès & van Tiggelen (2002) for elastic wave scattering and transport in heterogeneous media, except for the adding term P T 1 2δ µ(r) 2 ǫ jk between square brackets in the middle of the right column of Eq. (19). It arises because of the additional term that appears in the Jiang-Liu formulation of the elastic stress [see Eq. (10)] compared to the traditional expression given by Eq. (12). It is this remarkable difference along with the stress-dependent moduli that allow for a theoretical description of granular features such as volume dilatancy, mechanical yield, and anisotropy in the stress distribution, which are always absent in a pure elastic medium under deformation.

Multiple scattering, radiative transport and diffusion approximation
In the previous section we have presented the main steps to build up a theory for the propagation of elastic waves in disordered granular packings. Now we proceed to develop the rigorous basis to modeling the multiple scattering and the diffusive wave motion in granular media by employing the same mathematical framework used to describe the vibrational properties of heterogeneous materials (Frisch (1968); Karal & Keller (1964); Ryzhik et al., (1996); Sheng (2006); Weaver (1990)). The inclusion of spatially-varying constitutive relations (i.e., Eqns. (4)-(6)) to capture local disorder in the nonlinear granular elastic theory and the formulation of elastic wave equation in terms of a vector-field formalism, Eq. (17), are both important steps to build up a theory of diffusivity of ultrasound in granular media. In this section, we derive and analyze a radiative transport equation for the energy density of waves in a granular medium. Then, we derive the related diffusion equation and calculate the transmitted intensity by a plane-wave pulse.

Radiative transport and quantum field theory formalism
The theory of radiative transport provides a mathematical framework for studying the propagation of energy throughout a medium under the effects of absorption, emission and scattering processes (e.g., ; Weaver (1990)). The formulation we present here is well known, but most closely follows Frisch (1968); Ryzhik et al., (1996); Trégourès & van Tiggelen (2002); Weaver (1990). As the starting point, we take the Laplace transform of Eq. (17) to find the solution where Im(z) > 0, with z = ω + iǫ and ǫ ∼ 0 in order to ensure analyticity for all values of the frequency ω. The operator G(z) is the Green's function G(z) := [z − K] −1 , defined by the equation [z − K] G(z)=Iδ(r − r ′ ), where I is the identity tensor. Physically, it represents the response of the system to the force field for a range of frequencies ω and defines the source for waves at t = 0. A clear introduction to Green's function and notation used here is given in the book by Economou (2006). We shall be mainly interested in two average Green's functions: (i) the configurational averaged Green's function, related to the mean field; (ii) the covariance between two Green's function, related to the ensemble-averages intensity. Mathematical problems of this kind arise in the application of the methods of quantum field theory (QFT) to the statistical theory of waves in random media (Frisch (1968)).
In what follows, we derive a multiple scattering formalism for the mean Green's function (analogous to the Dyson equation), and the covariance of the Green's function (analogous to the Bethe-Salpeter equation). The covariance is found to obey an equation of radiative transfer for which a diffusion limit is taken and then compared with the experiments.

Configuration-specific acoustic transmission
A deterministic description of the transmitted signal through a granular medium is almost impossible, and would also be of little interest. For example, a fundamental difference between the coherent E and incoherent S signals lies in their sensitivity to changes in packing configurations. This appears when comparing a first signal measured under a static load P with that detected after performing a "loading cycle", i.e., complete unloading, then reloading to the same P level. As illustrated in Fig. 4 S is highly non reproducible, i.e., configuration sensitive. This kind of phenomenon arises in almost every branch of physics that is concerned with systems having a large number of degrees of freedom, such as the many-body problem. It usually does not matter, because only average quantities are of interest. In order to obtain such average equation, one must use a statistical description of both the medium and the wave. To calculate the response of the granular packing to wave propagation we first perform a configurational averaging over random realizations of the disorder contained in the constitutive relations for the elastic moduli and their local fluctuations (see subsection 3.1.1). As the fluctuations in the Lamé coefficients λ(r) and µ(r) can be expressed in terms of the fluctuating local compression (see Eq. (6)), then the operator K (Eq. (19)) is a stochastic operator. The mathematical formulation of the problem leads to a partial differential equation whose coefficients are random functions of space. Due to the well-known difficulty to obtaining exact solutions, our goal is to construct a perturbative solution for the ensemble averaged quantities based on the smallness of the random fluctuations of the system. For simplicity, we shall ignore variations of the density and assume that ρ(r) ≈ ρ 0 , where ρ 0 is a constant reference density. This latter assumption represents a good approximation for systems under strong compression, which is the case for the experiments analyzed here. We then introduce the disorder perturbation as a small fluctuation δK of operator (19) so that where K 0 is the unperturbed time-evolution operator in the "homogeneous" Jiang-Liu nonlinear elasticity. Using Eq. (19) along with Eqs. (4)- (6), we obtain after some algebraic manipulations the perturbation operator where Δ 1 = 1 + 5Δ(r)/(4δ 0 ) and Δ 2 = 1 − Δ(r)/δ 0 .

The Dyson equation and mode conversion
We may now write the ensemble average Green's function as where G 0 (ω)=[ω + iǫ − K 0 ] −1 is the "retarded" (outgoing) Green's function for the bare medium, i.e., the solution to (20) when Δ(r)=0. The second equality is the Dyson equation and Σ denotes the "self-energy" or "mass" operator, in deference to its original definition in the context of quantum field theory (Das (2008)). This equation is exact. An approximation is, however, necessary for the evaluation of Σ. The lowest order contribution is calculated under the closure hypothesis of local independence using the method of smoothing perturbation (Frisch (1968)). The expression for Σ is The Green's function is calculated by means of a standard expansion in an orthonormal and complete set of its eigenmodes Ψ n , each with a natural frequency ω n (Economou (2006)). If the perturbation is weak, we can use first-order perturbation theory (Frisch (1968)) and write the expanded Green's function as with The eigenmodes obey the orthonormality condition Ψ n |Ψ m = d 3 rΨ * n · Ψ m = δ nm .A straightforward calculation, employing integration by parts, leads to the mode conversion effective cross-section We may now derive an expression for the scattering mean free-time from Eqs. (26) and (27).
To do so we first recall that the extinction time of mode n is given by 1/τ n = −2ImΣ n (ω) and replace in Eq. (27) the integers n and m by ik i and jk j , respectively, where i and j are the branch indices obtained from the scattering relations that arise when we solve the eigenvalue problem for a homogeneous and isotropic elastic plate (Trégourès & van Tiggelen (2002)). In this way, mode n corresponds to the mode at frequency ω on the ith branch with wave vector k i . Similarly, mode m is the mode on the jth branch with wave vector k j . With the above replacements, the sum Σ m on the right-hand side of Eq. (26) becomes ∑ i A d 2k i /(2π) 2 . Finally, if we use Eq. (27) into Eq. (26) with the above provisions, we obtain the expression for the scattering mean free-time, or extinction time where is the mode scattering cross-section and n i (ω) := k i (ω)/v i is the spectral weight per unit surface of mode i at frequency ω in phase space. In Eq.(29) we have made use of the dyadic strain tensor S = 1/2[∇u +(∇u) T ].

The Bethe-Salpeter equation
To track the wave transport behavior after phase coherence is destroyed by disordered scatterings, we must consider the energy density of a pulse which is injected into the granular medium. We start by noting that the wave energy density is proportional to the Green's function squared. Moreover, the evaluation of the ensemble average of two Green's functions requires an equation that relates it to the effect of scattering. The main observable is given by the ensemble-average intensity Green's function G(ω + ) ⊗ G * (ω − ) , where ⊗ denotes the outer product, ω ± = ω ± Ω/2, where Ω is a slowly varying envelope frequency, and G(ω + ), G * (ω − ) are, respectively, the retarded and the advanced Green's functions. The covariance between these two Green's functions is given by The above equation is known as the Bethe-Salpeter equation and is the analog of the Dyson equation for G(ω + ). It defines the irreducible vertex function U, which is analogous to the self-energy operator Σ. This equation can be expanded in the complete base Ψ n of the homogeneous case. In this base, we find that G(ω + ) ⊗ G * (ω − ) = L nn ′ mm ′ (ω, Ω), which defines the object that determines the exact microscopic space-time behavior of the disturbance, where G(ω + ) ⊗ G * (ω − )=G n (ω + )G * n ′ (ω − )δ nm δ n ′ m ′ . The Bethe-Salpeter equation for this object reads

Radiative transport equation
Equation (32) is formally exact and contains all the information required to derive the radiative transport equation (RTE), but approximations are required for the operator U. Using the method of smoothing perturbation, we have that U nn ′ ll ′ (ω, Ω) ≈ Ψ n |δK|Ψ l Ψ n ′ |δK|Ψ l ′ . In most cases ω >> Ω. Therefore, we may neglect Ω in any functional dependence on frequency. The integer index n consists of one discrete branch index j, with the discrete contribution of k becoming continuous in the limit when A → ∞. In the quasi-two-dimensional approximation we can also neglect all overlaps between the different branches (Trégourès & van Tiggelen (2002)) and use the equivalence . As a next step, we need to introduce the following definition for the specific intensity L nk (q, Ω) of mode jk j at frequency ω where S m (ω) is the source of radiation at frequency ω, which in mode representation can be written as and q = k − k ′ is the scattering wave vector. If we now multiply Eq. (32) by S m S * m ′ and sum over the integer indices m and m ′ , we obtain that According to Eq. (33), we may then write Substituting the above relation into Eq. (35) and performing the summations over the indices n ′ , m, m ′ , and l ′ , we obtain after some algebraic manipulations that Eq. (35) reduces to Note that since the imaginary part of ω n is small, we can drop the term 2Im(ω n ) on the left-hand side of Eq. (58). Moreover, since Ω ≪ ω then ω + ≈ ω − = ω and hence ΔΣ nn = −i/τ nk n (ω). On the other hand, recalling that U nnll = | Ψ n |δK|Ψ l | 2 , replacing n by j, and making the equivalence Finally, if we turn out to the real space-time domain by taking the inverse Fourier transform of Eq. (37), it follows that This is the desired RTE. The first two terms between brackets on the left-hand side of Eq. (38) define the mobile operator d/dt = ∂ t + v j ·∇, where ∂ t is the Lagrangian time derivative and v j ·∇is a hydrodynamic convective flow term, while the 1/τ jk j -term comes from the average amplitude and represents the loss of energy (extinction). The second term on the right-hand side of Eq. (38) contains crucial new information. It represents the scattered intensity from all directions k ′ into the direction k. The object W(jk j , j ′ k j ′ ) is the rigorous theoretical microscopic building block for scattering processes in the granular medium. The first term is a source term that shows up from the initial value problem. The physical interpretation of Eq. (38) can therefore be summarized in the following statement: which mathematically describes the phenomenon of multiple scattering of elastic waves in granular media. This completes our derivation of the transport equation for the propagation of elastic waves in these systems. Remark: For granular media the contribution to the loss of energy due to absorption must be included in the extinction time 1/τ j . We refer the reader to Brunet et al. (2008b) for a recent discussion on the mechanisms for wave absorption. Whereas in the context of the nonlinear elastic theory employed in the present analysis intrinsic attenuation is not explicitly considered (similar to the "classical" elastic theory), its effects can be easily accounted for by letting the total extinction time be the sum of two terms: 1/τ j = 1/τ s j + 1/τ a j , where 1/τ a j is the extinction-time due to absorption. A rigorous calculation of this term would demand modifying the scattering cross-section (Papanicolaou et el. (1996)), implying that the non-linear elastic theory should be extended to account for inelastic contributions. In this chapter we do not go further on this way and keep the inclusion of the extinction time due to absorption at a heuristic level.

Diffusion equation
Now we derive the form of Eq. (38) in the diffusion limit and solve it to study the diffusive behavior of elastic wave propagation in granular media. Integrating Eq. (38) overk and performing some rearrangements we obtain the equation where is the spectral energy density (or fluence rate) U i , is the current density (or energy flux) J i , and is the mode conversion matrix C ij . The diffusion approximation is basically a first-order approximation to Eq. (38) with respect to the angular dependence. This approximation assumes that wave propagation occurs in a medium in which very few absorption events take place compared to the number of scattering events and therefore the radiance will be nearly isotropic. Under these assumptions the fractional change of the current density remains small and the radiance can be approximated by the series expansion L ik (q, Ω) ≃ 1 n i U i (q, Ω)+ 2 where the zeroth-order term contains the spectral energy density and the first-order one involves the dot product between the flow velocity and the current density; the latter quantity being the vector counterpart of the fluence rate pointing in the prevalent direction of the energy flow. Replacing this series approximation into Eq. (38) produces the equation From the above assumptions we can make the following approximations: ∂ t U i → 0 and Moreover, we can also neglect the contribution of 1/τ a ik i . The absorption term modifies the solution of the scattering cross-section making it to decay exponentially, with a decay rate that vanishes when τ a ik i → ∞ (Papanicolaou et el. (1996)). Furthermore, noting that in the diffusive regime U j /n j ≈ U i /n i , the above equation can be manipulated and put into the more convenient form It is evident from this equation that we can define the diffusion matrix as which allows us to express the current density J i as a generalized Fick's Law: A generalized diffusion equation then follows by combining the continuity-like equation (40) with the Fick's law (47), which reads where the source S i (ω) is defined by the integral S i (ω)=n i d 2k 2π |S ik i (ω)| 2 . At this point it is a simple matter to derive the diffusion equation for the total energy density U = ∑ i U i . Summing all terms in Eq. (48) over the index i, introducing the definitions: S(ω)=∑ i S i (ω) for the total source along with D(ω) := ∑ ij D ij (ω)n j / ∑ j n j , for the total diffusion coefficient, and ξ := 1 τ a = ∑ i n i τ ik a i / ∑ i n i , for the total absorption rate, and noting that where we have made use of Eq. (43), we finally obtain the time-dependent equation which describes the diffusive propagation of elastic waves.

Transmitted intensity
In section 2.3, Fig.3 we showed that the averaged transmitted intensity I(t) decays exponentially at long times. This picture is reminiscent of the diffusively transmitted pulses of classical waves through strongly scattering random media (Sheng (2006); Snieder & Page (2007); Tourin et al. (2000)). This is the main result of the present work, which stimulated the construction of the theory for elastic wave propagation in granular media presented above. Now we conclude our analysis with the derivation of the mathematical formula for the transmitted intensity I(t), corroborating that it fix very well with the experimental data.
In the experiment the perturbation source and the measuring transducer were placed at the axisymmetric surfaces and the energy density was measured on the axis of the cylinder. We can make use of Eq. (49) to calculate the analytical expression for the transmitted flux. In order to keep the problem mathematically tractable we assume that the horizontal spatial domain is of infinite extent (i.e., −∞ < x < ∞ and −∞ < y < ∞), while in the z-direction the spatial domain is limited by the interval (0 < z < L). The former assumption is valid for not too long time scales and for a depth smaller than half of the container diameter. With the use of Cartesian coordinates, a solution to Eq. (49) can be readily found by separation of variables with appropriate boundary conditions at the bottom (z = 0) and top of the cylinder (z = L). The separation of variables is obtained by guessing a solution of the form U(x, y, z, t)= U x (x, t)U z (z, t). It is not difficult to show that if the surface of the cylinder is brought to infinity, Eq. (49) satisfies the solution for an infinite medium It is well known that for vanishing or total internal reflection the Dirichlet or the Neumann boundary conditions apply, respectively, for any function obeying a diffusion equation with open boundaries. In the case of granular packings we need to take into account the internal reflections. In this way, there will be some incoming flux due to the reflection at the boundaries and appropriate boundary conditions will require introducing a reflection coefficient R, which is defined as the ratio of the incoming flux to the outgoing flux at the boundaries (Sheng (2006)). Mixed boundary conditions are implemented for the z-coordinate, which in terms of the mean free path l * are simply (Sheng (2006)): where the coefficient c ≡ 2l *

1+R
1−R . Therefore, the solution for U z reads where Z n (z)= sin (α n z) + κβ n cos (α n z) with α n = β n /L, κ = c/L, and the discrete values of β n determined by the roots of tan β n = 2β n κ/(β 2 n κ 2 − 1). Finally, using Eqs. (50) and (53) we can ensemble the solution for the total energy density where The total transmitted flux at the top wall of the cylinder can be readily calculated by taking the z-derivative of E as defined by Eq. (55) and by evaluating the result at z = L to give If, as mentioned by Jia (2004), the reflectivity of the wall is very high, then R ≈ 1. In the full reflection limit the following limits can be verified: κ → ∞, tan β → 0 =⇒ β n = nπ for all n = 0, 1, 2, . . . , lim κ→∞ C n = 0, and lim κ→∞ κC n =( 2/β n L) cos(α n z 0 ). For a plane-wave source we need to integrate Eq. (57) over x =(x, y) to obtain where v is the energy transport velocity and z 0 ≈ l * . This equation tells us that the flux transmitted to the detector behaves as I(t)=vU/4, when R ≈ 1. This result provides the theoretical interpretation of the acoustic coda in the context of the present radiative transport theory and assesses the validity of the diffusion approximation for a high-albedo (predominantly scattering) medium as may be the case of granular packings.

Energy partitioning
In section 3.3 we have shown that the total energy E T is given by Eq. (14), and that the Cartesian scalar inner product of the vector field Ψ is exactly to the total energy. In the diffusion limit, the conversion between compressional E P and shear E S energies equilibrates in a universal way, independent of the details of the scattering processes and of the nature of the excitation source. The energy ratio is governed by the equipartition of energy law, K = E S /E P = 2(c P /c S ) 3 , where the factor 2 is due to the polarization of the shear waves (Jia et al (2009);Papanicolaou et el. (1996); Ryzhik et al., (1996); Weaver (1990)). For typical values of c P /c S ≥ √ 3, the equipartition law predicts the energy ratio K ≥ 10. This shows that in the diffusive regime the shear waves dominate in the scattering wave field, which is observed in seismological data (Hennino et al. (2001); Papanicolaou et el. (1996)). The diffusion coefficient D is a weighted mean of the individual diffusion coefficients of the compressional wave D P and shear wave D S : D =( D P + D S )/(1 + K). With the weights K ≥ 10 the diffusion coefficient is approximated to D ≈ D S = c S l * S /3. This demonstration confirms the applicability of the diffusion equation for describing the multiple scattering of elastic waves (Jia (2004)).

Conclusions
In summary, the experiments presented in this chapter permit one to bridge between two apparently disconnected approaches to acoustic propagation in granular media, namely, the effective medium approach (Duffy & Mindlin (1957); Goddard (1990)), and the extreme configuration sensitive effects (Liu & Nagel (1992)). This unified picture is evidenced in fig.2 with the coexistence of a coherent ballistic pulse E P and a multiply scattered signal S. The coherent signal was shown to be independent of the packing topological configuration, whereas the coda-like portion of the signal behaves like a fingerprint of the topological configuration as showed in fig.2 b and fig.4. The experimental confirmation of the applicability of the diffusion approximation to describe the multiple scattering of elastic waves through a compressed granular medium, was decisive to guide the construction of the theoretical model for elastic waves propagation. We have shown that the nonlinear elastic theory proposed by Jiang & Liu (2007) can be used to derive a time-evolution equation for the displacement field. Introducing spatial variations into the elastic coefficients λ and µ, we were able to describe the disorder due to the inhomogeneous force networks. The link between the local disorder expressed through the constitutive relations, and the continuum granular elastic theory, permit us to put together within a single theoretical framework the micro-macro description of a granular packing. The mathematical formulation of the problem leads to a vector-field theoretic formalism analogous to the analytical structure of a quantum field theory, in which the total energy satisfies a Schrödinger-like equation. Then, introducing the disorder perturbation as a small fluctuation of the time-evolution operator associated to the Schrödinger-like equation, the RTE and the related diffusion equation have been constructed. We have shown that the temporal evolution for the averaged transmitted intensity I(t), Eq.58, fits very well with the experimental data presented in fig.3, providing the theoretical interpretation of the intensity of scattered waves propagating through a granular packing. This opens new theoretical perspectives in this interdisciplinary field, where useful concepts coming from different areas of physics (quantum field theory, statistical mechanics, and condensed-matter physics) are now merging together as an organic outgrowth of an attempt to describe wave motion and classical fields of a stochastic character. As perspectives for future research let us mention the study of the evolution of the wave transport behavior in a more tenuous granular network when the applied stress is decreased. The disordered nature of a granular packing has a strong effect on the displacements and forces of in individual realizations, which depends on the intensity of the external loads. A hot topic is the study of acoustic probing to the jamming transition in granular media (Vitelli et al. (2010)). This is related to anisotropic effects and the emergence of non-affine deformations of the granular packing. It is necessary a systematic study of the transport properties of elastic waves between the different regimes of external load: strong compression ↔ weak compression ↔ zero compression. The last one is related to the behavior of waves at the free surface of the granular packing (Bonneau et al. (2007;2008) ;Gusev et al. (2006)). The propagation of sound at the surface of sand is related to the localization of preys by scorpions (Brownell (1977)) and the spontaneous emission of sound by sand avalanches (the so-called song of dunes) (Bonneau et al. (2007)). We believe that the experiments presented in this chapter point out to the considerable interest in acoustic probing as a tool for studying of the mechanical properties of confined granular media. Clearly, before this can be undertaken, one should study in detail the sensitivity of the acoustic response to configurational variations. On the other hand, the present theory represents a powerful tool to understand complex granular media as, for example, sedimentary rocks whose geometrical configuration is affected by deposition ambients, sediments, accommodation phase, lithostatic overburden, etc. This explains why anisotropy is always present and characterization is so difficult. Therefore, the study of acoustic waves in such complex media gives useful information to sedimentologists. It can also be applied to important oil industry issues such as hole stability in wells. Important geotechnical applications involve accurate seismic migration, seismo-creep motions, and friction dynamics. Finally, let us mention the similarity between the scattering of elastic waves in granular media with the seismic wave propagation in the crust of Earth and Moon (Dainty & Toksöz (1981); Hennino et al. (2001); Snieder & Page (2007)). In particular, the late-arriving coda waves in the lunar seismograms bear a striking resemblance to the multiple scattering of elastic waves in the dry granular packing. Some features of the laboratory experiments may be used to explain some seismic observations in the high-frequency coda of local earthquakes in rocky soils and the granular medium may be useful as model system for the characterization of seismic sources. Acoustics is an discipline that deals with many types of fields wave phenomena. Originally the field of Acoustics was consecrated to the sound, that is, the study of small pressure waves in air detected by the human ear. The scope of this field of physics has been extended to higher and lower frequencies and to higher intensity levels. Moreover, structural vibrations are also included in acoustics as a wave phenomena produced by elastic waves. This book is focused on acoustic waves in fluid media and elastic perturbations in heterogeneous media. Many different systems are analyzed in this book like layered media, solitons, piezoelectric substrates, crystalline systems, granular materials, interface waves, phononic crystals, acoustic levitation and soft media. Numerical methods are also presented as a fourth-order Runge-Kutta method and an inverse scattering method.