JAMP  Vol.7 No.11 , November 2019
Mathematical Modeling of Porous Medium for Sound Absorption Simulations II: Wave Propagation and Interface Conditions
Abstract: The application of porous medium has a myriad of applications in different industries: automotive, aerospace, civil (commercial, residential), environmental noise control, and biomedical. In the past, design questions involving porous material were addressed with seat-of-the-pants decisions that led to multiple/iterative prototypes and experiments that were costly and time consuming. The objective, in this series of publications pertaining to porous medium, is to establish tools that will lead to effective and accurate simulations involving porous medium. In this third installment of this series the focus is on establishing the constitutive equations using tensors and then applying Transfer Matrix Method (TMM) to calculate diffuse field Transmission Loss (TL) across structures that comprises of layers of different porous medium. The constitutive equations are obtained by relating information regarding the micro-structure make-up to macro level properties. In order to apply the TMM, the equations for wave propagation across different mediums need to be developed and in turn represent these propagation properties in a matrix format. Additionally, the boundary condition between each layer type is defined in order to ensure numerical stability. The author’s current research effort is running simulations for the automotive industry to predict NVH environments. Therefore, TL calculations pertaining to the materials that are utilized in the interior of automobiles are used, in this paper, as a test bed for the developed analytical tools. Case in point, the TL for a multi-layered material consisting of one panel and two different layers of foam is calculated and compared to experimental data. Future publication goals will be to apply these tools in the biomedical field; an example will be to model and run simulations of different organs like the liver and lungs that are po-rous in nature.

1. Introduction

This paper is the third installment from a series of publications pertaining to modeling of porous medium. The first paper, Teagle et al. [1], derived a coupled set of fluid/structure equations for a porous medium applying asymptotic and homogenization techniques. In [1] it is shown that there are mainly 3 modes of energy transformation: 1) The first mode is through the connection between the micro and macro structural framework of the porous skeleton, 2) The second is via the viscous boundary layer, and 3) The third interaction is through thermal (entropy) boundary layer. The combination of the viscous boundary layer and how tortuous the porous material is, results in the encapsulation of the fluid medium. This, in turn, changes the apparent mass of the structural medium. Details pertaining to these encapsulating phenomena can be found in the work by Johnson, et al. [2]. In [2] the concepts of tortuosity, viscous length, and viscous permeability are explained. The thermal interchange is described in Teagle, et al. [3]. In [3], the mathematical description of the thermal energy dissipated by the thermal boundary layer is explained along with the relaxation process. This thermal exchange changes the acoustic bulk modulus of the porous medium and thus the speed of sound inside the porous layer. These thermal phenomena are then represented by the parameters of thermal length and thermal permeability.

This paper is a continuation of [1]. The equations of motion derived in the aforementioned publication will be presented in a form that is used for calculation of TL for porous material used in the automobile industry. The combination of 1) the change of density due to the viscous effect and 2) the change of acoustic bulk modulus due the thermal exchange will be applied. Additionally, a matrix representation of the wave propagation (for both forward and reflected propagating wave) that incorporates the aforementioned viscous and thermal effects is developed. Each porous layer makes contact with another type of porous medium, elastic panel, or air. This study will establish the correct boundary conditions between mediums in order to run numerical simulations that lead to stable and unique solutions. In order to translate these boundary conditions to the numerical model, interface matrices are developed. In continuation, the global matrix that represents the multi-layered material is assembled and applying the definition of impedance the diffuse field TL is calculated. In the TMM formulation it is assumed that each layer is of infinite extent. A correction applying Green’s function technique is applied.

2. Formulation

2.1. Basic Tensor Calculus and Notation

Tensor calculus and concepts from differential geometry are used extensively in this paper. This notation affords a level of abstraction that leads to an efficient explanation of the stresses, strains, and their relationships. Here, the usual notation and represent contravariant and covariant vectors (tensors of order 2). Physical tensors like stress, that are neither contravariant or covariant, are designated here as or (with over bar and void of any ij subscript/superscript or ij subscripts in parenthesis). The relationship between physical tensors, like stress, and their cotravariant and covariant counterparts are the following


