# Sudden Pressurization of a Spherical Cavity in a Poroelastic Medium

Authors: Mehmet Ozyazicioglu
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Governing equations of poroelastodynamics in time and frequency domain are derived. The continuity equation complements the momentum balance equations. After reduction for spherical symmetry (geometry and loading), the governing equations in frequency domain are solved by introducing wave potentials. The wave propagation velocities are obtained as the real parts of the characteristic equation of the coupled ODE system. Time domain solution for Dirac type boundary pressure is obtained through numerical inversion of transformed solutions. The results are compared to the solution in classical elasticity theory found in the literature.


---

## Body

## 1. Introduction

In wave mechanics, cavity expansion problems (circular and spherical) have received considerable attention since early 1930s [1]. These problems generally are amenable to exact solution, and the analytic solutions cast new light to the nature of wave propagation in solid media. The problem has practical applications in geophysics, seismology, and tunnel and mining engineering, as earthquake sources, underground detonation, and seismic probing. The exact solutions of such simple problems serve to understand more complex wave motions; these solutions can also be used as benchmark problems to assess the accuracy of numerical methods (FEM, BEM, FDM, etc.).The literature on cavity problems in elastodynamics abounds [2–6], while corresponding work in poroelasticity has only recently been emerging [7]. An analytical solution for a dynamically loaded poroelastic column [8, 9] is also available.This work concerns frequency domain analytical solution of dynamic pressurization of a spherical cavity in an infinite poroelastic medium. Time domain solution for Dirac type boundary pressure (pressure suddenly applied and removed) is also investigated using FFT (fast Fourier transform) technique. The developed solutions have parallel applications as the classical elastic cavity problem and clearly show the existence of a second pressure wave, the so-called slow wave.In the sequel, a short review of Biot’s poroelasticity theory is given. The governing equations are reduced for spherical symmetry. The analytical solution in frequency domain is obtained for a permeable boundary with arbitrary but spherically symmetric pressure inside the cavity. Time domain solution is obtained via a numerical inversion using FFT.


### 1.1. Biot’s Theory of Poroelasticity

Unlike the classical elasticity, Biot’s theory of poroelasticity is a coupled defromation flow theory of a porous solid matrix with interstitial fluid. Biot introduced his linear quasistatic theory in 1941 [10] and later extended it to cover the dynamic range [11, 12]. An extensive review of quasi-static poroelasticity can be found in [13, 14].The constitutive equations of linear-isotropic poroelasticity are given as (summation convention applies)(1a)τij=2μεij+λδijεkk-αδijp,(1b)θ=αεkk+1Qp,where εij=(1/2)(ui,j+uj,i) are the strains in the solid ui, and τij are components of solid displacement vector and total stress tensor, p is the fluid pressure, θ is the variation of fluid volume per unit reference volume, and δij is the Kronecker delta. Here, tensile τij and εij are positive, while pore pressure p is positive when it is compressive. Four elastic constants,λ: drained Lame’s modulus, as defined in classical elasticity (dimension ≡ F/L2);μ: drained shear modulus, as defined in classical elasticity (dimension ≡ F/L2);α: Biot’s effective stress coefficient (dimensionless);Q: Biot deformation modulus, corresponds to the reciprocal of constrained storage coefficient in hydrogeology (dimension ≡ F/L2),characterize the poroelastic material.The balance of linear momentum for the solid and fluid parts gives, respectively,(2)τji,j+fi=ρu¨i+nρfw¨i,(3)qi=-κ[p,i+ρfu¨i+nρf+ρanw¨i-gif].
Here, wi are relative fluid displacements, and qi=n·wi are the components of the fluid flux vector, the three densities ρ, ρf, and ρa are the bulk density, fluid density, and added fluid mass density, respectively. n is the porosity of the porous matrix. κ is the coefficient of permeability in units L3T/M. fi, gi are the body forces per unit volume of solid and fluid, respectively. The third equation is sometimes unconventionally referred to as generalized Darcy’s law. The last balance law is the continuity equation for the fluid part
(4)∂θ∂t+qk,k=a(t),
where a(t) is fluid generation rate per unit volume.In the absence of body forces and fluid generation, Fourier transforms of (1a), (1b), (2), (3), and (4) with respect to time give, when combined, the 3D governing equations in frequency domain
(5)μui,jjF+μuj,ijF+λuj,jiF-(α+β)p,iF=-ω2(ρ+βρf)uiF,-β1ω2ρfp,kkF+(α+β)uk,kF+1QpF=0.
It has been assumed that the medium is initially at rest. Here, pF, uF represent the Fourier transformed variables, that is,
(6)uF(ω)=∫-∞∞u(t)e-iωtdt⟷u(t)=12π∫-∞∞uF(ω)eiωtdω,
and β is a coefficient defined as
(7)β=n2κρfω2iωn2-ω2κ(ρa+nρf).
In dyadic form (5) reads
(8)μ∇2u_F+(λ+μ)∇∇·u_F-(α+β)∇pF+ω2(ρ+βρf)u_F=0_,-βω2ρf∇2pF+(α+β)∇·u_F+1QpF=0,
where u_ is the solid displacement vector. Analytical solutions of the 3D governing equations are scarse if not impossible. There are the point load solutions for infinite poroelastic medium, the fundamental solutions, for both quasi-static and dynamic poroelasticity [15–18] and poroelastic Lamb’s problem [19] available in the literature. Also a good compendium of anlytical solutions in the quasi-static range can be found in [14].


