Poroelastic acoustic medium in frequency domain

Stiffness properties

Consider a homogeneous isotropic linear poroelastic material with porosity $ϕ$ defined as the ratio of void volume to the whole volume of the sample. The void volume is fully saturated with fluid; the fluid is initially at rest. The Biot (1956) stress-strain constitutive relations are

\begin{matrix} (1) & σ_{i j}^{s} = [(P - 2 G) Θ^{s} + Q Θ^{f}] δ_{i j} + 2 G ε_{i j}^{s}, \end{matrix}

\begin{matrix} (2) & σ_{i j}^{f} = - ϕ p δ_{i j} = [Q Θ^{s} + R Θ^{f}] δ_{i j} . \end{matrix}

Here

$σ_{i j}^{s}$: is the structural stress,
$σ_{i j}^{f}$: is the fluid stress,
$p$: is the fluid pore pressure,
$ε_{i j}^{s}$: is the structural strain,
$Θ^{s} = ε_{i i}^{s}$: is the structural dilatation,
$Θ^{f}$: is the fluid dilatation, and
$δ_{i j}$: is the Kronecker delta.

P, G, Q, and R are the material coefficients.

Three static physical experiments determine the P, G, Q, and R values, as outlined in Allard and Atalla (2009).

The first experiment is a pure shear test, $Θ^{s} = Θ^{f} = 0$ and $σ_{i j}^{s} = 2 G ε_{i j}^{s}, σ_{i j}^{f} = 0$ (fluid cannot resist shear in statics), so G is the shear modulus of the material determined in this experiment.

In the second experiment the material is wrapped in a flexible membrane jacket that is impervious to the fluid and is subjected to an external hydrostatic pressure $p_{2}$ while zero pressure is maintained inside ( $σ_{i j}^{f} = 0$ ).

From this experiment the skeleton "in vacuo" bulk modulus is determined as $K_{b} = - p_{2} / Θ_{2}^{s},$ where $Θ_{2}^{s}$ is the structural dilatation measured in the experiment. The constitutive relations become

\begin{matrix} (3) & - p_{2} = (P - \frac{4}{3} G) Θ_{2}^{s} + Q Θ_{2}^{f}, \end{matrix}

\begin{matrix} (4) & 0 = Q Θ_{2}^{s} + R Θ_{2}^{f} . \end{matrix}

The porosity changes during this experiment as far as the material is compressed. For Abaqus you indirectly specify the G and $K_{b}$ through the elastic skeleton Young's modulus and Poisson's ratio.

In the third experiment the unwrapped material sample is placed inside fluid under pressure $p_{3};$ therefore, the structural stress tensor is

σ_{i j}^{s} = - p_{3} {(1 - ϕ) δ}_{i j} .

The constitutive equations are rewritten as

\begin{matrix} (5) & - p_{3} {(1 - ϕ) δ}_{i j} = (P - \frac{4}{3} G) Θ_{3}^{s} + Q Θ_{3}^{f}, \end{matrix}

\begin{matrix} (6) & - ϕ p_{3} = Q Θ_{3}^{s} + R Θ_{3}^{f} . \end{matrix}

The skeleton material bulk modulus is then $K_{s} = - p_{3} / Θ_{3}^{s}$ , and the fluid bulk modulus is $K_{f} = - p_{3} / Θ_{3}^{f}$ . There is no change in porosity in this experiment, and the dilatation of the sample is the same as if the material is not porous. For steady-state dynamic analyses you specify the real-valued $K_{s}$ directly, whereas the generally complex $K_{f}$ is specified either directly in the Biot-Atalla material model or is computed based on other parameters in the Biot-Johnson model.

From Equation 6

\begin{matrix} (7) & ϕ = - Q \frac{Θ_{3}^{s}}{p_{3}} - R \frac{Θ_{3}^{f}}{p_{3}} = \frac{Q}{K_{s}} + \frac{R}{K_{f}} . \end{matrix}

From Equation 5

\begin{matrix} (8) & 1 - ϕ = - \frac{(P - \frac{4}{3} G) Θ_{3}^{s}}{p_{3}} + Q (- \frac{Θ_{3}^{f}}{p_{3}}) = \frac{(P - \frac{4}{3} G)}{K_{s}} + \frac{Q}{K_{f}} . \end{matrix}

