Coherent wave propagation in viscoelastic media with mode conversions and pair-correlated scatterers

Wave Motion 72 (2017) 244–259

Contents lists available at ScienceDirect

Wave Motion journal homepage: www.elsevier.com/locate/wavemoti

Coherent wave propagation in viscoelastic media with mode conversions and pair-correlated scatterers Francine Luppé a,∗ , Tony Valier-Brasier b , Jean-Marc Conoir c , Pascal Pareige a a

Laboratoire Ondes et Milieux Complexes, Normandie Univ, UNIHAVRE, CNRS, LOMC, 76600 Le Havre, France

b

Sorbonne Universités, UPMC Univ Paris 06, UMR 7190, Institut Jean Le Rond d’Alembert, F-75005, Paris, France

c

CNRS, UMR 7190, Institut Jean Le Rond d’Alembert, F-75005, Paris, France

highlights • • • • •

The effective wavenumber is expanded up to order 3 in concentration (c). The effect of correlation between scatterers is studied. It is of order 3 in concentration but of order 2 in scattering. It can modify greatly the wave velocity around a sub wavelength resonance. It increases the attenuation by a small constant amount at higher frequency.

article

info

Article history: Received 14 June 2016 Received in revised form 14 December 2016 Accepted 12 March 2017 Available online 22 March 2017 Keywords: Scattering Correlation Coherent wave Ultrasound

abstract The influence of correlation between scatterers on coherent waves propagation is studied in the case of a viscoelastic medium hosting a random configuration of either spherical or cylindrical scatterers. A distinction is made between the hole correction and the additional disturbances to the pair correlation function beyond the excluded volume via a radial and concentration dependent Ursell function. The effect of the Ursell function on the effective wavenumber is shown to be of order 3 in concentration and order 2 in scattering, and the corresponding formulas generalize those of Caleap et al. (2012) for an ideal fluid host medium. The whole order 3 in concentration is calculated; its other part is of order 3 in scattering. Both parts of the order 3 in concentration are the sum of two terms, one related to mode conversions, the other not. The numerical study is performed mostly for aluminum spheres in epoxy, which is a rather illustrative situation of the different phenomena that participate to the coherent propagation. The Ursell function effect is enhanced at low frequency, while counteracted partly at higher frequency, by the other term of order 3 in concentration. The most visible effects of both terms are on the attenuation. The Ursell term related to mode conversions is larger than the one with no mode conversions included in the low frequency regime. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Multiple scattering of waves in random media is of practical interest in many fields, ranging from ultrasonic monitoring of particulate mixtures [1] to the detection of flaws in industrial [2], biological [3,4], or other natural media [5]. Metamaterials

∗

Corresponding author. E-mail addresses: [email protected] (F. Luppé), [email protected] (T. Valier-Brasier), [email protected] (J.-M. Conoir), [email protected] (P. Pareige). http://dx.doi.org/10.1016/j.wavemoti.2017.03.002 0165-2125/© 2017 Elsevier B.V. All rights reserved.

F. Luppé et al. / Wave Motion 72 (2017) 244–259

245

are often built by the periodic arrangement of scatterers in a host matrix [6–8], but it is also possible to get interesting properties from random ones [9]. In such cases, the concentration of scatterers is rather large (approximately 20%) and it becomes important to have models well adapted to such concentrations. While the well-known effective-field theories [10,11] provide pretty good approximations of the wavenumber of the coherent wave in sparse concentrations of scatterers in an ideal host fluid, they fail at larger concentrations [12], and people usually either turn to other theories, such as the Coherent Potential Approximation [13], or introduce a pair distribution function [14,15] that takes into account the average microstructure of the random medium [16]. Self-consistent theories have been developed also for elastic media [17] (see references therein), but their main drawback is that they lead to a characteristic equation that is to be solved numerically to obtain the coherent wavenumbers. This is the reason why we rather consider here the effect of a pair distribution function on the formulas for the wavenumber obtained in the frame of an effective field-theory. Caleap’s study [15] was conducted for spheres in an ideal fluid, and the aim of this paper is to conduct a similar one for host media in which more than one type of wave may propagate: elastic, porous, thermo-viscoelastic media. . . . Section 2 is thus devoted to the description of the multiple scattering medium, consisting in a random distribution of either identical spheres (3D case) or identical infinitely long cylinders (2D case) in an isotropic host medium in which P different types of waves may propagate. The pair distribution function is the sum of a hole correction term and a radial concentration dependent Ursell function. The dispersion equation of the P coherent waves is established in Section 3, in which the 3D case only is detailed, the corresponding useful equations of the 2D case being given in Appendix A. The analytic expansion of the effective wavenumbers is performed in Section 4 up to power 3 in concentration, and it is shown in Section 5 that the term due to the Ursell function is one order less in scattering than in concentration. Finally, numerical simulations are performed in Section 6 in order to compare the relative importance, on the coherent wave properties, of the hole correction and the Ursell function. 2. Statistics of the medium We consider N = n0 V identical spheres of radius a randomly distributed in a given region of volume V with a concentration c equal to (4/3) n0 π a3 . The probability density of finding one sphere centered at r2 if one is known to be at r1 is p (r2 |r1 ) =

1 N −1

n (r2 |r1 ) .

(1)

The simplest choice for the conditional number density n (r2 |r1 ), n (r2 |r1 ) =



n0 0

for ∥r2 − r1 ∥ > b otherwise,

(2)

is known as the ‘‘Hole Correction’’. Usually, b is assumed to satisfy b > 2a so that the spheres are not allowed to overlap. However, it is known that the hole correction, whilst physically reasonable for sparse concentration, becomes a poor approximation with the increase of concentration. More generally we can use [16] n (r2 |r1 ) =



n0 0

(1 + U (r , n0 ))

for ∥r2 − r1 ∥ > b otherwise.

(3)

The U function, in addition to decaying rapidly to zero as r = ∥r2 − r1 ∥ tends to infinity, obeys lim U (r , n0 ) = 0.

n0 →0

(4)

The pair distribution function 1 + U (r , n0 ) of radially symmetric scatterers such as those considered here is also known as the ‘‘radial distribution function’’, and U (r , n0 ) is called the (two-particle) Ursell function [18,19]. There are many approximations of the radial distribution function for natural random distributions of solid spherical particles at equilibrium in an isotropic medium, depending on the number and type (thermo dynamical, or mechanical) of interaction forces between particles and surrounding medium. The most well known is the Percus–Yevick approximation, which is a solution of the exact Percus–Yevick equation for the description of a classical fluid, and which is shown in Ref. [15] to describe the average over a large number of realizations of artificial random distributions as well. The Ursell function is obviously negligible for sparse concentration, as shown with Eq. (4). Studying the effect of the pair distribution function term on the effective wavenumber expansion into powers of the concentration has been done in Ref. [15] for an ideal fluid host medium. The next sections use Refs. [15,20,21] in order to extend this study to host media in which mode conversions may occur during scattering. The different steps of the calculation are detailed for the 3D case only, but the final useful formulas are given in both cases in the Appendices.

246

