Poroelastic Response of a Borehole (FLAC2D)

Problem Statement

Note

The project file for this example is available to be viewed/run in FLAC2D.[1] The main data files used are shown at the end of this example. The remaining data files can be found in the project.

A borehole (see Figure 1) is excavated in a saturated porous rock subject to an anisotropic in-situ stress field with isotropic component \(P_0\) and deviator \(S_0\). The borehole boundary is free to drain and is exposed to atmospheric pressure. The initial pore pressure field is \(p_0\). The problem is analyzed assuming plane-strain conditions and instantaneous drilling of the borehole. The objective of the FLAC2D simulation is to capture the poroelastic effects taking place during the short time response of the system. This problem provides a validation test for the FLAC2D simulation of two-dimensional coupled hydraulic-mechanical processes.

../../../../../_images/problem.png

Figure 1: Problem definition.

For this problem, the initial conditions are characterized by an in-situ pore pressure (\(p_0\)) = 1 MPa, an in-situ isotropic stress (\(P_0\)) = 3 MPa, and an in-situ deviatoric stress (\(S_0\)) = 1 MPa. The following material properties are prescribed.

shear modulus (\(G\))

1.5 × 109 Pa

bulk modulus (\(K\))

1.5 × 109 Pa

porosity (\(n\))

0.3

Biot coefficient (\(α\))

0.65

Biot modulus, fluid (\(M\))

2.0 × 109 Pa

permeability (\(k\))

10-12 (m/sec)/(Pa/m)

Closed-Form Solution

The two-dimensional poroelastic solution for a borehole in a non-hydrostatic stress field can be found in Detournay and Cheng (1988). The general solution is derived in the Laplace transform domain. Results in the time domain can be obtained using a numerical inversion technique. The short-time asymptotic solution for the region near the borehole is also presented. This solution is formulated by superposition of asymptotic solutions for three loading modes: (1) a far-field isotropic stress; (2) an initial pore-pressure distribution; and (3) a far-field stress deviator. The stresses (\(σ_{rr}, , σ_{θθ} , σ_{rθ}\)), pore pressure (\(p\)) and displacements (\(u_r,u_θ\)) induced by each loading mode are defined in the equations below. As shown in Figure 1, \(r\) and \(θ\) are the polar coordinates, and a is the borehole radius.

(1) Loading Mode 1: Far-Field Isotropic Stress

This loading mode produces the classical Lamé solution in elasticity. The rock deformation is entirely associated with the deviatoric strain, and there is no mechanism for pore pressure generation. Also, the stress field and displacements are independent of the bulk modulus of the rock and of time.

\[\eqalignno{ \sigma_{rr} &=\ - P_0 (1 - {{a^{2}}\over{r^{2}}}) \cr \sigma_{\theta \theta} &=\ - P_0 (1 + {{a^{2}}\over{r^{2}}}) \cr \cr \sigma_{r \theta} &=\ 0 \cr \ \ \ & &\ \ \cr u_r &=\ - {{P_0}\over{2G}} {{a^{2}}\over{r}} \cr \cr u_{\theta} &=\ 0 \cr \cr p &=\ 0 }\]

(2) Loading Mode 2: Initial Pore Pressure Distribution

The short-time asymptotic solution for stresses, pore pressure and displacements are given:

