# Acoustic wave propagation in a macroscopically inhomogeneous porous medium saturated by a fluid

Laurent de Ryck, Jean-Philippe Groby, Philippe Leclaire, Walter Lauriks,

Armand Wirgin, Zine El Abiddine Fellah, Claude Depollier

► To cite this version:

Laurent de Ryck, Jean-Philippe Groby, Philippe Leclaire, Walter Lauriks, Armand Wirgin, et al. Acoustic wave propagation in a macroscopically inhomogeneous porous medium saturated by a fluid. Applied Physics Letters, 2007, 10.1063/1.2431570. [hal-00077652v4](https://hal.science/hal-00077652v4)

HAL Id: hal-00077652

https://hal.science/hal-00077652v4

Submitted on 1 Jun 2016

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Acoustic wave propagation in a macroscopically inhomogeneous porous medium saturated by a fluid

L. De Ryck, J.-P. Groby, P. Leclaire, and W. Lauriks

Laboratorium voor Akoestiek en Thermische Fysica, K. U. Leuven, Celestijnenlaan 200D, 3001 Heverlee, Belgium

A. Wirgin and Z. E. A. Fellah

Laboratoire de Mécanique et d'Acoustique, UPR 7051 CNRS, 31 chemin Joseph Aiguier, 13009 Marseille, France

C. Depollier

Laboratoire d'Acoustique de l'Université du Maine, Avenue Olivier Messiaen, 72000 Le Mans, France

The equations of motion in a macroscopically inhomogeneous porous medium saturated by a fluid are derived. As a verification of their validity, these equations are reduced and solved for rigid frame porous systems. The reflection and transmission coefficients and their corresponding time signals are calculated numerically using wave splitting Green's function approach for a two-layer porous system considered as one single porous layer with a sudden change in physical properties. The results are compared to experimental results and to those of the classical transfer matrix method for materials saturated by air in the ultrasonic frequency range.

An inhomogeneous medium is one with properties that vary with position. Inhomogeneous and layered materials can be found in many fields of physics, for instance, in optics and electromagnetism¹,² or in acoustics.³,⁴ Other examples are geophysical⁵ or granular⁶ media. The study of the acoustic wave propagation in inhomogeneous porous media is of great interest in building and civil engineering and in petroleum prospection.

In this letter, acoustic propagation in fluid-saturated macroscopically inhomogeneous porous materials is studied. These media have received far less attention than homogeneous porous media.⁷⁻¹⁰ It is assumed that the wavelengths are greater than the average heterogeneity size at the pore scale so that the physical properties are homogenized. However, these properties can vary with the observation point within the material at the macroscopic scale of the specimen. The equations of motion are derived from Biot's alternative formulation of 1962 (Ref. 11) in which the total stress tensor, the fluid pressure, the solid displacement, and the fluid/solid relative displacement are used. It was briefly stated by Biot¹¹ and confirmed¹² that these variables should be employed for porous media with inhomogeneous porosity. The equations of motion for an inhomogeneous porous medium with elastic skeleton are derived. A verification of their validity is the study of a porous material saturated by air in the rigid frame approximation, i.e., when the fluid is light and the solid skeleton immobile.¹⁰ In this case, a wave equation is derived and solved numerically for a two-layer porous system treated as one single porous medium with a sudden but continuous change in physical properties. This provides an excellent means of comparing the proposed method: wave splitting Green's function approach (WS-GF),² which is applicable to any depth-dependent inhomogeneous system, to the results of the well established transfer matrix method (TMM) devel-

oped to calculate the acoustical properties of multilayer porous systems.¹⁰ Experimental results obtained at ultrasonic frequencies in a two-layer material saturated by air are compared to the simulations.

The constitutive linear stress-strain relations in an initially stress-free, isotropic porous medium are¹¹

$$ \sigma_{ij} = 2\mu\epsilon_{ij} + \delta_{ij}(\lambda_c\theta - \alpha M\zeta), \quad (1) $$

$$ p = M(-\alpha\theta + \zeta), \quad (2) $$

where $\sigma_{ij}$ is the total stress tensor and $p$ the fluid pressure in the pores; $\delta_{ij}$ denotes the Kronecker symbol (the summation on repeated indices is implied); $\theta=\nabla \cdot \mathbf{u}$ and $\zeta=-\nabla \cdot \mathbf{w}$ are, respectively, the dilatation of the solid and the variation of fluid content where $\mathbf{u}$ is the solid displacement and $\mathbf{w} = \phi(\mathbf{U}-\mathbf{u})$ the fluid/solid relative displacement ($\mathbf{U}$ is the fluid displacement); $\phi$ is the porosity; $\epsilon_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})$ the strain tensor of the solid (the comma denotes spatial partial derivatives); $\lambda_c=\lambda+\alpha^2 M$, where $\lambda, \mu, M$ are elastic constants and $\alpha$ a coefficient of elastic coupling. These parameters were defined by Biot and Willis.¹³ Applying the momentum conservation law in the absence of body forces, the equations of motion are written