F. Luppé et al. / Wave Motion 72 (2017) 244–259

3. Dispersion equation This section is written so as to avoid useless and cumbersome calculations. Starting from Ref. [22] we highlight the most important formulas and show how they can be modified or supplemented by the introduction of the Ursell function of Eqs. (3), (4). The notations followed are mostly the same as in Ref. [22]. The host medium allows the propagation of P different types of waves (P = 1 for ideal fluids, P = 2 for elastic solids, P = 3 for thermo-viscous fluids or Biot-porous solid,. . . ), and each is associated to a different coherent wave. All effective theories follow the same procedure. An incident plane wave is supposed to enter the volume hosting the scatterers, giving rise to a multiple scattering process. The multiple scattering equations express the incident fields on a given scatterer as the sum of an incident plane wave plus all the waves scattered by the other scatterers, as well as scattered waves in terms of incident waves on a given scatterer via the use of the transition matrix of that scatterer. The conditional ensemble average, over the random positions of the scatterers, of the incident field on one given scatterer at a given position, leads to a hierarchy of equations that is closed usually with the quasi-crystalline approximation. The average field incident on one given scatterer, the coherent field, is then expanded into spherical (or cylindrical, depending on the geometry of the scatterers) harmonics, in the same way as if it were a plane wave, propagating with a complex wavenumber, the effective wave number. Equating the coefficients of the exponentials related to the effective wavenumber ξ to zero provides the so(p) called Lorentz–Lorenz equation to which the partial amplitudes An in the expansion of a coherent wave of type p must obey. It is given in Eq. (18) in Ref. [22] as: A(np) −

P +∞ n+ν iεπ b   

yp

q=1 ν=0 ℓ=|n−ν|

(−1)ν+n (2ν + 1) Tν(qp) A(νq) Nℓ( ) (ξ b) G (0, ν|0, n|ℓ) = 0, p

(5)

with yp = ξ 2 − k2p

and ε = −4in0 .

(6)

In these equations, kp is the wave number of a p-type wave in the absence of scatterers, T (qp) is the diagonal scattering matrix of a wave of type q into one of type p, and the Gaunt coefficients G (0, ν|0, n|ℓ) are defined by Pnm (cos θ ) Pνµ (cos θ) = (p)

The Nℓ

 ℓ

m+µ

G (m, n|µ, ν|ℓ) Pℓ

(cos θ) .

(7)

function, defined as

(p)

(1)

Nℓ (ξ ) = ξ bj′ℓ (ξ b) hℓ

(1)′

kp b − kp bjℓ (ξ b) hℓ







kp b ,



(8)

is a direct consequence of imposing the hole correction Eq. (2). The generalization to an arbitrary pair distribution function as in Eq. (3) is done in Ref. [15]; it leads to the introduction of an additional term to Eq. (5), which becomes (p)

An −

P +∞ n+ν iεπ b   

yp

 (−1)

ν+n

q=1 ν=0 ℓ=|n−ν|

(2ν +

1) Tν(qp) A(νq)

(p) Nℓ



1 yp (p) L (ξ ) G (0, ν|0, n|ℓ) = 0, (ξ ) + kp b k2p ℓ

(9)

with L(mp) (ξ ) =

+∞



jm (ξ r ) h(m1) kp r U (r , n0 ) k3p r 2 dr .





(10)

b

The procedure to recover the effective wavenumbers ξ from Eq. (10) is exactly the same as in Ref. [22], which followed Ref. [21]. It consists in writing the system of equations obtained for all values of p and n in Eq. (9) in a matrix form with block matrices whose coefficients are matrices related to wave types (p, q), while the coefficients of these latter matrices are related to orders (n, ν ) of the afore-mentioned expansions in spherical harmonics. The trick is then to obtain the dispersion equation by setting to zero the determinant of a (P , P) block matrix, Det [Y (ξ ) − ε M (ξ )] = 0,

(11)

where Y (ξ ) and M (ξ ) are matrices of rank P, Y (ξ ) is diagonal with Ypp (ξ ) = yp . The Mqp (ξ ) elements of M (ξ ) are given by (misprints in Eqs. (27, 28) of Ref. [22])

    +∞   n   Mqp (ξ ) = eq |T | ep + eq  T Q T  ep ϵ n ,   n =1 

(0) Mqp (ξ ) = Mqp +



+∞  n =1

(n) n Mqp ε .

(12a)

(12b)

F. Luppé et al. / Wave Motion 72 (2017) 244–259

247

   ep are block vectors whose components other than the pth, which is equal to π /kp |e⟩, are zero vectors, T is a block matrix   of coefficients (T )pq = T (qp) , and the nth component of |e⟩ is (−1)n . The Q block matrix is diagonal, with Q pp = Q (p) , and, as we start here with Eq. (9), in place of Eq. (5), Eq. (19) in Ref. [22] is to be replaced by the following: (p) Q nν



π

(ξ ) =

−1 + ikp b

kp yp



n+ν 

ℓ=|n−ν|

(p) Nℓ





1 yp (p) L (ξ ) G (0, ν|0, n|ℓ) (−1)ν+n . (ξ ) + kp b k2p ℓ

(13)

As expected, we recover the characteristic equation of Ref. [15] when P = 1. We look now for the expansion of the solutions of the dispersion equation into powers of the concentration of scatterers. 4. Wavenumbers of the coherent waves For a given type of wave p, we assume the corresponding effective wavenumber to obey the following formal asymptotic expansion [20,21] up to order No in concentration,

ξp2 = k2p +

N0 

ε n y(pn) ,

(14)

n=1

(n)

and we look for the yp s by inserting Eq. (14) into the dispersion equation, Eq. (11). In order to do so, we need to expand (p)

as well the Ursell and Qnν functions that depend on the concentration, but, because of Eqs. (4), (12)–(13), only up to order No − 2:

  j  ∂ ∂j U r , n = n U r , n + · · · = n0 u ( r ) + · · · , (15) ( ) ( ) 0 0 0 ∂ nr0 ∂ nr0 n0 =0 n0 =0 j =1    N j (q) o −2  Q ξ d (ε) p n ν (q)   (q)    Qnν ξp = Qnν kp + εj  dε j j =1 ε=0 (q)   (q,1)   = Qnν kp + ε qnν kp + · · · . (16)   (q) (q) Qnν kp is the zero order term in ε = −4in0 of Qnν (ξ ), and, as such, it is found from replacing n0 with 0 in Eq. (13), so that N o −2

U (r , n0 ) =



j



n0

its expression remains the same as in Ref. [22], (q)

 

(p)

 

Qnm kp

  n +m    π n +m (q)  (−1) Nν kp G (0, m|0, n|ν) , 1 − ikq b =−  2 kq kp − k2q ν=|n−m| π

n +m



(−1)n+m

(q ̸= p) ,

(17a)

D(νp,0) kp G (0, m|0, n|ν) ,

 

(17b)

    ν (ν + 1) − k b2 j k b h(1) k b  p ν p q ν = ikq b .         −k bk bj′ k b h(1)′ k b − k bj′ k b h(1) k b  p q ν p q p ν p q ν ν