\[ \begin{align}\begin{aligned}\eqalign{\sigma_{rr} &= -2 \eta p_0 \left\{{\left({a}\over{r}\right)}^{{3}\over{2}} \left[ \sqrt{{4 c t}\over{a^2 \pi}} \exp \left( -{{{(r-a)}^2}\over{4 c t}} \right) - \left( {{r}\over{a}} - 1 \right) \hbox{erfc} \left( {{r - a}\over{\sqrt{4 c t}}} \right) \right] - {{a^2}\over{r^2}} \sqrt{{4 c t}\over{a^2 \pi}} \right\} \cr\\\sigma_{\theta \theta } &= 2 \eta p_0 \left\{{\left({a}\over{r}\right)}^{{3}\over{2}} \left[ \sqrt{{4 c t}\over{a^2 \pi}} \exp \left( -{{{(r-a)}^2}\over{4 c t}} \right) + \hbox{erfc} \left( {{r - a}\over{\sqrt{4 c t}}} \right) \right] - {{a^2}\over{r^2}} \sqrt{{4 c t}\over{a^2 \pi}} \right\}\\&+\ \ \ \ \ \ \ 2 \eta p_0 \left\{{{{1}\over{8}} \sqrt{{{a}\over{r}}} \left( 1 - {{a}\over{r}} \right) \left[ \sqrt{{4 c t}\over{a^2 \pi}} \ \exp \left( -{{{(r-a)}^2}\over{4 c t}} \right) - \left( {{r}\over{a}} - 1 \right) \hbox{erfc} \left( {{r - a}\over{\sqrt{4 c t}}} \right) \right] } \right\} \cr\\\sigma_{r \theta} &= 0 \cr\\u_r &= 2 \eta {{a p_0}\over{2G}} \left\{\sqrt{{a}\over{r}} \left[ \sqrt{{4 c t}\over{a^2 \pi}} \ \exp \left( -{{{(r-a)}^2}\over{4 c t}} \right) - \left( {{r}\over{a}} - 1 \right) \hbox{erfc} \left( {{r - a}\over{\sqrt{4 c t}}} \right) \right] + {{a}\over{r}} \sqrt{{4 c t}\over{a^2 \pi}} \right\} \cr\\u_{\theta} &= 0 \cr\\p\over{p_0} &= 1 - \sqrt{{{a}\over{r}}} \ \hbox{erfc} \left( {{r - a}\over{\sqrt{4 c t}}} \right) - {{1}\over{8}} \sqrt{{{a}\over{r}}} \left( 1 - {{a}\over{r}} \right)\\&\times\ \ \ \ \ \ \ \left[ \sqrt{{4 c t}\over{a^2 \pi}} \ \exp \left( -{{{(r-a)}^2}\over{4 c t}} \right) - \left( {{r}\over{a}} - 1 \right) \hbox{erfc} \left( {{r - a}\over{\sqrt{4 c t}}} \right) \right]\\}\end{aligned}\end{align} \]

The poroelastic coefficient, \(η\), and the diffusivity, \(c\), involved in the above equations are defined as

\[\eta = {{\alpha (1 - 2 \nu_u)} \over {2 (1 - \nu_u)}}\]

and

\[c = k / ({1 \over M} + {{\alpha^2}\over{{K + {{4}\over{3}} G} }})\]

where \(k\) is the mobility coefficient (FLAC2D permeability), \(K\) and \(G\) are the drained bulk and shear moduli of the medium, \(M\) is the Biot modulus, \(ν\) is the drained Poisson’s ratio, and \(ν_u\) is the undrained Poisson’s ratio. The undrained Poisson’s ratio, \(ν_u\), is defined as

\[\nu_u = {{3K_u-2G}\over{2(3K_u+G)}}\]

where the undrained bulk modulus, \(K_u\), is

\[K_u = \alpha^2 M + K\]

The mode 2 solution indicates that a maximum tensile stress \(σ_θθ=2ƞp_0\) is reached at the borehole wall.

(3) Loading Mode 3: Far-Field Deviator

The mode 3 solution is for deviatoric loading. The following short-time asymptotic equations for stresses, pore pressure and displacements are obtained for this mode:

\[\eqalignno{ \sigma_{rr} &=\ S_0 \left[ 1 - 4 {{a^2}\over{r^2}} +3 {{a^4}\over{r^4}} \right] \cos 2 \theta \cr \cr \cr \sigma_{\theta \theta} &=\ S_0 \left[ -1 + 4 {{\nu_u - \nu}\over{1 - \nu}} \sqrt{{{a}\over{r}}} \hbox{erfc} \left( {{r - a}\over{\sqrt{4 c t}}} \right) -3 {{a^4}\over{r^4}} \right] \cos 2 \theta \cr \cr \cr \sigma_{r \theta} &=\ S_0 \left[ -1 - 2 {{a^2}\over{r^2}} +3 {{a^4}\over{r^4}} \right] \sin 2 \theta \cr \ \ \ & &\ \ \cr u_r &=\ {{a S_0}\over{2G}} \left[ 4(1-\nu_u) {{a}\over{r}} - {{a^3}\over{r^3}} \right] \cos 2 \theta \cr \cr \cr u_{\theta} &=\ {{a S_0}\over{2G}} \left[ -2(1-2 \nu_u) {{a}\over{r}} - {{a^3}\over{r^3}} \right] \sin 2 \theta \cr \cr \cr p &=\ S_0 {{4}\over{3}} B(1+\nu_u) \left[ -\sqrt{{{a}\over{r}}} \hbox{erfc} \left( {{r - a}\over{\sqrt{4 c t}}} \right) + {{a^2}\over{r^2}} \right] \cos 2 \theta \cr }\]

