# Propagation of acoustic waves in a one-dimensional macroscopically inhomogeneous poroelastic material

Guillaume Gautier, Luc Kelders, Jean-Philippe Groby, Olivier Dazel, Laurent de Ryck, Philippe Leclaire

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

## To cite this version:

Guillaume Gautier, Luc Kelders, Jean-Philippe Groby, Olivier Dazel, Laurent de Ryck, et al. Propagation of acoustic waves in a one-dimensional macroscopically inhomogeneous poroelastic material. Journal of the Acoustical Society of America, 2011, 130 (3), pp.1390-1398. 10.1121/1.3605530. hal-01323683

HAL Id: hal-01323683

https://hal.science/hal-01323683

Submitted on 31 May 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# Propagation of acoustic waves in a one-dimensional macroscopically inhomogeneous poroelastic material

G. Gautier and L. Kelders

Laboratorium voor Akoestiek en Thermische Fysica, KULeuven, Celestijnenlaan 200D, B-3001 Leuven, Belgium

J. P. Groby and O. Dazel

Laboratoire d'Acoustique de l'Université du Maine, UMR6613 CNRS/Université du Maine, F-72085 Le Mans Cedex 9, France

L. De Ryck

LMS international NV, 3001 Leuven, Belgium

P. Leclaire

Laboratoire de Recherche en Mécanique et Acoustique, ISAT/Université de Bourgogne, 49 rue Mademoiselle Bourgeois, B.P. 31-58027 Nevers cedex, France

Wave propagation in macroscopically inhomogeneous porous materials has received much attention in recent years. The wave equation, derived from the alternative formulation of Biot's theory of 1962, was reduced and solved recently in the case of rigid frame inhomogeneous porous materials. This paper focuses on the solution of the full wave equation in which the acoustic and the elastic properties of the poroelastic material vary in one-dimension. The reflection coefficient of a one-dimensional macroscopically inhomogeneous porous material on a rigid backing is obtained numerically using the state vector (or the so-called Stroh) formalism and Peano series. This coefficient can then be used to straightforwardly calculate the scattered field. To validate the method of resolution, results obtained by the present method are compared to those calculated by the classical transfer matrix method at both normal and oblique incidence and to experimental measurements at normal incidence for a known two-layers porous material, considered as a single inhomogeneous layer. Finally, discussion about the absorption coefficient for various inhomogeneity profiles gives further perspectives.

## I. INTRODUCTION

The study of wave propagation in macroscopically inhomogeneous porous media was initially motivated by (1) the design of sound absorbing porous materials with optimal material and geometrical property profiles¹ and (2) the retrieval of the spatially varying material parameters of industrial foams.² These, and other inverse problems, are of great importance in connection with the characterization of the mechanical properties of naturally occurring macroscopically inhomogeneous porous materials, such as bones. The literature on inhomogeneous media is extensive in several fields of physics, from optics and electromagnetism,³,⁴ to acoustics,⁵,⁶ and geophysics⁷ and granular media.⁸ Many natural and man-made materials are porous, and therefore heterogeneous at a microscopic scale. The wave equation in macroscopically inhomogeneous porous media was derived from the alternative formulation of Biot's theory⁹ in Ref. 10 and solved in the case of rigid frame inhomogeneous porous materials via the Wave Splitting method and “transmission” Green's functions approach or via an iterative Born approximation procedure based on the specific Green's function of

the configuration.¹¹ The recovery of several profiles of spatially varying material parameters by means of an optimization approach, was then achieved in Ref. 2 still in the rigid frame approximation.

When saturated by a light fluid such as air, the frame is moving below the solid/fluid decoupling frequency. Moreover, in many applications porous materials are saturated by a heavy fluid such as water or bone marrow. They can also be excited mechanically. In these cases, the rigid frame approximation is not valid and the full macroscopically inhomogeneous poroelastic model should be used and solved.

It is assumed that the wavelengths are larger 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. Macroscopically inhomogeneous poroelastic materials imply that both acoustic and elastic properties are space-dependent at the macroscopic scale.

First, the constitutive linear stress-strain relations and the momentum conservation law in the absence of body forces are recalled for an inhomogeneous poroelastic material. These equations are then solved for a one-dimensional macroscopically inhomogeneous poroelastic material via the
state vector formalism or the so-called Stroh formalism¹² together with Peano series.¹³,¹⁴ The Stroh formalism is largely used to model acoustic wave propagation in stratified anisotropic elastic materials.¹⁵,¹⁶ Similar methods (transmission matrix method) are used in electromagnetism to model the propagation of electromagnetic waves in anisotropic or gyrotropic stratified materials.¹⁷,¹⁸ In the frequency domain, it consists in the rewriting of the constitutive linear stress-strain relations and the momentum conservation law in terms of the state vector, whose components can be chosen arbitrarily as long as they are continuous along the inhomogeneity of the material. This leads after spatial Fourier transform to a first-order ordinary differential system of equations whose unknown is the state vector. The spatial dependence of the materials properties are accounted for through the spatially dependent matrix components. The solution of the system appears in terms of Peano series, which are well fitted for wave propagation problems in functionally graded materials.¹⁹ These series apply to continuously varying poroelastic properties and avoids problems related to a lack of discretization when the materials are approximated by stratified ones. Numerical results obtained with this method are compared to calculations of the classical transfer matrix method performed with the validated MAINE3A code²⁰ for a known two-layers porous material, considered as a single inhomogeneous layer. The transfer matrix method is particularly suitable to solve problems involving a layered configuration.

Numerical results are also compared to experimental measurements at normal incidence. The experiment consists in recording waves reflected by the chosen two-layered porous medium laid on the floor of a semi-anechoic room when excited at normal incidence by a dipolar source.

Finally, applications in material design for engineering applications are treated, by comparing the absorption coefficient of a macroscopically inhomogeneous porous plate with various inhomogeneity profiles that are either continuous or discontinuous.

## II. EQUATIONS OF MACROSCOPICALLY INHOMOGENEOUS POROUS MATERIALS

As pointed out by several authors,⁹,²¹,²² the generalized formulation of the Biot theory⁹ is suitable to macroscopically inhomogeneous porous media and also to take into account anisotropy and viscoelastic frames. Recently, another formulation was proposed in Ref. 23 that is also suitable to macroscopically inhomogeneous porous media. The article focuses on the alternative Biot's formulation, which is largely employed in acoustics and geophysics.

