Kinetic model of elastic waves in anisotropic media

Kinetic modeling of multiple scattering of elastic waves in heterogeneous anisotropic media

Ibrahim Baydoun Laboratoire MSS-Mat, École Centrale Paris, UMR 8579 CNRS, France Éric Savin ONERA–The French Aerospace Lab, France Régis Cottereau Laboratoire MSS-Mat, École Centrale Paris, UMR 8579 CNRS, France Didier Clouteau Laboratoire MSS-Mat, École Centrale Paris, UMR 8579 CNRS, France  and  Johann Guilleminot Laboratoire Modélisation et Simulation Multi-Échelle, Université Paris-Est, UMR 8208 CNRS, France
July 5, 2019

In this paper we develop a multiple scattering model for elastic waves in random anisotropic media. It relies on a kinetic approach of wave propagation phenomena pertaining to the situation whereby the wavelength is comparable to the correlation length of the weak random inhomogeneities–the so-called weak coupling limit. The waves are described in terms of their associated energy densities in the phase space position  wave vector. They satisfy radiative transfer equations in this scaling, characterized by collision operators depending on the correlation structure of the heterogeneities. The derivation is based on a multi-scale asymptotic analysis using spatio-temporal Wigner transforms and their interpretation in terms of semiclassical operators, along the same lines as Bal [Wave Motion 43, 132-157 (2005)]. The model accounts for all possible polarizations of waves in anisotropic elastic media and their interactions, as well as for the degeneracy directions of propagation when two phase speeds possibly coincide. Thus it embodies isotropic elasticity which was considered in several previous publications. Some particular anisotropic cases of engineering interest are derived in detail.