(18)

Qnm kp =

2k3p

ν=|n−m|

with

  (q,0)

Dν

kp

(q,j)

In the following, we extend the results of Ref. [22] up to order No = 3. We thus need the qnν (q,1)

qnm

  kp

=

π 2

 

kp terms for j = 1:

n+m



(−1)n+m

G (0, m|0, n|ν)

ν=|n−m|  (0) 2



kp − k2q

Mpp

   2 2 2  ×  kq kp − kq  1 (q)   − 3 Λν kp

k2p



  (q,0)



  (q)

kp + 2 1 − ikq bNν

Dν

kp

    (q ̸= p) , 

(19)

2kq

(p,1)

qnm

 

kp = −

π 4



n+m

(−1)

n +m

 ν=|n−m|

G (0, m|0, n|ν)

(0)     Mpp  (p,0)   1 2Dν kp + dν kp + 3 Λ(νp) kp , 5 2kp kp



(20)

248

F. Luppé et al. / Wave Motion 72 (2017) 244–259

with

      2 + ν (ν + 1) jν kp b h(ν1) kp b + kp b  3      2     − ν (ν + 1) − i kp b j′ν kp b hν(1)′ kp b − i kp b j′ν kp b h(ν1) kp b ,

 

= ikp b

dν kp

π

(0) Mpp =

+∞ 

kp n=0



kp b

2

(21)

(2n + 1) Tn(pp) ,

(22)

and

Λ(νq) (ξ ) =

+∞



jν (ξ r ) h(ν1) kq r u (r ) k3q r 2 dr .



b



(23)

(q,1)

  kp terms, which are related to order 3 in concentration, contain not only the effect of the   (q)

It is worth noticing that the qnν

Ursell function through the Λν kp integrals, as usual [14–18], but that also of the hole correction through the functions evaluated at kp b in Eqs. (19), (20). This point had never been clearly discussed in the literature as far as we know. Collecting together Eqs. (11), (12), (14)–(23) provides the same orders 1 and 2 of the asymptotic expansion in Eq. (14) as in the case of the hole correction, i.e. (0) y(p1) = Mpp ,

(24a) (0) (0) P  Mpq Mqp , k2p − k2q q=1

(1) y(p2) = Mpp +

(24b)

q ̸= p

with

π

(0) Mqp =

+∞ 

1/2 1/2 kp kq n=0

π

(1) Mqp =

(2n + 1) Tn(pq) ,

P  

1/2 1/2 kq kp r =1

n

(25a)

(r )   (−1)n+m (2n + 1) (2m + 1) Tn(rq) Tm(pr ) Q nm kp .

(25b)

m

The third order term is found as (2) y(3) = Mpp +

(0) (0) (0) (0) (0) (0) P  Mqq Mpq Mqp − Mpp Mpq Mqp  2 k2p − k2q q=1 q ̸= p

P

+

(0)

(0)

P 



(0)

Mpq Mqr Mrp



k2p − k2q

q=1 r =1 q ̸= p r ̸= p, r ̸= q



+ 2

k2p − kr

(1) (0) (0) (1) P  Mpq Mqp + Mpq Mqp , k2p − k2q q=1

(26)

q ̸= p

with (2) = Mpp

P  π 

kp q=1

n

(−1)n+m (2n + 1) (2m + 1)

m

 (qp) (pq) (q,1)

× Tn

Tm qnm

 

kp +

P   r =1

(2ℓ +

ℓ

(qr ) (r ) 1) Tℓ Tn(rp) Tm(pq) Qnℓ

(q) kp Qℓm

 

   kp

.

(27)

All second and third order terms may be collected together and Eq. (14) can be written in the more comprehensive and practical form

ξp2 k2p









P P P 3    ε2   ε  (pq)  = 1 + 3 δ1(pp) + 6 δ2(pp) + δ2(pq)  + 9 δ3(pp) + δ3(pq) + δ3Ursell  + ···,

ε

kp

kp

q=1 q ̸= p

kp

q=1 q ̸= p

q =1

with the δ terms given in Appendix B, and the order 3 terms being the main result of this paper.

(28)

F. Luppé et al. / Wave Motion 72 (2017) 244–259

249

5. On the correlation Many authors consider that the correlation as a whole can be regarded as second order in concentration [14,15,23,24]. The asymptotic expansion, Eq. (28), shows that while the hole correction effect is indeed, the additional Ursell function term in the pair distribution function is third order only in concentration. In order to understand better the difference of order between those two correlation-related terms, let us consider the correlation term given for spherical scatterers in a fluid (P = 1, kp = k) in Refs. [14,24] through the following integral, Iν (ξ ) = ε 2

+∞

 b

[g (r ) − 1] r 2 h(ν1) (kr ) jν (ξ r ) dr ,

(29)

with the pair correlation (or distribution) function g defined, with our notations, as (cf. Eq. (15)) g (r ) = 1 + U (r , n0 ) = 1 + n0 u (r ) + · · · .

(30)

The key point here is that the pair distribution function g does not depend on r only, but on the concentration as well. Moreover, it tends to unity when the concentration vanishes, as discussed in Section 2 of this paper as well as in the first paper of Ref. [11]. Considering for example the Percus–Yevick’s pair correlation function, its Virial expansion, up to third order in concentration, is given in Ref. [15] as g (r ) =

 1

for r ≥ 2b

 1 + n0

4π 3

b3

 1−

3r 4b

+

1 r3 16



b3

(31)

for b ≤ r ≤ 2b,

which is clearly the same form as the right hand side of Eq. (30), and, used in Eq. (29), provides a third order term in concentration: Iν(q)

(ξ ) = ε n0 2

4π 3

b

3

2b



 1−

b

3r 4b

+

1 r3 16 b3



r 2 h(ν1) kq r jν (ξ r ) dr .





(32)

It follows also that, while correctly taken into account in Ref. [15], the effect of the Ursell function should not be a priori separated from the (other part of) order 3 in concentration that the authors had found too long to warrant including in their (pp) (pq) paper, and we focus now on the third order in concentration. The δ3 and δ3 coefficients include the hole correction, and (pq)

