Using a fully coupled flow and geomechanical simulator to model injection into heavy oil reservoirs, Sci ...

[ Pobierz całość w formacie PDF ]
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN FLUIDS
Int. J. Numer. Meth. Fluids
(2012)
Published online in Wiley Online Library (wileyonlinelibrary.com/journal/nmf). DOI: 10.1002/fld.3679
Using a fully coupled flow and geomechanical simulator to model
injection into heavy oil reservoirs
Hao Huang
1,
*
,†
, R. Chick Wattenbarger
2
, Xiuli Gai
1
, William P. Brown
3
,
Owen J. Hehmeyer
1
, Jianlin Wang
1
and Ted A. Long
1
1
ExxonMobil Upstream Research Company, 3319 Mercer Street, Houston, TX 77027, USA
2
ExxonMobil Production Company, 800 Bell Street, Houston, TX 77002, USA
3
ExxonMobil Technical Computing Company, 13401 North Freeway, Houston, TX 77060, USA
SUMMARY
In this paper, the geomechanical factors that may affect injection processes in heavy oil recovery are
investigated. To accurately capture the geomechanical effects, we employed a numerical formulation that
allows fully coupling of nonlinear geomechanical deformation and multicomponent porous media flows.
Two salient features of this new coupling formulation are the following: (1) all flow and geomechanical
equations are solved implicitly in one single matrix equation, and (2) it allows reuse of matrices from both
a traditional fully implicit multicomponent reservoir simulator and a nonlinear geomechanics simulator.
The former feature ensures stable coupling between the reservoir flow and geomechanics, and the latter
significantly reduces the programming work. Numerical examples are given to demonstrate the accuracy
and convergence performance of the new formulation.
The proposed formulation is then applied to model injection into heavy oil reservoirs. The numerical
investigation revealed that geomechanical factors, such as
in situ
stress anisotropy and the uneven deforma-
tion of reservoir rock and attached impermeable rock, can result in skewed or nonuniform plastic strain and,
hence, alter the sweep of the injected fluid. Coupled geomechanics simulation also gives rather different
transient pressure response from that of uncoupled simulation. Copyright © 2012 John Wiley & Sons, Ltd.
Received 20 December 2011; Revised 29 February 2012; Accepted 20 March 2012
KEY WORDS
:
coupled geomechanics; heavy oil recovery; reservoir simulation
1. INTRODUCTION
Injection is involved in many heavy oil recovery processes, such as polymer or water flooding, cyclic
steam stimulation, steam-assisted gravity drainage, and solids-stabilized emulsions [1]. During
injection, pore pressure usually dissipates relatively slow because of low mobility of the heavy oil.
To admit more fluid, rock has to increase the pore volume by deforming itself. The rock deforma-
tion is determined by its stresses and mechanical constitutive relation and can only be accurately
computed by solving the force balance equation over the domain of interest. Additionally, the
process of solving the force balance equation will inevitably involve the consideration of other
geomechanical factors such as
in situ
stress anisotropy, thermal stress, and mechanical boundary
conditions. Consequently, these factors may also play a role in influencing reservoir flow.
Many reservoir simulators capture rock deformation by using variable compressibility models
and solve the porous media flow equation without coupling with the force balance equation. Such
a model has been used in cyclic steam simulation and achieved success in history matching [2].
However, the uncoupled simulation cannot easily account for geomechanical factors such as shear
*Correspondence to: Hao Huang, ExxonMobil Upstream Research Company, 3319 Mercer Street, Houston, TX
77027, USA.
†
E-mail: hao.huang@exxonmobil.com
Copyright © 2012 John Wiley & Sons, Ltd.
H. HUANG
ET AL.
stress,
in situ
stress anisotropy, and stress arching effects of overburden [3]. This limitation prevents
us from performing simulation studies to gain more understanding of the role of geomechanics in
injection processes.
The iterative coupling scheme has been proposed [4–6] to couple reservoir flow and
geomechanics. However, this scheme suffers several drawbacks in the simulation of heavy oil
reservoirs. First, the convergence for iterative coupling scheme may be challenging if both reservoir
flow and geomechanical problems are nonlinear. In fact, the scheme is tantamount to solving the
coupled equations without using consistent tangent matrix. It has been shown by Simo and Taylor
[7] that an inconsistent tangent matrix may lead to a many-fold increase in the number of iterations
and failure to converge in moderately complex problems. Second, as discussed in [8], iteratively,
coupling may lead to ‘negative compressibility’ in shear dilation. To circumvent the iterative con-
vergence issue, some simulators allow the user to set the maximum number of iterations and simply
let the simulation advance once this number is reached, regardless of convergence status. The prac-
tice may significantly affect the accuracy of the solution by leaving either flow or geomechanical
equation unbalanced when exiting an iteration.
In this paper, we present a fully coupled formulation that allows force balance and porous media
flow equations to be solved in a single matrix. The resulting matrix is the consistent tangent
matrix of the coupled system and enables fast numerical convergence. The convergence perfor-
mance of the proposed formulation is demonstrated through a numerical example. The proposed
method reuses matrices generated by both a reservoir simulator using control volume method
(CVM) and those generated by a nonlinear geomechanics solver using the finite element method
(FEM). The meshes employed for the two discretization methods do not have to coincide. After
the numerical formulation description, simulation results are shown to demonstrate the effects of
geomechanics on injection precesses of heavy oil reservoirs. Specifically, the numerical investiga-
tion revealed that geomechanical factors, such as
in situ
stress anisotropy and uneven deformation
of reservoir rock and attached impermeable rock, can alter the sweep of the injected fluid. Coupled
geomechanics simulation also gives a rather different transient pressure response from that of
uncoupled simulation.
2. NUMERICAL FORMULATION
2.1. Governing equations
The coupled pore fluid and solid system under study is governed by the mass conservation equation
of the pore fluid and the force balance equation of the rock–fluid system. The compositional
description of the mass conservation equation for the component i (ignoring source and sink) is
given by
@. N
i
/
@t
Cr
F
i
D 0,
(1)
where is the rock porosity in deformed configuration. N
i
and
F
i
are the mass density and flux of
the component i , respectively. They may be expressed as
N
p
N
p
X
X
k
j
j
K
r.p
j
gh/,
N
i
D
C
ij
S
j
j
, d
F
i
D
j
C
ij
(2)
j D1
j D1
where N
p
is the number of phases in the fluid and C
ij
is the concentration of the component i in
phase j .
j
, S
j
, k
j
,and
j
are the density, saturation, relative permeability, and viscosity of phase
j , respectively.
K
, p,andh are the absolute permeability tensor, pore pressure, and depth, respec-
tively. To simplify the presentation, capillary pressure is not considered, and p is assumed to be the
same for all phases.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Fluids
(2012)
DOI: 10.1002/fld
 COUPLED FLOW AND GEOMECHANICS SIMULATION