where are scaling factors that satisfy the following


It should be noted that is the components of the metric tensor. The Einstein summation convention is used in which covariant index followed by the identical contravariant index is implicitly summed over, thus contracting the order of the tensor. In this paper, if 2 tensors of different order are shown next to each other, this contraction rule is followed. To obtain expression for strains, the gradient of the first term of the asymptotic expansion of the displacement vector, , is required. The expression for the components of the tensor of order 2 is the following


and the expression for strain is given as


The elastic constant, , is a fourth rank tensor. Due to tensor contraction rules the stress tensor, , is of second rank. Additionally, the divergence of the stress tensor, , results in the force field. From Equation (1), the contravariant component of the physical stress tensor, , is


Thus, the physical component of the force field in the jth direction is given by


The use of orthogonal coordinates results in


is a tensor operator that transforms the stress data into a hydrostatic form.

2.2. Fluid-Structure Interaction: Dynamic Equations

This paper applies the definitions and results obtained in [1]. In that publication, the fluid/structure interaction equations are derived via asymptotic and homogenization techniques. A slight modification is applied to the averaged relative displacement of the fluid with respect to the skeleton,. In this paper this equation is represented as


In [1], Equations (54) and (55) are derived and represent the structural stress, , and internal fluid pressure,. These are rewritten below as Equations (9) and (10)



where and are elastic constant tensor (fourth rank) per their explanations given in [1].


is the interstitial fluid bulk modulus, details on its derivation can be found in Teagle et al. [4]. Additionally, in [4] it is shown that


They act as scalar multiples of the identity operator.

Set and therefore


Analyzing the derivations in [1], represents the proportion of fluid pressure that produces the same strains as the total stress. Additionally, it is important to note that Equations (12a, b) (see Equations (21, 22)) are conditions that make it possible for there to be a strain-energy term for the porous material. This justifies Biot’s assumption of the existence of a potential when he derived his landmark equation, pertaining to porous rocks, for oil exploration.

Equations (44) and (45) from [1] represent the inertial forces for the structural and fluid components and are rewritten below


where, and


Subtracting Equation (13) from (12c) and substituting the definition of, Equation (8), the following expression for the structural portion is obtained







and applying the definition for, an efficient representation of Equation (13) is


Substituting for, in Equation (9), the expression in (10), and in turn taking the divergence the dynamical equation for the structural portion of the foam can be expressed as


In order to derive (17), the definition for, Equation (12b), was used. When, it is easy to verify that and are new Lame’ constants of the elastic portion when the porous material is drained.

Additionally, taking the gradient of Equation (10) followed by a multiplication by, the following expression dealing with pressure is obtained


Combining Equations (14) and (17), the dynamic equations pertaining to the structural/skeleton portion of the porous medium is obtained


Similarly, combining Equations (16) and (18) results in the macro level equations describing the fluid motion


Equation (19) can be expressed in a compact form


and by applying the coupling definition of Q, Equation (20) becomes


2.3. Wave Participation Factor: Eigenvalue Problem

In this section, the study adapts the concepts from Brouard, et al. [5]. The motion of the structural/skeleton and the fluid portion are described by introducing the potential scalar functions and and the potential vector functions, and, such that



Here, the basic definition from mechanics is used where the gradient represent portion of the displacement vector that is purely dilatational or in compression. The curl of and represent any shear motion associated with the displacement. Substitute these definitions into Equations (21) and (22)




Recognize that the operator commutes with i.e. and for any general vector.

Using these relationships, Equations (25) and (26) are rewritten in the following form



Gathering terms corresponding to the gradient operator results in



These equations will give the formulation to solve for the compression waves.

Equations (29) and (30) are represented in matrix form





and (34)

Applying (34), Equation (33) is turned into an eigenvalue problem


The eigenvalues produce the complex wave number for the compression waves and they will have the following expression




Such that


