Propagation of Acoustic Waves in a One-dimensional Macroscopically Inhomogeneous Poroelastic Material

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 1 and (2) the retrieval of the spatially varying material parameters of industrial foams. 2 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, 3,4 to acoustics, 5,6 and geophysics 7 and granular media. 8 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 9 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. 11 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 12 together with Peano series. 13,14 The Stroh formalism is largely used to model acoustic wave propagation in stratified anisotropic elastic materials. 15,16 Similar methods (transmission matrix method) are used in electromagnetism to model the propagation of electromagnetic waves in anisotropic or gyrotropic stratified materials. 17,18 In the frequency domain, it consists in the rewriting of the constitutive linear stressstrain 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. 19 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 20 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, 9,21,22 the generalized formulation of the Biot theory 9 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 sðx; tÞ [with x ¼ (x 1 ,x 2 )], we prefer to deal with the s(x,x), related to sðx; tÞ by the Fourier transform sðx; tÞ ¼ Ð 1 À1 sðx; xÞe ixt dx, wherein x ¼ 2p is the angular frequency, with v the frequency. Henceforth, we drop the x in s(x,x) so that it is written s(x).
From the alternative formulation of Biot, 9 the stress strain relations in an initially stress free, statistically isotropic poroelastic material take the form where in d ij denotes the Kronecker symbol. The components of the total stress tensor are r ij , the fluid pressure in the pores is p, and the components of the strain tensor are e ij ¼ 1/2(u i,j þ u ji ), with the solid displacement u, h ¼ u i,i , and f ¼ -w i , i with the fluid/solid relative displacement w ¼ / (U -u) (U being the fluid displacement and / the porosity). The Einstein summation notation is implicit in the expressions of h and f. The material properties 24 are the bulk modulus of the closed porosity system, i.e., in which the pore volume is sealed, k c ¼ k þ a 2 M, the Lamé coefficients of the elastic frame k and N, an additional elastic parameter M, and an elastic coupling coefficient a. These mechanical coefficients are related to the more commonly used in acoustics P, Q, and R coefficients from the original formulation 25 through The expressions of P, Q, and R, in case of air saturated materials, reduce to 26 wherein the bulk modulus of the skeleton is The correction function G(xPr), introduced in Ref. 27 to account for the thermal losses, is c the specific heat ratio, Pr the Prandtl number, s 1 the tortuosity, R 0 t the "thermal resistivity," and K 0 the thermal characteristic length. The thermal resistivity is related to the thermal characteristic length 27 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: wherein the bulk density of the porous medium is q, such that q ¼ ð1 À /Þq s þ /q f with q s the density of the solid and q f the density of the saturating fluid, and the mass parameterq eq (Refs. 9 and 28) is The correction function F(x), introduced in Ref. 28 and which accounts for the viscous losses, is given by with the Biot frequency x c ¼ R f /=q f s 1 , the flow resistivity R f and the viscous characteristic length K. The mass parameterq eq is the complex frequency dependent equivalent density used in the rigid frame approximation, while the complex frequency dependent equivalent bulk modulusK eq used in this approximation is related toK f throughK f ¼ /K eq . For most foams, the elastic coupling coefficient a reduces to 1, M ¼K eq and k c reduces to In the previous equations, N, k (and v), k c , a, M, /, s 1 ,

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 coordinate x 3 . A cross-sectional x 1 -x 2 plane view of the 2D scattering problem is shown in Fig. 1.
The upper and lower boundaries of the layer are flat and parallel. They are designated by C L and C 0 and their x 2 coordinates are L and 0. The porous material M [1] occupies the domain X [1] . The inhomogeneity of the plate occurs along the x 2 direction, i.e., N, k (and v), k c , a, M, /,s 1 , K, K 0 , R 0 t , R f are x 2 -dependent. An isotropic macroscopically inhomogeneous porous material is considered. The material can be viewed as a functionally graded material. The surrounding and saturating fluid is the air medium (density q f ¼ 1.213 kg m À3 , atmospheric pressure P 0 ¼ 1.01325 Â 10 5 Pa, and viscosity g ¼ 1.839 Â 10 5 kg m À1 s À1 ). The inhomogeneous porous layer is backed by a rigid plate at C 0 .
We denote the total pressure, wavevector and density, respectively, with p [0] , k [0] and q [0] in X [0] and the total stress tensor, the fluid pressure in the pores, the solid displacement, and solid/fluid relative displacement respectively with r [1] , p [1] , u [1] and w [1] in X [1] .
The wavevector k i of the incident plane wave lies in the sagittal plane and the angle of incidence is h i measured counterclockwise from the positive x 1 axis. The incident wave, initially propagating in X [0] , is expressed by p i ðxÞ The uniqueness of the solution to the forward-scattering problem is ensured by the radiation condition: p ½0 ðxÞ À p i ðxÞ $ outgoing waves; x j j ! 1; x 2 > L: (9) B. The state vector formalism and Peano series The spatial Fourier transform of the Eqs. (1) and (6) is first performed. The transformŝðx 2 ; k 1 Þ of s(x) is also introduced and can be written in the form sðxÞ ¼ŝð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ŝðx 2 ; k 1 Þ so that it is writtenŝðx 2 Þ. The geometry 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(x) or toŝðx 2 Þ. Inside the domain X [1] the normal components of the total stress tensor,r 1 ½ 12 and r 1 ½ 22 , the pressure p [1] , the solid displacements, u 1 ½ 1 and u 1 ½ 2 , and the normal component of the solid/fluid relative displacement w 1 ½ 2 are continuous along the x 2 -axis. It seems natural to choose these 6 parameters 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 C L , the normal components of the stress tensor ðr 2 are continuous at any x 2 in X [1] , the normal component of the velocity Àixðw  11 , 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: and 2   6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  4   3   7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  5 : The matrices B and D are x 2 -dependent, because k c (x 2 ), a(x 2 ), M(x 2 ), N(x 2 ),q eq x 2 ð Þ, and q(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 wherein M is the so-called matricant, 14 which relates the value of the state vectorŴð0Þ, at x 2 ¼ 0, to the value of the state vectorŴð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 , 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. 13,14,19,29 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 (14) and the evaluation of M is performed via the iterative scheme M f0g ¼ I; such that lim n!1 M n ¼ M.

C. The boundary problem
The application of the boundary conditions at both interfaces C L ½r vectorsŴðLÞ andŴð0Þ, which are required to solve the problem. From the separation of variables, the radiation conditions, and the spatial Fourier transform, the pressure field in X [0] can be written as where R is the reflection coefficient. The state vectors become and wherein S accounts for the excitation of the system by the plane incident wave, L L relates the unknowns R;û Finally, the introduction of (17) and (18) in (13) gives rise to the final system of equations, whose solution contain the reflection coefficient R This system is solved for each frequency and directly provide the reflection coefficient associated with the plane incident wave. The pressure field in X [0] can then be calculated through p ½0 ¼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: wherein C is the step value, which is different for each parameter in Table I, erf is the error function, x 0 2 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 0 2 ¼ 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 (h) are also plotted 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, 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 s 2 ¼ 83 cm above the ground. A microphone is placed successively at x ð1Þ 2 ¼ 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 guaranties 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 ðjÞ 2 þ A i Re ik ½0 x ðjÞ 2 , j ¼ 1, 2 at normal incidence. It directly equals the spatial Fourier transformp ð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 This formula is commonly used when the measuring technique, sometimes referred to as the Tamura method 33 or sometimes referred to as PP-method, 34 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)  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 0 2 ¼ 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: 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. The expanded form of the spatial Fourier transform of (1) and (6) iŝ The solid/fluid relative displacementŵ ½1 2 is then replaced bŷ V q eq x 2 aMp ½1 : (A2) The introduction of (A2) in the six other equations of the system (A1) leads to six first-order linear differential equations in terms ofr