each of them is a sum of terms that involve the product of three scattering coefficients. The last coefficient, δ3Ursell , which vanishes with the Ursell function, involves the product of two scattering coefficients only. In other words, while all other terms of a given order in concentration are of the same order in scattering, those due to the Ursell function are of one order (pp) (pq) (pq) less in scattering than in concentration. It might thus be that the δ3 and δ3 terms are negligible, compared to δ3Ursell . The relative importance of all correlation terms, in both the second and third orders in concentration, is studied numerically in the next section. 6. Numerical study We focus here on the coherent wave whose wavenumber is closest to that of a longitudinal wave in the absence of scatterers (p = 1). The aim of this section is to investigate numerically the influence of the different terms of the expansion of its wavenumber on both its velocity and its attenuation. Silica nano-spheres in (viscous) water have been studied in Ref. [25] in the sub-wavelength domain for the compression wave and concentrations in the range 5%–20%; it has shown that agreement between theory and experiment could not be achieved if one neglected the compression/viscous—shear mode conversions, and that, for the highest values of the concentration, third order terms in concentration were needed as well. Ref. [26] was about the sub-wavelength multiple scattering by Tungsten-Carbide spheres in an epoxy resin, with concentrations ranging from 2% to 10%; again, the influence of the mode conversions was clearly demonstrated from the comparison with experiment, along with the need of taking into account the hole correction (as opposed to the low frequency k1 b → 0 assumption). In both references, however, discrepancy between theory and experiment appeared as the concentration increased: at concentrations larger than 25% in Ref. [25], in which third order terms in concentration were included but the Ursell function was neglected, and at concentrations larger than 5% in Ref. [26], in which there were no third order terms. The Tungsten spheres in Ref. [26] exhibited a strong sub-wavelength resonant behavior that was not present in the case of Ref. [25], and we believe it might be the reason why the model failed at relatively low concentrations in Ref. [26]. In the following, we avoid this very peculiar case where the mass density ratio between the scatterers and the matrix is very large, which may lead to negative attenuation, as discussed later in this section. We consider here at first different concentrations of aluminum spheres in an epoxy matrix. The aluminum beads are then replaced by much softer spheres, made of Polydimethylsiloxane (PDMS), whose mass density, contrary to that of aluminum, is smaller than that of epoxy. Both Epoxy and PDMS are viscoelastic isotropic media, modeled

250

F. Luppé et al. / Wave Motion 72 (2017) 244–259 Table 1 Mass density ρ , first (λ0 ) and second (µ0 ) elastic Lamé coefficients of the different media in the numerical study.

Aluminum Epoxy PDMS

ρ (kg/m3 )

λ0 (GPa)

µ0 (GPa)

2700 1100 965

45 3.6 0.93

25 1.05 2.4 10−3

Table 2 Relaxation times τn used in the calculation of the elastic Lamé coefficients λ and µ of the epoxy in the numerical study.

λ µ

τ1 (µs)

τ2 (µs)

τ3 (ns)

τ4 (ns)

5 5

0.5 0.5

50 50

8 14

with frequency dependent complex Lamé coefficients. There are thus P = 2 different types of waves in the epoxy host medium: longitudinal (p = 1) and shear (p = 2) waves. The Zener model used for epoxy [27], C = C0 +

4 

Cn

n =1

(ωτn )2 − iωτn , 1 + (ωτn )2

(33)

with C standing for either λ or µ, is derived from the experiment and described in Ref. [26], while the Kelvin–Voigt model is used for the PDMS,



2



λ = λ0 − iω ηB − ηS , 3

µ = µ0 − iωηS ,

(34)

with ηB = 1.3 Pa s and ηS = 0.6 Pa s. The elastic Lamé coefficients λ0 and µ0 , along with the mass density, are given in Table 1 for (elastic) aluminum, epoxy and PDMS, while the relaxation times in Eq. (33) are given in Table 2. These values lead to a practically real wavenumber k1 of the longitudinal wave in epoxy on the frequency range we studied. Most figures in the two following sections will show the dispersion curves (velocity and attenuation) versus frequency obtained from different approximations of the effective wavenumber. The (green) curves with no marks correspond to k1 b → 0, i.e. no hole correction is taken into account, contrary to the (black) curves with triangles for which the exclusion radius is equal to the diameter of the scatterers. Both sets of curves correspond to the asymptotic expansion of the wavenumber limited to order two in concentration and scattering. The additional correlation term due to the Ursell function is added to provide the (blue) curves with stars. The pair distribution function considered is that of Percus–Yevick; its Virial expansion, Eq. (31), is used with the exclusion radius b equal to twice the radius a of the cylinders. Finally, the (magenta) curves with squares have been drawn with the whole order three in concentration taken into account. 6.1. Aluminum spheres in epoxy We consider here an epoxy matrix with aluminum spheres within. As can be seen from Table 1, aluminum is ten to twenty times stiffer than epoxy; its resonant behavior, thus is governed mainly by the ratio of its density to that of epoxy. The spheres have a 0.5 mm radius a. Figs. 1, 2 show respectively the speed and the attenuation of the longitudinal coherent wave for c = 15% on the reduced frequency range 0.01 ≤ k1 a ≤ 10, and Figs. 3, 4 the same at low frequency, 0.001 ≤ k1 a ≤ 1. The influence of the concentration is observed through the comparison of the latest with Figs. 5, 6, which are their equivalent at higher concentration, c = 25%. Same way as in Ref. [26], Figs. 1–4 show clearly a peculiar behavior of the dispersion curves obtained at order 2 when the hole correction is not taken into account (green curves, no marks), when compared to the other curves. Numerically, this is due to the fact that k2 b is five times greater than k1 a, so that the Hankel functions of argument k2 b cannot be replaced by their asymptotic expression at low argument as was done in the low frequency section of Ref. [22]: Eq. (36) of that reference, in fact, is rather a quasi-static approximation than a low frequency one. Physically, this underlines the importance of taking into account the size of the scatterers, even at low frequency. However, it may happen that forgetting it provides nonetheless a reasonable agreement with the experiment at larger than unity k1 a values and non overlapping scatterers [23]. The Lloyd–Berry effective wavenumber [10,11] corresponds to the dashed green curve with no marks in Figs. 3, 4. It is recovered from the order 2 with no hole correction by neglecting the mode conversions. The difference between those two curves with no marks clearly shows the importance of mode conversions at low frequency. We focus now on the other curves of Figs. 1–4, obtained with the hole correction accounted for and b = 2a. All three curves with marks in Fig. 1 provide practically the same velocity for high frequencies f ≥ 3 MHz (k1 a around 4 and larger). By contrast, Fig. 2 shows that the differences between the attenuation curves are negligible only at rather low frequencies,

F. Luppé et al. / Wave Motion 72 (2017) 244–259

251

Fig. 1. (Color online) Aluminum spheres in epoxy, a = 500 µm, c = 15%. Effective phase velocity ω/ [ℜe (ξ1 )] versus frequency, 0.01 ≤ k1 a ≤ 10, 0.025 ≤ k2 a ≤ 25. Green line, no marks: from the expansion of ξ1 up to order 2 in concentration, no hole correction (b → 0). Black line + triangles: from the expansion of ξ1 up to order 2, hole correction included (b = 2a). Blue line + stars: third order Ursell term added. Magenta line and squares: from the expansion of ξ1 up to whole order 3 in concentration.

Fig. 2. (Color online) Aluminum spheres in epoxy, a = 500 µm, c = 15%. Effective attenuation ℑm (ξ1 ) versus frequency, 0.01 ≤ k1 a ≤ 10, 0.025 ≤ k2 a ≤ 25. Green line, no marks: from the expansion of ξ1 up to order 2 in concentration, no hole correction (b → 0). Black line + triangles: from the expansion of ξ1 up to order 2, hole correction included (b = 2a). Blue line + stars: third order Ursell term added. Magenta line and squares: from the expansion of ξ1 up to whole order 3 in concentration.