Two sets of eigenvectors are generated the relationship is known as the participation factor. Substituting the eigenvectors into Equation (31) ones gets

for (38)

These results indicate that there exist two compression waves traveling the porous medium, one is fast and the other one is slower. is important since the number will indicate which wave, whether acoustic or solid, has the most contribution at that particular frequency. For the shear wave we accumulate all terms in Equation (27) and (28) that contain the curl operator



Equation (40) results in, insert this in (39)


Setting, the complex wave number for the shear wave can be given as


and the participation factor is


2.4. Wave Propagation: Porous Medium

This portion will analyze wave propagation in a semi infinite porous medium. Consider Figure 1.

The disturbance wave is traveling down at an angle of incidence of θ in the general xz plane and will impinge the porous medium at the z = 0 level. A portion of this wave will be reflected and some of it will be transmitted to travel through the porous medium until it impinges the next layer, at z = L, of different impedance. A portion of the wave will again be reflected back into the porous medium and the rest will be transmitted to the next layer.

Since there exist 3 wave types, 2 compression (dilatation) waves and 1 shear wave, and since each wave type has 2 waves (forward moving and the reflected wave) there will be 6 variables. These variables are:

structural velocity in the x and z direction

acoustic velocity in the x direction and z direction

= structural stress component perpendicular to z face in the z direction

structural shear component

acoustic stress component

An array of these components will be represented by . The continuity of these variables at the interface between 2 layers will become the boundary conditions that the traveling wave has to satisfy. Using the definition of the participation factor above, the potential functions can be expressed as



Figure 1. Porous medium element: Thickness L, angle of incidence θ.



forward wave

reflected wave

where, indicates the matching of acoustic and porous waves along the x direction (k is the acoustic wave number and k sin(θ) is the projected acoustic wave number in the x direction). projects the porous wave in the z direction. Note: 1) if is real, there is a propagating wave, 2) if is complex there is a decaying wave in the z direction. The following notations will be used, and . An expression for can now be derived


The six variables can be written as a function of the potential functions and using the above notation:




The following three equations are essential to define the stress components




Utilizing the stress definitions established in section 2.2, the stress components can be written as




Applying the corresponding expressions in (44) into Equations (45)-(54) a set of equations involving, and complex trigonometry equations will be obtained. For example if (44) is applied to (45) the following expression in sues


Similar operations are done to derive expressions for and. The final results are represented in matrix form





The idea behind a transfer matrix is to relate the 6 variable, V(z), at z = L to the conditions at z = 0, i.e. relate the matrix with respect V(L) . The conditions at z = 0 can be written as


Here, is known as the transfer matrix (58)

2.5. Interface Conditions: Uniqueness

A general form for the required boundary conditions between the layers of multi-layered materials is developed. These boundary conditions are developed to ensure uniqueness and numerical stability. Special focus is on deriving interface conditions for the porous material. This section puts concepts introduced by Deresiewicz, et al. [6], under the context and language set forth in this paper.

The kinetic energy, per unit volume for the two phase system is expressed as


, , and are, , void the portion with tortuosity (without viscous effect). The dynamic equations in Equation (17) assume that the porous material is statistically isotropic. Recall that is the stress in the skeleton and is the stress in the acoustic medium, , is the macro-level identity (hydrostatic) tensor. The strain energy, per unit volume is given by

, = Strain tensor for liquid (solid) (60)

Taking the time derivative of the kinetic and strain energies results in


The following tensorial relationship holds


Combining Equations (61, 62) and the divergence theorem the expression for power is obtained


The portion inside the volume integral represents the traction and inertial work applied to the skeleton and acoustic medium respectively. Equations (21) and (22) are written in shorthand form



The term is the viscous force term due to the interstitial fluid, Β is a viscous transfer function replicating the viscous effects. When (64) and (65) are plugged into Equation (63), the volume portion of the equation becomes


is the dissipation function per unit volume. Taking (66) into consideration the power expression is expressed as