Key words and phrases:
Anisotropic elasticity, Elastic waves, Kinetic model, Transport equation, Radiative transfer
Corresponding author: É. Savin, ONERA–The French Aerospace Lab, 29 avenue de la Division Leclerc, F-92322 Châtillon cedex, France (

1. Introduction and summary

1.1. Modeling of wave propagation phenomena in random media

The study of multiply-scattered elastic waves in heterogeneous, anisotropic media has relevance to non destructive evaluation of materials and structures, seismic waves characterization, acoustic emission and backscattered echo analyses, with possible applications in geophysical prospection, biomedical imaging, or structural health monitoring, among others. In this respect, the use of ultrasound to infer the microstructure of polycrystalline materials has been widely considered in the past since the earlier work of Mason & McSkimin [37]. The nondestructive techniques elaborated afterwards are based on the measurement of exponential rates of spatial decay (attenuations) and speeds of averaged plane waves, i.e. coherent fields. The difficulty raised by this approach is the impossibility to distinguish the various sources of potential decays between scattering, geometrical spreading, internal absorption, or the influence of the reflections at the opposite faces of the sample. These shortcomings have prompted the development of probing techniques based on the measurement of the evolution of the incoherent part of ultrasonic waves, i.e. multiply scattered, possibly diffusive fields [20]. The earlier attempts in this direction can be tracked back to the works of Guo et al. [22] and Weaver [54]. This alternative approach has received a considerable attention in the last decade since it was observed that the empirical cross-correlations of such diffuse fields could be directly related to the Green function of the propagation medium [14, 55, 11]. The numerous developments and applications in geophysics which have followed are described in e.g. [12] and references therein.

These advances call for accurate analytical models of ultrasonic wave propagation phenomena in unstructured or structured heterogeneous media. Iterative perturbation expansions for weakly random polycrystalline materials were considered in [46, 25, 26, 53] in the spirit of the seminal developments of Karal & Keller [29]. Other approaches are based upon diagrammatic expansions in which the mean response is governed by a Dyson equation, and the mean square response is governed by a Bethe-Salpeter equation [18, 52]. Both are analytically intractable unless low-order truncations are enforced, typically a so-called first-order smoothing approximation (FOSA) for the Dyson equation, and a so-called ladder approximation for the Bethe-Salpeter equation. These approximations have been used in [54, 51, 50, 48, 49] for the derivation of (i) scattering-based attenuation coefficients on one hand, and of (ii) radiative transfer equations for the (renormalized) mean square wave fields in the limit of small wavelengths (high frequencies) with respect to the macroscopic features of the medium on the other hand. The works evoked above have been mainly confined to untextured or textured aggregates of cubic-symmetry crystallites. Besides, radiative transfer models of high-frequency wave propagation in heterogeneous media have a broad range of validity and were derived in many fields from either phenomenological principles [13, 27, 44, 41] or, more recently, systematic formal multiple-scale asymptotic expansions [40, 23, 10, 6, 7].

Transport and radiative transfer equations describe the mesoscopic regime of wave propagation when the wavelength is comparable to the characteristic length of the heterogeneities, typically a correlation length in a random medium (hereafter referred to as the fast scale). It corresponds to a situation of strong interaction between waves and random heterogeneities which cannot be addressed by usual homogenization and multi-scale techniques. Also it considers large propagation distances compared to the wavelength, and weak amplitudes of the random perturbations of the material parameters with respect to a bare, possibly heterogeneous, background medium varying at a length scale (the slow scale) one order of magnitude larger than the wave/correlation lengths. This corresponds to the so-called weak coupling limit as defined in the dedicated literature, whereby an explicit separation of scales can be invoked. The analysis developed in [40, 23, 10, 6, 7] is based on the use of a Wigner transform of the wave field, of which high-frequency, non-negative limit captures the angularly resolved energy density in time and space. It can be made mathematically rigorous as in [32, 19, 1] ignoring however the influence of random inhomogeneities, except for some particular situations [15, 34].

The purpose of the research presented in this paper is to assess the influence of material full anisotropy on the radiative transfer regime of elastic waves in randomly heterogeneous media. Anisotropy is considered at two levels. The first one is related to the constitutive law of random materials. The second level is related to the correlation structure of these random materials, referred to as anisomery in the dedicated literature [35]. More specifically we have developed formal models for the consideration of anisotropy in the collision kernels of the radiative transfer equations pertaining to multiply-scattered elastic waves. These models describe the evolution of their energy density in the phase space position wave vector in terms of the Wigner measure of the wave fields. The analysis follows to a great extent the techniques used by Bal [6] and Akian [1] in that it handles a second-order wave equation and introduces the spatio-temporal Wigner transform of the elastic waves. However, as opposed to [6] it considers vector wave fields, and as opposed to [1] it considers the influence of random perturbations in the weak coupling regime. Therefore our derivation generalizes those previous works and the classical reference [40] to fully anisotropic bare elastic media with fully anisotropic random perturbations.

1.2. Summary of the main results

We now summarize our main results. We aim at describing elastic waves in a random medium taking into account the non-uniformity of the background medium and their scattering by random inhomogeneities, with due consideration of the effects of coupling between their different polarizations. We also consider the regime where the leading wavelength is comparable to the (small) correlation length of the heterogeneities, in order to ensure maximum interactions with the waves. This is a necessary condition if we want to probe the medium and its fluctuations. It defines the high-frequency range terminology we shall use throughout the paper. At last, our fundamentally new results are that this objective is achieved for arbitrary anisotropy of the random medium. As a scalar wave propagates in a random medium with an incident wave vector , it can be scattered at any time and position into any direction and wave vector (such that ). Therefore it is relevant to consider an angularly resolved scalar energy density for this wave, defined for all positions in phase space. In [40, 6] it is shown that energy conservation takes the form of a scalar radiative transfer equation:


where is the frequency of the waves at with wave vector , and is the rate of conversion of energy with wave vector into energy with wave vector at position –the so-called scattering cross-section. The total scattering cross-section is:

such that the transport equation is conservative because the former relationship yields:

for all times. The scattering cross-section is explicitly determined by the power spectral density of the inhomogeneities [40, 6]. The transport equation (1) also holds when the waves are scattered by randomly distributed discrete scatterers, in which case the scattering cross-section is the cross-section of a single scatterer multiplied by their density. Here we only consider continuous random inhomogeneities.

For vector waves we must in addition keep track of their state of polarization. In a three-dimensional anisotropic elastic medium three orthogonal polarization directions exist, corresponding to at most three different directionally-dependent phase velocities: one for quasi-longitudinal compressional wave, and two for quasi-transverse shear waves. Labeling the polarization states by or , each one has its own energy density but it may be converted to any other state and any other direction at any position when scattered by the inhomogeneities. Conservation of energy is now expressed in terms of coupled radiative transfer equations for the energy densities of the different polarizations:


where is the frequency of the waves at with wave vector and polarization state , and is the rate of conversion of energy with wave vector and polarization state into energy with wave vector and polarization state , at position . The total scattering cross-sections are:

such that the coupled transport equations are conservative since:

in the vector case. The scattering cross-sections are explicitly determined in terms of the power spectral density of the inhomogeneities–which may be characterized by up to coefficients in a triclinic medium. They were originally derived in [40] for isotropic media solely, such that the quasi-transverse shear velocities are the same everywhere and independent of the propagation direction. Here we obtain generalized formulas for them accounting for all possible symmetry classes of anisotropic media, namely Eq. (68) together with the definitions of Eq. (50) and Eq. (39). In this latter equation the elasticity tensor of the background medium and the elasticity tensor of the random inhomogeneities are formally allowed to belong to different symmetry classes, and to vary at different scales. This implies that they have different contributions to the wave dynamics:

  • The variations of at the slow scale contribute only to the left-hand side of the radiative transfer equations (1) and (2). They basically characterize the phase velocities and polarizations in this regime.

  • The variations of contribute only to the right-hand side of these radiative transfer equations in terms of the power spectral densities of its elasticity coefficients. It describes how high-frequency waves are continuously scattered by the material inhomogeneities at the fast scale, which is also their (small) wavelength. These collision operators account for both the elastic (change of direction without changing polarization) and inelastic (with a change of polarization) processes.

In Eq. (2) the energy densities are scalars as long as the phase velocities are all distinct. Our derivation considers the most general case when they may coincide for some modes, though. This is the case for isotropic media for example, for which Eq. (2) becomes a matrix system. This situation is fully addressed in the subsequent analyses.

1.3. Outline

The rest of the paper is organized as follows. In Sect. 2 we introduce the basic framework and notations used throughout. We more particularly focus on the characterization of anisotropy in terms of the acoustic, or Christoffel tensor of the bare medium, and the relevance of using a Wigner transform and its non-negative limit measure for the analysis of multiple scattering phenomena in the high-frequency range. This limit is simply the energy density introduced above for each polarization state . The corresponding transport model is derived in detail in Sect. 3 ignoring the influence of random inhomogeneities in a first step. The spatio-temporal Wigner transform and the formal mathematical tools used for the subsequent analyses are also introduced there. The main contribution of the paper is Sect. 4 which outlines the extension of the previous transport model to account for anisotropic random inhomogeneities. Matrix radiative transfer equations are obtained in the most general case. Their collision operators are explicitly described in terms of the correlation structure of the heterogeneities for all possible symmetries arising in elastic constitutive relations. In this respect, it should be noted that the proposed theory requires a full characterization of the power spectral densities of the tensorial random fluctuations (assumed to be statistically homogeneous at the small wavelength scale), but no other statistical information. These data may be obtained from the random matrix theory of the elasticity tensor developed recently (see [47, 21] and references therein) for example. We will however not pursue the analysis presented here in that direction, considering that simplified correlation models as classically encountered in the literature are sufficient for the purpose and scope of this paper. Using these models, the collision kernels (the scattering cross-sections) for some selected material symmetries are plotted in Sect. 5 in order to illustrate our results. Some conclusions and perspectives are finally drawn in Sect. 6.

2. Elastic bulk wave propagation in a high-frequency setting

In this section we establish the high-frequency setting we are interested in for the derivation of the multiple-scattering kinetic (transport) model that will be detailed in the subsequent parts of this paper. The material below is essentially adapted from [43]. The primary objective is to introduce the main notations that will be used throughout the paper.

2.1. Elastic wave equation

We first recall how the vector wave equation arises in an elastic medium occupying an open domain , where stands for the usual three-dimensional Euclidean space. That medium is constituted by a heterogeneous linear viscoelastic material, of which density is denoted by and fourth-order relaxation, or elasticity tensor is denoted by , for . Its displacement field about a quasi-static equilibrium considered as the reference configuration is denoted by , and its second-order Cauchy stress tensor is denoted by . The subscript stands for the (small) spatial scale of variation of the initial conditions that will be imposed to the materials and will be propagated to the displacement and stress fields at later times by hyperbolicity. Then the balance of momentum in a fixed reference frame ignoring the action of body forces reads:


Here the divergence of a second-order tensor is defined by for any constant vector , and is the gradient vector with respect to . At last stands for the matrix transpose. The initial conditions for Eq. (3) read:


They are parameterized by the small parameter , which quantifies the rate of change of and with respect to the dimensions of or the propagation/observation distances. Since high-frequency waves will be generated by an initial vibrational energy oscillating at a scale , the functions and shall be considered as strongly -oscillating functions in the sense of Gérard et al[19]. The plane waves and for a given wave vector and , typically fulfill this condition.

In addition, the stress field is given as a function of the linearized strain tensor by the material constitutive equation:


Both and are parameterized by the scale of the applied loads since also is. Here is the symmetrised tensor product of two vectors . The elastic wave equation for is derived plugging this relation into Eq. (3) and thus reads:


Indeed, since and are symmetric the elasticity tensor shall satisfy the minor symmetries . Invoking a thermodynamical reversibility argument, it also satisfies the major symmetry .

2.2. The Christoffel tensor

Let us define the matrix by:


and the symmetric, positive semi-definite matrix which is constituted by the independent coefficients of the elasticity tensor . More precisely, is a block matrix of which block is the matrix with elements . Then the second-order acoustic, or Christoffel tensor of the propagation medium is defined by:


where stands for the conjugate transpose matrix. It is symmetric, real and positive definite in . So for a given direction on the unit sphere of it has three real positive (possibly multiple) eigenvalues for or , and the associated eigenvectors can be chosen real and orthogonal. They can also be normalized such that:


where is the identity matrix of . The eigenvalues and eigenvectors of the Christoffel tensor correspond to the phase velocities and polarizations of plane waves propagating along the direction in the medium . Here we do not assume any particular ordering of the eigenvalues with indices , and . The polarization with the closest direction to is called the quasi-longitudinal wave, the other two components being the quasi-transversal waves. Also in view of Eq. (8), the eigenvalues have the form where the ’s have dimension of celerities. One has in addition . The above spectral expansion of the Christoffel tensor is valid for all directions but the so-called acoustic axes [4, 8, 39], along which two eigenvalues may coincide. For such a direction of degeneracy where, say , the expansion reduces to:


It follows from Eq. (10) that any vector which is orthogonal to is an eigenvector of the Christoffel tensor , i.e. it is an allowed polarization for a wave propagating along with a wave celerity . A common example is elastic isotropy, when where and are the Lamé parameters and , for standing for the symmetric part of a square matrix . Then the Christoffel tensor reads:


where and are the velocities for compressive and shear waves, respectively. The various properties of the Christoffel tensor for anisotropic media are discussed in e.g. [4, 3, 9, 28] and references therein.

2.3. High-frequency setting

The high-frequency limit in the previous setting shall be considered for quadratic observables of the wave displacement field as explained now. We rely on the simple following example which was already discussed in [42, 43]. The real function oscillating with an amplitude about its mean :


has no strong limit when , although the functions and vary slowly. However for any smooth function with compact support on , one can obtain a vague limit as:


Thus the purpose of the observation function is to compute a smoothened version of the ”energy” of as given by . It allows to estimate the deviation of oscillations with amplitude about the mean at any point selected by the support of . This feature is illustrated on Fig. 1.

Figure 1. The energy limit of a strongly oscillating sequence: the oscillating function (thin line), the mean function (thick line), and the square-root limit (thick dashed line). After [43].

A mathematical generalization of this idea is possibly given by the notion of Wigner measure [40, 19, 36, 7, 1]. Let us now consider that the observable function is an arbitrary , compactly supported matrix function of both the space variable and the wave vector . For a vector field , the functional space of -valued, square integrable functions equipped with the scalar product , consider the (semiclassical) operator:


for . This parameter defines the so-called quantization of the operator. The case corresponds to the standard quantization. It is simply denoted by such that:


where stands for the Fourier transform of . The case corresponds to the Weyl quantization, which is usually denoted by . Then for a sequence uniformly bounded in , there exists a positive, Hermitian measure such that, up to extracting a subsequence if need be:


independently of the quantization . is the so-called Wigner measure of the sequence because it can also be interpreted as the weak limit of its Wigner transform . Indeed, if the latter is defined for -valued temperate distributions by:


then one has:

Thus describes the limit energy of the sequence in the phase space . As in Eq. (13), the matrix function is used to select any quadratic observable or quantity of interest associated to this energy: the kinetic energy, or the free energy, or the power flow, etc. For example, the high-frequency strain energy in may be estimated with :


up to some possible boundary effects on . Similarly, the kinetic energy is estimated by:


The vibrational energy density does not solve a closed-form equation in the high-frequency limit . However the Wigner measure, which provides a decomposition of these quantities in the phase space, does so as explained in the subsequent derivations. This is another reason why we shall now focus on such limit measure rather than directly or quadratic quantities of .

3. Wigner measure of high-frequency elastic waves

In this section we show how to obtain explicitly the Wigner measure (16) of the high-frequency solutions of the elastic wave equation (20) in the setting outlined in the foregoing section. Indeed, it is needed in (18) and (19) for the computation of the evolution of the strain and kinetic energies, for example. The objectives are also to outline its main properties for a slowly varying medium, as well as the (formal) mathematical tools used for its derivation. Both will prove useful in the subsequent Sect. 4 where elastic waves in a rapidly varying random medium with correlation lengths comparable to the small wavelength are considered. The analysis presented here is derived from [19], where first-order hyperbolic systems with constant and slowly varying coefficients are addressed, and [1], where arbitrary order hyperbolic systems with slowly varying coefficients are addressed. The dispersion properties of the elastic Wigner measure are derived in Sect. 3.4, and its evolution properties are derived in Sect. 3.5. Here it is shown that its components in the eigenspaces of the Christoffel tensor, the so-called specific intensities, satisfy transport equations of the Liouville type. The latter states that the energy densities in these eigenspaces are transported in phase space with celerities corresponding to the associated eigenvalues. Before doing so, we shall need some formal mathematical tools in order to compute the Wigner transform and its limit for high-frequency elastic waves. They are introduced in Sect. 3.2 and Sect. 3.3 below. Now in order to hopefully clarify the subsequent derivations, we start by reformulating the elastic wave equation in a form that is adapted to the analysis developed in the remaining of the paper.

3.1. Elastic wave equation as a semiclassical operator

Here we write Eq. (6) in a more convenient form for the derivation of the high-frequency regime . We shall first consider a slowly fluctuating medium characterized by an elastic tensor which is independent of the small parameter . The corresponding Christoffel tensor being given by Eq. (8) as where has been defined in Eq. (7), the elastic wave equation (6) then reads [1]:


with , , and the operator . The latter may be expanded as:




with the notation for two matrices and . It should be observed that is a first-order partial differential operator (independent of time), and that in a homogeneous medium.

3.2. Some formal rules of pseudo-differential calculus

Let be a smooth matrix-valued observable defined on , we recall the notation of Eq. (15) for the homogeneous semiclassical operator in the standard quantization. We then have (see e.g. [6]):


Here we assume that the differential operator within the observable acts on so that, for instance, should be interpreted as the component-wise inverse Fourier transform of the matrix . These relations can be extended for a non homogeneous semiclassical operator (as in Eq. (15) again) in the form [19, 6]:




since . Here stands for the usual Poisson bracket such that .

3.3. Spatio-temporal Wigner transform

In the sequel, the Wigner measure of the solutions of Eq. (20) shall be obtained using a spatio-temporal Wigner transform of that equation and its high-frequency limit as . This is because the spatial and temporal scales in the wave equation (20) play a symmetric role, and their oscillations should be accounted for altogether. Therefore a larger phase space than the one considered in the definition (17) has to be introduced. For two sequences and of square-integrable functions, the spatio-temporal Wigner transform is defined as:


Note that if it is applied to and , will be denoted by as implicitly done in Eq. (16). The spatio-temporal Wigner transform (26) is related to the (time-dependent) spatial Wigner transform (17) by , where arises as the dual variable of . Besides, one should observe that the rules expounded in Sect. 3.2 above can be extended further on to a time-dependent observable or a space-time observable . It suffices, for example, to redefine the Poisson bracket in Eq. (24) as a differential operator with respect to the space-time variable and the impulse variable : .

3.4. Dispersion properties

The pseudo-differential calculus and the spatio-temporal Wigner transform are now used for the wave equation (20). Computing the space-time Wigner transform (26) of and , yields:

But the partial derivative with respect to time reads:


where . Thus invoking the rules of calculus of the previous section, we get:


Considering the leading-order term, one obtains:


for the Wigner measure of the sequence given by Eq. (16). Owing to the properties (9) of the Christoffel tensor, one thus has:


where and:

Here is the number of different eigenvalues of the Christoffel tensor , and is the family of eigenvectors associated to the positive eigenvalue of which order of algebraic multiplicity is . It is assumed in the remaining that these multiplicities remain constant in phase space. We shall however see in Sect. 5, following [8, 39], that this is not always true for the classes of symmetry considered there and that the eigenvalues of the Christoffel tensor may coincide at some points or lines. The analysis of such crossings in terms of Wigner measures is a difficult task which has been addressed mathematically in [16, 17]. An explicit transition coefficient can be obtained in terms of the gap between the eigenvalues thanks to a proper rescaling of the crossing process. However, the extension of the results presented in these works to elasticity is not straightforward and requires further analyses out of the scope of the present publication. Going back to Eq. (30), the families of eigenvectors , , form an orthonormal basis of : , where is the identity matrix and stands for the Kronecker symbol. Then the ’s are projectors, implying that each term in the sum above cancels. Thus:

Likewise, being hermitian, on and consequently one has on for all . Taking the difference of these equalities yields:

But on as soon as , and in this case. Expanding on on the basis such that :

the expansion on finally reduces to the diagonal terms solely:


where . So the Wigner measure of the -oscillating elastic wave fields is expanded into the finite sum of its orthogonal projections onto the different energy paths of the propagation operator with symbol . These paths are determined by the equation , , in phase space. They correspond to the rays for the medium arising in classical Hamiltonian dynamics, as shown in the subsequent section.

3.5. Evolution properties

On the other hand, considering the Wigner transform of the adjoint wave equation and we have:

since the wave operator is actually formally self-adjoint, . This yields:


or, upon substracting Eq. (28) and Eq. (32), multiplying by , and observing that :


Here stands for the Lie bracket. One should observe that the matrix above is skew-symmetric. Following [40], we then introduce the matrices for such that with , and compute the projection of Eq. (33) on the eigen directions . Thus multiplying Eq. (33) by on the left side and by on the right side, one first obtains:

where . Indeed, the normalization condition yields and a similar relation for the partial derivatives