The strong form of the force balance equation is
r CO
g
D 0,
(3)
where is the total Cauchy stress and may be defined as
D
0
˛p
I
,
(4)
where ˛ is the Biot’s constant,
0
is the effective stress felt by the rock skeleton, and
I
is the identity
matrix. The O
in Equation (3) is the density of the coupled rock–fluid system and is defined as
N
i
X
O
D
N
i
C .1 /
rock
,
(5)
i D1
where N
i
is the total number of components of the fluid system and
rock
is the density of rock grains.
The interdependence of the reservoir flow and geomechanics is summarized in the following.
The discussion of the assumptions and derivations of the listed relations may be found in the
literature [8].
1. The effect of pore pressure on the force balance is captured by Equation (4).
2. The rock deformation (as a result of solving the force balance equation—Equation (3)) will
affect the rock porosity and permeability. The derivatives of and
K
with respect to rock strain
are
D
I
3K
g
I
W
D
@
@
1
(6)
1
3
v
I
,
@
K
@
@
K
@
v
I
C
1
dev
@
K
@
dev
D
(7)
where
D
and K
g
are the rock’s constitutive matrix and the grain compressive modulus, respec-
tively. Note that
K
is assumed to depend on the volumetric strain
v
and the deviatoric strain
dev
, which are defined as
v
D tr./
(8)
s
1
3
v
I
W
1
3
v
I
T
dev
D
.
(9)
3. The pore pressure may also cause grain deformation and, hence, a change of porosity. The
derivative of porosity with respect to pressure is
@p
D
"
1
#
@
I
W
D
W
I
9K
g
.
(10)
K
g
2.2. Discretized formulation
Denote
f
and
m
to be the spatial domain occupied by the pore fluid and the rock, respectively. In
a typical coupled simulation,
f
m
. Traditionally, the CVM is used to solve the mass balance
equation (Equation (1)) and the FEM is used to solve the force balance equation (Equation (3)).
These two discretization schemes will be adopted here and are illustrated in Figure 1. In the CVM,
f
is discretized into cells and unknowns are sought at cell centers. The domain occupied by the
union of the cells is denoted by
f
. Unknowns for the flow equation at a cell c are denoted by
a vector l . One choice of unknowns may be pressure, N
p
1 saturations, and N
i
N
p
concen-
trations. In FEM,
m
is discretized by elements and unknowns are usually located at nodes. To
facilitate computations involved in incremental constitutive equations, we choose the displacement
increments over the time step,
u
, as unknowns.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Fluids
(2012)
DOI: 10.1002/fld
 H. HUANG