Rather than dealing directly with the arbitrary field $\bar{s}(x, t)$ [with $x = (x_1, x_2)$], we prefer to deal with the $s(x, \omega)$, related to $\bar{s}(x, t)$ by the Fourier transform $\bar{s}(x, t) = \int_{-\infty}^{\infty} s(x, \omega) e^{i\omega t} d\omega$, wherein $\omega = 2\pi\nu$ is the angular frequency, with $\nu$ the frequency. Henceforth, we drop the $\omega$ in $s(x, \omega)$ so that it is written $s(x)$.

From the alternative formulation of Biot,⁹ the stress strain relations in an initially stress free, statistically isotropic poroelastic material take the form

$$ \begin{cases} \sigma_{ij} = 2N\epsilon_{ij} + (\lambda_c\theta - \alpha M\zeta)\delta_{ij}, \\ p = M(\zeta - \alpha\theta), \end{cases} \quad (1) $$

where in $\delta_{ij}$ denotes the Kronecker symbol. The components of the total stress tensor are $\sigma_{ij}$, the fluid pressure in the pores is $p$, and the components of the strain tensor are $\epsilon_{ij} = 1/2(u_{ij} + u_{ji})$, with the solid displacement $\mathbf{u}$, $\theta = u_{i,i}$, and $\zeta = -w_{i,i}$ with the fluid/solid relative displacement $\mathbf{w} = \phi(\mathbf{U} - \mathbf{u})$ ($\mathbf{U}$ being the fluid displacement and $\phi$ the porosity). The Einstein summation notation is implicit in the expressions of $\theta$ and $\zeta$. The material properties²⁴ are the bulk modulus of the closed porosity system, i.e., in which the pore volume is sealed, $\lambda_c = \lambda + \alpha^2 M$, the Lamé coefficients of the elastic frame $\lambda$ and $N$, an additional elastic parameter $M$, and an elastic coupling coefficient $\alpha$.

These mechanical coefficients are related to the more commonly used in acoustics $P$, $Q$, and $R$ coefficients from the original formulation²⁵ through

$$ \begin{aligned} \alpha &= \phi(Q + R)/R, \quad M = R/\phi^2, \\ \lambda_c &= P - Q^2/R - 2N + \alpha^2 M. \end{aligned} \quad (2) $$

The expressions of $P$, $Q$, and $R$, in case of air saturated materials, reduce to²⁶

$$ \begin{aligned} R &= \phi \tilde{K}_f, \quad Q = \tilde{K}_f(1 - \phi), \\ P &= K_b + 4N/3 + (1 - \phi)^2 \tilde{K}_f / \phi, \end{aligned} \quad (3) $$

wherein the bulk modulus of the skeleton is $K_b = 2N(1+\nu)/3(1-2\nu)$, the Poisson coefficient (of the skeleton) is $\nu$, and

$$ \tilde{K}_f = \frac{\gamma P_0}{\gamma - (\gamma - 1) \left[ 1 + i \frac{\omega_c}{\omega \text{Pr}} G(\omega \text{Pr}) \right]^{-1}} \quad (4) $$

The correction function $G(\omega \text{Pr})$, introduced in Ref. 27 to account for the thermal losses, is

$$ G(\omega \text{Pr}) = \sqrt{1 - i\eta\rho_f\omega\text{Pr}\left(\frac{2\tau_\infty}{R'_c\phi\Lambda'}\right)^2}, \quad (5) $$

with $\omega_c' = R'_c\phi/\rho_f\tau_\infty$, $\gamma$ the specific heat ratio, Pr the Prandtl number, $\tau_\infty$ the tortuosity, $R'_c$ the “thermal resistivity,” and $\Lambda'$ the thermal characteristic length. The thermal resistivity is related to the thermal characteristic length²⁷ through $R'_c = 8\tau_\infty\eta/\phi\Lambda'^2$.

In the absence of body forces, the conservation of momentum and the generalized Darcy's law lead to the following equations in the frequency domain:

$$ \begin{cases} \omega^2 \rho_f w_i + \omega^2 \rho u_i = -\sigma_{ij,j}, \\ \omega^2 \rho_f u_i + \omega^2 \tilde{\rho}_{eq} w_i = p_{,i}, \end{cases} \quad (6) $$

wherein the bulk density of the porous medium is $\rho$, such that $\rho = (1 - \phi)\rho_s + \phi\rho_f$ with $\rho_s$ the density of the solid
and $\rho_f$ the density of the saturating fluid, and the mass pa- rameter $\tilde{\rho}_{eq}$ (Refs. 9 and 28) is

$$ \tilde{\rho}_{\text{eq}} = \frac{\rho_f \tau_{\infty}}{\phi} \left[ 1 + i \frac{\omega_c}{\omega} F(\omega) \right]. \quad (7) $$

The correction function $F(\omega)$, introduced in Ref. 28 and which accounts for the viscous losses, is given by

$$ F(\omega) = \sqrt{1 - i\eta\rho_f\omega \left(\frac{2\tau_{\infty}}{R_f\phi\Lambda}\right)^2}, \quad (8) $$

with the Biot frequency $\omega_c = R_f\phi/\rho_f\tau_\infty$, the flow resistivity $R_f$ and the viscous characteristic length $\Lambda$. The mass parame- ter $\tilde{\rho}_{eq}$ is the complex frequency dependent equivalent density used in the rigid frame approximation, while the complex fre- quency dependent equivalent bulk modulus $\tilde{K}_{eq}$ used in this approximation is related to $\tilde{K}_f$ through $\tilde{K}_f = \phi\tilde{K}_{eq}$. For most foams, the elastic coupling coefficient $\alpha$ reduces to 1, $M = \tilde{K}_{eq}$ and $\lambda_c$ reduces to $\lambda_c = K_b - 2/3N + \alpha^2M$.

In the previous equations, $N$, $\lambda$ (and $\nu$), $\lambda_c$, $\alpha$, $M$, $\phi$, $\tau_\infty$, $\Lambda$, $\Lambda'$, $R'_t$, $R_f$ are x dependent.

### III. NUMERICAL EVALUATION OF THE PRESSURE FIELD

#### A. Description of the configuration

Both the incident plane acoustic wave and the plate are assumed to be invariant with respect to the Cartesian coordi- nate $x_3$. A cross-sectional $x_1-x_2$ plane view of the 2D scatter- ing problem is shown in Fig. 1.

The upper and lower boundaries of the layer are flat and parallel. They are designated by $\Gamma_L$ and $\Gamma_0$ and their $x_2$ coor- dinates are $L$ and $0$. The porous material $M^{[1]}$ occupies the domain $\Omega^{[1]}$. The inhomogeneity of the plate occurs along the $x_2$ direction, i.e., $N$, $\lambda$ (and $\nu$), $\lambda_c$, $\alpha$, $M$, $\phi$, $\tau_\infty$, $\Lambda$, $\Lambda'$, $R'_t$, $R_f$ are $x_2$-dependent. An isotropic macroscopically inhomoso- geneous porous material is considered. The material can be viewed as a functionally graded material. The surrounding and saturating fluid is the air medium (density $\rho_f=1.213$ kg m$^{-3}$, atmospheric pressure $P_0=1.01325 \times 10^5$ Pa, and vis- cosity $\eta=1.839 \times 10^5$ kg m$^{-1}$ s$^{-1}$). The inhomogeneous po- rous layer is backed by a rigid plate at $\Gamma_0$.

