<figure><img src="image_1.png" /></figure>

HAL
open science

# Homogenization model for a dense elastic medium of cylindrical scatterers and based on a realistic pair correlation function

Nouhayla Khalid, Michel Darmon, Jean-François Chaix

► To cite this version:

Nouhayla Khalid, Michel Darmon, Jean-François Chaix. Homogenization model for a dense elastic medium of cylindrical scatterers and based on a realistic pair correlation function. Journal of Sound and Vibration, 2026, 620, pp.119442. 10.1016/j.jsv.2025.119442. cea-05371793

HAL Id: cea-05371793

https://cea.hal.science/cea-05371793v1

Submitted on 18 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

ELSEVIER

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

# Homogenization model for a dense elastic medium of cylindrical scatterers and based on a realistic pair correlation function

N. Khalid<sup>a ID</sup>, M. Darmon<sup>a ID</sup>,*, J.-F. Chaix<sup>b ID</sup>

## ARTICLE INFO

### Keywords:

- Analytical modeling
- Ultrasounds
- Simulation
- Concrete
- Multiple scattering

## ABSTRACT

Multiple scattering effects due to a random distribution of identical cylindrical inclusions in an elastic medium are investigated. The approach is based on the analysis proposed by Fikioris and Waterman. The solution of the developed modal equations yields the effective wavenumbers of elastic coherent waves. Attention is made in this paper to a more realistic description of scatterers' distribution in the host medium: this distribution is taken into account in modeling by introducing the notion of pair correlation function. The existing Conoir and Norris approach has been established by using the Hole Correction as pair correlation function, which simplifies the description of scatterers' distribution and may lead to unphysical results for dense media. New formulas for the effective wavenumbers are proposed here by generalization of the latter theory and enable the use of any pair correlation function for possible application to dense media. The new generalized analytical model coupled to a realistic pair correlation function is compared to the previous approach in different materials configurations and validated for concrete structures by comparison to numerical simulations.

## 1. Introduction

The propagation of elastic waves through a heterogeneous medium containing random distributions of inhomogeneities is accompanied by scattering events. The presence of numerous scatterers in a medium leads to interaction with the scattered fields, giving rise to waves that are scattered several times. This phenomenon is known as multiple scattering. The study of multiple scattering of elastic waves is of huge importance for characterization of microstructures in a wide range of disciplines, ranging from seismology, biomedical ultrasound to ultrasonic nondestructive testing (NDT). It is very important to take account of multiple scattering phenomenon for the success of ultrasound applications. Nevertheless, the random distribution of scatterers in heterogeneous media makes theoretical analysis of the multiple scattering phenomenon extremely challenging. The heterogeneous medium is therefore replaced by an effective homogeneous medium in which coherent waves propagate. The coherent waves correspond to the average, over all possible positions of the scatterers, of the multi-scattered waves.

Over the years, many researchers have been interested in explaining the rich phenomenon by developing analytical models and theories expressing the coherent field and therefore describing wave propagation in a heterogeneous medium. Various methods for analyzing multiple scattering of elastic waves have been considered in the literature, including perturbations methods [1] and the classical Stanke and Kino's model for polycrystalline materials [2]; self consistent methods, which are based on an iterative scheme in which, at each iteration, the characteristics of the matrix are replaced by those of the effective medium calculated at the

<footnote><sup>a</sup> Université Paris-Saclay, CEA, List, F-91120 Palaiseau, France</footnote>

<footnote><sup>b</sup> LMA, Aix-Marseille Université, CNRS, Centrale Méditerranée, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: michel.darmon@cea.fr (M. Darmon).</footnote>
previous iteration, until the results converge. One talks about Coherent potential approximation (CPA) [3,4] and homogenization methods derived from it [5-7], or about Generalized Self Consistent model (GSCM) [8]. There are also the effective-field theories, these are based on the theory of homogenization and provide a description of propagation in heterogeneous media in a simplified mathematical form. They investigate the interaction between the field and the different randomly distributed scatterers, and they can be classified into two major categories: static models whose starting equations are the equilibrium equations, and dynamic models that are based on the wave propagation equation. In the former, the notion of frequency does not come into play and the models are only valid at very long wavelengths when scattering is negligible, whereas the dynamic models are frequency dependant models and they propose as effective characteristics the longitudinal and transverse wave numbers which define the phase velocity and attenuation of the wave under consideration. Therefore, they can be classified according to their frequency range of validity. Here we give a brief overview of the most popular dynamic models and approximations. The most famous one is the Born approximation or the simple scattering approximation, which consists of taking into account the waves scattered only once. The Born approximation only works at very low frequencies: at very long wavelengths, the presence of scatterers has little impact on the coherent wave. The Rytov approximation [9], the geometric approximation [10] and the parabolic approximation [11] enable the analysis of scattering effects in the high frequency regime. However, they are still limited to weak scattering problems (media with very low concentration of scatterers or with very low contrasts between the matrix and the scatterers). Dynamic multiple scattering homogenization models can also be classified according to whether they deal with scalar or elastic waves. For scalar waves, Foldy [12] studied the multiple scattering of isotropic point-like scatterers. He initiated the use of the T-matrix formalism [13]. The Independent Scattering Approximation (ISA) model [14] was obtained using another formalism by assuming that each scatterer can only interact once with the wave during a scattering experiment. Its formulation is similar to Foldy's multiple scattering first-order terms in concentration. Both Foldy's and the Independent Scattering Approximation (ISA) are not sufficient to characterize the coherent waves when the density of diffusers increases. Several classical works have proposed improvements of Foldy's theory. Lax [15,16] introduced the Quasi Crystalline Approximation (QCA) for anisotropic scatterers. This approximation, which is more precise than Foldy's, incorporates information related to the correlation between pairs of diffusers, which enables the description of the aggregates' distribution in the microstructure by defining the pair correlation function. Several widely used models were then derived from the QCA. Waterman and Truell (WT) [17] considered both forward and backscattering in their multiple scattering formula, which corrects Foldy's formulation by a second-order term in scatterers concentration. Fikioris and Waterman [18] introduced in the QCA, the “Hole Correction” as a pair correlation function between scatterers. The Hole-Correction is a suitable pair correlation function for sparse densities of inclusions. Vander Meulen [19] introduced in QCA the Percus–Yevick equation (QCA-PY) as a more sophisticated spatial pair correlation function between scatterers. This latter has shown its validity for modeling strong densities of scatterers. Caleap [20] proposed a comparison between three different correlation functions: The Hole Correction, the analytical solution of Percus–Yevick equation for hard spheres, and a Virial expansion of the Percus–Yevick's solution. Linton and Martin [21,22] mathematically proved and also revisited for both circular cylinders and spheres the Lloyd and Berry's model [23] developed initially for spherical scatterers, including the Hole Correction. For elastic waves, Conoir and Norris (CN) [24] developed a model with waves diffracted and converted in all directions, generalizing the Linton and Martin' model (LM) [21] in the acoustic 2D case (matrix containing cylindrical scatterers). Subsequently, Luppé, Conoir and Norris [25] developed a model for the case of elastic scattering of longitudinal and transverse waves from spherical scatterers. In 2017, Luppé, Valier-Brasier, Conoir and Pareige [26] generalized Caleap's model (2012) initially developed for spherical scatterers (liquids or gas) immersed in an ideal fluid host medium. They proposed an expression of the effective wavenumber expanded up to order 3 in concentration and then studied the effect of correlation between scatterers. According to the authors, the effect of correlation between scatterers is of order 3 in concentration and order 2 in scattering. The study was questioned by the same authors most recently; They say that the effect of pair correlation between scatterers appears not only at third order in concentration, but also at second order. Therefore, Luppé, Conoir and Valier-Brasier (LCVB) [27] developed another model taking into account more rigorously both longitudinal and transverse waves in a medium embedded with spheres. They proposed the use of both the scalar and vector addition theorems, which is not considered in the work of Luppé, Conoir and Norris (LCN) [25]. They presented the LCVB model as a corrected version of the LCN model.