## 1.1. Biot’s Theory of Poroelasticity

Unlike the classical elasticity, Biot’s theory of poroelasticity is a coupled defromation flow theory of a porous solid matrix with interstitial fluid. Biot introduced his linear quasistatic theory in 1941 [10] and later extended it to cover the dynamic range [11, 12]. An extensive review of quasi-static poroelasticity can be found in [13, 14].The constitutive equations of linear-isotropic poroelasticity are given as (summation convention applies)(1a)τij=2μεij+λδijεkk-αδijp,(1b)θ=αεkk+1Qp,where εij=(1/2)(ui,j+uj,i) are the strains in the solid ui, and τij are components of solid displacement vector and total stress tensor, p is the fluid pressure, θ is the variation of fluid volume per unit reference volume, and δij is the Kronecker delta. Here, tensile τij and εij are positive, while pore pressure p is positive when it is compressive. Four elastic constants,λ: drained Lame’s modulus, as defined in classical elasticity (dimension ≡ F/L2);μ: drained shear modulus, as defined in classical elasticity (dimension ≡ F/L2);α: Biot’s effective stress coefficient (dimensionless);Q: Biot deformation modulus, corresponds to the reciprocal of constrained storage coefficient in hydrogeology (dimension ≡ F/L2),characterize the poroelastic material.The balance of linear momentum for the solid and fluid parts gives, respectively,(2)τji,j+fi=ρu¨i+nρfw¨i,(3)qi=-κ[p,i+ρfu¨i+nρf+ρanw¨i-gif].
Here, wi are relative fluid displacements, and qi=n·wi are the components of the fluid flux vector, the three densities ρ, ρf, and ρa are the bulk density, fluid density, and added fluid mass density, respectively. n is the porosity of the porous matrix. κ is the coefficient of permeability in units L3T/M. fi, gi are the body forces per unit volume of solid and fluid, respectively. The third equation is sometimes unconventionally referred to as generalized Darcy’s law. The last balance law is the continuity equation for the fluid part
(4)∂θ∂t+qk,k=a(t),
where a(t) is fluid generation rate per unit volume.In the absence of body forces and fluid generation, Fourier transforms of (1a), (1b), (2), (3), and (4) with respect to time give, when combined, the 3D governing equations in frequency domain
(5)μui,jjF+μuj,ijF+λuj,jiF-(α+β)p,iF=-ω2(ρ+βρf)uiF,-β1ω2ρfp,kkF+(α+β)uk,kF+1QpF=0.
It has been assumed that the medium is initially at rest. Here, pF, uF represent the Fourier transformed variables, that is,
(6)uF(ω)=∫-∞∞u(t)e-iωtdt⟷u(t)=12π∫-∞∞uF(ω)eiωtdω,
and β is a coefficient defined as
(7)β=n2κρfω2iωn2-ω2κ(ρa+nρf).
In dyadic form (5) reads
(8)μ∇2u_F+(λ+μ)∇∇·u_F-(α+β)∇pF+ω2(ρ+βρf)u_F=0_,-βω2ρf∇2pF+(α+β)∇·u_F+1QpF=0,
where u_ is the solid displacement vector. Analytical solutions of the 3D governing equations are scarse if not impossible. There are the point load solutions for infinite poroelastic medium, the fundamental solutions, for both quasi-static and dynamic poroelasticity [15–18] and poroelastic Lamb’s problem [19] available in the literature. Also a good compendium of anlytical solutions in the quasi-static range can be found in [14].