From Equation 4 $Θ_{2}^{f} = - \frac{Q}{R} Θ_{2}^{s}$ . Divide Equation 3 by $(- p_{2})$ to obtain

\begin{matrix} (9) & 1 = (P - \frac{4}{3} G) (- \frac{Θ_{2}^{s}}{p_{2}}) - \frac{Q^{2}}{R} (- \frac{Θ_{2}^{s}}{p_{2}}) = \frac{(P - \frac{4}{3} G - \frac{Q^{2}}{R})}{K_{b}} . \end{matrix}

From Equation 7, Equation 8, and Equation 9 the three elastic coefficients P, Q, and R are obtained as

P = \frac{(1 - ϕ) [1 - ϕ - \frac{K_{b}}{K_{s}}] K_{s} + \frac{ϕ K_{s}}{K_{f}} K_{b}}{1 - ϕ - \frac{K_{b}}{K_{s}} + \frac{ϕ K_{s}}{K_{f}}} + \frac{4}{3} G,

Q = \frac{[1 - ϕ - \frac{K_{b}}{K_{s}}] ϕ K_{s}}{1 - ϕ - \frac{K_{b}}{K_{s}} + \frac{ϕ K_{s}}{K_{f}}},

R = \frac{ϕ^{2} K_{s}}{1 - ϕ - \frac{K_{b}}{K_{s}} + \frac{ϕ K_{s}}{K_{f}}} .

For a typical poroelastic acoustic material with $K_{b} ≪ K_{s}$ and $K_{f} ≪ K_{s}$ , approximately $R = ϕ K_{f}, Q = (1 - ϕ) K_{f}, P - 2 G = λ + \frac{{(1 - ϕ)}^{2}}{ϕ} K_{f}$ , where $λ$ is Lame coefficient.

In a steady-state dynamic analysis the $K_{f}$ is usually complex and frequency dependent; therefore, other special physical experiments are required to determine $K_{f}$ while $K_{b}$ , $K_{s}$ , and G are purely real and frequency independent. Thus, the P, Q, and R are also generally complex and frequency dependent. The structural damping is accounted for by using only the structural stiffness part.

Inertial and damping dynamic properties and equilibrium equations

Instead of the usual structural-only equilibrium equations in the absence of body forces and damping

σ_{i j, j}^{s} = ρ \ddot{u_{i}},

Biot (1956) introduced

\begin{matrix} (10) & σ_{i j, j}^{s} = ρ_{11} {\ddot{u}}_{i}^{s} + ρ_{12} {\ddot{u}}_{i}^{f} + \tilde{b} ({\dot{u}}_{i}^{s} - {\dot{u}}_{i}^{f}), \end{matrix}

\begin{matrix} (11) & σ_{i j, j}^{f} = ρ_{12} {\ddot{u}}_{i}^{s} + ρ_{22} {\ddot{u}}_{i}^{f} - \tilde{b} ({\dot{u}}_{i}^{s} - {\dot{u}}_{i}^{f}) . \end{matrix}