The right hand side of (67a) expresses the rate at which work is done on the material by surface forces, is the normal vector to the surface. Now consider 2 different porous medium with volumes and, boundaries and and where represents the common boundary between porous medium 1 and 2. To establish uniqueness, for each of the two mediums the field quantities will be replaced by a difference representing the possibility of two different solutions but with the same boundary conditions, making the right hand side of Equation (67a) equal to zero. is positive definite, therefore the only way that the difference version will satisfy the zero condition is that everywhere, indicating that the solution is unique. The surface integral in Equation (67a) is partitioned in the following way,

is the element number (67b)

In the common portion,. Combining the two surface integrals results in


For the non-intersecting boundaries (S1 and S2) the boundary conditions will have to be given. At the Sc boundary, continuity is required across the interface to be able to maintain uniqueness. If the skeleton phase is to remain in contact with each other and to maintain the principle of conservation of mass of the acoustic medium it is continuity of the normal relative velocity of the acoustic medium with respect to the skeleton, i.e.


Applying (69) into (68), the Sc integrand becomes


To assure continuity of condition (70), continuity in the following quantities has to be maintained


The non-alignment of the pores can produce a pressure drop across the interface so the continuity condition for pressure p is modeled as, is a coefficient of resistance. The materials used in this study are highly porous so was set equal to zero. This also seemed to be the coefficient that gave the best results after a quick parameter study. If the intersection or interface is between a porous medium and a plate then, or. The continuity conditions for this case (porous-plate) will be the following



The over-bar in conditions (72) and (73) represent quantities in the plate. In case the interface is between porous medium and air, for the layer pertaining of air will be set equal to one, the stresses and will be set equal to zero and. The interface conditions for porous and air becomes



The continuity of was applied to obtain Equation (75).

3. Results and Conclusion

A Simple application of TMM is applied to the simple layer configuration shown in Figure 2. The figure shows a multi layer system where there is a plate that is glued to foam 1 which in turn is connected to foam 2. The parameters pertaining to plate 1, foam 1, and foam 2 are listed in Table 1. This simulation also applied the following parameters that are not listed in the table: Viscous length = 40 μm, Thermal length = 80 μm. The goal of this simulation is to calculate Random Incidence Transmission Loss (TL). This is achieved by applying Finite Size correction Transfer Matrix Method (FTMM) [7] [8]. The Finite Size Correction is achieved by incorporating a Green’s Function integration to the radiation efficiency. The dynamic equations pertaining to the skeleton and fluid along with the boundary conditions depicted in Equations (71)-(75) are applied in order to simulate how the structural and fluid stresses will change as the acoustic wave travels through the multi-layered material.

Equations (9) and (10) represent the structural equations in its most general form. A detailed Finite Element (FE) can be developed in order to analyze the interplay between macro and micro levels and in turn obtain expressions for P, Q, R, and N; parameters required in Equations (21) and (22). In this study there is not enough information to construct an FE model, the only structural information available are the experimentally obtained bulk modulus, , and the

Figure 2. Multi layer configuration.

Table 1. Parameters pertaining to plate and foam.

shear modulus, N, both pertaining to the porous material. Materials used in the acoustic field usually consist of a highly porous material whose structural portions are considerably stiffer. Due to one of Biot’s [9] experiment, a jacketed porous medium experiences a hydrostatic pressure pj, but the air inside the porous medium experiences no change in pressure. The definition for porosity

before and after deformation respectively are and.

Assuming highly stiff solid frame will lead to small changes in solid volume, therefore and, is jacketed frame dilation. This information leads to the following equation for


Similarly, a relationship for the fluid dilatation is obtained,. Combining this last expression with that of Equation (76) the following important relationship which relates the fluid and frame dilatation is obtained


Analyzing Equations (17) and (18), the following expressions for structural stress is obtained



For the porous material with relatively very stiff frame, Equation (11) shows that M can be approximated by. In turn, taking the trace of Equation

(78), dividing that result by 3, and setting in Equation (79), the jacketed experiment is simulated

= porous bulk modulus (80)


Equation (77) is applied to obtain (80) and (81). Firstly, from Equation (81), it is deduced that