$$ \nabla \cdot \boldsymbol{\sigma} = \rho \ddot{\mathbf{u}} + \rho_f \ddot{\mathbf{w}}, \quad (3) $$

$$ -\nabla p = \rho_f \ddot{\mathbf{u}} + m \ddot{\mathbf{w}} + \frac{\eta}{\kappa} F \dot{\mathbf{w}}, \quad (4) $$

where the dot and double dot notations refer to the first and second order time derivatives, respectively; $\rho_f$ is the density of the fluid in the (interconnected) pores, $\rho$ the bulk density of the porous medium, such that $\rho=(1-\phi)\rho_s+\phi\rho_f$, where $\rho_s$ is the density of the solid; $m=\rho_f\tau_\infty/\phi$ is a mass parameter defined by Biot,¹¹ $\tau_\infty$ is the tortuosity, $\eta$ the viscosity of the fluid, $\kappa$ the permeability, and $F$ the viscosity correction function. For an inhomogeneous porous layer or a half space
<figure><img src="image_2.png" /><figcaption>FIG. 1. Slab of inhomogeneous porous material and example of a simulated inhomogeneity profile for $\phi(x)$.</figcaption></figure>

whose properties vary along the depth $x$, the following parameters in the above equations are now dependent on $x$: $\lambda$, $\mu$, $\lambda_c$, $\alpha$, $M$, $\phi$, $\rho$, $m$, $\tau_\infty$, $\kappa$, and $F$. Inserting Eqs. (1) and (2) into Eqs. (3) and (4) yields the equations of motion in terms of the displacements

$$ \nabla[(\lambda_c + 2\mu) \nabla \cdot \mathbf{u} + \alpha M \nabla \cdot \mathbf{w}] - \nabla \wedge [\mu \nabla \wedge \mathbf{u}] - 2 \nabla \mu \nabla \cdot \mathbf{u} + 2 \nabla \mu \wedge (\nabla \wedge \mathbf{u}) + 2[\nabla \mu \cdot \nabla]\mathbf{u} = \rho\ddot{\mathbf{u}} + \rho_f \ddot{\mathbf{w}}, $$

$$ \nabla[M \nabla \cdot \mathbf{w} + \alpha M \nabla \cdot \mathbf{u}] = \rho_f \ddot{\mathbf{u}} + m \ddot{\mathbf{w}} + \frac{\eta}{\kappa} F \dot{\mathbf{w}}, \quad (5) $$

where the $x$ dependence of the constitutive parameters has been removed to simplify the notations.

Under the assumption of a rigid frame, $\mathbf{u}=0$, Eqs. (1) and (3) vanish. The system of equations reduces to Eqs. (2) and (4) and can be written in the frequency domain in a more suitable form for acoustical applications as

$$ -j\omega p = K_e(x, \omega) \nabla \cdot [\phi(x) \dot{\mathbf{U}}], \quad (6) $$

$$ -\nabla p = j\omega\rho_e(x, \omega)\phi(x)\dot{\mathbf{U}}, \quad (7) $$

where $\rho_e(x, \omega)$ and $K_e(x, \omega)$ are, respectively, the effective density and bulk modulus of the inhomogeneous equivalent fluid. Their expressions are

$$ \rho_e(x, \omega) = \rho_f \frac{\tau_\infty(x)}{\phi(x)} \left[ 1 - j \frac{R_f(x) \phi(x)}{\omega \rho_f \tau_\infty} F(x, \omega) \right], \quad (8) $$

$$ K_e(x, \omega) = \frac{\gamma P_0 / \phi(x)}{\gamma - (\gamma - 1) [1 - j R_f(x) \phi(x) G(x, B^2 \omega) / B^2 \omega \rho_f \tau_\infty]^{-1}}, \quad (9) $$

where $\gamma$ is the specific heat ratio, $P_0$ the atmospheric pressure, and $B^2$ the Prandtl number. The correction functions $F(x, \omega)$ and $G(x, B^2\omega)$ depend on $x$ for inhomogeneous media. They are well-defined functions incorporating, respectively, the viscous characteristic length $\Lambda$ of Johnson et al.⁹ and the thermal characteristic length $\Lambda'$ of Champoux and Allard.¹⁴ The effective density and bulk modulus of the inhomogeneous equivalent fluid are functions of the frequency-independent parameters $\phi(x)$, $\tau_\infty(x)$, $\Lambda(x)$, $\Lambda'(x)$, and of the flow resistivity $R_f(x) = \eta/\kappa(x)$. The wave equation in $p$ can be obtained by combining Eqs. (6) and (7).