In the mode 3 equations, the Skempton coefficient, \(B\), is expressed as

\[B = (K_u - K)/ (\alpha K_u)\]

The superposition of the three loading mode solutions is performed in three FISH functions “stor_pp”, “stor_sigt” and “stor_ur” to provide solutions for pore pressure, tangential stress and radial displacement, respectively, for comparison to the FLAC2D results. Note that the FISH function “ERFC.FIS” is used to calculate the complementary error function.

FLAC2D Model

The FLAC2D grid used in the simulation is shown in Figure 2. The model takes advantage of the problem quarter symmetry: it corresponds to a quarter of a ring with inner radius a equal to 1 m, and outer radius b equal to 50 m. The origin of the system-of-reference axes is at the borehole center; the \(x\)- and \(y\)-axes are in the direction of the minimum and maximum compressive in-situ principal stresses, respectively. The boundary conditions correspond to roller boundaries along the symmetry lines, zero pore pressure and no traction at the borehole, and fixed displacements at the far boundary. The grid has 50 zones in the radial direction and 16 zones on the circumference, The radial extent of the model is selected such that the stress state at the outer boundary, estimated using the Kirsch solution, compares with good accuracy to the in-situ stress field.

../../../../../_images/grid1.png

Figure 2: FLAC2D grid.

To qualify as a short-time simulation, the total simulation time must correspond to a value of the normalized time \(t^*=ct/a^2\) smaller than \(10^{-2}\) (see Detournay and Cheng 1988). The selected target time, \(t^*\) , of approximately \(10^{-3}\) is well within the range of applicability of the short-time asymptotic solution. With a diffusivity, \(c\), for this problem of approximately \(2 × 10^{-3}\) m2/s, and a borehole radius of 1 m, this time translates to roughly 0.5 second. The magnitude of the radius of influence of the borehole at that time is estimated, using \(R =\sqrt{4ct}\), at roughly 0.1 m. The grid is graded in the radial direction with a ratio of 1.18, so that 15 zones cover the radius of influence. (Note that the total number of zones in the radial and tangential directions must be large enough for the zone aspect ratio to remain below approximately 5:1 in order to avoid numerical inaccuracies.) The FLA2D model is run to three times: \(t\) = 0.003 s, 0.03 s and 0.3 s. These times correspond to normalized times \(t^*=0.7×10^{-5}\), \(0.7×10^{-4}\) and \(0.7×10^{-3}\) respectively. The analytical and FLAC2D results for pore pressure, tangential stress and radial displacement are stored in tables.

The simultion is run in two different ways. First, fully coupled, where a single fluid timestep is taken and mechanical solves to equilibrium for each fluid step. This is shown in the data file borehole.dat. The second approach, shown in borehole_uncoupled.dat in an uncoupled approach where fluid solves to a specified time, followed by solving to mechanical equilibrium. In this case, fluid modulus needs to be adjusted to account for the lack of storativity in the fluid-only solution (see zone fluid modulus-scale). Also, the generation of pore pressure during the mechanical solution is turned off with zone fluid property pore-pressure-generation off.

Results and Discussion

The results and discussion of the FLAC2D simulation closely follow those in Detournay and Cheng (1988). Results of the coupled and uncoupled simulations are very similar. The main difference between the two solutions is that the uncoupled approach solves more quickly. Figures below show results from the coupled simulation.

Isochrones of pore pressure variation with radius are compared to the analytical solution in Figure 3 for \(t\) = 0.003 s, 0.03 s, 0.3 s, and the direction \(θ\) = 0. The steep radial gradient of pore pressure developing at the borehole wall is illustrated in this figure. This gradient is associated with a rapid drainage of fluid near the borehole, which influences the apparent stiffness properties of the rock.

../../../../../_images/pore-presssure-comparison.png