## 2. Governing Equations for Spherical Symmetry

Here, a spherical cavity of radius “a” in an infinite PE medium (Figure 1) will be considered. By substituting the gradient and divergence in (8), governing equations in spherical coordinates can be obtained.Figure 1
Description of problem, pressurized spherical cavity in infinite poroelastic medium.Letur, uϕ, uθ be the spherical components of solid displacements, because the cavity is spherically symmetric, and we assume a time varying but spherically symmetric pressure inside; waves emanating from such a source will have spherical symmetry, that is,
(9)ur=u(r,t)≠0,uθ=uϕ=∂∂θ(·)=∂∂ϕ(·)=0.
Hence, the governing equations of 3D poroelasticity in FTS, in this case, reduce to the following:
(10)(λ+2μ)[1r2ddr(r2duFdr)-2r2uF]-(α+β)dpFdr+ω2(ρ+βρf)uF=0,-βω2ρf1r2ddr(r2dpFdr)+(α+β)1r2ddr(r2uF)+pFQ=0.
Here, u denotes the only nonzero displacement for simplicity, that is, radial component ur.The stress components are related to the radial displacement as(11)τrr=(λ+2μ)∂u∂r+2λur-αp,τθθ=τϕϕ=(λ+2μ)ur+λ[∂u∂r+ur]-αp.
Other τij=0.


## 3. General Solution of the Governing Equations

If one introduces a displacement potentialΦ such that ur=∂Φ/∂r then (10) and (11) become
(12)d2dr2(ΨF)-α+βλ+2μ(χF)-λ42(ΨF)=0,d2dr2(χF)-ω2ρf(α+β)βd2dr2(ΨF)-ω2ρfβQ(χF)=0,
where
(13)Ψ=rΦ,χ=rp,λ42=-ω2(ρ+βρf)λ+2μ.
Now, we assume solutions of the form
(14)ΨF=Ψ0eimωr,χF=χ0eimωr,
where m is a frequency dependent parameter. Substituting (14) in (12), one obtains the eigenvalue problem:
(15)[-m2ω2-λ42-α+βλ+2μm2ω4ρf(α+β)β-ω2ρfβQ-m2ω2][Ψ0χ0]=[00].
Hence, ω are the roots of the characteristic equation
(16)(m2ω2+λ42)(ω2ρfβQ+m2ω2)+α+βλ+2μm2ω4ρf(α+β)β=0.
With the following abbreviations:
(17)B=ρ+βρfλ+2μ-ρfβQ-ρf(α+β)2β(λ+2μ),C=ρf(ρ+βρf)βQ(λ+2μ),
the roots of the characteristic equation are
(18)m1=12[B-B2+4C],m2=12[B+B2+4C],m3=-m1,m4=-m2,
and the solution becomes
(19)ΨF=Ψ1e-im1ωr+Ψ2e-im2ωr+Ψ3eim1ωr+Ψ4eim2ωr,χF=χ1e-im1ωr+χ2e-im2ωr+χ3eim1ωr+χ4eim2ωr.
It is clear that the terms eimkωr represent waves propagating inwards. In an infinite medium, no waves propagate from infinity (Sommerfeld’s radiation condition [1]); discarding these terms, (19) becomes
(20)ΨF=Ψ1e-im1ωr+Ψ2e-im2ωr,χF=χ1e-im1ωr+χ2e-im2ωr.
Moreover, the amplitudes Ψk, χk are related through (15) as
(21)χk=-(λ+2μ)(mk2ω2+λ42)α+βΨk;
introducing (21) in (20)
(22)ΨF=Ψ1e-im1ωr+Ψ2e-im2ωr,χF=-λ+2μα+β[Ψ1(m12ω2+λ42)e-im1ωr+Ψ2(m22ω2+λ42)e-im2ωr].
Furthermore, since Ψ=rΦ and ur=∂Φ/∂r(23)ΦF=Ψ1e-im1ωrr+Ψ2e-im2ωrr.
Finally, the general solution for radial displacement and pore pressure becomes
(24)uF=-Ψ1[1r+iωm1]e-iωm1rr-Ψ2[1r+iωm2]e-iωm2rr,pF=-λ+2μα+β[Ψ1(m12ω2+λ42)e-iωm1rr+Ψ2(m22ω2+λ42)e-iωm2rr].
The integration constants Ψ1 and Ψ2 are to be evaluated from boundary conditions at r=a.