f ≤ 2 MHz (k1 a around 2.8 and smaller). Other calculations, not shown here, for titanium beads have led us to think that such a difference between velocity and attenuation behavior is quite a general result for scatterers that are heavier than the matrix. In addition, one can check from the comparison of Figs. 5, 6, drawn for c = 25%, to Figs. 3, 4, drawn for c = 15%, that the frequency range in which there are noticeable differences between the different approximations increases with concentration. Let us concentrate now on the influence of the Ursell function, which can be asserted from the comparison between the (black) curves with triangles and the (blue) ones with stars. The Ursell function has no influence on the velocity at low and high frequency, f ≤ 0.05 MHz and f ≥ 3 MHz (k1 a around 4 and larger), while the situation is quite different for the attenuation, which is sensitive to this additional correlation term on the whole frequency range 2 ≤ f ≤ 7 MHz (k1 a roughly between 2.8 and 10). Again, it seems to be a relatively general result that the influence of the Ursell function on the attenuation takes place on a larger frequency band than it does on the velocity: we observed it for the other bead materials tested, as well as the authors of Ref. [23] for steel rods in water and of Ref. [24] for bubbles in PDMS. The sub-wavelength resonance in Ref. [26] was shown to be due to the dipole mode n = 1. It is also the case here for the resonance at 340 kHz (from Eq. (3) in Ref. [26]) that is more clearly seen in Fig. 5, because of the higher concentration, than in Figs. 1–4. In order to investigate its effect, we have plotted in Fig. 7 the velocity obtained from the wavenumber expansion up to order 2 in concentration, with no hole correction taken into account (b → 0), for different values of the concentration. This figure clearly shows a curvature inversion around the resonance frequency (marked by the vertical dashed line) that is

252

F. Luppé et al. / Wave Motion 72 (2017) 244–259

Fig. 3. (Color online) Aluminum spheres in epoxy, a = 500 µm, c = 15%. Effective phase velocity ω/ [ℜe (ξ1 )] versus frequency. Zoom of Fig. 1 at low frequency, 0.001 ≤ k1 a ≤ 1, 0.0025 ≤ k2 a ≤ 2.5. Solid green line, no marks: from the expansion of ξ1 up to order 2 in concentration, no hole correction (b → 0). Dashed green line, no marks: order 2, no hole correction (b → 0) and no mode conversions. Black line + triangles: from the expansion of ξ1 up to order 2, hole correction included (b = 2a). Blue line + stars: third order Ursell term added. Magenta line and squares: from the expansion of ξ1 up to whole order 3 in concentration.

Fig. 4. (Color online) Aluminum spheres in epoxy, a = 500 µm, c = 15%. Effective attenuation ℑm (ξ1 ) versus frequency. Zoom of Fig. 2 at low frequency, 0.001 ≤ k1 a ≤ 1, 0.0025 ≤ k2 a ≤ 2.5. Solid green line, no marks: from the expansion of ξ1 up to order 2 in concentration, no hole correction (b → 0). Dashed green line, no marks: order 2, no hole correction (b → 0) and no mode conversions. Black line + triangles: from the expansion of ξ1 up to order 2, hole correction included (b = 2a). Blue line + stars: third order Ursell term added. Magenta line and squares: from the expansion of ξ1 up to whole order 3 in concentration.

enhanced by the increase of the concentration. Turning then back again to Fig. 5, corresponding to c = 25%, one can see that the main effect of the Ursell function is a shift towards high frequency of the whole velocity curve: for example, the [0.2 MHz, 0.4 MHz] frequency interval over which the phase velocity decreases is shifted to [0.3 MHz, 0.5 MHz] by the introduction of the Ursell term. A similar shift, due in that case to the hole correction only, can be observed from the comparison of Figs. 5–7. The hole correction is usually seen only as a mean to proscribe inter penetration of scatterers. In the light of the similarity of its effect with that of the Ursell function, it can also be seen as the low concentration part of the correlation term, as described indeed by Eq. (31). The influence of mode conversions on the terms related to the Ursell function is depicted by Fig. 8, which clearly shows that the amplitude of the Ursell term related to mode conversions can be larger than that with no mode conversion. The same holds for both the real and imaginary parts of those terms: mode conversions must be taken into account at low frequency (k1 a < 1). As for the remaining part of order three in concentration, it seems to enhance the Ursell function effect at low frequency (the blue curves with stars are in between the black ones with triangles and the magenta ones with squares in Figs. 3, 4) and to counteract it partly at higher frequency (the black curves with triangles are between the blue one with stars and the magenta one with squares in Figs. 1, 2). This concurs with the need of its introduction in the low frequency study of Ref. [25].

F. Luppé et al. / Wave Motion 72 (2017) 244–259

253

Fig. 5. (Color online) Aluminum spheres in epoxy, a = 500 µm, c = 25%. Effective phase velocity ω/ [ℜe (ξ1 )] versus frequency, 0.001 ≤ k1 a ≤ 1, 0.0025 ≤ k2 a ≤ 2.5. Black line + triangles: from the expansion of ξ1 up to order 2 in concentration, hole correction included (b = 2a). Blue line + stars: third order Ursell term added. Magenta line and squares: from the expansion of ξ1 up to whole order 3 in concentration.

Fig. 6. (Color online) Aluminum spheres in epoxy, a = 500 µm, c = 25%. Effective attenuation ℑm (ξ1 ) versus frequency, 0.001 ≤ k1 a ≤ 1, 0.0025 ≤ k2 a ≤ 2.5. Black line + triangles: from the expansion of ξ1 up to order 2 in concentration, hole correction included (b = 2a). Blue line + stars: third order Ursell term added. Magenta line and squares: from the expansion of ξ1 up to whole order 3 in concentration.

A last remark should be made before concluding this section: while hardly noticeable in Fig. 6, the attenuation at order 2 in concentration (black curve with triangles) is negative around 0.28 MHz, while it is always positive at the lower concentration of Fig. 4. Such a phenomenon had already been noticed by others [15,28,29]. Introduction of the correlation gives back a positive attenuation in Fig. 6, as it was shown to do in Ref. [15]. However, other studies we have conducted (not shown here) have shown that the correlation term merely increases the critical concentration at which a negative attenuation may occur. This critical concentration, in fact, appears mostly in the vicinity of a sub-wavelength resonance, and depends on the scatterers mass density: the higher the latter, the smaller the critical concentration [26]. Now, the question is: is a negative attenuation merely a manifestation of the limit of the multiple scattering model itself, or a truly physical metamaterial like property, as is, for example, a negative index [9]? Answering this question is beyond our ability and aim at the present time. 6.2. PDMS spheres in epoxy PDMS is both lighter and much softer than epoxy (see Table 1), and, in contrast with aluminum, its resonant behavior involves more than one parameter (the mass density ratio with the matrix). Its high compressibility gives rise to a large number of Mie resonances [30] on the whole frequency range that cause the modulations observed in Figs. 9, 10. The subwavelength resonance is not as strong as that of the aluminum spheres, hence the phase velocity variations (f < 0.5 MHz, k1 a smaller than 1.4) in Fig. 9 are smaller than those in Fig. 5 for the same k1 a reduced frequency range. A much stronger dispersion of the phase velocity is observed in Fig. 9 than in Fig. 7 for the aluminum spheres and the same concentration