In this section the tilde indicates that the physical property is complex and frequency dependent. The $u_{i}^{s}$ and $u_{i}^{f}$ are solid and fluid macroscopic displacements; that is, the average volume displacements per unit cross-section area. The $\tilde{b}$ terms represent viscous damping. The densities $ρ_{11}, ρ_{22},$ and $ρ_{12}$ are related to the skeleton density $ρ_{s}$ and fluid density $ρ_{f}$ by the relations (see Allard and Atalla (2009),

\begin{matrix} (12) & ρ_{11} = ρ_{s} - ρ_{12}, ρ_{22} = {ϕ ρ}_{f} - ρ_{12} . \end{matrix}

The skeleton density $ρ_{s}$ is $(1 - ϕ)$ times the skeleton material density. All these densities do not depend on frequency.

For harmonic motion $e^{j ω t}$ ,

j: is the imaginary unit,
$ω$: is the frequency,
t: is the time;

Atalla et al. (1998) define

\begin{matrix} (13) & {\tilde{ρ}}_{11} = ρ_{11} + \frac{\tilde{b}}{j ω}, \end{matrix}

\begin{matrix} (14) & {\tilde{ρ}}_{22} = ρ_{22} + \frac{\tilde{b}}{j ω}, \end{matrix}

\begin{matrix} (15) & {\tilde{ρ}}_{12} = ρ_{12} - \frac{\tilde{b}}{j ω}, \end{matrix}

so Equation 10 and Equation 11 become

\begin{matrix} (16) & σ_{i j, j}^{s} + ω^{2} {\tilde{ρ}}_{11} u_{i}^{s} + ω^{2} {\tilde{ρ}}_{12} u_{i}^{f} = 0, \end{matrix}

\begin{matrix} (17) & - ϕ p_{, i} + ω^{2} {\tilde{ρ}}_{12} u_{i}^{s} + ω^{2} {\tilde{ρ}}_{22} u_{i}^{f} = 0. \end{matrix}

From Equation 17

\begin{matrix} (18) & u_{i}^{f} = \frac{ϕ p_{, i}}{ω^{2} {\tilde{ρ}}_{22}} - \frac{{\tilde{ρ}}_{12}}{{\tilde{ρ}}_{22}} u_{i}^{s} . \end{matrix}

With Equation 18, Equation 16 becomes

\begin{matrix} (19) & σ_{i j, j}^{s} + ω^{2} \tilde{ρ} u_{i}^{s} + ϕ \frac{{\tilde{ρ}}_{12}}{{\tilde{ρ}}_{22}} p_{, i} = 0, \end{matrix}

\begin{matrix} (20) & \tilde{ρ} = {\tilde{ρ}}_{11} - \frac{{({\tilde{ρ}}_{12})}^{2}}{{\tilde{ρ}}_{22}} . \end{matrix}

To obtain a formulation with just structural displacements, $u_{i}^{s},$ and fluid pore pressure, $p,$ as primary variables, the fluid displacement, $u_{i}^{f},$ needs to be eliminated, but it is still present in the constitutive equations Equation 1 and Equation 2 through fluid dilatation $Θ^{f} = u_{k, k}^{f}$ :

\begin{matrix} (21) & σ_{i j}^{s} = [(\tilde{P} - 2 \tilde{G}) u_{k, k}^{s} + \tilde{Q} u_{k, k}^{f}] δ_{i j} + 2 \tilde{G} ε_{i j}^{s}, \end{matrix}

\begin{matrix} (22) & σ_{i j}^{f} = - ϕ p δ_{i j} = [\tilde{Q} u_{k, k}^{s} + \tilde{R} u_{k, k}^{f}] δ_{i j} . \end{matrix}

In Equation 21 and Equation 22 the tildes are properly used for clarity; the tildes are not shown in Equation 1 and Equation 2. This leads to

\begin{matrix} (23) & σ_{i j}^{s} = (\tilde{P} - 2 \tilde{G} - \frac{{\tilde{Q}}^{2}}{\tilde{R}}) u_{k, k}^{s} δ_{i j} + 2 \tilde{G} ε_{i j}^{s} - ϕ \frac{\tilde{Q}}{\tilde{R}} p δ_{i j} . \end{matrix}

Defining

{\hat{σ}}_{i j}^{s} = (\tilde{P} - 2 \tilde{G} - \frac{{\tilde{Q}}^{2}}{\tilde{R}}) u_{k, k}^{s} δ_{i j} + 2 \tilde{G} ε_{i j}^{s}

leads to

σ_{i j}^{s} = {\hat{σ}}_{i j}^{s} - ϕ \frac{\tilde{Q}}{\tilde{R}} p δ_{i j} .

Substituting this into Equation 19 gives

\begin{matrix} (24) & {\hat{σ}}_{i j, j}^{s} + ω^{2} \tilde{ρ} u_{i}^{s} + \tilde{γ} p_{, i} = 0, \end{matrix}

\tilde{γ} = ϕ (\frac{{\tilde{ρ}}_{12}}{{\tilde{ρ}}_{22}} - \frac{\tilde{Q}}{\tilde{R}}) .

It can be verified that

{\hat{σ}}_{i j}^{s} = (K_{b} - \frac{2}{3} G) u_{k, k}^{s} δ_{i j} + 2 G ε_{i j}^{s},

which is the classical constitutive equation for an elastic solid (with no tildes).

From Equation 18

u_{i, i}^{f} = \frac{ϕ Δ p}{ω^{2} {\tilde{ρ}}_{22}} - \frac{{\tilde{ρ}}_{11}}{{\tilde{ρ}}_{22}} Θ^{s}, Δ p = p_{, i i} .

Combining it with Equation 22 produces

\begin{matrix} (25) & Δ p + \frac{{\tilde{ρ}}_{22}}{\tilde{R}} ω^{2} p - \frac{{\tilde{ρ}}_{22}}{ϕ^{2}} \tilde{γ} ω^{2} Θ^{s} = 0. \end{matrix}

Without the coupling $\tilde{γ}$ terms, Equation 24 and Equation 25 are classical structural and fluid equilibrium equations. The coupling is volumetric.

According to Johnson et al. (1987), the ${\tilde{ρ}}_{22}$ is an effective fluid density in a rigid skeleton and can be represented as

\begin{matrix} (26) & {\tilde{ρ}}_{22} = ρ_{f} ϕ \tilde{α}, \end{matrix}

where the $\tilde{α}$ is frequency-dependent dynamic tortuosity. It is related to $\tilde{γ}$ :

\tilde{α} = \frac{{\tilde{ρ}}_{22}}{{ϕ ρ}_{f}} = \frac{{\tilde{ρ}}_{22}}{ρ_{22} + ρ_{12}} = \frac{{\tilde{ρ}}_{22}}{{\tilde{ρ}}_{22} + {\tilde{ρ}}_{12}},

(use Equation 12, Equation 14, and Equation 15)

\frac{ϕ}{\tilde{α}} = ϕ (1 + \frac{{\tilde{ρ}}_{12}}{{\tilde{ρ}}_{22}}) = ϕ (\frac{{\tilde{ρ}}_{12}}{{\tilde{ρ}}_{22}} - \frac{\tilde{Q}}{\tilde{R}}) + ϕ (1 + \frac{\tilde{Q}}{\tilde{R}}) = \tilde{γ} + ϕ (1 + \frac{\tilde{Q}}{\tilde{R}}) .

The limit $ω \to + \infty$ relation $ρ_{22} = ρ_{f} ϕ α_{\infty}$ , where $α_{\infty} \geq 1$ is the real-valued tortuosity at infinite frequency, is also valid. The $ρ_{f}, ϕ$ , and $α_{\infty}$ are specified by the user. Knowing the ${\tilde{ρ}}_{22}$ and $ρ_{22},$ the $\tilde{b}$ can be determined from Equation 14. With Equation 12, all coefficients for Equation 13, Equation 14, and Equation 15 are now determined.

In the Biot-Atalla material model the ${\tilde{ρ}}_{22}$ is a user-defined frequency-dependent complex density of acoustic medium, so the dynamic tortuosity $\tilde{α}$ can be obtained using Equation 26. The acoustic medium complex bulk modulus ${\tilde{K}}_{f}$ is similarly a user-defined property.

In the Biot-Johnson material model the $\tilde{α}$ and ${\tilde{K}}_{f}$ are computed as

\tilde{α} (ω) = α_{\infty} [1 + \frac{σ ϕ}{j ω ρ_{f} α_{\infty}} G_{J} (ω)], G_{J} (ω) = {(1 + \frac{4 j ω α_{\infty}^{2} η ρ_{f}}{σ^{2} Λ^{2} ϕ^{2}})}^{\frac{1}{2}},

{\tilde{K}}_{f} (ω) = \frac{γ P_{0}}{γ - (γ - 1) {[α' (ω)]}^{- 1}}, α' (ω) = 1 + \frac{8 η}{j ω \Pr {Λ'}^{2} ρ_{f}} {[1 + \frac{j ω \Pr {Λ'}^{2} ρ_{f}}{16 η}]}^{\frac{1}{2}} .

The following user-defined material properties are used:

$σ$: is the static flow resistivity,
$Λ$: is the viscous characteristic length,
$Λ'$: is the thermal characteristic length,
$η$: is the fluid dynamic viscosity,
$P_{0}$: is the ambient fluid pressure,
$γ$: is the heat capacity ratio, the dimensionless ratio of the specific heat per unit mass at constant pressure over the specific heat per unit mass at constant volume, and
$\Pr$: is the Prandtl number, the dimensionless ratio of viscous diffusion rate over thermal diffusion rate.

See Allard and Atalla (2009) for details.

The description above provides all required coefficients for the finite element weak integral $(u, p)$ formulation. Such a formulation, together with the compatibility boundary conditions (poroelastic-to-poroelastic, poroelastic-to-elastic, poroelastic-to-acoustic media), is presented in Atalla et al. (2001).