### 3.1. Boundary Conditions

Among various possible combinations of traction, displacement, pore pressure, and fluid flux, only the permeable boundary condition will be considered here, that is,(25)τrr|r=a=-S0,p=0.
The stress boundary condition corresponds to a loading of σn=S0δ(t) in time domain, where S0 is the strength of loading and δ(t) is the Dirac delta function (or distribution). When (25) are introduced in (24) one obtains
(26)Ψ1=-aS0Δ(m22ω2+λ42)eiωm1a,Ψ2=aS0Δ(m12ω2+λ42)eiωm2a,
where
(27)Δ=[4μa(iωm1+1a)-(λ+2μ)ω2m1](m22ω2+λ42)-[4μa(iωm2+1a)-(λ+2μ)ω2m2](m12ω2+λ42).


### 3.2. Analytical Solution in Frequency Domain

Thus, the complete analytical solutions in frequency domain for Dirac type pulse are(28)uF=aS0Δ[(m22ω2+λ42)[1r+iωm1]e-iωm1(r-a)r-(m12ω2+λ42)[1r+iωm2]e-iωm2(r-a)r],pF=λ+2μα+βaS0Δ[(m22ω2+λ42)(m12ω2+λ42)e-iωm1(r-a)r-Ψ2(m12ω2+λ42)e-iωm2(r-a)r].
Stresses can be found by inserting (28) in (11). The two terms in each of (28) represent two p-phases with different propagation characteristics. The reciprocal of real part of m1 and m2 corresponds to fast and slow (second) longitudinal wave velocities. The second (slow) p-wave is highly dispersive and highly damped.


### 3.3. Solution in the Classical Theory of Elasticity

The elastodynamic solution in frequency domain for suddenly pressurized spherical cavity can be derived following a similar outline and is found in [1, 20]; here we restate the form of the solution, which complies with the Fourier transform definition in (6):
(29)ueF(ω,r)=a3S04μr(1+i(ω/c1)a+(i(ω/c2)(a/2))2)×(1r+iωc1)e-i(ω/c1)(r-a),
where ue is the displacement, c1=(λ+2μ)/ρ, c2=μ/ρ are the p- and s-wave velocities in an idealized elastic medium, and other parameters are as defined previously.


## 3.1. Boundary Conditions

Among various possible combinations of traction, displacement, pore pressure, and fluid flux, only the permeable boundary condition will be considered here, that is,(25)τrr|r=a=-S0,p=0.
The stress boundary condition corresponds to a loading of σn=S0δ(t) in time domain, where S0 is the strength of loading and δ(t) is the Dirac delta function (or distribution). When (25) are introduced in (24) one obtains
(26)Ψ1=-aS0Δ(m22ω2+λ42)eiωm1a,Ψ2=aS0Δ(m12ω2+λ42)eiωm2a,
where
(27)Δ=[4μa(iωm1+1a)-(λ+2μ)ω2m1](m22ω2+λ42)-[4μa(iωm2+1a)-(λ+2μ)ω2m2](m12ω2+λ42).


## 3.2. Analytical Solution in Frequency Domain

Thus, the complete analytical solutions in frequency domain for Dirac type pulse are(28)uF=aS0Δ[(m22ω2+λ42)[1r+iωm1]e-iωm1(r-a)r-(m12ω2+λ42)[1r+iωm2]e-iωm2(r-a)r],pF=λ+2μα+βaS0Δ[(m22ω2+λ42)(m12ω2+λ42)e-iωm1(r-a)r-Ψ2(m12ω2+λ42)e-iωm2(r-a)r].
Stresses can be found by inserting (28) in (11). The two terms in each of (28) represent two p-phases with different propagation characteristics. The reciprocal of real part of m1 and m2 corresponds to fast and slow (second) longitudinal wave velocities. The second (slow) p-wave is highly dispersive and highly damped.


## 3.3. Solution in the Classical Theory of Elasticity

The elastodynamic solution in frequency domain for suddenly pressurized spherical cavity can be derived following a similar outline and is found in [1, 20]; here we restate the form of the solution, which complies with the Fourier transform definition in (6):
(29)ueF(ω,r)=a3S04μr(1+i(ω/c1)a+(i(ω/c2)(a/2))2)×(1r+iωc1)e-i(ω/c1)(r-a),
where ue is the displacement, c1=(λ+2μ)/ρ, c2=μ/ρ are the p- and s-wave velocities in an idealized elastic medium, and other parameters are as defined previously.