ET AL.
m
f
(a)
h
m
f
(b)
(c)
Figure 1. (a) The flow and mechanical domain; (b) a discretization of the flow domain by using control
volume method; (c) a discretization of mechanical domain by using finite element method.
By using a backward Euler method to integrate the time, the numerical formulation to solve
Equation (1) in a cell c by using CVM is
V
c
C
X
s
. N
i
/ j
t
nC1
. N
i
/ j
t
n
t
E
f ,i
D
F
cs
i
j
t
nC1
D 0,
(11)
where E
f ,i
denotes the mass balance of component i in cell c, V
c
is the volume of the cell c,and
F
cs
i
is the flux from cell c to s. The calculation of F
c
i
depends on the specific flux calculation
method. This type of equation needs to be solved for all cells and components.
The weak formulation of Equation (3) over
m
is
Z
Wr d C
Z
g
d D
0
(12)
m
m
for all admissible test functions ,thatis, W
m
! R
3
with zero trace on the displacement
boundaries. Following standard finite element techniques, Equation (12) results in
Z
X
N
e
B
T
d C
f
ext
D
0
,
E
m
D
(13)
m,e
eD1
where N
e
is the number of elements,
m,e
is the domain occupied by element e,
B
is the strain-
displacement matrix [9], and
f
ext
is the external nodal forces arising from the loading terms in
Equation (3). E
m
computes the force balance of all the nodes in the mesh.
In summary, the discretized equations that need to be solved to get a numerical approximation to
the balance equations (Equations (1) and (3)) are
°
E
f ,1
, E
f ,2
, :::, E
f ,N
i
, E
f ,1
, :::, E
N
c
f ,N
i
, E
m
±
T
D
0
,
(14)
where N
c
is the number of cells in the CVM.
2.3. Computation of Jacobian matrix
The Newton–Raphson method is employed to solve the discretized equations (Equation (14)). At
each time step, one needs to iteratively compute corrections to unknowns. This process continues
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Fluids
(2012)
DOI: 10.1002/fld
COUPLED FLOW AND GEOMECHANICS SIMULATION
until a solution converges. Inside each iteration, Equation (14) is first linearized, which will lead to
a matrix equation. Then, the matrix equation is solved to obtain corrections to the unknowns. The
matrix equation has the form
2
4
3
5
<
=
<
=
:::
:::
:::
:::
:::
C
l
c
:::
C
u
:::
R
f ,i
:::
R
m
@E
f,i
@l
c
@E
f,i
@
u
:::
:::
D
,
(15)
:::
:::
:::
:::
:
;
:
;
@E
m
@l
c
@E
m
@
u
:::
:::
where
C
l
c
is the correction for the flow unknowns at cell c and
C
u
is the correction to displacement
increment. R
f
and R
m
are the residuals of the flow and geomechanical equations, respectively.
The computation of each term in the matrix of Equation (15) is detailed as follows.
1. The pure flow contribution, @E
f ,i
=@l
c
, captures the effect of flow unknowns l on the mass
balance. It can be directly requested from the existing reservoir simulator.
2. The pure geomechanical contribution, @E
m
=@
u
, captures the effect of displacement unknowns
on the force balance. It can be directly requested from the existing geomechanical simulator.
3. The mechanical-to-flow coupling term, @E
f ,i
=@
u
, captures the effect of deformation on the
fluid flow and may be expressed as
C
X
s
@E
f ,i
@
u
DN
i
@.V
c
/
@
@.
x
c
/
@
u
@F
cs
i
@
K
@.
x
m
/
@
u
@
K
@
,
(16)
where
x
c
is the location of the center of the cell c and
x
m
is the location where
K
is evalu-
ated. For example, if an upstream weighting scheme is used,
x
m
is located at the center of the
upstream cell. N
i
and @F
c
i
=@
K
may be exported from the reservoir simulator. @. V
c
/=@ and
@
K
=@ may be evaluated using Equations (6) and (7), respectively. The derivative of with
respect to the nodal displacement at a node k can be evaluated using
"
#
@.
x
/
@
u
k
@
N
k
./
@.
x
t
C
u
=2/
D sym
,
(17)
where symΠis the symmetric operator,
N
k
is the interpolation functions for node k,
x
t
is the
location of
x
at time t,and is the natural coordinates of
x
.
4. The flow-to-mechanical coupling term, @E
m
=@l
c
, accounts for the effect of pore pressure on
the force balance equation. This term involves integration over elements, and Gaussian quadra-
ture is used to evaluate the integration. If meshes for CVM and FEM do not coincide, the pore
pressure at each quadrature point needs to be interpolated from values at neighboring cell
centers. The moving least square method [10], which was designed for interpolating values
between two particle systems in large deformation simulation, may be employed as an inter-
polant to compute pore pressure at quadrature points. The relation for the interpolation may
be expressed as
p.
x
q
/ D
X
r2N
r
T
qr
p
r
,
(18)
where N
r
is the set of cell centers in the neighborhood of
x
q
, p
r
is the pore pressure at the cell
center r,andT
qr
is the weight associated with p
r
. T
qr
needs to satisfy at least
P
r2N
r
T
qr
D 1
to reproduce a constant function. Specific formulations to reproduce higher order functions
may be found in [10].
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Fluids
(2012)
DOI: 10.1002/fld
  [ Pobierz całość w formacie PDF ]

  • zanotowane.pl
  • doc.pisz.pl
  • pdf.pisz.pl
  • emaginacja.xlx.pl

    •