254

F. Luppé et al. / Wave Motion 72 (2017) 244–259

Fig. 7. (Color online) Aluminum spheres in epoxy, a = 500 µm. Effective phase velocity ω/ [ℜe (ξ1 )] versus frequency, 0.001 ≤ k1 a ≤ 1, 0.0025 ≤ k2 a ≤ 2.5, for different values of the concentration. From the expansion of ξ1 up to order 2 in concentration, no hole correction (b → 0). The vertical dashed line is situated at the resonance frequency of the spheres, f = 340 kHz.

Fig. 8. (Color online) Aluminum spheres in epoxy, a = 500 µm. Comparison of the magnitude of the Ursell related  parts of the third order term in

  3  ε (11)  δ , blue dashed line: 3 3Ursell  k1 c 

concentration in the expansion of ξ12 /k21 versus frequency 0.01 ≤ k1 a ≤ 10, 0.025 ≤ k2 a ≤ 25. Solid red line: 

   3    ε (12)   k3 c δ3Ursell , from Eqs. (B.1), (B.3.6). This figure does not depend on the concentration, due to the division par c 3 .  1 

(25%). We have, however tried to identify neither the vibration mode associated to that behavior, nor the resonances observed at higher frequency or the origin of the small variations superposed to those resonances in the attenuation curves of Fig. 10, as it is out of the scope of this paper. The main thing to be retained, however, is that the previous conclusions made for heavy and ‘‘rigid’’ particles are still valid here: the attenuation is more sensitive to correlation than the velocity, with the attenuation decrease due to the Ursell function counteracted by the remaining terms of order 3 in concentration, so that the attenuation, when all terms at order 3 in concentration are taken into account, is larger than the expansion limited to order 2 in concentration provides. 7. Conclusion We have studied the properties of coherent waves in elastic media through the expansion of the effective wavenumber up to the third power of the concentration. Order three contains two different contributions. The first one is a consequence of correlation between scatterers that we modeled by the introduction of a pair distribution function. While of order three

F. Luppé et al. / Wave Motion 72 (2017) 244–259

255

Fig. 9. (Color online) PDMS spheres in epoxy, a = 1 mm, c = 25%. Effective phase velocity ω/ [ℜe (ξ1 )] versus frequency, 0.01 ≤ k1 a ≤ 10, 0.025 ≤ k2 a ≤ 25. Green line, no marks: from the expansion of ξ1 up to order 2 in concentration, no hole correction (b → 0). Black line + triangles: from the expansion of ξ1 up to order 2, hole correction included (b = 2a). Blue line + stars: third order Ursell term added. Magenta line and squares: from the expansion of ξ1 up to whole order 3 in concentration.

Fig. 10. (Color online) PDMS spheres in epoxy, a = 1 mm, c = 25%. Effective attenuation ℑm (ξ1 ) versus frequency, 0.01 ≤ k1 a ≤ 10, 0.025 ≤ k2 a ≤ 25. Green line, no marks: from the expansion of ξ1 up to order 2 in concentration, no hole correction (b → 0). Black line + triangles: from the expansion of ξ1 up to order 2, hole correction included (b = 2a). Blue line + stars: third order Ursell term added. Magenta line and squares: from the expansion of ξ1 up to whole order 3 in concentration.

in concentration, it is of order two only in scattering. The second one is, as all other terms of the expansion, of the same order in scattering as in concentration. When compared to the influence of the correlation on the properties of the coherent wave, that of the other part of the third order is no different on a qualitative point of view. It needs, however, to be taken into account for quantitative purposes. The effect of correlation depends strongly on frequency, through the ratio of the longitudinal wavelength to the radius of the scatterers. At low frequency, it may modify the phase velocity rather greatly, as the frequency range over which the sub-wavelength dipolar resonance induces a high velocity dispersion is shifted towards higher frequencies. As a general rule, while the correlation has a rather small effect on the phase velocity, its influence on the attenuation is not negligible, except at very low frequency (k1 a < 1). Adding only the term related to the Ursell function to the expansion up to order two in concentration decreases the attenuation, but adding the whole order three in concentration increases it by a small and practically constant amount. The authors wish to thank one of the reviewers who has helped them clarifying the original manuscript by pointing out that the general term of correlation included both the hole correction and the Ursell function.

256

F. Luppé et al. / Wave Motion 72 (2017) 244–259

Appendix A. The 2D case The expressions needed for the 2D case are given here +∞ 

(0) Mqp =

Tn(pq)

(A.1.1)

n=−∞

+∞ +∞  P  

(1) Mqp =

(r )

Tn(rq) Tm(pr ) Q nm kp

 

(A.1.2)

n=−∞ m=−∞ r =1

+∞ +∞  P  

(2) Mpp =

(q,1)

Tn(qp) Tm(pq) qnm

n=−∞ m=−∞ q=1

+∞ +∞ +∞  P  P   

+

n=−∞ m=−∞ ℓ=−∞ q=1 r =1

(q)

(q)

i π2 Nn−m kp − 1

 

Qnm kp = (p)

1

 

(p)

k2p

Tn

(r )

Tm Qnℓ

 

kp Qℓm kp

(A.1.3)

(q ̸= p)

(A.1.4)

 

Dn−m kp

(A.1.5)

(q) Nm (ξ ) = ξ bJm′ (ξ b) Hm(1) kq b − kq bJm (ξ b) Hm(1)′ kq b



(q)

(q)

 

 

k2p − k2q

Qnm kp =

(qr ) (rp) (pq)

Tℓ







(A.1.6)

 

Dm

kp

= (q,1)

qnm

2    (1)     (1)′   iπ  2  ′ m − kp b Jm kp b Hm kq b − kp bkq bJm kp b Hm kq b 4 (0)

  kp

Mpp

=

k2p

− (p,1)

qnm

π

(0)

kp =

Λ(mq) (ξ ) =

k4p +∞



k2p − k2q

(q) Λn−m 2

2

  kp

8kq

Mpp

 



1



(q,0)

k2p − k2q Dn−m kp − k2p





 



iπ 2

(q)

(A.1.7)



 

Nn−m kp − 1

(q ̸= p)

(A.1.8)

      1 (p,0)   π − Dn−m kp + dn−m kp − 2 Λ(np−) m kp , 2

(A.1.9)

8kp

(1) Jm (ξ r ) Hm kq r u (r ) k2q rdr ,





(A.1.10)

b

 

dm kp =

1 8

m2 − kp b



2

 2     − iπ kp b Jm kp b Hm(1) kp b .

(A.1.11)

Appendix B. Final expressions The most practical form of the wavenumber expansion, up to order 3 in concentration, is:

ξ

 2 p

k2p

≃ 1+

ε kdp

δ1(pp) +







 (pq)  (pq)  ε  (pp)  (pq)  ε  δ2  + 3d δ3(pp) + δ3 + δ3Ursell  , δ2 + 2d 2