Figure 3: Pore pressure variation with radius at \(θ\) =0° (\(t\) =0.003s, 0.03s and 0.3s).

The drained modulus characterizes the rock in the vicinity of the borehole, and the stiffer undrained modulus characterizes the rock farther away. The shielding effect on the stress concentration near the borehole caused by this stiffness contrast is shown in Figure 4. In this figure, the isochrones of tangential stress variation with radius are compared to the analytical solution for the three values of time considered earlier and along the \(x\)-direction. As this figure shows, at the very small times considered in the simulation, the peak tangential stress is actually located inside the rock. This mechanism explains field observations that indicate that failure can be initiated at a small distance inside the rock, rather than at the borehole wall as predicted by an elastic analysis.

../../../../../_images/tangential-stress-comparison.png

Figure 4: Tangential stress variation with radius at \(θ\) =0° (\(t\) =0.003s, 0.03s and 0.3s).

The radial displacement versus θ is compared to the analytical solution in Figure 5 for a simulation time of 0.3 second. The results of the FLAC2D simulation show good agreement with the analytical predictions.

../../../../../_images/radial-displacement.png

Figure 5: Radial displacement variation with \(θ\) at \(r=a\) (\(t\) =0.3s).

Reference

Detournay, E., and A. H.-D. Cheng. “Poroelastic Response of a Borehole in a Non-Hydrostatic Stress Field,” Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 25(3), 171-182 (1988).

Data Files

borehole.dat

;---------------------------------------------------------------------
;;; Poroelastic Response of a Borehole
;---------------------------------------------------------------------
model new
model large-strain off
fish automatic-create off
model deterministic on
model title "Poroelastic Response of a Borehole"
model configure fluid-flow

; from FLAC 8.1
;zone import 'radial'

zone create2d annular-sector point 0 0,0 point 1 50,0 point 2 0,50 ...
  point 3 1,0 point 4 0,1 size 50 16 ratio 1.18 1

program call 'Erfc.fis' suppress
program call 'fishfunctions.dat' suppress
zone face skin
zone cmodel assign elastic

; --- properties ---
zone property bulk [c_bu] shear [c_sh] density 2500

;fluid properties
zone fluid property porosity [c_n] mobility-coefficient [c_k] 
zone fluid property biot-coefficient [c_bc] biot-modulus [c_bm]
zone fluid property fluid-density 1000 
zone fluid property pore-pressure-generation on
zone fluid unsaturated cutoff 

; boundary conditions
zone gridpoint fix velocity-x range group 'West2'
zone gridpoint fix velocity-x range group 'East'
zone gridpoint fix velocity-y range group 'East'
zone gridpoint fix velocity-y range group 'Bottom

; initial conditions
zone initialize stress xx [sxx0] yy [syy0] zz [sxx0]
zone gridpoint initialize pore-pressure [pp0]

; --- undrained response ---
model solve-static
model save 'bh_a'

; --- drained response ---
zone face apply pore-pressure 0 range group 'West1'
model fluid timestep fix 3e-4

model solve-fluid-coupled fluid time-total 3e-3
[stor_pp('pp t=0.003s')]
[stor_sigt('sigt t=0.003s')]
model save 'bh1'

model solve-fluid-coupled fluid time-total 3e-2
[stor_pp('pp t=0.03s')]
[stor_sigt('sigt t=0.03s')]
model save 'bh2'

model solve-fluid-coupled fluid time-total 3e-1
[stor_pp('pp t=0.3s')]
[stor_sigt('sigt t=0.3s')]
[stor_ur('ur t=0.3s')]
model save 'bh3'

borehole_uncoupled.dat

;---------------------------------------------------------------------
;;; Poroelastic Response of a Borehole
;---------------------------------------------------------------------
model new
model large-strain off
fish automatic-create off
model deterministic on
model title "Poroelastic Response of a Borehole"
model configure fluid-flow

; from FLAC 8.1
;zone import 'radial'

zone create2d annular-sector point 0 0,0 point 1 50,0 point 2 0,50 ...
  point 3 1,0 point 4 0,1 size 50 16 ratio 1.18 1

program call 'Erfc.fis' suppress
program call 'fishfunctions.dat' suppress
zone face skin
zone cmodel assign elastic

; --- properties ---
zone property bulk [c_bu] shear [c_sh] density 2500