We denote the total pressure, wavevector and density, respectively, with $p^{[0]}$, $k^{[0]}$ and $\rho^{[0]}$ in $\Omega^{[0]}$ and the total stress tensor, the fluid pressure in the pores, the solid displacement, and solid/fluid relative displacement respectively with $\sigma^{[1]}$, $p^{[1]}$, $\mathbf{u}^{[1]}$ and $\mathbf{w}^{[1]}$ in $\Omega^{[1]}$.

The wavevector $\mathbf{k}^i$ of the incident plane wave lies in the sagittal plane and the angle of incidence is $\theta^i$ measured counterclockwise from the positive $x_1$ axis. The incident wave, initially propagating in $\Omega^{[0]}$, is expressed by $p^i(\mathbf{x}) = A^i e^{i(k_1 x_1 - k_2^0 (x_2 - L))}$, wherein $k_1 = -k^{[0]} \cos \theta^i$, $k_2^0 = k^{[0]} \sin \theta^i$ and $A^i = A^i(\omega)$ is the signal spectrum.

The uniqueness of the solution to the forward-scattering problem is ensured by the radiation condition:

$$ p^{[0]}(x) - p^i(x) \sim \text{outgoing waves; } |x| \to \infty, x_2 > L. \quad (9) $$

#### B. The state vector formalism and Peano series

The spatial Fourier transform of the Eqs. (1) and (6) is first performed. The transform $\hat{s}(x_2, k_1)$ of $s(\mathbf{x})$ is also introduced and can be written in the form $s(\mathbf{x}) = \hat{s}(x_2, k_1)e^{ik_1 x_1}$ because of the plane wave nature of the excitation. Henceforth, we drop the $k_1$ in $\hat{s}(x_2, k_1)$ so that it is written $\hat{s}(x_2)$. The geome- try of the configuration being planar and infinite along the $x_2$-axis, the fields components that are continuous along the inhomogeneity and the boundary conditions apply either to $s(\mathbf{x})$ or to $\hat{s}(x_2)$.

Inside the domain $\Omega^{[1]}$ the normal components of the total stress tensor, $\sigma_{12}^{[1]}$ and $\sigma_{22}^{[1]}$, the pressure $p^{[1]}$, the solid displacements, $u_1^{[1]}$ and $u_2^{[1]}$, and the normal component of the solid/fluid relative displacement $w_2^{[1]}$ are continuous along the $x_2$-axis. It seems natural to choose these 6 parame- ters as components of the state vector. Nevertheless, it seems better to adapt these components to the considered boundary problem. On one hand, at the interface $\Gamma_L$, the normal components of the stress tensor ($\sigma_{12}^{[1]} = 0$, $\sigma_{22}^{[1]} = -p^{[0]}$), the pressure ($p^{[1]}=p^{[0]}$), and normal component of the velocity $-i\omega(w_2^{[1]} + u_2^{[1]}) = V_2^{[0]}$ with $V_2^{[0]}$ the normal component of the velocity in $\Omega^{[0]}$ are continu- ous. On the other hand, the solid displacement and the normal component of the solid/fluid relative velocity vanish at the interface $\Gamma_0$, i.e., $u_1^{[1]} = 0, u_2^{[1]} = 0$ and $-i\omega(w_2^{[1]} + u_2^{[1]}) = 0$. Because $u_2^{[1]}$ and $w_2^{[1]}$ are continuous at any $x_2$ in $\Omega^{[1]}$, the normal component of the velocity $-i\omega(w_2^{[1]} + u_2^{[1]})$ is also continuous. The normal component of the solid/fluid relative displacement $w_2$ is also replaced by $V_2^{[0]} = -i\omega(w_2^{[1]} + u_2^{[1]})$ in Eqs. (1) and (6). This parame- ter is preferred to $w_2^{[1]}$ to define the wave vector.

After spatial transform, the Eqs. (1) and (6) split into two systems of equations (see the Appendix): one set of six first order differential equations only depending on the components of the column state vector $\hat{W} = \langle \hat{\sigma}_{12}^{[1]}, \hat{\sigma}_{22}^{[1]}, \hat{p}^{[1]} \rangle$, $\hat{u}_1^{[1]}, \hat{u}_2^{[1]}, \hat{V}_2^{[1]} \rangle$, and one set of two equations that linked the two last unknowns, $\hat{\sigma}_{11}^{[1]}$ and $\hat{w}_{11}^{[1]}$, to the components of the state vector. The problem also reduces to the solution of the first order differential matrix system composed of the first six first order differential equations:

$$ \frac{\partial}{\partial x_2} \hat{W} - A \hat{W} = 0, \quad (10) $$

with $A = -B^{-1}D$, wherein

$$ B = \begin{bmatrix} 1 & 0 & 0 & 0 & ik_1(\lambda_c - \alpha M) & -\frac{k_1 \alpha M}{\omega} \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & N & 0 & 0 \\ 0 & 0 & 0 & 0 & \lambda_c + 2N - \alpha M & \frac{\alpha M}{\omega} \\ 0 & 0 & 0 & 0 & (\alpha - 1)M & \frac{M}{\omega} \end{bmatrix}, \quad (11) $$
and

$$
D = \begin{bmatrix}
0 & 0 & ik_1 \frac{\rho_f}{\tilde{\rho}_{eq}} - \frac{k_1^2 \alpha M}{\tilde{\rho}_{eq} \omega^2} & \rho \omega^2 - \frac{\rho_f^2}{\tilde{\rho}_{eq}} \omega^2 - k_1^2 \lambda_c + 2N - \frac{\rho_f}{\tilde{\rho}_{eq}} \alpha M & 0 & 0 \\
ik_1 & 0 & 0 & 0 & (\rho - \rho_f) \omega^2 & i \rho_f \omega \\
0 & 0 & 0 & 0 & (\rho_f - \tilde{\rho}_{eq}) \omega^2 & i \tilde{\rho}_{eq} \omega \\
-1 & 0 & 0 & 0 & ik_1 N & 0 \\
0 & -1 & -\frac{k_1^2 \alpha M}{\tilde{\rho}_{eq} \omega^2} & ik_1 \lambda_c - \frac{\rho_f}{\tilde{\rho}_{eq}} \alpha M & 0 & 0 \\
0 & 0 & 1 - \frac{k_1^2 M}{\tilde{\rho}_{eq} \omega^2} & ik_1 M & 0 & 0
\end{bmatrix} . \quad (12)
$$