The second order differential operator of the wave equation in a homogeneous fluid can be factorized and this yields a system of two coupled first order differential equations. This is the wave splitting description, which was

<figure><img src="image_3.png" /><figcaption>FIG. 2. Signal (a) reflected and (b) transmitted by a two-layer porous system. The incidence is normal.</figcaption></figure>

mainly used in scattering problems in the time domain in electromagnetism¹⁵ and then adapted to the frequency domain.² An inhomogeneous porous slab on which impinges an incident wave is shown in Fig. 1. Applied to the wave equation in the ambient fluid, the so-called “vacuum wave splitting transformation”² is $p^\pm=[p\pm Z_0\phi(x)\dot{\mathbf{U}}\cdot\mathbf{n}]/2$, where $Z_0=\rho_f c_0$ is the characteristic impedance of the fluid surrounding the slab [used instead of $Z_e(x, \omega) = \sqrt{\rho_e(x, \omega)K_e(x, \omega)}$ of the fluid in the slab] and $\mathbf{n}$ the unit normal vector (Fig. 1). A system of linear first order coupled differential equations is obtained from Eqs. (6) and (7):

$$ \partial_x p^+ = A^+(x, \omega)p^+ + A^-(x, \omega)p^-, $$

$$ \partial_x p^- = -A^-(x, \omega)p^+ - A^+(x, \omega)p^-, \quad (10) $$

with $A^\pm(x, \omega)=j\omega/2[(Z_0/K_e(x, \omega))\pm(\rho_e(x, \omega)/Z_0)]$.

The computation principle is the following: $p^\pm$ are first calculated in the surrounding homogeneous fluid at $x=L$. An infinitely thin homogeneous layer of thickness $dx$ is inserted at $x=L-dx$ with the corresponding values of $\rho_e$ and $K_e$. At $x=L$ a new set of $p^\pm$ is determined with the help of Eq. (10). A new layer is added at $x=L-2dx$. Using the updated values of $p^\pm$ at $x=L$, the pressure subfields $p^\pm(L-dx, \omega)$ are calculated. The operation is repeated until the last layer is added at $x=0$. For each new layer, the continuity conditions on $p$ and $\phi(x)\dot{\mathbf{U}}\cdot\mathbf{n}$ are implicitly accounted for on both sides of the cumulated slab. The initialization of the procedure requires that $p^\pm$ must be determined at $x=L$. To avoid this calculation, Green's function approach²,¹⁵ is used. Two Green's functions $G^\pm$ are defined by $\forall x \in [0, L]$, $p^\pm(x, \omega)=G^\pm(x, \omega)p^\pm(L, \omega)$. Green's functions are characteristic of the sole material properties and describe the internal field. The boundary conditions at $x=L$ are known and are $G^+(L, \omega)=1$ and $G^-(L, \omega)=0$. The system of coupled first order linear differential equations in $G^\pm$ obtained by inserting Green's functions in Eq.

<table><caption>TABLE I. Properties of the two-layer medium studied.</caption><thead><tr><th></th><th>φ</th><th>τ<sub>∞</sub></th><th>Λ (μm)</th><th>Λ' (μm)</th><th>R<sub>f</sub> (Ns m<sup>-4</sup>)</th><th>Thickness (mm)</th></tr></thead><tbody><tr><td>Layer 1</td><td>0.96</td><td>1.07</td><td>273</td><td>672</td><td>2843</td><td>7.1</td></tr><tr><td>Layer 2</td><td>0.99</td><td>1.001</td><td>230</td><td>250</td><td>12 000</td><td>10.0</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Normal incidence reflection and transmission coefficients R(ω) and T(ω).</figcaption></figure>

(10) can be solved numerically using a Runge-Kutta routine. The reflection and transmission coefficients R(ω) and T(ω) are deduced from p^±,

$$p^{-}(0, \omega) = R(\omega)p^{+}(0, \omega), \qquad (11)$$

$$p^{+}(L, \omega) = T(\omega)p^{+}(0, \omega). \qquad (12)$$

In the numerical simulations of the reflected and transmitted waves, the physical properties are multiplied by a function I(x) to create the change of their values with depth. This function can be called “inhomogeneity function” and is given by $I(x)=1+C(1-\text{erf}(-(x-x_0)/r))$, where C is a constant (different for each property modeled), erf the error function, $x_0$ the position of the jump, and r a steepness factor. The steeper is the jump the finer the stepping dx must be for better accuracy. 400 points were chosen to discretize the total slab and $dx=17.1/400=0.0428$ mm. The value chosen for r was $r=0.1dx$. Smaller values of r had little effect on the computed results. Smoothing the jump by taking $r=10dx$ resulted in an important reduction of the signal reflected at the interface between the two layers.

