Sudden Pressurization of a Spherical Cavity in a Poroelastic Medium

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.

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)𝜏_𝑖𝑗=2𝜇𝜀_𝑖𝑗+𝜆𝛿_𝑖𝑗𝜀_𝑘𝑘−𝛼𝛿_𝑖𝑗𝑝,𝜃=𝛼𝜀_𝑘𝑘+1𝑄𝑝, (1a)(1b)where 𝜀_𝑖𝑗=(1/2)(𝑢_𝑖,𝑗+𝑢_𝑗,𝑖) are the strains in the solid 𝑢_𝑖, and 𝜏_𝑖𝑗 are components of solid displacement vector and total stress tensor, 𝑝 is the fluid pressure, 𝜃 is the variation of fluid volume per unit reference volume, and 𝛿_𝑖𝑗 is the Kronecker delta. Here, tensile 𝜏_𝑖𝑗 and 𝜀_𝑖𝑗 are positive, while pore pressure 𝑝 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); 𝑄: 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)(3) Here, 𝑤_𝑖 are relative fluid displacements, and 𝑞_𝑖=𝑛⋅𝑤_𝑖 are the components of the fluid flux vector, the three densities 𝜌, 𝜌_𝑓, and 𝜌_𝑎 are the bulk density, fluid density, and added fluid mass density, respectively. 𝑛 is the porosity of the porous matrix. 𝜅 is the coefficient of permeability in units L³T/M. 𝑓_𝑖, 𝑔_𝑖 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) where 𝑎(𝑡) 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 𝜇𝑢^𝐹_{𝑖,𝑗𝑗}+𝜇𝑢^𝐹_{𝑗,𝑖𝑗}+𝜆𝑢^𝐹_{𝑗,𝑗𝑖}−𝛼+𝛽𝑝^𝐹_,𝑖=−𝜔²𝜌+𝛽𝜌_𝑓𝑢^𝐹_𝑖,1−𝛽𝜔²𝜌_𝑓𝑝^𝐹_,𝑘𝑘+𝛼+𝛽𝑢^𝐹_𝑘,𝑘+1𝑄𝑝^𝐹=0. (5) It has been assumed that the medium is initially at rest. Here, 𝑝^𝐹, 𝑢^𝐹 represent the Fourier transformed variables, that is, 𝑢^𝐹(𝜔)=^∞_−∞𝑢(𝑡)𝑒^{−𝑖𝜔𝑡}(𝑡)=1𝑑𝑡⟷𝑢2𝜋^∞_−∞𝑢^𝐹(𝜔)𝑒^𝑖𝜔𝑡𝑑𝜔, (6) and 𝛽 is a coefficient defined as 𝑛𝛽=²𝜅𝜌_𝑓𝜔²𝑖𝜔𝑛²−𝜔²𝜅𝜌_𝑎+𝑛𝜌_𝑓. (7) In dyadic form (5) reads 𝜇∇²𝐮^𝐹+𝜆+𝜇∇∇⋅𝐮^𝐹−𝛼+𝛽∇𝑝^𝐹+𝜔²𝜌+𝛽𝜌_𝑓𝐮^𝐹=𝟎,−𝛽𝜔²𝜌_𝑓∇²𝑝^𝐹+𝛼+𝛽∇⋅𝐮^𝐹+1𝑄𝑝^𝐹=0, (8) where 𝐮 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 “𝑎” 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.

Let 𝑢_𝑟, 𝑢_𝜙, 𝑢_𝜃 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, 𝑢_𝑟()𝑢=𝑢𝑟,𝑡≠0,_𝜃=𝑢_𝜙=𝜕(⋅)=𝜕𝜕𝜃(⋅)𝜕𝜙=0. (9) Hence, the governing equations of 3D poroelasticity in FTS, in this case, reduce to the following: 1𝜆+2𝜇𝑟²𝑑𝑟𝑑𝑟²𝑑𝑢^𝐹−2𝑑𝑟𝑟²𝑢^𝐹−𝛼+𝛽𝑑𝑝^𝐹𝑑𝑟+𝜔²𝜌+𝛽𝜌_𝑓𝑢^𝐹−𝛽=0,𝜔²𝜌_𝑓1𝑟²𝑑𝑟𝑑𝑟²𝑑𝑝^𝐹+1𝑑𝑟𝛼+𝛽𝑟²𝑑𝑟𝑑𝑟²𝑢^𝐹+𝑝^𝐹𝑄=0. (10) Here, 𝑢 denotes the only nonzero displacement for simplicity, that is, radial component 𝑢_𝑟.

The stress components are related to the radial displacement as 𝜏_𝑟𝑟=𝜆+2𝜇𝜕𝑢𝑢𝜕𝑟+2𝜆𝑟𝜏−𝛼𝑝,_𝜃𝜃=𝜏_𝜙𝜙=𝑢𝜆+2𝜇𝑟+𝜆𝜕𝑢+𝑢𝜕𝑟𝑟−𝛼𝑝. (11) Other 𝜏_𝑖𝑗=0.