Plugging Equation (82a, b) into (80) and solving for P results in


The results incorporating Equations (82abc) (indicated as “case 9”) into the simulation are shown in Figure 3. The results are compared to measured

Figure 3. Top: TL comparison (dB); Bottom: TL difference (dB) measured data is the reference.

(“meas”) results and also to the empirical formulas of Delaney and Bazley (“DB”) [10] [11]. The graph shows that the calculated results come within 0.7 dB.


tr = trace of a tensor, = averaged at the microscopic level

v = Fluid Velocity, p = Fluid Pressure

u = Structural displacement, U = Fluid displacement

= Porosity, n = unit normal pointing into the solid

= Viscosity, = Second Viscosity

= Domain Occupied by Fluid, = Domain Occupied by Structure

= Coefficient of thermal Conductivity, Pr = Prandtl Number =

= Stress Tensor in the Fluid, = Stress Tensor in the structure

= density of the Fluid, = density of the Structure

= Lame Parameters, n = unit normal pointing into the solid

T = Temperature deviation, = Coefficient of thermal Conductivity

= Dynamic Viscous Permeability, = Dynamic Thermal Permeability

Cp = specific heat at constant pressure, Cv = Specific heat at constant volume

= Elastic fourth ranked Contra variant tensor, operates on

= Acoustic bulk Modulus

= asterisk superscript means scaled variable (dimensionless)

Cite this paper: Teagle-Hernandez, A. , Ohtmer, O. , Nguyen, D. (2019) Mathematical Modeling of Porous Medium for Sound Absorption Simulations II: Wave Propagation and Interface Conditions. Journal of Applied Mathematics and Physics, 7, 2780-2795. doi: 10.4236/jamp.2019.711191.

[1]   Teagle-Hernandez, A., Ohtmer, O. and Nguyen, D. (2018) Mathematical Modeling of Porous Medium for Sound Absorption Simulations: Application of Multi-Scales and Homogenization. Journal of Applied Mathematics and Physcis, 6, 2705-2717.

[2]   Johnson, D. and Koplik, J. (1987) Theory of Dynamic Permeability and Totuosity in Fluid-Saturated Porous Media. Journal of Fluid Mechanics, 176, 379-402.

[3]   Teagle-Hernandez, A., Ohtmer, O. and Nguyen, D. (2019) Model for Frequency Dependence of Thermal Permeability in Order to Quantify the Effects of Thermal Exchange on Wave Propagation in Multi Layered Porous Medium. Journal of Applied Mathematics and Physics, 7, 1440-1462.

[4]   Burridge, R. and Keller, J. (1981) Poroelasticity Equations Derived from Microstructure. J. Acoust. Soc. Am., 70, 1140-1146.

[5]   Brouard, B., Lafarge, D. and Allard, J.F. (1995) A General Method of Modeling Sound Propagation in Layered Media. J. Sound Vib, 183, 129-142.

[6]   Deresiewicz, H. and Skalak, R. (1963) On Uniqueness in Dynamic Poroelasticity. Bulletin of the Seismological Society of America, 53, 409-416.

[7]   Nguyen, D. and Teagle-Hernandez, A. (2019) Application of Transfer Matrix Method, Green’s Functions, and Variational Techniques to Predict Diffuse Absorption Coefficients. CalTech-SOCAMS, Pasadena, CA, 27 April 2019.

[8]   Nguyen, D. and Teagle-Hernandez, A. (2019) Greens Functions and Variational Techniques for the Transfer Matrix Correction. Acoustical Society of America.

[9]   Biot, M. (1962) Mechanics of Deformation and Acoustic Propagation in Porous Media. Journal of Applied Physiscs, 34, 1482-1498.

[10]   Delany, M. and Bazley, E. (1970) Acoustical Properties of Fibrous Materials. Applied Acoustics, 3, 105-116.

[11]   Miki, Y. (1990) Acoustical Properties of Porous Materials—Modification of Delany-Bazley Models. Journal of the Acoustical Society of Japan, 11, 19-24.