kp

P

q=1 q ̸= p

3

kp

P

P

q=1 q ̸= p

q =1

(B.1)

with d = 2 in the 2D case and d = 3 in the 3D one, and the different delta terms given hereafter for both cases. The 2D case

δ1(pp) =

+∞  n=−∞

Tn(pp)

(B.2.1)

F. Luppé et al. / Wave Motion 72 (2017) 244–259

+∞ +∞  

δ2(pp) =

(p,0)

257

Tn(pp) Tm(pp) Dn−m kp

 

(B.2.2)

n=−∞ m=−∞

δ2(pq) = i

π 2

δ3(pp) = −

k2p k2p

+∞ +∞  

k2q n=−∞ m=−∞

−

+∞ 1 

+∞ +∞  

2 n=−∞ m=−∞ ℓ=−∞ +∞ +∞ +∞   

+

n=−∞ m=−∞ ℓ=−∞

δ3(pq)

(q)

Tn(qp) Tm(pq) Nm−n kp

 

(pp) (p,0) Dn−m

Tn(pp) Tm(pp) Tℓ Tn

 

kp +

+∞ +∞ +∞    n=−∞ m=−∞ ℓ=−∞

(p,0)  

(pp) (pp) (pp)

Tℓ

(B.2.3)

(p,0)

(pp)

Tn(pp) Tm(pp) Tℓ

 

dn−m kp

 

Tm Dn−ℓ kp Dℓ−m kp