The experimental principle and the inhomogeneity function are shown in Fig. 1, where an airborne ultrasonic wave is generated and detected by specially designed (ULTRAN) transducers in a frequency range between 150 and 250 kHz. The incident wave is partly reflected, partly transmitted, and partly absorbed by layers of highly porous polyurethane foams put in contact, not glued. The physical parameters (Table I) of each layer were measured<fnref n="16" /> and, within the bounds of the experimental error, those of layer 2 were adjusted to best fit the results for the layered system. The ther-

mal length can be interpreted as a measure of the average pore size (although the “pore” is not straightforwardly defined) while the viscous characteristic length would correspond to an average size of the “constrictions” in the porous medium. The two materials considered are structurally different as one has a fairly high thermal length, resulting in a fairly low flow resistivity. The reflection experiment was carried out with a single transducer used in the pulse-echo mode.<fnref n="17" /> Another transducer was required on the other side of the specimen to measure the transmission. A particularly good agreement between the experimental and the simulated results is found for the reflected and transmitted signals [Figs. 2(a) and 2(b)]. The waveforms were calculated from the experimental incident wave form, which was recorded and introduced into the simulation routines. The agreement is also very good for the reflection and transmission coefficients (Fig. 3) in the frequency range of 170–230 kHz. In all figures, the WS-GF curves cannot be distinguished from the TMM curves. The discrepancies below 170 and above 230 kHz in Fig. 3 are attributed to the limited bandwidth of the incident signal, making the experimental results more sensitive to noise outside the useful bandwidth.

The authors are grateful to O. Matsuda for his useful comments on earlier versions of this letter.

<footnote n="1">¹O. Matsuda and O. B. Wright, J. Opt. Soc. Am. B **19**, 3028 (2002).</footnote>

<footnote n="2">²O. Forslund and S. He, Prog. Electromagn. Res. **19**, 147 (1998).</footnote>

<footnote n="3">³J. L. Buchanan, R. P. Gilbert, A. Wirgin, and Y. S. Xu, *Marine Acoustics: Direct and Inverse Problems* (SIAM, Philadelphia, 2004), Vol. 3, Section 1.5.</footnote>

<footnote n="4">⁴D. J. Van Manen, J. O. A. Robertsson, and A. Curtis, Phys. Rev. Lett. **94**, 164301 (2005).</footnote>

<footnote n="5">⁵J. G. Berryman and R. R. Greene, Geophysics **45**, 213 (1980).</footnote>

<footnote n="6">⁶V. Tournat, V. Zaitsev, V. Gusev, V. Nazarov, P. Béquin, and B. Castagnède, Phys. Rev. Lett. **92**, 085502 (2004).</footnote>

<footnote n="7">⁷M. A. Biot, J. Acoust. Soc. Am. **28**, 168 (1956).</footnote>

<footnote n="8">⁸K. Attenborough, Phys. Rep. **82**, 181 (1982).</footnote>

<footnote n="9">⁹D. L. Johnson, J. Koplik, and R. Dashen, J. Fluid Mech. **176**, 379 (1987).</footnote>

<footnote n="10">¹⁰J. F. Allard, *Propagation of Sound in Porous Media: Modeling Sound Absorbing Materials* (Chapman and Hall, London, 1993), pp. 79 and 145.</footnote>

<footnote n="11">¹¹M. A. Biot, J. Appl. Phys. **33**, 1482 (1962).</footnote>

<footnote n="12">¹²J. M. Carcione, *Wavefield in Real Media: Wave Propagation in Anisotropic, Anelastic and Porous Media* (Pergamon, Amsterdam, 2001), Vol. 31, p. 261.</footnote>

<footnote n="13">¹³M. A. Biot and D. G. Willis, J. Appl. Mech. **24**, 594 (1957).</footnote>

<footnote n="14">¹⁴Y. Champoux and J. F. Allard, J. Appl. Phys. **70**, 1975 (1991).</footnote>

<footnote n="15">¹⁵S. He, J. Math. Phys. **34**, 4628 (1993).</footnote>

<footnote n="16">¹⁶P. Leclaire, L. Kelders, W. Lauriks, N. R. Brown, M. Melon, and B. Castagnède, J. Appl. Phys. **80**, 2009 (1996).</footnote>

<footnote n="17">¹⁷Z. E. A. Fellah, S. Berger, W. Lauriks, C. Depollier, P. Trompette, and J. Y. Chapelon, J. Appl. Phys. **93**, 9352 (2003).</footnote>