;fluid properties
zone fluid property porosity [c_n] mobility-coefficient [c_k] 
zone fluid property biot-coefficient [c_bc] biot-modulus [c_bm]
zone fluid property fluid-density 1000 
zone fluid property pore-pressure-generation on
zone fluid unsaturated cutoff 

; boundary conditions
zone gridpoint fix velocity-x range group 'West2'
zone gridpoint fix velocity-x range group 'East'
zone gridpoint fix velocity-y range group 'East'
zone gridpoint fix velocity-y range group 'Bottom

; initial conditions
zone initialize stress xx [sxx0] yy [syy0] zz [sxx0]
zone gridpoint initialize pore-pressure [pp0]

; --- undrained response ---
model solve-static
model save 'bh_a'

; --- drained response ---
zone face apply pore-pressure 0 range group 'West1'
model fluid timestep fix 3e-4

model solve-fluid-coupled fluid time-total 3e-3
[stor_pp('pp t=0.003s')]
[stor_sigt('sigt t=0.003s')]
model save 'bh1'

model solve-fluid-coupled fluid time-total 3e-2
[stor_pp('pp t=0.03s')]
[stor_sigt('sigt t=0.03s')]
model save 'bh2'

model solve-fluid-coupled fluid time-total 3e-1
[stor_pp('pp t=0.3s')]
[stor_sigt('sigt t=0.3s')]
[stor_ur('ur t=0.3s')]
model save 'bh3'

fishfunctions.dat

fish define ini_bh(rmin)
   global P0     = 0.;3.e6
   global S0     = 1.e6
   global pp0    = 1.e6
   global sxx0   = -P0 + S0
   global syy0   = -P0 - S0
   global c_bu   = 1.5e9
   global c_sh   = 1.5e9
   global c_k    = 1e-12
   global c_n    = 0.3
;
   global a      = rmin
;  Biot coef:
   global c_bc   = 0.65
;  Biot modulus: M
   global c_bm   = 2e9
; undrained bulk
   global K_u   = c_bu + c_bc*c_bc*c_bm
; drained poisson ratio
   global aux    = c_bu/c_sh
   global nu     = (3.*aux-2.)/(6.*aux+2.)
; undrained poisson ratio
   global aux    = K_u/c_sh
   global nu_u   = (3.*aux-2.)/(6.*aux+2.)
; Skempton coefficient
   global B=(K_u-c_bu)/(c_bc*K_u)
;  diffusion coef
   global stor   = 1./c_bm + c_bc*c_bc/(c_bu+4.*c_sh/3.)
   global diff   = c_k / stor
;  eta
   global eta    = c_bc*(1.-2.*nu)/(2.*(1.-nu))
   global c_bma  = 1./stor
;
   global b0_ii  = 0
   global ip     = 1
   global jp     = 1
   ;global tabin  = 1
   ;global tabout = 1
end
[ini_bh(1)]

fish define stor_pp(tab_in)
   local tab_out = tab_in + ' Analytical'
   local tab1  = tab_in + ' FLAC2D'
   local c0    = 1./(a*math.sqrt(math.pi))
   local c3    = S0*4.*B*(1.+nu_u)/3.
   global n_num1 = 0
;[pnt = list(gp.list)(gp.isgroup(::gp.list,"Bottom"))]
   loop foreach local pnt list(gp.list)(gp.isgroup(::gp.list,"Bottom"))
      n_num1 = n_num1 + 1
      local xval  = gp.pos.x(pnt)
      local yval  = gp.pos.y(pnt)
      local rad   = math.sqrt(xval*xval+yval*yval)
      local theta = math.atan2(yval,xval)
      local c1    = a/rad
      local c2    = math.sqrt(c1)
      local c_t   = zone.fluid.time.total
      local c4    = math.sqrt(4.*diff*c_t)
      local e_val = (rad - a)/c4
      local verfc = erfc(e_val)
      local cmod1 = 0.
      local aux   = c4*c0*math.exp(-e_val*e_val) - (rad/a-1.)*verfc
      local cmod2 = ((1.-c2*verfc)-0.125*(c2*(1.-c1))*aux)*pp0
      local cmod3 = c3*(-c2*verfc+c1*c1)*math.cos(2.*theta)
      table.x(tab_out,n_num1) = rad - a
      table.y(tab_out,n_num1) = cmod1+cmod2+cmod3
      table.x(tab1,n_num1)   = rad - a
      table.y(tab1,n_num1)   = gp.pp(pnt)
   end_loop