(B.2.4)

  +∞ +∞ +∞    k2p (pp) (q,0)    Tn(qp) Tm(pq) Tℓ Dn−m kp + 2   k − k2 p q n=−∞ m=−∞ ℓ=−∞     4 +∞ +∞ +∞      kp     iπ (q) (qp) (pq) (pp)  −  T T T N k p  n − m n m ℓ   2 k2 − k2 2 n=−∞ m=−∞ ℓ=−∞   p q   2 +∞ +∞ +∞      kp π  = qp) (pp) (pq) pq) (pp) (qp) ( (  +i   Tn Tm Tℓ + Tn Tm Tℓ   2 k2 − k2 p q n=−∞ m=−∞ ℓ=−∞         (p,0) (q)  Dℓ−m kp Nℓ−n kp   P +∞ +∞ +∞      (q)    π 2 k4p  1 qr ) (rp) (pq) (r ) ( − T T T N k N k p p  n m ℓ−n m−ℓ 4 k2p − k2q r = 1 k2p − k2r n=−∞ m=−∞ ℓ=−∞ ℓ

(B.2.5)

r ̸= p

(pq)

δ3Ursell = −

π2 

kp a

8

+∞ +∞  2 

(qp) (pq)

Tn

n=−∞ m=−∞



+∞

Tm

(1)

Jn−m kp r Hn−m kq r u0 (r ) k2p rdr ,









(B.2.6)

b

with u0 the non dimensional function defined from function u, u0 (r ) =



∂U ∂c

 = n 0 =0

1

π a2

u (r ) .

(B.2.7)

The 3D case

δ1(pp) = π

+∞ 

(2n + 1) Tn(pp)

(B.3.1)

n =0

δ2(pp) =

+∞  +∞ π2 

2

k2p k2p

n+m

 ν=|n−m|

n=0 m=0

  δ2(pq) = iπ 2 kp b 

δ3(pp)

(2n + 1) (2m + 1) Tn(pp) Tm(pp)

−

+∞  +∞ 

k2q

 n=0 m=0

D(νp,0) kp G (0, m|0, n|ν)

 

(2n + 1) (2m + 1) Tn(qp) Tm(pq)

n +m

 ν=|n−m|

Nν(q) kp G (0, m|0, n|ν)

 

   +∞  +∞  +∞ n +m    pp −2 D(νp,0) kp G (0, m|0, n|ν) (2n + 1) (2m + 1) (2ℓ + 1) Tn(pp) Tm(pp) Tℓ( )   n=0 m=0 ℓ=0 ν=|n−m|    +∞ +∞ +∞  n +m       (pp) (pp) (pp) −  dν kp G (0, m|0, n|ν) (2n + 1) (2m + 1) (2ℓ + 1) Tn Tm Tℓ     n = 0 m = 0 | ℓ= 0 ν=| n − m  π3     +∞  +∞  +∞ =   pp ( ) 8 +2  2n + 1) (2m + 1) (2ℓ + 1) Tℓ Tn(pp) Tm(pp) (     n=0 m=0 ℓ=0   n+ℓ ℓ+m       (p,0)   (p,0)  ×  D k D k G 0 , ℓ| 0 , n |ν) ( p p   ν µ ν=|n−ℓ| ν=|ℓ−m|

G (0, m|0, ℓ|µ)

(B.3.2)

(B.3.3)

(B.3.4)

258

F. Luppé et al. / Wave Motion 72 (2017) 244–259

k3p



+∞  +∞  +∞ 



δ3(pq)



(2n + 1) (2m + 1) (2ℓ + 1) Tn(qp) Tm(pq) Tℓ(

pp)



 2k k2 − k2   q p  q n=0 m=0 ℓ=0   n + m       ×  D(νq,0) kp G (0, m|0, n|ν)     ν=|n−m|   4 +∞  +∞  +∞      k p −i k b  (qp) (pq) (pp) (2n + 1) (2m + 1) (2ℓ + 1) Tn Tm Tℓ p 2    2 − k2   k n=0 m=0 ℓ=0 p q   n +m       ×  (q) Nν kp G (0, m|0, n|ν)     ν=|n−m|   3 2 +∞ +∞ +∞    =π    kp  (B.3.5) pq qp +i kp b   (2n + 1) (2m + 1) (2ℓ + 1) Tn(qp) Tm(pp) Tℓ( ) + Tn(pq) Tm(pp) Tℓ( )    2 2 2 kp − kq n=0 m=0 ℓ=0     m+ℓ n+ℓ          ×  (q) (p,0) Dν kp Nµ kp G (0, m|0, ℓ|ν) G (0, ℓ|0, n|µ)     ν=|m−ℓ| µ=|n−ℓ|     4 P +∞ +∞ +∞     2  kp 1 (rp) (pq) (qr ) − kp b      2n + 1 2m + 1 2 ℓ + 1 T T T ( ) ( ) ( ) n m ℓ   2 2 2 2 kp − kq r = 1 kp − kr n=0 m=0 ℓ=0     r ̸= p     m +ℓ n +ℓ         × Nν(r ) kp Nµ(q) kp G (0, ℓ|0, n|ν) G (0, m|0, ℓ|µ) ν=|n−ℓ| µ=|m−ℓ|

(pq) δ3Ursell = −

π3 

kp a

3

+∞  +∞ 3 

n+m

 

×

(2n + 1) (2m + 1) Tn(qp) Tm(pq)

n=0 m=0

ν=|n−m| b

+∞

jν kp r h(ν1) kq r u0 (r ) k3p r 2 drG (0, m|0, n|ν)









(B.3.6)

with u0 the non dimensional function defined from function u, u0 ( r ) =



∂U ∂c

 = n 0 =0

3 1 4 π a3

u (r ) .

(B.3.7)

References [1] R.E. Challis, M.J.W. Povey, M.L. Mather, A.K. Holmes, Ultrasound techniques for characterizing colloidal dispersions, Rep. Progr. Phys. 68 (2005) 1541–1637. [2] S. Shahjahan, F. Rupin, A. Aubry, B. Chassignole, F. Fouquet, A. Derode, Comparison between experimental and 2-D numerical studies of multiple R by means of array probes, Ultrasonics 54 (2014) 358–367. scattering in Inconel600⃝ [3] F. Mézière, M. Muller, E. Bossy, A. Derode, Measurements of ultrasound velocity and attenuation in numerical anisotropic porous media compared to Biot’s and multiple scattering models, Ultrasonics 54 (2014) 1146–1154. [4] G. Haïat, A. Lhémery, F. Renaud, F. Padilla, P. Laugier, S. Naili, Velocity dispersion in trabecular bone: Influence of multiple scattering and of absorption, J. Acoust. Soc. Am. 124 (2008) 4047–4058. [5] J. Groenenboom, R. Snieder, Attenuation, dispersion, and anisotropy by multiple scattering of transmitted waves through distributions of scatterers, J. Acoust. Soc. Am. 98 (1995) 3482–3492. [6] Z. Liu, X. Zhang, Y. Mao, Y.Y. Zhu, Z. Yang, C.T. Chan, P. Sheng, Resonant sonic materials, Science 289 (2000) 1734–1736. [7] D.R. Smith, J.B. Pendry, M.C.K. Wiltshire, Metamaterials and negative refractive index, Science 305 (2004) 788–792. [8] Z. Yang, J. Mei, M. Yang, N.H. Chan, P. Sheng, Membrane-type acoustic metamaterial with negative dynamic mass, Phys. Rev. Lett. 101 (2008) 204301 1-4. [9] T. Brunet, A. Merlin, B. Mascaro, K. Zimny, J. Leng, O. Poncelet, C. Aristégui, O. Mondain-Monval, Soft 3D acoustic metamaterial with negative index, Nature Mater. 14 (2015) 384–388. [10] P. Lloyd, M.V. Berry, Wave propagation through an assembly of spheres. IV Relations between different multiple scattering theories, Proc. Phys. Soc. 91 (1967) 678–688. [11] C.M. Linton, P.A. Martin, Multiple scattering by random configurations of circular cylinders: second-order corrections for the effective wavenumber, J. Acoust. Soc. Am. 117 (2005) 3413–3423. and; Multiple scattering by multiple spheres: a new proof of the Lloyd-Berry formula for the effective wavenumber, SIAM J. Appl. Math. 66 (2006) 1649–1668. [12] R.E. Challis, V.J. Pinfield, Ultrasonic wave propagation in concentrated slurries - The modelling problem, Ultrasonics 54 (2014) 1737–1744. [13] P. Sheng, Wave scattering and the Effective medium, in: P. Sheng (Ed.), Introduction to Wave Scattering, Localization, and Mesoscopic Phenomena, Academic Press, New York, 1995, pp. 49–113. [14] F. Vander Meulen, G. Feuillard, O. Bou Matar, F. Levassort, M. Lethiecq, Theoretical and experimental study of the influence of the particle size distribution on acoustic wave properties of strongly inhomogeneous media, J. Acoust. Soc. Am. 110 (2001) 2301–2307. [15] M. Caleap, B.W. Drinkwater, P.D. Wilcox, Coherent acoustic wave propagation in media with pair-correlated spheres, J. Acoust. Soc. Am. 131 (2012) 2036–2047.

F. Luppé et al. / Wave Motion 72 (2017) 244–259

259

[16] L. Tsang, J.A. Kong, K.-H. Ding, Characteristics of discrete scatterers and rough surfaces, in: J.A. Kong (Ed.), Scattering of Electromagnetic Waves. Theories and Applications, Wiley-Interscience, 2000, pp. 173–176. [17] S.K. Kanaun, V.M. Levin, Propagation of longitudinal elastic waves in composites with a random set of spherical inclusions (effective field approach), Arch. Appl. Mech. 77 (2007) 627–651. [18] L. Tsang, J.A. Kong, K.-H. Ding, C.O. Ao, Particle positions for dense media characterizations and simulations, in: Jin Au Kong (Ed.), Scattering of Electromagnetic Waves. Numerical Simulations, Wiley-Interscience, 2001, pp. 403–450. [19] G. Stell, The Percus-Yevick equation for the radial distribution function of a fluid, Physica 29 (1963) 517–534. [20] A.N. Norris, J.M. Conoir, Multiple scattering by cylinders immersed in fluid: high order approximations for the effective wavenumbers, J. Acoust. Soc. Am. 129 (2011) 104–113. [21] J.M. Conoir, A.N. Norris, Effective wavenumbers and reflection coefficients for an elastic medium containing random configurations of cylindrical scatterers, Wave Motion 47 (2010) 183–197. [22] F. Luppé, J.M. Conoir, A.N. Norris, Effective wavenumbers for thermo-viscoelastic media containing random configurations of spherical scatterers, J. Acoust. Soc. Am. 131 (2012) 1113–1120. [23] A. Derode, V. Mamou, A. Tourin, Influence of correlations between scatterers on the attenuation of the coherent wave in a random medium, Phys. Rev. E 74 (2006) 036606.1–036606.9. [24] V. Leroy, A. Strybulevych, J.H. Page, M.G. Scanlon, Influence of positional correlations on the propagation of waves in a complex medium with polydisperse resonant scatterers, Phys. Rev. E 83 (2011) 046605 1-12. [25] D.M. Forrester, J. Huang, V.J. Pinfield, F. Luppé, Experimental verification of nanofluid shear-wave reconversion in ultrasonic fields, Nanoscale 8 (2016) 5497–5506. [26] M. Duranteau, T. Valier-Brasier, J.M. Conoir, R. Wunenburger, Random acoustic metamaterial with a subwavelength dipolar resonance, J. Acoust. Soc. Am. 139 (2016) 3341–3352. [27] R.E. Challis, F. Blarel, M.E. Unwin, X. Guo, On the modelling of ultrasonic bulk wave propagation in epoxies, J. Phys. Conf. Ser. 269 (1) (2011) 012001. [28] M. Chekroun, L. Le Marrec, B. Lombard, J. Piraux, Time-domain numerical simulations of multiple scattering to extract elastic effective wavenumbers, Wave Random Complex 22 (2012) 398–422. [29] S. Kanaun, V. Levin, Propagation of shear elastic waves in composites with a random set of spherical inclusions (effective field approach), Int. J. Solids Struct. 42 (2005) 3971–3997. [30] T. Brunet, K. Zimny, B. Mascaro, O. Sandre, O. Poncelet, C. Aristégui, O. Mondain-Monval, Tuning Mie scattering resonances in soft materials with magnetic fields, Phys. Rev. Lett. 111 (2013) 264301.