## 4. Results in Frequency and Time Domain

In all of the computations material constants (Table1) for Berea Sandstone are used [21].Table 1
Material data for Berea sand stone (poroelastic medium).
n
α
Q (Pa)
μ (Pa)
ν
κ (m4/N/s)
ρ (kg/m3)
ρ
f  (kg/m3)
ρ
a  (kg/m3)
0.19
0.778
1.353 * 1010
6 * 109
0.2
1.9 * 10−10
2458
1000
125.4The material parameters for corresponding elastic medium are as in Table2.Table 2
Material data for elastic medium.
μ (Pa)
λ (Pa)
ν
ρ  (kg/m3)
6 * 109
4 * 109
0.2
2458The real and complex parts of displacement and pore pressure solutions in frequency domain are plotted in Figures2 and 3, respectively, note that the elastic theory doesn’t involve pore pressure.Frequency domain solution for nondimensional displacements,E represents μ.
(a)
(b)Frequency domain solution for nondimensional pore pressure.
(a)
(b)The following (Figure4) are the dispersion curves for the fast and slow wave. As seen in the figure, the wave front of fast p-wave travels with a speed equal to 3137.2 m/s where as the slow wave front velocity is 1037 m/s. While, both waves display dispersion, the slow wave is rather more dispersive than the fast wave.Dispersion curves for fast and slow waves (real parts ofm1, m2).
(a)
(b)


### 4.1. Time Domain Solution Using FFT Technique

The frequency domain analytical solutions (28) are extremely complex (if not impossible) to be inverse transformed (by analytical means) into time domain. Thus, time domain solutions are obtained using numerically inverting (28) by FFT (fast Fourier transform) technique.Since Dirac impulse is highly transient, a high frequency resolution is required, thus the following parameters are chosen for FFT:T=50 seconds (analysis interval), M=17 (exponent in FFT),which give a frequency resolution of 0.628 rad/s and a Nyquist frequency of 41177 rad/sec. The results are plotted in Figure5.Figure 5
Displacement wave propagating from a spherical cavity surface (Dirac type pressure).As seen from Figure5, the fast primary wave in poroelastic medium travels faster (arrives earlier) than p-wave in elastic medium, with the amplitude of elastic wave being slightly larger. The second slow wave in poroelastic wave cannot have been captured, probably because of its small amplitude and high damping and dispersion.


## 4.1. Time Domain Solution Using FFT Technique

The frequency domain analytical solutions (28) are extremely complex (if not impossible) to be inverse transformed (by analytical means) into time domain. Thus, time domain solutions are obtained using numerically inverting (28) by FFT (fast Fourier transform) technique.Since Dirac impulse is highly transient, a high frequency resolution is required, thus the following parameters are chosen for FFT:T=50 seconds (analysis interval), M=17 (exponent in FFT),which give a frequency resolution of 0.628 rad/s and a Nyquist frequency of 41177 rad/sec. The results are plotted in Figure5.Figure 5
Displacement wave propagating from a spherical cavity surface (Dirac type pressure).As seen from Figure5, the fast primary wave in poroelastic medium travels faster (arrives earlier) than p-wave in elastic medium, with the amplitude of elastic wave being slightly larger. The second slow wave in poroelastic wave cannot have been captured, probably because of its small amplitude and high damping and dispersion.


## 5. Conclusions

Analytical solution of waves propagating from the surface of a suddenly pressurized poroelastic spherical cavity with permeable boundary is obtained in frequency domain. The problem is one of the classical problems in wave propagation in elastic solids [1, 20]. Among many other possibilities of boundary conditions, this paper, which employs permeable boundary conditions, is one of the first attempts to the solution of a variety of other problems related to spherical wave propagation from cavities in poroelastic media.Time domain solution is obtained with inverse FFT. Since, Dirac type pulse is a highly transient excitation, a high frequency resolution is necessary. The results clearly show thatp-waves in poroelastic media propagate faster, with slightly lower amplitude, as compared to elastic media with similar material properties. The slow wave could not have been visualized in the graphics, because of small amplitude, high damping, and dispersion characteristics of this phase.



---

*Source: 632634-2013-03-28.xml*