end

fish define stor_sigt(tab_in)
   local tab_out = tab_in + ' Analytical'
   local tab1  = tab_in + ' FLAC2D'
   
   local c0 = 1./(a*math.sqrt(math.pi))
   local c3    = P0*a/(2.*c_sh)
   local c5    = eta*a*pp0/c_sh
   local c6    = S0*a/(2.*c_sh)
   local c8    = 4.*(nu_u-nu)/(1.-nu)
   global n_num2 = 0
   local count=1
   loop foreach local zpnt zone.list
    if count <50 then
    if zone.id(zpnt) = count  then
      count = count + 2
      n_num2 = n_num2 + 1
      local xc  = zone.pos.x(zpnt)
      local yc  = zone.pos.y(zpnt)
      local rad = math.sqrt(xc*xc+yc*yc)
      local theta = math.atan2(yc,xc)
      local sn  = math.sin(theta)
      local cn  = math.cos(theta)
      local sigx = zone.stress.xx(zpnt)
      local sigy = zone.stress.yy(zpnt)
      local sigxy= zone.stress.xy(zpnt)
      local sigt = sigx*sn*sn - 2.0*sigxy*sn*cn + sigy*cn*cn
      local c1    = a/rad
      local c2    = math.sqrt(c1)
      local c7    = c8*c2
      local cval  = c1*c1
      local c_t   = zone.fluid.time.total
      local c4    = math.sqrt(4.*diff*c_t)
      local e_val = (rad - a)/c4
      local verfc = erfc(e_val)
;     tangential
      local aux   = c4*c0*math.exp(-e_val*e_val) - (rad/a-1.)*verfc
      local cmod1  = (-cval-1.)*P0
      local aux1  = c4*c0*math.exp(-e_val*e_val) + verfc
      local cmod2 = 2.*eta*(c1*c2*aux1-cval*c4*c0+0.125*(c2*(1.-c1))*aux)*pp0
      local cmod3 = (-1.+c7*verfc-3.*cval*cval)*math.cos(2.*theta)*S0
      table.x(tab_out,n_num2) = rad - a
      table.y(tab_out,n_num2) = cmod1+cmod2+cmod3
      table.x(tab1,n_num2)   = rad - a
      table.y(tab1,n_num2)   = sigt
    endif
    endif
    endloop
end

fish define stor_ur(tab_in)
   local tab_out = tab_in + ' Analytical'
   local tab1  = tab_in + ' FLAC2D'
; --- radial displacement ---
   global n_num3 = 0
   local c0    = 1./(a*math.sqrt(math.pi))
   local c3    = P0*a/(2.*c_sh)
   local c5    = eta*a*pp0/c_sh
   local c6    = S0*a/(2.*c_sh)
   loop foreach local pnt list(gp.list)(gp.isgroup(::gp.list,"West1"))
      n_num3 = n_num3 + 1
      local xval  = gp.pos.x(pnt)
      local yval  = gp.pos.y(pnt)
      local rad   = math.sqrt(xval*xval+yval*yval)
      local theta = math.atan2(yval,xval)
      local thetad= theta / math.degrad
      local c1    = a/rad
      local c2    = math.sqrt(c1)
      local c_t   = zone.fluid.time.total
      local c4    = math.sqrt(4.*diff*c_t)
      local e_val = (rad - a)/c4
      local verfc = erfc(e_val)
      local cmod1 = -c1*c3
      local aux   = c4*c0*math.exp(-e_val*e_val) - (rad/a-1.)*verfc
      local cmod2 = (c2*aux+c1*c4*c0)*c5
      local cmod3 = c1*(4.*(1.-nu_u)-c1*c1)*math.cos(2.*theta)*c6
      table.x(tab_out,n_num3) = thetad
      table.y(tab_out,n_num3) = cmod1+cmod2+cmod3
      table.x(tab1,n_num3)    = thetad
      table.y(tab1,n_num3)    = ...
         gp.disp.x(pnt)*math.cos(theta)+gp.disp.y(pnt)*math.sin(theta)
   end_loop
end