3. General Solution of the Governing Equations

If one introduces a displacement potential Φ such that 𝑢_𝑟=𝜕Φ/𝜕𝑟 then (10) and (11) become 𝑑²𝑑𝑟²Ψ^𝐹−𝛼+𝛽𝜒𝜆+2𝜇^𝐹−𝜆²₄Ψ^𝐹𝑑=0,₂𝑑𝑟²𝜒^𝐹−𝜔²𝜌_𝑓𝛼+𝛽𝛽𝑑²𝑑𝑟²Ψ^𝐹−𝜔²𝜌_𝑓𝜒𝛽𝑄^𝐹=0, (12) where 𝜆Ψ=𝑟Φ,𝜒=𝑟𝑝,₂₄𝜔=−²𝜌+𝛽𝜌_𝑓.𝜆+2𝜇 (13) Now, we assume solutions of the form Ψ^𝐹=Ψ₀𝑒^{𝑖𝑚𝜔𝑟},𝜒^𝐹=𝜒₀𝑒^{𝑖𝑚𝜔𝑟}, (14) where 𝑚 is a frequency dependent parameter. Substituting (14) in (12), one obtains the eigenvalue problem: −𝑚²𝜔²−𝜆²₄−𝛼+𝛽𝑚𝜆+2𝜇₂𝜔⁴𝜌_𝑓𝛼+𝛽𝛽−𝜔²𝜌_𝑓𝛽𝑄−𝑚²𝜔²Ψ₀𝜒₀=00. (15) Hence, 𝜔 are the roots of the characteristic equation 𝑚²𝜔²+𝜆²₄𝜔²𝜌_𝑓𝛽𝑄+𝑚²𝜔²+𝛼+𝛽𝑚𝜆+2𝜇²𝜔⁴𝜌_𝑓𝛼+𝛽𝛽=0. (16) With the following abbreviations: 𝐵=𝜌+𝛽𝜌_𝑓−𝜌𝜆+2𝜇^𝑓−𝜌𝛽𝑄^𝑓𝛼+𝛽²𝛽,𝜌𝜆+2𝜇𝐶=^𝑓𝜌+𝛽𝜌_𝑓,𝛽𝑄𝜆+2𝜇 (17) the roots of the characteristic equation are 𝑚₁=12𝐵−𝐵²,𝑚+4𝐶₂=12𝐵+𝐵²,𝑚+4𝐶₃=−𝑚₁,𝑚₄=−𝑚₂, (18) and the solution becomes Ψ^𝐹=Ψ₁𝑒^{−𝑖𝑚1𝜔𝑟}+Ψ₂𝑒^{−𝑖𝑚2𝜔𝑟}+Ψ₃𝑒^{𝑖𝑚1𝜔𝑟}+Ψ₄𝑒^{𝑖𝑚2𝜔𝑟},𝜒^𝐹=𝜒₁𝑒^{−𝑖𝑚1𝜔𝑟}+𝜒₂𝑒^{−𝑖𝑚2𝜔𝑟}+𝜒₃𝑒^{𝑖𝑚1𝜔𝑟}+𝜒₄𝑒^{𝑖𝑚2𝜔𝑟}. (19) It is clear that the terms 𝑒^{𝑖𝑚𝑘𝜔𝑟} represent waves propagating inwards. In an infinite medium, no waves propagate from infinity (Sommerfeld’s radiation condition [1]); discarding these terms, (19) becomes Ψ^𝐹=Ψ₁𝑒^{−𝑖𝑚1𝜔𝑟}+Ψ₂𝑒^{−𝑖𝑚2𝜔𝑟},𝜒^𝐹=𝜒₁𝑒^{−𝑖𝑚1𝜔𝑟}+𝜒₂𝑒^{−𝑖𝑚2𝜔𝑟}. (20) Moreover, the amplitudes Ψ_𝑘, 𝜒_𝑘 are related through (15) as 𝜒_𝑘=−𝜆+2𝜇𝑚²_𝑘𝜔²+𝜆²₄Ψ𝛼+𝛽^𝑘; (21) introducing (21) in (20) Ψ^𝐹=Ψ₁𝑒^{−𝑖𝑚1𝜔𝑟}+Ψ₂𝑒^{−𝑖𝑚2𝜔𝑟},𝜒^𝐹=−𝜆+2𝜇Ψ𝛼+𝛽¹𝑚²₁𝜔²+𝜆²₄𝑒^{−𝑖𝑚1𝜔𝑟}+Ψ₂𝑚²₂𝜔²+𝜆²₄𝑒^{−𝑖𝑚2𝜔𝑟}. (22) Furthermore, since Ψ=𝑟Φ and 𝑢_𝑟=𝜕Φ/𝜕𝑟Φ^𝐹=Ψ₁𝑒^{−𝑖𝑚1𝜔𝑟}𝑟+Ψ₂𝑒^{−𝑖𝑚2𝜔𝑟}𝑟. (23) Finally, the general solution for radial displacement and pore pressure becomes 𝑢^𝐹=−Ψ₁1𝑟+𝑖𝜔𝑚₁𝑒^{−𝑖𝜔𝑚1𝑟}𝑟−Ψ₂1𝑟+𝑖𝜔𝑚₂𝑒^{−𝑖𝜔𝑚2𝑟}𝑟,𝑝^𝐹=−𝜆+2𝜇Ψ𝛼+𝛽¹𝑚²₁𝜔²+𝜆²₄𝑒^{−𝑖𝜔𝑚1𝑟}𝑟+Ψ₂𝑚²₂𝜔²+𝜆²₄𝑒^{−𝑖𝜔𝑚2𝑟}𝑟. (24) The integration constants Ψ₁ and Ψ₂ are to be evaluated from boundary conditions at 𝑟=𝑎.

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, 𝜏_{𝑟𝑟𝑟=𝑎}=−𝑆₀,𝑝=0. (25) The stress boundary condition corresponds to a loading of 𝜎_𝑛=𝑆₀𝛿(𝑡) in time domain, where 𝑆₀ is the strength of loading and 𝛿(𝑡) is the Dirac delta function (or distribution). When (25) are introduced in (24) one obtains Ψ₁=−𝑎𝑆₀Δ𝑚²₂𝜔²+𝜆²₄𝑒^{𝑖𝜔𝑚1𝑎},Ψ₂=𝑎𝑆₀Δ𝑚²₁𝜔²+𝜆²₄𝑒^{𝑖𝜔𝑚2𝑎}, (26) where Δ=4𝜇𝑎𝑖𝜔𝑚₁+1𝑎−𝜆+2𝜇𝜔²𝑚₁𝑚²₂𝜔²+𝜆²₄−4𝜇𝑎𝑖𝜔𝑚₂+1𝑎−𝜆+2𝜇𝜔²𝑚₂𝑚²₁𝜔²+𝜆²₄. (27)