The matrices **B** and **D** are $x_2$-dependent, because $\lambda_c(x_2)$, $\alpha(x_2)$, $M(x_2)$, $N(x_2)$, $\tilde{\rho}_{eq}(x_2)$, and $\rho(x_2)$ are $x_2$-dependent. This dependence has not been detailed in (11) and (12) for conciseness.

The solution of system (10) takes the form

$$
\hat{W}(L) = M\hat{W}(0), \quad (13)
$$

wherein **M** is the so-called matricant,¹⁴ which relates the value of the state vector $\hat{W}(0)$, at $x_2=0$, to the value of the state vector $\hat{W}(L)$, at $x_2=L$. Since A is $x_2$ dependent (i.e., the plate is not homogeneous or piecewise constant) and A does not commute for different values of $x_2$, i.e., $[A(x_2), A(x'_2)] = A(x_2)A(x'_2) - A(x'_2)A(x_2) \neq 0$, $\forall (x_2, x'_2) \in [0, L]^2$, $x_2 \neq x'_2$, the matricant **M** does not contain matrix exponentials or multiplications of matrix exponentials. The matricant is rather defined by the so-called multiplicative integral satisfying the Peano expansion.¹³,¹⁴,¹⁹,²⁹ This avoids any problem related to lack of discretization when the continuously varying material is approximated by a piecewise constant material. The Peano series reads as

$$
M = I + \int_{0}^{L} A(x_2) dx_2 + \int_{0}^{L} A(x_2) \left[ \int_{0}^{x_2} A(\zeta) d\zeta \right] dx_2 + \dots \quad (14)
$$

and the evaluation of **M** is performed via the iterative
scheme

<figure><img src="image_3.png" /><figcaption>FIG. 1. Cross-sectional plane view of the configuration.</figcaption></figure>

$$
\left\{
\begin{aligned}
\mathbf{M}^{\{0\}} &= \mathbf{I}, \\
\mathbf{M}^{\{n\}} &= \mathbf{I} + \int_{0}^{L} A(x_2) \mathbf{M}^{\{n-1\}}(x_2) dx_2,
\end{aligned}
\right.
\qquad (15)
$$

such that $\lim_{n \to \infty} \mathbf{M}^n = \mathbf{M}$.

**C. The boundary problem**

The application of the boundary conditions at both inter-
faces $\Gamma_L$ [$\hat{\sigma}_{12}^{[1]} = 0$, $\hat{\sigma}_{22}^{[1]} = -\hat{p}^{[0]}$, $\hat{p}^{[1]} = \hat{p}^{[0]}$, $\hat{V}_2^{[1]} = \hat{V}_2^{[0]} = -i/\omega \rho \frac{\partial \hat{p}^{[0]}}{\partial x_2}$] and $\Gamma_0$ ($\hat{u}_1 = 0, \hat{u}_2 = 0, \hat{V}_2^{[1]} = 0$) yields the state
vectors $\hat{W}(L)$ and $\hat{W}(0)$, which are required to solve the
problem. From the separation of variables, the radiation con-
ditions, and the spatial Fourier transform, the pressure field
in $\Omega^{[0]}$ can be written as

$$
\hat{p}^{[0]} = A^i e^{-ik_2^{[0]}(x_2-L)} + A^i \mathrm{Re} i k_2^{[0]}(x_2-L), \quad (16)
$$

where *R* is the reflection coefficient. The state vectors become

$$
\hat{\mathbf{W}}(L) = \begin{bmatrix} \hat{\sigma}_{12}^{[1]}(L) \\ \hat{\sigma}_{22}^{[1]}(L) \\ \hat{p}^{[1]}(L) \\ \hat{u}_1^{[1]}(L) \\ \hat{u}_2^{[1]}(L) \\ \hat{u}_2^{[1]}(L) \end{bmatrix} = \begin{bmatrix} 0 \\ -\hat{p}^{[0]}(L) \\ \hat{p}^{[0]}(L) \\ \hat{u}_1^{[1]}(L) \\ \hat{u}_2^{[1]}(L) \\ -\frac{i}{\omega \rho} \frac{\partial \hat{p}^{[0]}}{\partial x_2} \Big|_{x_2=L} \end{bmatrix} = \mathbf{S} + \mathbf{L}^L \cdot \begin{bmatrix} A^i R \\ \hat{u}_1^{[1]}(L) \\ \hat{u}_2^{[1]}(L) \end{bmatrix}, \quad (17)
$$

and

$$
\hat{\mathbf{W}}(0) = \begin{bmatrix} \hat{\sigma}_{12}^{[1]}(0) \\ \hat{\sigma}_{22}^{[1]}(0) \\ \hat{p}^{[1]}(0) \\ \hat{u}_1^{[1]}(0) \\ \hat{u}_2^{[1]}(0) \\ \hat{V}_2^{[1]}(0) \end{bmatrix} = \begin{bmatrix} \hat{\sigma}_{12}^{[1]}(0) \\ \hat{\sigma}_{22}^{[1]}(0) \\ \hat{p}^{[1]}(0) \\ 0 \\ 0 \\ 0 \end{bmatrix} = \mathbf{L}^0 \cdot \begin{bmatrix} \hat{\sigma}_{12}^{[1]}(0) \\ \hat{\sigma}_{22}^{[1]}(0) \\ \hat{p}^{[1]}(0) \end{bmatrix}, \quad (18)
$$
wherein S accounts for the excitation of the system by the plane incident wave, $\mathbf{L}^L$ relates the unknowns $R, \hat{u}_1^{[1]}(L)$ and $\hat{u}_2^{[1]}(L)$ to the state vector $\hat{\mathbf{W}}(L)$, and $\mathbf{L}^0$ relates the unknowns $\hat{\sigma}_{12}^{[1]}(0), \hat{\sigma}_{22}^{[1]}(0)$ and $\hat{p}^{[1]}(0)$ to the state vector $\hat{\mathbf{W}}(0)$. Their expressions are

$$ \mathbf{S} = \begin{bmatrix} 0 \\ -A^i \\ A^i \\ 0 \\ 0 \\ -A^i k_2^{[0]} / \omega \rho^{[0]} \end{bmatrix}, \quad \mathbf{L}^0 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad \text{and} \\ \mathbf{L}^L = \begin{bmatrix} 0 \\ -1 \\ 1 \\ 0 \\ 0 \\ k_2^{[0]} / \omega \rho^{[0]} \end{bmatrix} $$

Finally, the introduction of (17) and (18) in (13) gives rise to the final system of equations, whose solution contain the reflection coefficient $R$

$$ [[\mathbf{L}^L] [-\mathbf{M} \cdot \mathbf{L}^0]] \cdot \begin{bmatrix} A^i R \\ \hat{u}_1^{[1]}(L) \\ \hat{u}_2^{[1]}(L) \\ \hat{\sigma}_{12}^{[1]}(0) \\ \hat{\sigma}_{22}^{[1]}(0) \\ \hat{p}^{[1]}(0) \end{bmatrix} = -\mathbf{S}. \quad (20) $$

This system is solved for each frequency and directly provide the reflection coefficient associated with the plane incident wave. The pressure field in $\Omega^{[0]}$ can then be calculated through $p^{[0]} = \hat{p}^{[0]} e^{ik_1 x_1}$. A similar system can be cast when the inhomogeneous porous plate is not backed and radiate in the air medium. In this case, a transmitted field is required and the transmission coefficient is also calculated.

## IV. VALIDATION ON A MULTILAYERED POROUS MEDIUM

In order to validate the present method, calculations are performed for a known two-layers poroelastic medium configuration considered as a single poroelastic plate. Each layer of the plate is a porous foam saturated by air. The characteristic properties of each layer have been determined by classical methods$^{26}$ and are given in Table I. The layer 1 is a Eurocell foam while layer 2 is a Fireflex (Recticel, Wetteren/East-Flanders, Belgium) foam.

The choice of this configuration is motivated by the fact that the results of the present method can be compared with

known results from the classical Transfer Matrix Method (TMM) provided by MAINE3A$^{20}$ and with experimental measurements.

To consider the slab as a single inhomogeneous material, the jump discontinuities in the two-layered system are smoothed by using the following analytical continuous and continuously differentiable function:

$$ I(x_2) = 1 + \frac{C}{2} \left[ 1 + \operatorname{erf} \left( \frac{x_2 - x_2^0}{r} \right) \right], \quad (21) $$

wherein $C$ is the step value, which is different for each parameter in Table I, erf is the error function, $x_2^0$ is the position of the jump, and $r$ corresponds to the steepness of the continuous jump such that the smaller $r$ is, the steeper is the jump.

The number of iterations required for the correct evaluation of the matricant increases with frequency. A change of variable like that proposed in Refs. 30 and 31 seems inaccurate, because the characteristic parameters are frequency dependent for porous materials. Moreover, a large number of iterations is required for the correct evaluation, because poroelastic materials are highly dissipative.

The configuration is discretized in with 1000 points and the number of iterations for the evaluation of the Peano series is 250. The jump position and its slope are fixed to $x_2^0 = 50$ mm and $r = 10^{-6}$, while the total thickness is $L = 69.8$ mm.

### A. Numerical validation

The results calculated with the present method are first compared with those from the classical TMM provided by MAINE3A at normal and oblique incidences. Figure 2 depicts the real and imaginary parts of the reflection coefficient at normal and oblique incidences by the two-layers medium studied. The curves cannot be distinguished one from another.

This provides a validation of the model and of the solving method.

The real and imaginary part of the reflection coefficient calculated with MAINE3A under the rigid frame approximation (○) and with the Limp model$^{32}$ (□) are also plotted in Fig. 2. Below 4000 Hz the rigid frame approximation is not valid at the Biot resonances. This test emphasizes the added value of the new method modeling the wave propagation in macroscopically inhomogeneous poroelastic foams at low frequencies, and validates the method in the high frequency range. The peaks at low frequencies are also attributed to Biot and not to Limp effects. The method is also stable and robust.

### B. Experimental validation

The principle of the experiment is shown in Fig. 3. A 1.58 m × 1.06 m plate of Eurocell of 5 cm thick is laid on the floor of the semi-anechoic room available in ATF,

<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> (N s m<sup>-4</sup>)</th><th>ρ (kg m<sup>-3</sup>)</th><th>N (Pa)</th><th>ν</th><th>Thickness (mm)</th></tr></thead><tbody><tr><td>Layer 1</td><td>0.95</td><td>1.1</td><td>15</td><td>45</td><td>42000</td><td>126</td><td>280000 (1-0.05i)</td><td>0.24</td><td>50</td></tr><tr><td>Layer 2</td><td>0.96</td><td>2.2</td><td>110</td><td>352</td><td>9000</td><td>37</td><td>60000 (1-0.05i)</td><td>0.21</td><td>19.8</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Real and imaginary part of the reflection coefficient at normal incidence (a) and (b) and at $\theta^i = \pi/6$ (c) and (d) calculated with the present method (—), with MAINE3A (—), with MAINE3A under the approximation of the equivalent fluid model (○), and with MAINE3A with the Limp model (□).</figcaption></figure>

KULeuven. A 1.5 m × 2 m plate of Fireflex of 1.98 cm thick was then laid above of the Eurocell foam. The two foams are in contact and not glued. The floor is supposed to be rigid.

A dipolar source, radiating a sweep from 100 Hz to 5000 Hz, is placed at the center of the Eurocell foam at $x_2^s = 83$ cm above the ground. A microphone is placed successively at $x_2^{(1)} = 8.5$ cm and $x_2^{(2)} = 9.5$ cm above the ground at the center of the source. The source center was determined by recording the sweep signal with the microphone at various positions along the $x_1$ and $x_3$ axis. The source center corresponds to the maximum of amplitude recorded by the microphone. This guarantees that the recordings are performed at normal incidence. The $x_2$ position of both the source and of the microphone ensure that, locally, the field can be approximated by a plane wave. The recorded signals are averaged over 30 waveforms. The pressures recorded at both positions take the form $p^{(j)} = A^i e^{-ik^{[0]}x_2^{(j)}} + A^i \text{Re}i^{[k^{[0]}]x_2^{(j)}}$, $j=1, 2$ at normal incidence. It directly equals

<figure><img src="image_5.png" /><figcaption>FIG. 3. Experimental setup.</figcaption></figure>

the spatial Fourier transform $\hat{p}^{(j)} = p^{(j)}$, because the pressure is recorded at normal incidence. From these two measurements, the reflection coefficient at normal incidence can be determined through

$$R = \frac{p^{(1)}e^{-ik^{[0]}x_2^{(2)}} - p^{(2)}e^{-ik^{[0]}x_2^{(1)}}}{p^{(2)}e^{ik^{[0]}x_2^{(1)}} - p^{(1)}e^{ik^{[0]}x_2^{(2)}}}. \quad (22)$$

This formula is commonly used when the measuring technique, sometimes referred to as the Tamura method³³ or sometimes referred to as PP-method,³⁴ is used. Figure 4 depicts the absolute value of the experimental reflection coefficient calculated with the present method and measured between 500 Hz to 5000 Hz. This frequency bandwidth is usually considered as being the one over which the reflection methods give accurate results. A fairly good match is found between the curves. In particular, the peak around 800 Hz, which corresponds to a modified Biot resonance of the Eurocell, is well recovered. The discrepancies between the two curves can be explained by the fact that (1) a plane incident wave is assumed, (2) the foams are of finite size, (3) the boundary conditions between the rigid backing and the Eurocell foam does not correspond to perfectly rigid, (4) the two foams are not perfectly in contact, and by (5) the spacing between the two positions of the microphone. In order to emphasize the pertinence of a macroscopically inhomogeneous porous material approach, the absolute value of the reflection coefficient calculated for an homogeneous plate of thickness $L = 6.98$ cm of Eurocell and of Fireflex is also
<figure><img src="image_6.png" /><figcaption>FIG. 4. Comparison between $|R|$ calculated with the present method (—) and recovered with the measurement (○) and $|R|$ for a L-thick homogeneous plate of Fireflex (●●●) and of Eurocell (—).</figcaption></figure>

plotted Fig. 4. Both results are quite far from the experimental curve and from the numerical calculation when the two-layers configuration is modeled as a macroscopically inhomogeneous material.

## V. ASSUMPTION ON OTHER PROPERTY PROFILES

Macroscopically inhomogeneous porous materials offer the possibility of wider applications in sound absorbing material design than their macroscopically homogeneous counterparts. In the following, examples of absorbing material design are presented.

First, the acoustical properties of the previously studied two-layered plate see (Table I) are considered but the jump occurs at $x_2^0 = L/2$, i.e., the thickness of both foams is equal to 3.49 cm. The slope is then varied from $r = 10^{-6}$ to $r = 10^{-2}$. Then, a Hanning profile is considered by using the following analytical continuous and continuously differentiable function:

$$H(x) = 1 + \frac{C}{2} \left[ 1 - \cos\left(\frac{x}{L}\right) \right], \quad (23)$$

wherein C is the relative difference in amplitude between the value at both sides and the value in the middle of the slab.

Figure 5 depicts the profile generated for $R_f$ as an example (the other parameters also vary with an identical profile)

for $r = 10^{-6}$, $r = 10^{-2}$ and a Hanning profile, together with the corresponding absorption coefficient $A = 1 - |R|^2$. The absorption coefficient is one of the most important parameters when designing a foam. These last years, this coefficient has received a large attention, particularly the attempt of increasing its value at low frequencies. Porous foams suffer from a weak absorption at low frequency when compared to their efficiency at higher frequencies. While a modification of the slope only induces a frequency shift at higher frequencies, the use of the Hanning profile enables to increase the absorption coefficient at low frequency. When the foam stack is reversed, the absorption coefficient differs in case of single slope profiles, while it remains the same with symmetrical profiles as Hanning profile.

## VI. CONCLUSION

A model of the acoustic response of macroscopically inhomogeneous elastic frame porous materials derived from the alternative Biot's theory of 1962 was solved. A fast and stable numerical method, derived from the state vector formalism together with Peano series was developed to solve the macroscopically inhomogeneous poroelastic wave equations. To our knowledge, these equations are solved and these methods are derived for the first time for poroelastic materials. A validation of this method was made on the example of a two-layers medium, by comparison to the exact solution obtained by the transfer matrix method (MAINE3A) at both normal and oblique incidence and by comparison to experimental results at normal incidence. In the numerical procedure, the jump of properties between the layers was accounted for in the form of a single continuous function. This result validates the present procedure. Finally, examples of the absorption coefficient are given for various property profiles. These examples illustrate the possibility of designing acoustically absorbing materials by controlling the gradient of parameters. This last point requires further investigation, in particular, how the properties vary when the foam is compressed. Another possible application of the equations and method of solution derived here is the development of an optimized inversion procedure to characterize macroscopically inhomogeneous porous materials.

<figure><img src="image_7.png" /><figcaption>FIG. 5. Example of the a profile considered (a) and the corresponding absorption coefficient (b): slope profiles with $r = 10^{-6}$ (—) and $r = 10^{-2}$ (---) and Hanning profile (●●●).</figcaption></figure>
APPENDIX: DERIVATION OF THE FIRST-ORDER
LINEAR ORDINARY DIFFERENTIAL SYSTEM

The expanded form of the spatial Fourier transform of
(1) and (6) is

$$
\begin{align*}
\hat{\sigma}_{11}^{[1]} &= ik_1(\lambda_c + 2N)\hat{u}_1^{[1]} + \lambda_c \frac{\partial \hat{u}_2^{[1]}}{\partial x_2} + ik_1 \alpha M \hat{w}_1 + \alpha M \frac{\partial \hat{w}_2^{[1]}}{\partial x_2}, \\
\hat{\sigma}_{12}^{[1]} &= ik_1 N \hat{u}_2^{[1]} + N \frac{\partial \hat{u}_1^{[1]}}{\partial x_2}, \\
\hat{\sigma}_{22}^{[1]} &= (\lambda_c + 2N) \frac{\partial \hat{u}_2^{[1]}}{\partial x_2} + ik_1 \lambda_c \hat{u}_1^{[1]} + ik_1 \alpha M \hat{w}_1 + \alpha M \frac{\partial \hat{w}_2^{[1]}}{\partial x_2}, \\
\hat{p}_{[1]} &= -\alpha M ik_1 \hat{u}_1^{[1]} - \alpha M \frac{\partial \hat{u}_2^{[1]}}{\partial x_2} - ik_1 M \hat{w}_1 - M \frac{\partial \hat{w}_2^{[1]}}{\partial x_2}, \\
\rho \omega^2 \hat{u}_1^{[1]} + \rho_f \omega^2 \hat{w}_1 &= -ik_1 \hat{\sigma}_{11}^{[1]} - \frac{\partial \hat{\sigma}_{12}^{[1]}}{\partial x_2}, \\
\rho \omega^2 \hat{u}_2^{[1]} + \rho_f \omega^2 \hat{w}_2^{[1]} &= -ik_1 \hat{\sigma}_{12}^{[1]} - \frac{\partial \hat{\sigma}_{22}^{[1]}}{\partial x_2}, \\
\rho_f \omega^2 \hat{u}_1^{[1]} + \tilde{\rho}_{eq} \omega^2 \hat{w}_1 &= ik_1 \tilde{p}_{[1]}, \\
\rho_f \omega^2 \hat{u}_2^{[1]} + \tilde{\rho}_{eq} \omega^2 \hat{w}_2^{[1]} &= \frac{\partial \tilde{p}_{[1]}}{\partial x_2}. \tag{A1}
\end{align*}
$$

The solid/fluid relative displacement $\hat{w}_2^{[1]}$ is then replaced by
$\hat{V}_2^{[1]} = -i\omega(\hat{w}_2^{[1]} + \hat{u}_2^{[1]})$. From the first and seventh equations
of (A1), $\hat{u}_1^{[1]}$ and $\hat{\sigma}_{11}^{[1]}$ can be expressed in terms of $\hat{p}^{[1]}$, $\hat{u}_1^{[1]}$,
$\hat{u}_2^{[1]}$, $\hat{V}_2^{[1]}$

$$
\begin{align}
\hat{w}_1^{[1]} &= \frac{1}{\tilde{\rho}_{\text{eq}} \omega^2} (ik_1 \hat{p}^{[1]} - \rho_f \omega^2 \hat{u}_1^{[1]}), \hat{\sigma}_{11}^{[1]} = ik_1 (\lambda_c + 2N) \hat{u}_1^{[1]} \nonumber \\
&\quad + (\lambda_c - \alpha M) \frac{\partial \hat{u}_2^{[1]}}{\partial x_2} + i \frac{\alpha M \partial \hat{V}_2^{[1]}}{\omega \partial x_2} \nonumber \\
&\quad + ik_1 \alpha M \hat{w}_1 = ik_1 (\lambda_c + 2N - \alpha M \frac{\rho_f}{\tilde{\rho}_{\text{eq}}} \hat{u}_1^{[1]}) \nonumber \\
&\quad + (\lambda_c - \alpha M) \frac{\partial \hat{u}_2^{[1]}}{\partial x_2} + \frac{\alpha M \partial \hat{V}_2^{[1]}}{\omega \partial x_2} - \frac{k_1^2}{\tilde{\rho}_{\text{eq}} \omega^2} \alpha M \hat{p}^{[1]}. \tag{A2}
\end{align}
$$

The introduction of (A2) in the six other equations of the
system (A1) leads to six first-order linear differential equa-
tions in terms of $\hat{\sigma}_{12}^{[1]}$, $\hat{\sigma}_{22}^{[1]}$, $\hat{p}^{[1]}$, $\hat{u}_1^{[1]}$, $\hat{u}_2^{[1]}$, and $\hat{V}_2^{[1]}$:

$$
\begin{align*}
& \frac{\partial \hat{\sigma}_{12}^{[1]}}{\partial x_2} + ik_1(\lambda_c - \alpha M) \frac{\partial \hat{u}_2^{[1]}}{\partial x_2} - \frac{k_1 \alpha M \partial \hat{V}_2^{[1]}}{\omega \partial x_2} + ik_1 \frac{\rho_f}{\tilde{\rho}_{eq}} - \frac{k_1^2}{\tilde{\rho}_{eq} \omega^2} \alpha M \hat{p}^{[1]} + \rho \omega^2 - \frac{\rho_f^2}{\tilde{\rho}_{eq}} \omega - k_1^2 \lambda_c + 2N - \alpha M \frac{\rho_f}{\tilde{\rho}_{eq}} \hat{p}^{[1]} \Biggr) \hat{u}_1^{[1]} = 0, \\
& \frac{\partial \hat{\sigma}_{22}^{[1]}}{\partial x_2} + ik_1 \hat{\sigma}_{12}^{[1]} + \omega^2 (\rho - \rho_f) \hat{u}_2^{[1]} + i \rho_f \omega \hat{V}_2^{[1]} = 0, \quad - \frac{\partial \hat{p}^{[1]}}{\partial x_2} + \omega^2 (\rho_f - \tilde{\rho}_{eq}) \hat{u}_2^{[1]} + i \tilde{\rho}_{eq} \omega \hat{V}_2^{[1]} = 0, \\
& N \frac{\partial \hat{u}_1^{[1]}}{\partial x_2} - \hat{\sigma}_{12}^{[1]} + ik_1 N \hat{u}_2^{[1]} = 0, \\
& (\lambda_c + 2N - \alpha M) \frac{\partial \hat{u}_2^{[1]}}{\partial x_2} + i \frac{\alpha M}{\omega} \frac{\partial \hat{V}_2^{[1]}}{\partial x_2} - \hat{\sigma}_{22}^{[1]} + ik_1 \frac{\rho_f}{\tilde{\rho}_{eq}} \alpha M \hat{u}_1^{[1]} - \frac{k_1^2 \alpha M}{\tilde{\rho}_{eq} \omega^2} \hat{p}^{[1]} = 0, \\
& M(\alpha - 1) \frac{\partial \hat{u}_2^{[1]}}{\partial x_2} + i \frac{M \partial \hat{V}_2^{[1]}}{\omega \partial x_2} + 1 - \frac{k_1^2}{\tilde{\rho}_{eq} \omega^2 M} \hat{p}^{[1]} + ik_1 M \alpha - \frac{\rho_f}{\tilde{\rho}_{eq}} \hat{u}_1^{[1]} = 0. \tag{A3}
\end{align*}
$$

Introducing the state vector $\hat{\mathbf{W}} = \langle \hat{\sigma}_{12}^{[1]}, \hat{\sigma}_{22}^{[1]}, \hat{p}^{[1]}, \hat{u}_1^{[1]} | \hat{u}_2^{[1]}, \hat{V}_2^{[1]} \rangle$, the first-order linear ordinary differential system can be written in the form

$$
\mathbf{B} \frac{\partial}{\partial x_2} \hat{\mathbf{W}} + \mathbf{D} \hat{\mathbf{W}} = \mathbf{0}, \quad (\text{A4})
$$

wherein **B** and **D** are defined in Eqs. (11) and (12).

<footnote n="1">¹L. De Ryck, W. Lauriks, Z. E. A. Fellah, A. Wirgin, J. P. Groby, P. Leclaire, and C. Depollier, “Acoustic wave propagation and internal fields in rigid frame macroscopically inhomogeneous porous media,” J. Appl. Phys. **102**, 024910 (2007).</footnote>

<footnote n="2">²L. De Ryck, W. Lauriks, P. Leclaire, J. P. Groby, A. Wirgin, and C. Depollier, “Reconstruction of material properties profiles in one-</footnote>

<footnote>dimensional macroscopically inhomogeneous rigid frame porous
media in the frequency domain," J. Acoust. Soc. Am. **124**, 1591–1606
(2008).</footnote>

<footnote n="3">³O. Matsuda and O. B. Wright, “Reflection and transmission of light in multilayers perturbed by picosecond strain pulse propagation,” J. Opt. Soc. Am. B, **19**, 3028–3041 (2002).</footnote>

<footnote n="4">⁴O. Forslund and S. He, “Electromagnetic scattering from an inhomogeneous grating using a wave-splitting approach,” Prog. Electromagn. Res. **19**, 147–171 (1998).</footnote>

<footnote n="5">⁵J. L. Buchanan, R. P. Gilbert, A. Wirgin, and Y. Xu, Marine Acoustics: Direct and Inverse Problems (SIAM, Philadelphia, 2004).</footnote>

<footnote n="6">⁶J. van Manen, J. O. A. Robertsson, and A. Curtis, “Modeling of wave propagation in inhomogeneous media,” Phys. Rev. Lett., **94**, 164301 (2005).</footnote>

<footnote n="7">⁷J. G. Berryman and R. R. Greene, “Discrete inverse methods for elastic waves in layered media,” Geophysics, **45**, 213–233 (1980).</footnote>

<footnote n="8">⁸V. Tournat, V. Zaitsev, V. Gusev, V. Nazarov, P. Bequin, and B. Castagnede, “Probing weak forces in granular media through nonlinear dynamic</footnote>
dilatancy: Clapping contacts and polarization anisotropy,” Phys. Rev. Lett. 92, 085502 (2004).

⁹M. A. Biot, “Mechanics of deformation and acoustic propagation in porous media,” J. Appl. Phys. 33, 1482–1498 (1962).

¹⁰L. De Ryck, J.-P. Groby, P. Leclaire, W. Lauriks, A. Wirgin, C. Depollier, and Z. E. A. Fellah, “Acoustic wave propagation in a macroscopically inhomogeneous porous medium saturated by a fluid,” Appl. Phys. Lett. 90, 181901 (2007).

¹¹J.-P. Groby, L. De Ryck, P. Leclaire, A. Wirgin, W. Lauriks, R. P. Gilbert, and Y. S. Xu, “Use of specific Green’s functions for solving direct problems involving a heterogeneous rigid frame porous medium slab solicited by acoustic waves,” Math. Meth. Appl. Sci. 30, 91–122 (2007).

¹²A. N. Stroh, “Steady state problems in anisotropic elasticity,” J. Math. Phys. 41, 77–103 (1962).

¹³G. Peano, “Integration par serie des equations differentielles lineaires (Integration of linear differential equations by series expansion),” Math. Annal. 32, 450–456 (1888).

¹⁴M. C. Pease, *Methods of Matrix Algebra* (Academic Press, New York, 1965), Chaps. 6 and 7.

¹⁵A. Shuvalov, O. Poncelet, and M. Deschamps, “General formalism for plane guided waves in transversely inhomogeneous anisotropic plates,” Wave Motion 40, 413–426 (2004).

¹⁶A. Shuvalov, O. Poncelet, M. Deschamps, and C. Baron, “Long-wave-length dispersion of acoustic waves in transversely inhomogeneous anisotropic plates,” Wave Motion 42, 367–382 (2005).

¹⁷J. L. Tsalamangas and N. K. Uzunoglu, “Radiation from a dipole in the proximity of a general anisotropic grounded layer,” IEEE Trans. Antennas Propagat. 33, 165–172 (1985).

¹⁸J. R. Bowler, L. D. Sabbagh, and H. A. Sabbagh, “A theoretical and computational model of eddy-current probes incorporating volume integral and conjugate gradient methods,” IEEE Trans. Magn. 25, 2650–2664 (1989).

¹⁹A. L. Shuvalov, E. Le Clezio, and G. Feuillard, “The state-vector formalism and the Peano series solution for modelling guided waves in functionally graded anisotropic piezoelectric plates,” Int. J. Eng. Sci. 46, 929–947 (2008).

²⁰B. Brouard, D. Lafarge, and J. F. Allard, “A general method of modeling sound propagation in layered media,” J. Sound Vib. 83, 129–142 (1995).

²¹M. Carcione, *Wavefield in Real Media: Wave Propagation in Anisotropic, Anelastic and Porous Media* (Pergamon, New York, 2001), Chap. 7.

²²S. R. Pride and J. G. Berryman, “Connecting theory to experiment in poroelasticity,” J. Mech. Phys. Solids 46, 719–747 (1998).

²³O. Dazel, B. Brouard, C. Depollier, and S. Griffiths, “An alternative Biot’s displacement formulation for porous materials,” J. Acoust. Soc. Am. 121, 3509–3516 (2007).

²⁴M. A. Biot and D. G. Willis, “The elastic coefficients of the theory of consolidation,” J. Appl. Mech. 24, 594–601 (1957).

²⁵M. A. Biot, “Theory of propagation of elastic waves in fluid-saturated porous solid,” J. Acoust. Soc. Am. 28, 168–178 (1956).

²⁶J.-F. Allard and N. Atalla, *Propagation of Sound in Porous Media: Modeling Sound Absorbing Materials* (John Wiley and Sons, Chichester, UK, 2006), Chaps. 5 and 6.

²⁷J.-F. Allard and Y. Champoux, “New empirical equations for sound propagation in rigid frame porous materials,” J. Acoust. Soc. Am. 91, 3346–3353 (1992).

²⁸D. J. Johnson, J. Koplik, and R. Dashen, “Theory of dynamic permeability and tortuosity in fluid-saturated porous media,” J. Fluid Mech. 176, 379–402 (1987).

²⁹C. Baron, “Le developpement en serie de Peano du matricant pour l’etude de la propagation des ondes elastiques en milieux a proprieties continue-ment variables (Matricant Peano series development to study elastic waves propagation in continuously varying properties materials),” Ph.D. thesis, Universite de Bordeaux I (2005).

³⁰N. A. Haskell, “The dispersion of surface waves on multilayered media,” Bull. Seismol. Soc. Am. 43, 377–393 (1953).

³¹W. T. Thomson, “Transmission of elastic wave through a stratified solid medium,” J. Appl. Phys. 21, 89–93 (1950).

³²O. Doutres, N. Dauchez, J. M. Genevaux, and O. Dazel, “Validity of the limp model for porous materials: A criterion based on the Biot theory,” J. Acoust. Soc. Am. 122, 2038–2048 (2007).

³³M. Tamura, *Spatial Fourier transform method of measuring reflection coefficients at oblique incice I: Theory and numerical examples, Technical report* (National Institute for Environmental Studies, Tsukuba, Japan, 1990).

³⁴J.-F. Allard and Y. Champoux, “In-situ two-microphone technique for the measurement of the acoustic surface impedance at oblique incidence and low frequencies,” Noise Control Eng. J. 32, 15–23 (1989).