To our knowledge, all the models developed for multiple scattering for elastic media embedded with either spherical or cylindrical scatterers use only the simple Hole Correction pair correlation function for their proposed applications, and so do the most recent models of Luppé, Conoir and Valier-Brasier (LCVB) [27] developed for spherical scatterers, and Conoir and Norris (CN) [24] developed for cylindrical scatterers, which makes them valid for modeling the coherent waves in media with a sparse concentration of diffusers. Therefore, these elastic analytical homogenization models show limitations in terms of scatterers concentration. They are limited to the presence of heterogeneities at low concentration. Preliminary evaluation of the CN analytical model for concrete structures has been done in Ref. [28]. Notably unphysical negative values of attenuation are then predicted by the CN model for realistic high densities of aggregates in mortar. This phenomenon has also been reported in different dense media by other authors [20,29] when the Hole Correction is taken as the pair correlation function.

This study deals with modeling the propagation of elastic waves in an isotropic solid medium containing a random distribution of cylindrical scatterers. A generalized homogenization elastic model using any pair correlation function is proposed. The developed model is finally applied to simulate ultrasonic waves propagation in concrete for NDT application [30,31].

This paper is organized as follows: Section 2 is dedicated to the description of the host medium. The Fikioris and Watermann theory is presented, then focusing on the case of cylindrical scatterers in an elastic medium, and expanded to any pair correlation function. The proposed new generalized elastic homogenization model is developed. The simulations obtained with the generalized model are discussed in Section 3, in which the new analytical results are compared to those obtained with numerical tools in the case of cylindrical aggregates in a mortar matrix and to the Conoir and Norris model for different materials and impedance contrasts.
<figure><img src="image_6.png" /><figcaption>Fig. 1. Description of the scattering problem from an elastic cylindrical inclusion of radius a in an homogeneous infinite elastic matrix.</figcaption></figure>

## 2. The generalized Conoir and Norris model (GCN)

In this section, first, the multiple scattering theory of Fikioris and Waterman is presented. The formalism developed for the elastic case by Conoir and Norris is revisited, with the purpose of extending the theory to the case of using any pair correlation function. The aim is, therefore, to generalize the Conoir and Norris model, initially developed for the case of multiple scattering by a random and uniform distribution of elastic cylinders embedded in an elastic medium, with the Hole Correction pair correlation function as a simplified description of the correlation between scatterers.

### 2.1. Coherent waves in Fikioris and Waterman's theory's framework

We reuse here the same theoretical formalism as Conoir and Norris [24] - to which the reader can refer - and summarize it until the appearance of the pair correlation function.

An isotropic elastic medium in which L longitudinal and T transverse waves may propagate is considered. N identical cylindrical scatterers of radius a located at points $r_1, r_2, ..., r_N$ are immersed in the medium. The inclusions are elastic, homogeneous and isotropic. The structural parameters related to the inclusions are indicated by index equal 1, and those related to the matrix by index equal 0, as shown in Fig. 1. Respectively, the densities of the two media are $ρ^0$ and $ρ^1$, and $c_L^0$ and $c_L^1$ are the velocities of the compressional waves, the shear wave velocities are $c_T^0$ and $c_T^1$. Fig. 1 shows a circular elastic inclusion in an elastic homogeneous medium. The reference frame used to describe the system is that of cylindrical coordinates $(e_r, e_θ, z)$ and the medium is infinite in the z direction. The inclusion is insonified by a L or T plane incident wave of unit amplitude, propagating in the x-axis direction $(θ = 0)$. The time dependence of the fields in $e^{-iωt}$, where i is the imaginary unit, ω the angular frequency, and t time, will be omitted throughout the paper.

The potential $ψ_{a,inc}$ of the incident wave can be decomposed as follows:

$$ψ_{a,inc}(\vec{r}) = e^{ik_a x} = \sum_{n=-\infty}^{n=\infty} A_a i^n J_n(k_a r) e^{in\theta}, \quad (1)$$

with $(\vec{r}) = (x, y)$ and where $α ∈ \{L, T\}$ and $k_a$ are the incident mode and wave number, and the $J_n$ are Bessel functions of first kind.

Conoir and Norris introduced then at an observation point r the total potential field $ψ_a$, the sum of the incident field $ψ_{a,inc}$ and the field $ψ_{a,s}$ scattered by the inclusions. The scattered field $ψ_{a,s}(\vec{r}, \vec{r}_k)$ is also linked to the exciting field $ψ_{a,E}(\vec{r}, \vec{r}_k)$ acting on the kth scatterer positioned at $r_k$ by a linear scattering operator called the T-matrix [24] $\overline{\overline{T}(\vec{r}_k)} = \begin{pmatrix} T^{LL} & T^{LT} \\ T^{LT} & T^{TT} \end{pmatrix}$. The expression for the T-matrix coefficients for all the involved modes in the case of scattering by elastic cylindrical scatterers in an elastic medium is presented in [32]. Conoir and Norris expressed the multiple scattering equations as follows:

$$ψ_a(\vec{r}) = ψ_{a,inc}(\vec{r}) + \sum_{k=1}^{N} [T^{Lα}(\vec{r}_k)ψ_{L,E}(\vec{r}, \vec{r}_k) + T^{Tα}(\vec{r}_k)ψ_{T,E}(\vec{r}, \vec{r}_k)], \quad (2)$$

$$ψ_{a,E}(\vec{r}, \vec{r}_k) = ψ_{a,inc}(\vec{r}) + \sum_{\substack{j=1 \\ j \neq k}}^{N} [T^{Lα}(\vec{r}_j)ψ_{L,E}(\vec{r}, \vec{r}_j) + T^{Tα}(\vec{r}_j)ψ_{T,E}(\vec{r}, \vec{r}_j)]. \quad (3)$$

In order to derive the multiple scattering equations governing the coherent field, the configurational average, initially defined by Foldy [12], is applied under the Quasi-Crystalline approximation [15,16]:

$$\langle ψ_a(\vec{r}) \rangle = ψ_{a,inc}(\vec{r}) + \int d\vec{r}_k n(\vec{r}, \vec{r}_k) [T^{Lα}(\vec{r}_k) \langle ψ_{L,E}(\vec{r}, \vec{r}_k) \rangle_{r_k}$$
<figure><img src="image_7.png" /><figcaption>Fig. 2. The exclusion distance between two neighboring scatterers of radius a.</figcaption></figure>

$$+ T^{T\alpha}(\vec{r}_k) \langle \psi_{T,E}(\vec{r}, \vec{r}_k) \rangle_{r_k} \rangle, \quad (4)$$

$$\langle \psi_{\alpha,E}(\vec{r}, \vec{r}_k) \rangle_{r_k} = \psi_{\alpha,\text{inc}}(\vec{r}) + \int d\vec{r}_j n(\vec{r}, \vec{r}_j) [T^{L\alpha}(\vec{r}_j) \langle \psi_{L,E}(\vec{r}, \vec{r}_j) \rangle_{r_j} \\ + T^{T\alpha}(\vec{r}_j) \langle \psi_{T,E}(\vec{r}, \vec{r}_j) \rangle_{r_j}] \quad (5)$$

For a uniform distribution of identical cylinders, the density $n(\vec{r}, \vec{r}_j)$ is usually taken equal to a constant for distances beyond the exclusion distance b, referring to the Hole Correction pair correlation function:

$$n(\vec{r}, \vec{r}_j) = \begin{cases} n_0, & \text{for } |\vec{r} - \vec{r}_j| > b, \\ 0, & \text{otherwise.} \end{cases} \quad (6)$$

where the exclusion distance b or “the distance of closest approach” can be defined as the minimal distance allowed between the centers of two adjacent scatterers as shown in Fig. 2. It is bounded between $b_{min} = 2a < b < b_{max}$ [28,33].

The Hole Correction approximation reflects the mandatory non interpenetration condition required for solid scatterers. Physically, two solid scatterers are not allowed to overlap. The Hole Correction is well known to be a good description of the scatterers distribution for sparse scatterers concentration. However, it is a poor approximation for dense media. Let us note g the pair correlation function. The expression of the Hole Correction pair correlation function is as follows:

$$g(|\vec{r} - \vec{r}_j|) = \begin{cases} 1, & \text{for } |\vec{r} - \vec{r}_j| > b, \\ 0, & \text{otherwise.} \end{cases} \quad (7)$$

More generally, for any expression of the pair correlation function, one can express the scatterer density as follows:

$$n(\vec{r}, \vec{r}_j) = \begin{cases} n_0 g(|\vec{r} - \vec{r}_j|) = n_0(1 + h_2(|\vec{r} - \vec{r}_j|)), & \text{for } |\vec{r} - \vec{r}_j| > b, \\ 0, & \text{otherwise,} \end{cases} \quad (8)$$

where g is the pair correlation function also known as the “radial distribution function” represents the probability of finding a scatterer (located at $\vec{r}$) at a distance $\rho_j = |\vec{r} - \vec{r}_j|$ from another scatterer (located at $\vec{r}_j$), and is related to the “two particle” Ursell function $h_2(|\vec{r} - \vec{r}_j|) = g(|\vec{r} - \vec{r}_j|) - 1$.

Since $\langle \psi_{\alpha,E}(\vec{r}, \vec{r}_k) \rangle_{r_k}$ is satisfying the Helmholtz equation, it can take the following expression:

$$\langle \psi_{\alpha,E}(\vec{r}, \vec{r}_k) \rangle_{r_k} = \sum_{n=-\infty}^{+\infty} A_n^\alpha J_n(k_a \rho_k) e^{in\theta_k}, \quad (9)$$

where $\vec{\rho}_k = \vec{r} - \vec{r}_k$, $\rho_k = |\vec{\rho}_k|$ and $\theta_k = \arg(\vec{\rho}_k)$.

Now, we seek coherent waves in the equivalent homogeneous medium, as solutions of the multiple scattering equations:

$$A_n^\alpha(\vec{r}_j) = i^n A_n^\alpha e^{i\vec{\xi}'\vec{r}_j} + i^n B_n^\alpha e^{i\vec{\xi}'\vec{r}_j} = i^n A_n^\alpha e^{i\xi'x_j} + i^n B_n^\alpha e^{i\xi'x_j}, \quad (10)$$

$\vec{\xi}$ and $\vec{\xi}'$ are the effective wave vectors of coherent waves propagating in an effective medium considered as homogeneous elastic.

From the multiple scattering equations, Conoir and Norris finally derived a modal equation as follows:

$$A_n^\alpha e^{i\vec{\xi}'\vec{r}_1} + B_n^\alpha e^{i\vec{\xi}'\vec{r}_1} = A_n e^{ik_n x_1} + \sum_p i^{-p} [T_{n+p}^{L\alpha} A_{n+p}^L I_p^\alpha(\xi) + T_{n+p}^{L\alpha} B_{n+p}^L I_p^\alpha(\xi')] \\ + T_{n+p}^{T\alpha} A_{n+p}^T I_p^\alpha(\xi) + T_{n+p}^{T\alpha} B_{n+p}^T I_p^\alpha(\xi)], \quad (11)$$

with

$$I_p^\alpha(\xi) = \int d\vec{r}_j \vec{n}(\vec{r}_1, \vec{r}_j) H_p^{(1)}(k_a r_j) e^{ip\theta_j} e^{i\xi r_j}, \quad \xi \in \{\xi, \xi'\}, \alpha \in \{L, T\}. \quad (12)$$

In the previous work of Conoir and Norris [24], the integral $I_p^\alpha(\xi)$ of Eq. (12) was simply calculated by taking into account the Hole Correction assumption. Many works [20,26,29] have shown nonphysical results for high scatterers concentrations in the predictions of effective parameters due to this simplification of the scatterers distribution’s description.
Our goal in the remainder of this section is to develop a generalized model that considers any pair correlation function. To do so, Rohfritsch (in Appendix C of [33]) has calculated the integral of Eq. (12) for application to acoustic waves (taking into consideration the correlation between scatterers through any pair correlation function $g(r) = 1 + h_2(r)$):

$$I_p^\alpha(\xi) = \frac{2i^{\beta+1}}{k_\alpha(\xi - k_\alpha)} e^{ik_\alpha x_1} + \frac{2\pi n_0 e^{i\xi x_1}}{\xi^2 - k_\alpha^2} i^p \mathcal{N}_p^\alpha(\xi), \quad (13)$$

with

$$\mathcal{N}_p^\alpha(\xi) = N_p^\alpha(\xi) + (\xi^2 - k_\alpha^2) M_p(\xi b), \quad (14)$$

and

$$M_p^\alpha(\xi b) = \int_b^{+\infty} H_p^{(1)}(k_\alpha r) J_p(\xi r) h_2(r) dr, \quad (15)$$

$$N_p^\alpha(\xi) = \xi b J'_p(\xi b) H_p^{(1)}(k_\alpha b) - k_\alpha b J_p(\xi b) H_p^{(1)\prime}(k_\alpha b). \quad (16)$$

The Lorentz-Lorenz law.

The Eq. (11) is satisfied whatever the value of $x_1$ if the coefficients of $e^{i\xi x_1}$, $e^{i\xi' x_1}$, $e^{ik_L x_1}$ and $e^{ik_T x_1}$ are all zero. Equating the coefficients of $e^{i\xi' x_1}$ and $e^{i\xi' x_1}$ to zero gives what is known as the Lorentz-Lorenz law. We obtain, by writing this law, a generalization to any pair correlation function of the modal equations of Fikioris and Waterman obtained in the Hole-Correction approximation [18]:

$$A_n^\alpha - \frac{2\pi n_0}{\xi^2 - k_\alpha^2} \sum_p (T_p^{L\alpha} A_p^L + T_p^{T\alpha} A_p^T) \mathcal{N}_{n-p}^\alpha(\xi b) = 0, \quad \alpha \in L, T, \quad (17)$$

$$B_n^\alpha - \frac{2\pi n_0}{\xi^2 - k_\alpha^2} \sum_p (T_p^{L\alpha} B_p^L + T_p^{T\alpha} B_p^T) \mathcal{N}_{n-p}^\alpha(\xi' b) = 0, \quad \alpha \in L, T. \quad (18)$$

Solving the two identical homogeneous systems of equations (17) and (18) determines the effective wave numbers $\xi$ and $\xi'$. Coefficients $A_n^\alpha$ and $B_n^\alpha$ of the development in cylindrical harmonics of the coherent field, are then obtained using the effective wavenumber.

## 2.2. Matrix form of the generalized modal equations of Fikioris and Waterman

The modal equations Eqs. (17) and (18) are similar. We choose to solve the first system for longitudinal and transversal incident waves. For that, we begin by defining:

$$a_\alpha = (..., A_{-1}^\alpha, A_0^\alpha, A_1^\alpha, ..., )^t, \quad \alpha \in L, T, \quad (19)$$

$$\mathbf{e} = (..., 1, 1, 1, ..., )^t. \quad (20)$$

The (17) is then written in the matrix form:

$$a_\alpha - \frac{2\pi n_0}{\xi^2 - k_\alpha^2} \mathcal{N}^\alpha(T^{L\alpha} a_L + T^{T\alpha} a_T) = 0, \quad \alpha \in L, T, \quad (21)$$

where the infinite square matrices $\mathcal{N}^\alpha$ and $T^{\alpha\beta}$ are defined as follows:

$$\mathcal{N}_{mn}^{\alpha} = \mathcal{N}_{m-n}^{\alpha} \quad \text{and} \quad T_{mn}^{\alpha\beta} = T_{n-p}^{\alpha\beta} \delta_{mn}. \quad (22)$$

$\delta_{mn}$ is the Kronecker delta symbol defined as follows: $\delta_{mn} = \begin{cases} 0 & \text{if } m \neq n, \\ 1 & \text{otherwise.} \end{cases}$

Let us note $\overline{Q}_{mn}^\alpha$ the elements of the infinite square matrix $\overline{Q}^\alpha$:

$$\overline{Q}_{mn}^\alpha = \frac{i\pi}{2} \mathcal{N}_{m-n}^\alpha(\xi) - 1 \quad \xi^2 - k_\alpha^2. \quad (23)$$

By introducing the matrix $\overline{Q}^\alpha$, Eq. (21) can be expressed as:

$$a_\alpha - \epsilon(\overline{Q}^\alpha + \frac{ee'}{y_\alpha})(T^{L\alpha} a_L + T^{T\alpha} a_T) = 0, \quad \alpha \in L, T, \quad (24)$$

where the scalars $\epsilon$ and $y_\alpha$ are:

$$\epsilon = -4in_0, \quad (25)$$

$$y_\alpha = \xi^2 - k_\alpha^2, \quad \alpha \in L, T. \quad (26)$$
Knowing that Eq. (24) is valid for longitudinal and transverse incident waves, we combine these two cases, and we form a single equation involving doubly infinite vectors and matrices:

$$ \left\{ \tau - \epsilon \overline{QT} - \epsilon \frac{e_L e_L^t}{y_L} T - \epsilon \frac{e_T e_T^t}{y_T} T \right\} a = 0, \qquad (27) $$

where

$$ e_L = \begin{pmatrix} e \\ 0 \end{pmatrix}, e_T = \begin{pmatrix} 0 \\ e \end{pmatrix}, a = \begin{pmatrix} a_L \\ a_T \end{pmatrix}, e_L e_L^t = e_T e_T^t = 1, \qquad (28) $$

and

$$ \tau = \begin{pmatrix} I & 0 \\ 0 & I \end{pmatrix}, \bar{Q} = \begin{pmatrix} \overline{Q^L} & 0 \\ 0 & \overline{Q^T} \end{pmatrix}, T = \begin{pmatrix} T^{LL} & T^{TL} \\ T^{LT} & T^{TT} \end{pmatrix}. \qquad (29) $$

By multiplying the system given by Eq. (27) from the left by $T^{1/2}$, it is replaced by an equivalent system of simpler form:

$$ \left\{ \tau - \epsilon Q - \frac{\epsilon}{y_L} f_L f_L^t - \frac{\epsilon}{y_T} f_T f_T^t \right\} u = 0, \qquad (30) $$

where

$$ f_L = T^{1/2} e_L, \quad f_T = T^{1/2} e_T, \quad u = T^{1/2} a, \quad \text{and} \quad Q = T^{1/2} \overline{Q} T^{1/2}. \qquad (31) $$

We solve the infinite system of equations given in Eq. (30) in a similar way as in [24] to obtain the fundamental equation for determining the coherent wavenumbers $\xi$:

$$ (\xi^2 - k_L^2 - \epsilon M_{LL})(\xi^2 - k_T^2 - \epsilon M_{TT}) - \epsilon^2 M_{LT} M_{TL} = 0. \qquad (32) $$

with

$$ M_{\alpha\beta} = f_a^t (\tau - \epsilon Q)^{-1} f_\beta, \quad \alpha, \beta \in \{L, T\}. \qquad (33) $$

## 2.3. Asymptotic solutions of the fundamental wavenumber equation

Once we obtained the fundamental equation for determining the coherent wavenumbers, we aim in this section to solve this Eq. (32). One can choose to expand the solutions of the equation above about one of the two leading order solutions $y_\alpha = 0$, $\alpha \in L, T$. Here, the choice is to expand about the longitudinal incident wave root $y_L = 0$. Although, it is important to note that all the solutions can be applied by interchanging $L$ into $T$. At low concentration limit, the asymptotic expansion of $y_L$ in $\epsilon$ is:

$$ y_L = \epsilon y_L^{(1)} + \epsilon^2 y_L^{(2)} + \epsilon^3 y_L^{(3)} + \dots \qquad (34) $$

Substituting the expression of $y_L$ in Eq. (34) into Eq. (32) and replacing $y_T$ by its expression $y_T = y_L + k_L^2 - k_T^2$, one obtains:

$$ \begin{aligned} & \left( y_L^{(1)} - M_{LL} + \epsilon y_L^{(2)} + \epsilon^2 y_L^{(3)} + \mathcal{O}(\epsilon^3) \right) \\ & \times \left( k_L^2 - k_T^2 + \epsilon (y_L^{(1)} - M_{TT}) + \epsilon^2 y_L^{(2)} + \mathcal{O}(\epsilon^3) \right) \\ & - \epsilon M_{LT} M_{TL} = 0. \end{aligned} \qquad (35) $$

The coefficients in the expansion of Eq. (34) are obtained by considering derivatives of Eq. (35) with respect to $\epsilon$ at $\epsilon = 0$ (for small densities $n_0$).

We begin by writing the following expansion for small $\epsilon$:

$$ (\tau - \epsilon Q)^{-1} = \tau + \epsilon Q + \epsilon^2 Q^2 + \dots. \qquad (36) $$

Based on Eqs. (33) and (36), one obtains:

$$ M_{\alpha\beta} = e_a^t T e_\beta + \epsilon e_a^t T \overline{Q} T e_\beta + \epsilon^2 e_a^t T \overline{Q} T \overline{Q} T e_\beta + \mathcal{O}(\epsilon^3), \quad \alpha, \beta \in \{L, T\}. \qquad (37) $$

### Asymptotic expansion at $\mathcal{O}(1)$.

Considering $\epsilon$ very small, the approach consists on treating Eq. (35) as a regular perturbation problem, i.e. make the approximation $\epsilon = 0$, and hence find the solution to the problem. The leading order term yields using also Eq. (37):

$$ y_L^{(1)} = M_{LL}|(\epsilon = 0) = e^t T^{LL} e = \sum_{n=-\infty}^{n=+\infty} T_n^{LL}. \qquad (38) $$
Asymptotic expansion at $\mathcal{O}(\epsilon^1)$.

At the next order, Eq. (35) is equivalent to the following:

$$ (y_L^{(1)} - M_{LL} + \epsilon y_L^{(2)} + \mathcal{O}(\epsilon^2))(k_L^2 - k_T^2 + \epsilon (y_L^{(1)} - M_{TT}) + \mathcal{O}(\epsilon^2)) - \epsilon M_{LT} M_{TL} = 0. \quad (39) $$

In order of $\epsilon M_{LT} M_{TL}$ to be of order $\mathcal{O}(\epsilon^1)$, we stop at order $\mathcal{O}(1)$ for expressions of $M_{LT}$ and $M_{TL}$ in Eq. (37). Thus, $M_{a\beta} = e_a^t T e_\beta$ where $\alpha, \beta \in \{L, T\}$.

The derivative of Eq. (39) with respect to $\epsilon$ at $\epsilon = 0$ implies, using Eq. (38):

$$ \left( y_L^{(2)} - \left. \frac{d M_{LL}}{d \epsilon} \right|_{\epsilon=0} \right) (k_L^2 - k_T^2) + \left( y_L^{(1)} - M_{LL} \right) \left. \left( y_L^{(1)} - M_{TT} \right) \right|_{\epsilon=0} - M_{LT} M_{TL} = 0, \quad (40) $$

$$ (k_L^2 - k_T^2) \left( y_L^{(2)} - \frac{d M_{LL}}{d \epsilon} \right) (\epsilon = 0) - M_{LT} M_{TL} = 0. \quad (41) $$

Using the expansion in Eq. (37), the derivative of $M_{LL}$ with respect to $\epsilon$ at $\epsilon = 0$ is expressed as follows:

$$ \frac{d M_{LL}}{d \epsilon} (\epsilon = 0) = \frac{d}{d \epsilon} (e'_L T e_L + e'_L T \overline{Q} T e_L + \mathcal{O}(\epsilon^2)) (\epsilon = 0) = e'_L T \overline{Q} T e_L, \\ = e' (T^{LL} \overline{Q} T^{LL} + T^{TL} \overline{Q} T^{LT}) e, \quad (42) $$

Consequently, using Eq. (41) implies:

$$ y_L^{(2)} = \frac{d M_{LL}}{d \epsilon} (\epsilon = 0) - \frac{1}{k_T^2 - k_L^2} M_{LT} M_{TL}. \quad (43) $$

We now replace the derivative of $M_{LL}$ from Eq. (42) into Eq. (43), which yields

$$ y_L^{(2)} = e' \left( T^{LL} \overline{Q}^L \right)_{\xi=k_L} T^{LL} e + e' T^{TL} \left[ \overline{Q}^T \right]_{\xi=k_L} - \frac{e e'}{k_T^2 - k_L^2} T^{LT} e. \quad (44) $$

Using Eq. (23) and evaluating $\overline{Q}^L$ and $\overline{Q}^T$ at $\xi = k_L$, one obtains:

$$ y_L^{(2)} = e' T^{LL} \overline{Q}^L(k_L) T^{LL} e + e' T^{TL} \left( \frac{i \pi \mathcal{N}^T(k_L)}{2(k_L^2 - k_T^2)} \right) T^{LT} e. \quad (45) $$

$\overline{Q}^a(\xi)$ and $\mathcal{N}^a(\xi)$ are the square matrices whose elements are respectively $\overline{Q}_{mn}^a(\xi)$ and $\mathcal{N}_{mn}^a(\xi) = \mathcal{N}_{m-n}^a(\xi)$. Using an expansion of the function in Eq. (14) for small $(\xi - k_a)$, $N_{m-n}^a$ becomes:

$$ N_{m-n}^{a}(\xi) = \frac{2}{i\pi} + \frac{2(\xi^2 - k_a^2)}{i\pi k_a^2} D_{m-n}(\xi), \quad (46) $$

where

$$ D_p(\xi) = -\frac{i\pi}{4} \left( (\xi b)^2 - p^2 J_p(\xi b) H_p^{(1)}(\xi b) + (\xi b)^2 J'_p(\xi b) H_p^{(1)\prime}(\xi b) \right). \quad (47) $$

Based on Eqs. (23) and (46), $\overline{Q}_{mn}^a(\xi)$ is expressed as follows:

$$ \overline{Q}_{mn}^{a}(\xi) = \frac{i\pi/2(N_{m-n}^{a}(\xi) + (\xi^2 - k_a^2)M_{m-n}^{a}(\xi b)) - 1}{\xi^2 - k_a^2}, \quad (48) \\ = \frac{D_{m-n}(\xi)}{k_a^2} + i\frac{\pi}{2}M_{m-n}^{a}(\xi b). $$

Using Eq. (48), the expression of $y_L^{(2)}$ given in Eq. (45) becomes:

$$ y_L^{(2)} = \sum_{n=-\infty}^{+\infty} \sum_{m=-\infty}^{+\infty} T_{nn}^{LL} \left( \frac{D_{m-n}(k_L)}{k_L^2} + i \frac{\pi}{2} M_{m-n}^L(k_L b) \right) T_{mm}^{LL} \\ + \sum_{n=-\infty}^{+\infty} \sum_{m=-\infty}^{+\infty} T_{nn}^{TL} \left( \frac{i\pi}{2(k_L^2 - k_T^2)} \right) \\ \times (N_{m-n}^T(k_L) + (k_L^2 - k_T^2) M_{m-n}^T(k_L b)) T_{mm}^{LT}. \quad (49) $$

Finally, based on Eqs. (34), (38) and (49), the generalized version of Conoir and Norris (GCN) analytical model, enabling the use of any pair correlation function, is expressed as follows:

$$ k_{eff}^2 = k_L^2 + \delta_1 n_0 + \delta_2 n_0^2 + \delta_{correlation} n_0^2, \quad (50) $$
<table><caption>**Table 1**

Material parameters for composite pairs (matrix/inclusion), ordered by increasing longitudinal impedance contrast $\Delta Z/Z$.</caption><thead><tr><th>Matrix</th><th>ρ (kg/m<sup>-3</sup>)</th><th>c<sub>L</sub> (m/s)</th><th>c<sub>T</sub> (m/s)</th><th>Inclusion</th><th>ρ (kg/m<sup>-3</sup>)</th><th>c<sub>L</sub> (m/s)</th><th>c<sub>T</sub> (m/s)</th><th>ΔZ/Z</th></tr></thead><tbody><tr><td>Borosilicate</td><td>2230</td><td>5600</td><td>3400</td><td>Ordinary glass</td><td>2500</td><td>5400</td><td>3200</td><td>0.08</td></tr><tr><td>Epoxy</td><td>1150</td><td>2500</td><td>1200</td><td>Polystyrene</td><td>1050</td><td>2350</td><td>1120</td><td>0.15</td></tr><tr><td>Mortar</td><td>2050</td><td>3950</td><td>2250</td><td>Silica aggregate</td><td>2610</td><td>4300</td><td>2470</td><td>0.32</td></tr><tr><td>Epoxy</td><td>1150</td><td>2500</td><td>1200</td><td>Tungsten carbide</td><td>15500</td><td>6700</td><td>4100</td><td>2.80</td></tr></tbody></table>

where:

$$ \delta_1 = -4i \sum_n T_n^{LL}, \qquad (51) $$

$$ \delta_2 = -\frac{16}{k_L^2} \sum_{n,m} D_{m-n}(k_L) T_n^{LL} T_m^{LL} - \frac{8i\pi}{k_L^2 - k_T^2} \sum_{p,q} N_{p-q}^T(k_L) T_p^{TL} T_q^{LT}, \qquad (52) $$

$$ \delta_{\text{correlation}} = -8i\pi \sum_{n,m} M_{m-n}^L(k_L b) T_n^{LL} T_m^{LL} - 8i\pi \sum_{p,q} M_{p-q}^T(k_L b) T_p^{TL} T_q^{LT}, \qquad (53) $$

and $D_{m-n}$, $N_{m-n}^a$ and $M_{m-n}^a$ are explicit, respectively, in (15), (16) and (47) above.

The difference between the GCN and the classical Conoir and Norris model is the addition of the term $\delta_{\text{correlation}}$. The CN model considers only the Hole Correction pair correlation function defined in Eq. (8) in order to describe scatterers' distribution in the host medium. This description is shown to be poor for dense media. The GCN model enables the use of any pair correlation function.

### 3. Evaluation of elastic homogenization models

In the previous section, the 'Generalized Conoir and Norris model (GCN)' has been developed in order to handle the limitations of the CN model for high scatterers concentrations. The CN model is based on a simplified description of the correlation between scatterers, whereas the GCN model enables taking into account a more generalized description of the correlations between scatterers. In this section, the aim is to validate the new GCN model by comparison to numerical simulations and to the CN model.

#### 3.1. Comparison to numerical methods for concrete structures

The objective is to evaluate the validity of the GCN model using two different numerical codes: 2D CIVA FEM and 2D SPECFEM software.

CIVA is a NDT/SHM simulation and analysis platform, including simulation tools for different NDT techniques, notably ultrasonic models intensively validated [34,35]. System models [36] based on the reciprocity principle have been devised to simulate the ultrasonic response of flaws for various NDT configurations [34,37]. 2D CIVA FEM is a transient elastodynamic finite element code integrated into CIVA for various applications. Indeed, it allows simulations of crack-like flaws inspections [38] with a more precise prediction [39] of head waves than models [40–42] developed for canonical shapes; it is henceforth also adapted for concrete modeling [43]. 2D SPECFEM is a numerical code based on the spectral elements numerical method [44] that has also been employed for wave propagation simulation in concrete [45]. Both CIVA 2D FEM and SPECFEM 2D tools use high-order finite elements. Recently, another code (not used here but providing similar simulated results) of CIVA FEM simulation platform has been devised in 2D and 3D computations [46].

The Conoir and Norris approach calculates the effective wavenumber of a longitudinal or transverse incident wave propagating in an equivalent elastic medium. The pair correlation function (HC), used in this model, follows a law of equiprobability of scatterers' existence beyond an exclusion distance b between scatterers. The validity of existing homogenization theories have been widely discussed in a review notably devoted to concrete applications [47]. The evaluation of the Conoir and Norris model in the case of concrete structures has been recently studied in the following paper [28]. It has been shown that above 30% of aggregate concentration, the Conoir and Norris model exhibits a shift to negative unphysical values of the longitudinal effective attenuation in the low-frequency regime. Besides, an underestimation of those attenuation values takes place in the high-frequency regime compared to 2D CIVA and SPECFEM 2D numerical results. The reasons behind this disagreement between the analytical Conoir and Norris model and FEM numerical simulations may be due to two reasons: low concentration approximation and the simplified description of the distribution of aggregates in the medium. In order to confirm the hypotheses, we seek to evaluate the new model established above in Eq. (50), that takes into account the realistic pair correlation function.

In the following, elastic cylindrical aggregates of radius 6 mm are randomly distributed in a mortar matrix. The two materials are characterized by their mass density $\rho$ and velocities of the longitudinal and transverse waves, respectively $c_L$ and $c_T$. Their properties are given in Table 1 below.

The only difference between Eq. (61) of Ref. [24] and Eq. (50) of the current paper comes from the correlation term in $\delta_{\text{correlation}}$. The correlation term cancels out the simplification of the scatterer distribution description implemented by the Hole Correction function. The existence of a scatterer at a fixed distance from another scatterer is no longer considered equiprobable. A spatial correlation phenomena is introduced.
<figure><img src="image_8.png" /><figcaption>Fig. 3. 2D Pair distribution functions obtained with the Percus Yevick analytical solution for different concentrations of hard disks [50].</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 4. (a) Phase velocity and (b) attenuation as a function of frequency or of the nondimensional wave number (product of the wave number of the longitudinal wave in the matrix and the aggregate radius) for an incident longitudinal wave. Comparison between GCN (solid lines), CN (points), 2D CIVA FEM (dashed lines), and SPECFEM 2D (stars).</figcaption></figure>

Concrete is a dense, disordered system, however heterogeneities present in the concrete structure are likely to be spatially correlated. The introduction of correlations in the medium has then a strong impact on its effective properties.

The pair correlation function $g_2$ can then be used to provide information on the spatial correlation between scaterrers. A lot of work has been done to find out precisely the function $h_2(r)$. In 1958, Percus and Yevick developed a self-consistent method for accessing this quantity in three dimensions through an iterative process [48]. An integral formulation was proposed by Wertheim for the case of hard sphere interactions [49]. For even dimensions, Adda-Bedia proposed a semi-analytical solution for the integral equation of Percus and Yevick [50].

The use of the Percus–Yevick (PY) approximation in our present work follows widespread practice, acknowledging its historical role and analytical tractability. However, the Percus–Yevick pair correlation function is not considered as a universal “gold standard”, and its limitations particularly at very high densities are now recognized in recent comparative studies [51,52]. Nevertheless, in Chapter 4 of [52], the Percus–Yevick solution is shown to be fast to calculate and relatively valid in the tested configurations until around 60% concentration (maximal value for tests in the current paper).

As shown in Fig. 3, the visible oscillations on the Percus–Yevick pair correlation function characterize the correlation between the diffusers. They are particularly noticeable at high scatterer densities. For low concentrations, the correlation function converges to an equiprobable distribution (the Hole Correction pair correlation function). Consequently, incorporating the Percus Yevick pair correlation function into the new analytical model as a more realistic pair correlation function should theoretically lead to improved predictions of the effective parameters for strong densities.

In order to validate the new GCN model given by Eq. (50) for the case of the Percus Yevick pair correlation function, a simulation study is carried out for aggregate concentrations ranging from 12% to 42%. The exclusion distance is fixed to $b = 2a + d$, with $a$ being the radius of the scatterer, and $d$ is the minimum distance allowed between the edges of two scatterers. We suppose that it is unrealistic to have two scatterers set side by side in concrete, $d$ is then fixed to a non-zero value $d = 0.5$ mm in all the 2D simulations either analytical (GCN) or numerical (2D CIVA and SPECFEM 2D).
<figure><img src="image_10.png" /><figcaption>Fig. 5. Zoom of (a) phase velocity and (b) attenuation (shown in Fig. 4) for low frequencies (20 kHz–200 kHz). Comparison between GCN (solid lines) and CN (points) models.</figcaption></figure>

In Fig. 4, a good agreement between the CN analytical model and the two numerical ones for 12 and 24% concentration is observed. However, disagreement increases with the concentration. For 42% concentration, not only an unexpected shift to negative values of the effective attenuation in the low-frequency regime (20 kHz–70 kHz) takes place, as shown in the zoomed-in Fig. 5, but also an underestimation of those values in the high-frequency regime is noticed in Fig. 4(b).

Applying the Generalized Conoir and Norris (GCN) model with the Percus Yevick pair correlation function (GCN-PY), one can observe that there is a small (but increasing) effect for 12% and 24% of inclusion densities. The impact of integrating the pair correlation function into the GCN model is indeed visible for high densities. For 42% of scatterers concentration in the low-frequency regime, the shift to negative values of attenuation is not observed anymore, as shown in Fig. 5, and the overall prediction of the phase velocity is closer to the numerical one. The GCN model allows to better capture the effect of resonances on the attenuation: At low frequencies, mode conversions are significant, especially at the vicinity of dipolar resonance frequency (around 120 kHz) [see Chapter 7 of [52]] whose influence and related peak frequency are better predicted by GCN than CN. Nevertheless, the underestimation of the attenuation values in high frequencies is replaced by an overestimation of those values compared to the numerical results. On the other hand, a better overall agreement is noticed between the GCN-PY model and the numerical simulations.

## 3.2. Comparison between GCN and CN for different impedance contrasts

The aim of this part is to study the interest that the GCN can provide for modeling propagation in dense media. We focused on one possible interest: the ability to remove unphysical attenuations predicted by the previous CN model. From our experience, negative attenuation may occur: - At high volume fractions of inclusions (typically on average > 20%–30%) - At low frequencies (where many homogenization theories are expected to be valid) - For high contrasts in acoustic/elastic properties between the matrix and the inclusion. Therefore, different simulations with both the GCN and CN models assuming $b = 2a$ have been carried out for a 397 µm inclusion diameter for different pairs of materials (matrix/inclusion) as defined in Table 1 and only attenuations are plotted below. The cases of ordinary Glass in Borosilicate Glass (Fig. 6(a)) and of polystyrene (PS) in epoxy (Fig. 6(b)) correspond to low impedance contrasts (<0.2), whereas concrete – studied in the previous subsection – refers to a medium one (0.2 < contrast < 1) and Tungsten Carbide in Epoxy (Fig. 7) is a case of very high contrast (>1). For the low contrast cases, positive attenuations are well predicted by the GCN model for very high concentrations (42 and 60%) for which the CN model fails. For concrete (medium contrast), negative attenuation have also been removed by GCN for 42% concentration. For the very high contrast case of Tungsten Carbide in Epoxy (Fig. 7), the GCN improves the simulation until around 20% concentration; above this limit, it does not allow to avoid negative attenuations. Nevertheless, for all cases, the validity of the GCN model is much better than that of the CN model; it enables to handle much higher concentrations than the CN model.

## 4. Conclusion

The propagation of elastic waves in dense media has been studied. A generalized analytical elastic multiple scattering model has been developed for cylindrical scatterers. It enables the introduction of any pair correlation function (which is linked to the distribution of distances between two different scatterers, one being fixed). The solution of the Percus Yevick equation is used here to provide a realistic description of the pair correlation function for relatively high scatterers concentration. The new model has been applied to concrete structures and shows an overall good agreement by comparison to numerical tools: Unphysical results in
<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Attenuation for (a) ordinary Glass in Borosilicate Glass and (b) polystyrene (PS) in epoxy, as a function of frequency for an incident longitudinal wave and different concentrations (from 12 to 60%). Comparison between GCN (solid lines) and CN (points) models.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 7.** Attenuation for (a) Tungsten Carbide in Epoxy, as a function of frequency for an incident longitudinal wave and a 17.5% concentration. (b) Zoom of the attenuation for low frequencies (25 kHz–350 kHz). Comparison between GCN (solid lines) and CN (points) models.</figcaption></figure>

the low frequency regime initially noticed in the literature have been corrected. Further numerical and experimental validations of the new model can be performed, notably by studying the effect of pores in the mortar matrix. Other comparisons to a previous homogenization model have been carried out for materials leading to different contrast between inclusions and host medium: the new model leads to a more physical behavior than the existing model and extends its validity range to higher concentrations. It is undeniable that, as an approximate solution, the new model owns a validity which depends on the concentration and on the impedance contrast and is less guaranteed for very high values of the latter parameters. However, it could represent a very useful and fast tool – compared to numerical models – to simulate coherent waves in many elastic heterogeneous engineering materials as composites, cementitious,...

For spherical scatterers, Luppé, Conoir and Valier-Brasier [27] have suggested a corrected version of the Luppé Conoir and Norris model [25] and validated the proposed model for low concentrations under a simplified pair correlation function. Work is underway to evaluate the model using pair correlation functions realistic for dense media.

CRediT authorship contribution statement

**N. Khalid:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization.
**M. Darmon:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.
**J.-F. Chaix:** Writing – review & editing, Validation, Supervision, Software, Project administration, Methodology, Investigation, Formal analysis.
## Funding

This work was supported by the French Alternative Energies and Atomic Energy Commission (CEA), France.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

The authors do not have permission to share data.

## References

[1] F.C. Karal Jr., J.B. Keller, Elastic, electromagnetic, and other waves in a random medium, J. Math. Phys. 5 (4) (1964) 537-547, http://dx.doi.org/10.1063/1.1704145.

[2] F.E. Stanke, G.S. Kino, A unified theory for elastic wave propagation in polycrystalline materials, J. Acoust. Soc. Am. 75 (3) (1984) 665-681, http://dx.doi.org/10.1121/1.390577.

[3] A.J. Devaney, Multiple scattering theory for discrete, elastic, random media, J. Math. Phys. 21 (11) (1980) 2603-2611, http://dx.doi.org/10.1063/1.524373.

[4] J.-Y. Kim, J.-G. Ih, B.-H. Lee, Dispersion of elastic waves in random particulate composites, J. Acoust. Soc. Am. 97 (3) (1995) 1380-1388, http://dx.doi.org/10.1121/1.412080.

[5] F.J. Sabina, J.R. Willis, A simple self-consistent analysis of wave propagation in particulate composites, Wave Motion 10 (2) (1988) 127-142, http://dx.doi.org/10.1016/0165-2125(88)90038-8.

[6] P.G.J. Bussink, P.L. Iske, J. Oortwijn, G.L.M.M. Verbist, Self-consistent analysis of elastic wave propagation in two-dimensional matrix-inclusion composites, J. Mech. Phys. Solids 43 (10) (1995) 1673-1690, http://dx.doi.org/10.1016/0022-5096(95)00048-N.

[7] S.K. Kanaun, V.M. Levin, F.J. Sabina, Propagation of elastic waves in composites with random set of spherical inclusions (effective medium approach), Wave Motion 40 (1) (2004) 69-88, http://dx.doi.org/10.1016/j.wavemoti.2003.12.013.

[8] R.-B. Yang, A dynamic generalized self-consistent model for wave propagation in particulate composites, J. Appl. Mech. 70 (4) (2003) 575-582, http://dx.doi.org/10.1115/1.1576806.

[9] S.M. Rytov, Y.A. Kravtsov, V.I. Tatarskii, Principles of Statistical Radiophysics 4: Wave Propagation through Random Media, Springer Verlag, Berlin, 1989.

[10] H. Sato, M.C. Fehler, Seismic Wave Propagation and Scattering in the Heterogeneous Earth, Springer, Berlin, Heidelberg, 2009, http://dx.doi.org/10.1007/978-3-540-89623-4.

[11] C.E.R. Line, R.W. Hobbs, J.A. Hudson, D.B. Snyder, Statistical inversion of controlled-source seismic data using parabolic wave scattering theory, Geophys. J. Int. 132 (1) (1998) 61-78, http://dx.doi.org/10.1046/j.1365-246x.1998.00398.x.

[12] L.L. Foldy, The multiple scattering of waves. I. General theory of isotropic scattering by randomly distributed scatterers, Phys. Rev. 67 (3-4) (1945) 107-119, http://dx.doi.org/10.1103/PhysRev.67.107.

[13] V.K. Varadan, V.V. Varadan, Acoustic, Electromagnetic and Elastic Wave Scattering - Focus on the T-Matrix Approach; Proceedings of the International Symposium, Ohio State University, Columbus, Ohio, June 25-27, 1979, Pergamon Press, Oxford, 1980.

[14] A. Ishimaru, Wave Propagation and Scattering in Random Media. Volume 1 - Single Scattering and Transport Theory, vol. 1, Academic Press, 1978, http://dx.doi.org/10.1016/B978-0-12-374701-3.X5001-7.

[15] M. Lax, Multiple scattering of waves, Rev. Modern Phys. 23 (4) (1951) 287-310, http://dx.doi.org/10.1103/RevModPhys.23.287.

[16] M. Lax, Multiple scattering of waves. II. The effective field in dense systems, Phys. Rev. 85 (4) (1952) 621-629, http://dx.doi.org/10.1103/PhysRev.85.621.

[17] P.C. Waterman, R. Truell, Multiple scattering of waves, J. Math. Phys. 2 (4) (1961) 512-537, http://dx.doi.org/10.1063/1.1703737.

[18] J.G. Fikioris, P.C. Waterman, Multiple scattering of waves. II. "Hole corrections" in the scalar case, J. Math. Phys. 5 (10) (1964) 1413-1420, http://dx.doi.org/10.1063/1.1704077.

[19] F.V. Meulen, G. Feuillard, O.B. 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 (5) (2001) 2301-2307, http://dx.doi.org/10.1121/1.1404435.

[20] M. Caleap, B.W. Drinkwater, P.D. Wilcox, Coherent acoustic wave propagation in media with pair-correlated spheres, J. Acoust. Soc. Am. 131 (3) (2012) 2036-2047, http://dx.doi.org/10.1121/1.3675011.

[21] 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 (6) (2005) 3413-3423, http://dx.doi.org/10.1121/1.1904270.

[22] C.M. Linton, P.A. Martin, Multiple scattering by multiple spheres: A new proof of the lloyd-berry formula for the effective wavenumber, SIAM J. Appl. Math. 66 (5) (2006) 1649-1668, http://dx.doi.org/10.1137/050636401.

[23] C. Zhou, K. Li, J. Han, Characterizing the effect of compressive damage on transport properties of cracked concretes, Mater. Struct. 45 (3) (2012) 381-392, http://dx.doi.org/10.1617/s11527-011-9771-4.

[24] J.-M. Conoir, A.N. Norris, Effective wavenumbers and reflection coefficients for an elastic medium containing random configurations of cylindrical scatterers, Wave Motion 47 (3) (2010) 183-197, http://dx.doi.org/10.1016/j.wavemoti.2009.09.004.

[25] F. Luppé, J.-M. Conoir, A.N. Norris, Effective wave numbers for thermo-viscoelastic media containing random configurations of spherical scatterers, J. Acoust. Soc. Am. 131 (2) (2012) 1113-1120, http://dx.doi.org/10.1121/1.3672690.

[26] F. Luppé, T. Valier-Brasier, J.-M. Conoir, P. Pareige, Coherent wave propagation in viscoelastic media with mode conversions and pair-correlated scatterers, Wave Motion 72 (2017) 244-259, http://dx.doi.org/10.1016/j.wavemoti.2017.03.002.

[27] F. Luppé, J.-M. Conoir, T. Valier-Brasier, Longitudinal and transverse coherent waves in media containing randomly distributed spheres, Wave Motion 115 (2022) 103082, http://dx.doi.org/10.1016/j.wavemoti.2022.103082.

[28] N. Khalid, M. Darmon, J.-F. Chaix, 2D analytical modelling for ultrasonic inspection of concrete structures: Effects of scatterers position correlation, J. Phys.: Conf. Ser. 2768 (1) (2024) 012002, http://dx.doi.org/10.1088/1742-6596/2768/1/012002.

[29] M. Chekroun, L.L. Marrec, B. Lombard, J. Piraux, Time-domain numerical simulations of multiple scattering to extract elastic effective wavenumbers, Waves Random Complex Media 22 (3) (2012) 398-422, http://dx.doi.org/10.1080/17455030.2012.704432.

[30] M. Darmon (Ed.), Ultrasonic Modelling for Non-destructive Testing, MDPI, 2024, http://dx.doi.org/10.3390/books978-3-7258-0758-1.

[31] M. Darmon, Special issue on ultrasonic modeling for non-destructive testing, Appl. Sci. 14 (5) (2024) 2008, http://dx.doi.org/10.3390/app14052008.

[32] M. Chekroun, Caractérisation mécanique des premiers centimètres du béton avec des ondes de surface (Ph.D. thesis), Nantes Université, 2008.
[33] A. Rohfritsch, Diffusion multiple en milieux désordonnés : influence des propriétés macroscopiques et microscopiques sur les ondes cohérentes (Ph.D. thesis), Sorbonne Université, 2021.

[34] M. Darmon, A. Ferrand, V. Dorval, S. Chatillon, S. Lonné, Recent modelling advances for ultrasonic TOFD inspections, AIP Conf. Proc. 1650 (1) (2015) 1757-1765, http://dx.doi.org/10.1063/1.4914799.

[35] A. Ferrand, M. Darmon, S. Chatillon, M. Deschamps, Modeling of ray paths of head waves on irregular interfaces in TOFD inspection for NDE, Ultrasonics 54 (7) (2014) 1851-1860, http://dx.doi.org/10.1016/j.ultras.2013.12.007.

[36] V. Dorval, S. Chatillon, B. Lu, M. Darmon, S. Mahaut, A general Kirchhoff approximation for echo simulation in ultrasonic NDT, AIP Conf. Proc. 1430 (1) (2012) 193-200, http://dx.doi.org/10.1063/1.4716230.

[37] G. Toullelan, R. Raillon, S. Chatillon, V. Dorval, M. Darmon, S. Lonné, Results of the 2015 UT modeling benchmark obtained with models implemented in CIVA, AIP Conf. Proc. 1706 (1) (2016) 190001, http://dx.doi.org/10.1063/1.4940642.

[38] A. Imperiale, S. Chatillon, M. Darmon, N. Leymarie, E. Demaldent, UT simulation using a fully automated 3D hybrid model: Application to planar backwall breaking defects inspection, AIP Conf. Proc. 1949 (1) (2018) 050004, http://dx.doi.org/10.1063/1.5031546.

[39] M. Darmon, G. Toullelan, V. Dorval, An experimental and theoretical comparison of 3D models for ultrasonic non-destructive testing of cracks: Part I, embedded cracks, Appl. Sci. 12 (10) (2022) 5078, http://dx.doi.org/10.3390/app12105078.

[40] L.J. Fradkin, M. Darmon, S. Chatillon, P. Calmon, A semi-numerical model for near-critical angle scattering, J. Acoust. Soc. Am. 139 (1) (2016) 141, http://dx.doi.org/10.1121/1.4939494.

[41] V. Zernov, L. Fradkin, A. Gautesen, M. Darmon, P. Calmon, Wedge diffraction of a critically incident Gaussian beam, Wave Motion 50 (4) (2013) 708-722, http://dx.doi.org/10.1016/j.wavemoti.2013.01.004.

[42] L.J. Fradkin, A.K. Djakou, C. Prior, M. Darmon, S. Chatillon, P.-F. Calmon, The alternative Kirchhoff approximation in elastodynamics with applications in ultrasonic nondestructive testing, ANZIAM J. 62 (4) (2020) 406-422, http://dx.doi.org/10.1017/S1446181120000036.

[43] V. Dorval, A. Imperiale, M. Darmon, E. Demaldent, J.-M. Hénault, FEM-based simulation tools for ultrasonic concrete inspection, E-J. Nondestruct. Test. 27 (09) (2022) http://dx.doi.org/10.58286/27228.

[44] D. Komatitsch, J. Tromp, Introduction to the spectral element method for three-dimensional seismic wave propagation, Geophys. J. Int. 139 (3) (1999) 806-822, http://dx.doi.org/10.1046/j.1365-246x.1999.00967.x.

[45] T. Yu, J.-F. Chaix, L. Audibert, D. Komatitsch, V. Garnier, J.-M. Hénault, Simulations of ultrasonic wave propagation in concrete based on a two-dimensional numerical model validated analytically and experimentally, Ultrasonics 92 (2019) 21-34, http://dx.doi.org/10.1016/j.ultras.2018.07.018.

[46] V. Dorval, A. Imperiale, E. Demaldent, Determination of effective properties of concrete by finite element simulations for ultrasonic NDT, Acad. J. Civ. Eng. 41 (4) (2023) 117-124, http://dx.doi.org/10.26168/ajce.41.4.14.

[47] M. Darmon, N. Khalid, M. Ramaniraka, V. Dorval, J.-M. Henault, J.-F. Chaix, Modeling of ultrasonic velocity and attenuation for concrete inspection: a review, Ultrasonics 129 (14) (2025) 107690, http://dx.doi.org/10.1063/1.2991338.

[48] J.K. Percus, G.J. Yevick, Analysis of classical statistical mechanics by means of collective coordinates, Phys. Rev. 110 (1) (1958) 1-13, http://dx.doi.org/10.1103/PhysRev.110.1.

[49] M.S. Wertheim, Exact solution of the Percus-Yevick integral equation for hard spheres, Phys. Rev. Lett. 10 (8) (1963) 321-323, http://dx.doi.org/10.1103/PhysRevLett.10.321.

[50] M. Adda-Bedia, E. Katzav, D. Vella, Solution of the Percus-Yevick equation for hard hyperspheres in even dimensions, J. Chem. Phys. 129 (14) (2008) 144506, http://dx.doi.org/10.1063/1.2991338.

[51] I. Pihlajamaa, L.M. Janssen, Comparison of integral equation theories of the liquid state, Phys. Rev. E 110 (4) (2024) 044608, http://dx.doi.org/10.1103/PhysRevE.110.044608.

[52] N. Khalid, Advanced Semi-Analytical Modeling for Ultrasonic Inspection of Concrete Structure (Ph.D. thesis), Université Paris-Saclay, 2024.