3.2. Analytical Solution in Frequency Domain

Thus, the complete analytical solutions in frequency domain for Dirac type pulse are 𝑢^𝐹=𝑎𝑆₀Δ𝑚²₂𝜔²+𝜆²₄1𝑟+𝑖𝜔𝑚₁𝑒^{−𝑖𝜔𝑚1(𝑟−𝑎)}𝑟−𝑚²₁𝜔²+𝜆²₄1𝑟+𝑖𝜔𝑚₂𝑒^{−𝑖𝜔𝑚2(𝑟−𝑎)}𝑟,𝑝^𝐹=𝜆+2𝜇𝛼+𝛽𝑎𝑆₀Δ𝑚²₂𝜔²+𝜆²₄𝑚²₁𝜔²+𝜆²₄𝑒^{−𝑖𝜔𝑚1(𝑟−𝑎)}𝑟−Ψ₂𝑚²₁𝜔²+𝜆²₄𝑒^{−𝑖𝜔𝑚2(𝑟−𝑎)}𝑟. (28) Stresses can be found by inserting (28) in (11). The two terms in each of (28) represent two 𝑝-phases with different propagation characteristics. The reciprocal of real part of 𝑚₁ and 𝑚₂ corresponds to fast and slow (second) longitudinal wave velocities. The second (slow) 𝑝-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): 𝑢^𝐹_𝑒()=𝑎𝜔,𝑟³𝑆₀4𝜇𝑟1+𝑖𝜔/𝑐₁𝑎+𝑖𝜔/𝑐₂(𝑎/2)²×1𝑟𝜔+𝑖𝑐₁𝑒^{−𝑖(𝜔/𝑐1)(𝑟−𝑎)}, (29) where 𝑢_𝑒 is the displacement, 𝑐₁=(𝜆+2𝜇)/𝜌, 𝑐₂=𝜇/𝜌 are the 𝑝- and 𝑠-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 (Table 1) for Berea Sandstone are used [21].

Table 1 

Material data for Berea sand stone (poroelastic medium).

The material parameters for corresponding elastic medium are as in Table 2.

Table 2 

Material data for elastic medium.

The real and complex parts of displacement and pore pressure solutions in frequency domain are plotted in Figures 2 and 3, respectively, note that the elastic theory doesn’t involve pore pressure.

(a)

(b)

Figure 2 

Frequency domain solution for nondimensional displacements, 𝐸 represents 𝜇.

(a)

(b)

Figure 3 

Frequency domain solution for nondimensional pore pressure.

The following (Figure 4) are the dispersion curves for the fast and slow wave. As seen in the figure, the wave front of fast 𝑝-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.

(a)

(b)

Figure 4 

Dispersion curves for fast and slow waves (real parts of 𝑚1, 𝑚2).

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: 𝑇=50 seconds (analysis interval), 𝑀=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 Figure 5.

Figure 5 

Displacement wave propagating from a spherical cavity surface (Dirac type pressure).

As seen from Figure 5, the fast primary wave in poroelastic medium travels faster (arrives earlier) than 𝑝-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 that 𝑝-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.