Note: Descriptions are shown in the official language in which they were submitted.
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
OIL OR GAS PRODUCTION USING COMPUTER SIMULATION OF OIL OR GAS
FIELDS AND PRODUCTION FACILITIES
Qinghua Wang
Graham Fleming
Qin 1,u
BACKGROUND OF THE INVENTION
[0001] The present invention relates to oil and gas production using
computerized
simulation of oil or gas reservoirs and of production facilities.
[0002] Fig. 1 schematically illustrates production facilities for extracting
oil and/or gas
from a reservoir 110. Reservoir 110 consists of subsurface rock formations
below the
earth surface 112. A wellbore 120.1 has been drilled into the reservoir and
connected to
surface facilities 130 a system of pipes and pumps 130 to pump the fluid (oil
or gas,
possibly mixed with water or mud) out of the reservoir. A casing has been
inserted into
wellbore 120.1, and perforations 134 were made in the casing to allow the
wellbore to
communicate with the reservoir. Another wellbore 120.2, connected to surface
network
130 and having perforations 134, can optionally be used to inject water,
steam, carbon
dioxide or other chemicals into the reservoir to push the oil or gas towards
wellbore 120.1
or make sticky oil more mobile.
[0003] Design of production facilities involves choice of wellbore locations
and injected
chemicals, well shapes and dimensions, pipe and pump networks, drilling
schedules,
production schedules, etc.. The optimal choices of these and other parameters
depend on
physical properties of reservoir 110 such as pressures, saturations,
porosities, and
possibly others. Such choices are facilitated by computer simulation that
simulates
production for existing or hypothetical production facilities and production
schedules. A
suitable computer system 204 (Fig. 2) includes one or more computer processors
210
coupled to a memory 220. Typically, test measurements are made at selected
locations in
reservoir 110 or in existing production facilities to determine parameters of
interest. The
test data may include geological data from sources such as seismic analysis,
rock cores,
and well logs. These data may include data on rock porosities, permeabilities,
and other
information of interest, and allow determination of parameters of interest
such as
- 1 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
pressures or saturations at a number of locations. The test measurements are
stored in
memory 220 as shown at 230. A simulator program 240, stored in memory 220, is
executed by processors 210 to read the test data 230 and provide simulated
results which
can be stored in memory 220 and/or displayed on a printer or computer monitor
or other
output device, or provided to another computer program or another computer
over a
network or otherwise, for use by oil field designers or other personnel.
[0004] A discrete model of reservoir 110 and the production facilities is
created by
computer 204 in memory 220 for use by simulation program 240 as illustrated in
Fig. 3.
More particularly, an imaginary grid of cells (grid blocks) 320 is
superimposed on
reservoir 110, with data indicating the locations of wellbores 120 and
perforations 134. A
node (not shown) is defined in each cell (typically, but not necessarily, at
the center), and
the pressure and possibly other values in the cell will be modeled as the
pressure or other
values at the node. The wellbore portions above the perforations are shown as
120.1W,
120.2W and will be called just "wells" herein. Nodes 324 are defined along
each well
120W (i.e. 120.1W or 120.2W) in the computer model; these nodes represent
point
locations at which the pressures or other values will modeled. A well may have
one or
more nodes 324. Connections 330 represent the well portions between the nodes
324.
Similar nodes and connections (not shown) can represent the surface network
130. Other
nodes 325 may be defined to represent point locations in the perforated
sections of the
wellbores (below the wells). Connections 331 represent the wellbore portions
between
the nodes 325 and between the top wellbore node 325 and the first well node
324 (the
"bottom-hole" node). A flow of material between a perforated cell 320 (i.e. a
cell
perforated by a wellbore) and the wellbore will be modeled as the flow between
the
perforated cell and a node 325, or a flow between the cell and the
corresponding bottom-
hole node 324.
[0005] A coordinate system (XYZ) is chosen to represent various locations. A
time axis
304 with time points ti, t2, t3 ... is defined in computer 204.
[0006] A set of equations is chosen to model the oil or gas field, including
hydraulics
equations between nodes in the wells, wellbores, and the surface network;
Darcy and
non-Darcy flow between cells 320; perforation equations between cells 320 and
nodes
325; fluid component mass balance equations for nodes and cells; constraint
equations,
- 2 -
CA 02873406 2016-03-09
and possibly others. See U.S. patent no. 7,668,707 issued February 23, 2010 to
Watts,
III et at. The equations may include differential equations, and the
differential
equations are discretized using known techniques. The equations can be
symbolically
written as:
[0007] F(Vari k 1) = 0 (1)
,
[0008] where Varkl represents pressures, flow rates, or other variables at the
nodes
and connections at different times. The subscript "i" represents the time
point ti on
axis 304, and the subscripts j, k, I represent X, Y, and Z values
respectively. The
values Vari k I are treated as independent variables.
9.) 9
[0009] Since some physical parameters can be expressed through others, there
is
some latitude as to which parameters are chosen as the independent variables
(primary
variables) in equations (1). Further, the equations (1) depend on the
discretization
method, and in particular on whether any particular Var value is measured at a
cell
320 or a node 324 or 325, or a segment 330 or 331. Such choices will affect
how
efficiently the computer can solve the equations. The term "efficiently"
relates to
minimal use of the computer resources, and particularly the processor and
memory. It
is desirable to reduce the execution time and memory amount needed to solve
the
equations while providing a solution that serves as a good approximation to
the real
world.
[0010] There can be millions of different locations and time points, so the
equation
(1) may represent millions or billions of scalar equations. In a typical
procedure, the
equations (1) are solved for each of consecutive times ti, i.e. for time t1,
then for time
t2, and so on. The solution for time t1 depends on the test measurements
stored as
shown at 230 in Fig. 2. The solution for each subsequent time ti is found
using the
previously determined Var values for time ti_i and possibly earlier times and
possibly
using the test measurements 230.
[0011] Equations (1) are typically non-linear, and for each time ti these
equations can
be solved approximately using the Newton-Raphson method illustrated in Fig. 4.
This
is an iterative method in which the non-linear equations (1) are replaced by
linear
equations at each iteration. The successive iterations use different linear
equations to
provide
- 3 -
CA 02873406 2019-11-12
WO 2013/180709 PCT/US2012/040023
hopefully better and better approximate solutions at each iteration. To
simplify the
notation, let us re-write the equations (1) as:
[0012] F(x) = 0
[0013] Suppose xn is a current iterative solution in the Newton-Raphson
method. Then
the next iteration xn+1 is obtained as follows. The Jacobian matrix JF(xn) is
calculated;
this is the matrix of partial derivatives of F(x) evaluated at the current
iteration xn . Then
the following linear equation is solved for
[0014] JF(xn)( xn+i xn) = -F(xn)
[0015] This equation can be alternatively written as:
[0016] JF(xn)( xn+i) = JF(xn)( xn) - F(x)
(2)
[0017] or
[0018] F'L( xn+ 1) = b (3)
[0019] where F'L= JF , and b is the constant equal to the right-hand side of
equation (3).
[0020] This process is illustrated in Fig. 4. The symbol Var.NT denotes the
current
approximation xn . Since ti is fixed, the Vari j,k,i will sometimes be written
as simply
Vat) Jo. The values Var.NT of the last iteration are taken as the approximate
solution
Var,j,k,l of the equations (1) for the time t.
[0021] The initial iterated values Var.NT is chosen at step 410, and this can
be any
desired values, e.g. the values Vari_i ,k,1 found in the previous iteration i-
1 if i> 1, or to
some values obtained from the initial measurements 230 if i=1.
[0022] At step 414, the matrix JF(Var.NT) and the value
[0023] b= JF(Var.NT)Var.NT-F(Var.NT)
[0024] are determined for the equations (2)-(3). At step 418, these equations
are solved,
i.e. the following equation is solved:
[0025] F'L(Vari,m) = b (4)
[0026] Let us denote the solution as Var*; If this is the last iteration,
then the solution
J1 5
Var*; ki is taken as the solution of (1), i.e. the method of Fig. 4 terminates
for See step
J5 5
422. If other iterations are desired, then the solution Var*j,k J is used as
the next value of
- 4 -
CA 02873406 2019-11-12
WO 2013/180709 PCT/US2012/040023
Var.NT for the next iteration (step 426), and steps 414-418 are repeated. Any
number of
iterations can be used to obtain better and better approximations of the
solution of (1).
The iterations can be terminated when F(Var*; ki) is close to zero in some
sense.
J1
[0027] Equations (4) may include billions of equations with billions of
individual
variables shown collectively as Vari Jo. Therefore, computer system 240 solves
equations (4) approximately. An exemplary approximate method is a Schwarz
method
illustrated in Fig. 5. For simplicity, the function F'L will be denoted simply
as f, and
Vari,k,1 as x. Thus, equation (4) can be written as:
[0028] f(x) = b (4')
[0029] The function f represents many scalar linear functions, and the
variable x
represents just as many scalar components. The Schwarz method will be
illustrated on the
example in which each of x and f has 1000 scalar components, denoted as xl,
x2, ...,
x1000 for x, and denoted as fl, f2, f1000 for f.
[0030] At step 610, the values xl,..., x1000 are divided into two overlapping
groups of
values (more than two groups can also be used). In the example discussed
below, the two
groups are:
[0031] Gl: xl, ..., x600
[0032] G2: x400,... x1000
[0033] Also, the scalar functions fl,...,f1000 are divided into two
overlapping groups,
with the same number of functions in each group as there are values in Gl, G2.
For
example, the two groups of functions can be:
[0034] Gl': fl,...,f600
[0035] G2': f400,..., f1000
[0036] We will use the notation Gl'(x) to denote the values fl (x) through
f600(x); the
notation G2'(x) will denote the values f400(x) through f1 000(x).
[0037] The Schwarz method is an iterative method, and the iterative solution
will be
denoted as x.sc. The value x.sc output by the Schwarz method is the solution
x.
[0038] At step 620, the initial iterative value x.sc is determined using any
suitable
technique. For example, x.sc can be set to the value of Var.NT at step 418 of
the current
Newton-Raphson iteration of Fig. 4.
- 5 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[0039] At steps 630-640, the next iterated value is determined for x.sc as
follows.
[0040] At step 630, the x values from G1 \G2, i.e. the values x601 through
x1000, are set
to the corresponding values of the current iterative value x.sc, and the
following equation
is solved for xl through x600:
[0041] Gl'(x)= b 1 (5)
[0042] where bl is the "b" portion corresponding to G1'. In other words, the
following
equation is solved:
[0043] Gl'(xl,...,x600, x.sc601,...,x.se1000)= bl (5')
[0044] The solution will be denoted xl*,...,x600*.
[0045] The values xl* through x399* will be used as the corresponding values
of the
next iterative value of x.sc as explained below.
[0046] At step 640, the following equation is solved for x400* through x1000*:
[0047] G2'(x) = b2 (6)
[0048] where b2 is the b portion corresponding to G2. In this equation, the
values
xl,...,x399 are constant. In the additive Schwarz method, these values are set
to the
corresponding values of the current iterative value x.sc. In the
multiplicative Schwarz
method, these values are set to the corresponding values xl*,...,x399* found
at step 630.
The solution of (6) will be denoted x600*,...,xl 000*. The values xl* through
x1000*
will be collectively denoted as x*.
[0049] If this is the last iteration, then x*can be taken as the solution of
(4'), i.e. as the
value Var*; k I at step 418. See step 650. If other iterations are desired
(step 660), then
J,
the solution x* is used as the next iterated value of x.sc, and steps 630-640
are repeated.
Any number of iterations can be used to obtain better and better
approximations of the
solution of (4').
[0050] In some variations of the Schwarz method, iterations other than the
first iteration
are solved for the residual b-f(x) rather than x. The residual is updated at
each iteration.
In particular, the overlap values x400*,...,x600* found at step 630 are used
to update the
residual for step 640.
[0051] In a variation, the components x400*,...,x600* of x* are obtained for
step 650 by
combining the corresponding values obtained at steps 630 and 640, i.e. by
combining the
corresponding components of solutions of equations (5) and (6) so as to
minimize the
- 6 -
CA 02873406 2016-03-09
residual's normilb-f(x*)I . A known method for combining these components is
GMRES (Generalized Minimum Residual method).
[0052] Clearly, the computational complexity of solving the equations (1)
depends on
the choice of variables Var and the groups GI, G2, G1', G2'. Generally, the
reservoirs
are quite different from wells with respect to the physics models, so the
reservoirs and
the wells are usually in different groups. More particularly, each group GI,
G2 can be
defined through cells 320 or nodes 324 or 325: each group consists of
variables
pertaining to a group of cells 320 or nodes 324 or 325 or the related segments
330 or
331, and each function group G1', G2' consists of the equations for the
corresponding
cells and nodes. Thus, the reservoir cells can be placed in group G1 (i.e. in
the node
group defining G1), the wells in group 02, and the perforation cells 320 (the
cells
with perforations) and vvellbore nodes 325 can be the overlap between G1 and
G2.
This scheme is used in NEXUS Reservoir Simulation Software system available
from Halliburton of Houston, Texas. In NEXUS , the variables Var are pressures
and
component compositions at nodes 324 and 325, and pressures and component
masses
in the cells 320.
[0053] Background information can also be found in the following publications
(no
representation is being made as to whether or not these publications are prior
art
relative to the present disclosure):
[0054] - U.S. patent application publication 2009/0294122 Al (Dec. 3, 2009),
inventors Hansen et al.
[0055] ¨ U.S. patent no. 7,668,707 B2 (Feb. 23, 2010), inventors Watts. III et
al.
[0056] ¨ U.S. patent no. 7,516,056 B2 (April 7, 2009), inventors Wallis et al.
[0057] - U.S. patent application publication 2010/0286971 Al (Nov. 11. 2010),
inventors Middya et al.
[0058] ¨ UK patent application GB 2 434 235 A (July 18, 2007), applicant
Schlumberger Technology Corporation.
[0059] ¨ European patent application EP 2 192 507 Al (June 2, 2010), applicant
Maersk Olie OG Gas A/S.
[0060] ¨ PCT application WO 2006/020952 A2 (February 23, 2006), applicants
Saudi
Arabian Oil Company and ARAMCO Services Company.
- 7 -
[00611 ¨ H. Kazemi, "Numerical Simulation of Water-Oil Flow in Naturally
Fractured Reservoirs", Society of Petroleum Engineers Journal, December 1976,
pages 317-326.
[0062] ¨ Rafael Bru et al., "Overlapping additive and multiplicative Schwarz
iterations for H-matrices", Linear Algebra and its Applications, Elsevier,
2003.
SUMMARY
100631 This section summarizes some features of the invention. Other features
may be
described in the subsequent sections. The invention is defined by the appended
claims, which are incorporated into this section by reference.
[0064] The inventors have observed that it may be desirable to include the
wells
120W in the overlap to obtain a better coupling between the reservoir and the
surface
network. In some embodiments however, the "perforation" cells 320 (i.e. the
cells
with perforations 134) are removed from the overlap. Smaller overlap reduces
the
computational complexity. Other embodiments include the perforation cells in
the
overlap.
[0065} According to another aspect of the invention, the equations for the
well and
wellbore nodes are processed to eliminate some non-pressure variables.
Therefore,
when these variables are treated explicitly, the error caused by the explicit
treatment is
reduced, so the pressure solution is more accurate. In many cases, the
accuracy of the
pressure solution across the reservoir and surface network determines the
performance
of the global linear solver for equations (4), i.e. determines the ability to
get an
accurate solution for all the variables at a minimum computational cost.
[0065a1 According to another aspect there is provided a computer implemented
method for of determining a plurality of values for a computer simulation, the
plurality of values are determined for a facility comprising an oil and/or gas
reservoir
and a plurality of wellbores each of which has been made or is to be made,
each
wellbore providing one or more passages for fluid flow into or out of the
reservoir, the
method comprising: acquiring, by a computer system, a plurality of test
measurements
from at least one selected location of one of the reservoir and the facility;
obtaining,
by the computer system, a data model based on the plurality of test
measurements,
the data model representing the reservoir and each wellbore and comprising one
or
- 8 -
CA 2873406 2017-12-12
more nodes at said facility, the nodes comprising one or more sets of nodes,
each set
comprising a single node or a plurality of interconnected nodes, at least one
node in
each set being in a wellbore; obtaining, by the computer system, a system of
linear
algebraic equations interrelating said values, said values being represented
as
variables in the equations, wherein said values include a pressure for each
node in
each set and include flow rate values defined by flow rates into or out of the
nodes in
each set, and for each node the equations include a set of one or more
equations with
values at the node, wherein for at least one node of each set of nodes, the
set of one or
more equations comprises at least one equation with a pressure at the node and
one or
more flow rate values at the node; performing, by the computer system, a
linear
transformation of the system of equations, the linear transformation
comprising a
linear transformation of the equations of each set of equations to eliminate
at least one
flow rate value from at least one equation in each set of equations, the
linear
transformation providing a transformed system of equations, wherein
eliminating the
at least one flow rate value comprises: in each set of nodes, selecting a
reference
node; eliminating total mass flow rates for each set of nodes; in each set of
nodes,
adding mass balance equations for non-reference nodes to corresponding mass
balance equations for reference nodes, thereby eliminating off-diagonal
compositions
for the reference nodes; for each reference node in a set of nodes, adding
mass
balance equations for components to mass balance equations for total flow,
thereby
eliminating diagonal compositions for the reference nodes; and for each set of
nodes,
combining columns for pressures at the nodes in each set; and solving, by the
computer system, the transformed system of equations for said values; and
generating
the computer simulation of an oil or gas reservoir or production facility
using the
transformed system of equations.
[0066] The invention is not limited to the features and advantages described
above
except as defined by the appended claims.
BRIEF DESCRIPTION OF THE DRAWINGS
[0067] Fig. 1 shows a vertical a cross section of an oil or gas field with
production
facilities according to prior art.
[0068] Fig. 2 is a block diagram of a computer system that can be used for
simulation
of oil or gas fields and production facilities according to prior art and some
- 8a -
CA 2873406 2017-12-12
embodiments of the present invention.
[0069] Fig. 3 is a schematic diagram of a model of an oil or gas field and the
surface
network according to prior art
- 8b -
CA 2873406 2017-12-12
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[0070] Fig. 4 is a flow chart of a Newton-Raphson method used to simulate an
oil or gas
field and production facilities according to prior art.
[0071] Fig. 5 is a flow chart of a Schwarz method used to simulate an oil or
gas field and
production facilities according to prior art.
[0072] Fig. 6 is a schematic diagram of models of oil or gas fields and the
surface
networks which models are used to illustrate computer simulation according to
some
embodiments of the present invention.
[0073] Fig. 7 is an example of a Jacobian matrix at the wellbore and surface
network side
before elimination of perforation flow rates and total connection flow rates.
[0074] Fig. 8 is an example of a Jacobian matrix at the wellbore and surface
network side
after elimination of perforation flow rates and total connection flow rates.
[0075] Fig. 9 is an example of a Jacobian matrix at the wellbore and surface
network side
after adding equations for a perforation node (such as 325) to equations at a
bottom-hole
node serving as a reference node.
[0076] Fig. 10 is an example of a Jacobian matrix at the wellbore and surface
network
side after elimination of compositions in the total mass balance equation for
the bottom-
hole node serving as a reference node.
[0077] Fig. 11 is an example of a Jacobian matrix at the wellbore and surface
network
side for a gridded well model with a rate constraint.
[0078] Fig. 12 is an example of a coefficient matrix at a reference node after
the mass
balance equations for non-reference nodes have been added to the reference
node.
[0079] Fig. 13 is an example of a coefficient matrix at a reference node after
dropping of
hydraulics derivatives with respect to compositions.
[0080] Fig. 14 is an example of a coefficient matrix at a reference node after
elimination
of composition coupling terms at the reference node.
[0081] Fig. 15 illustrates an example of a global matrix structure for
implicit formulation
with implicit wells and surface network.
[0082] Fig. 16 illustrates an example of a global matrix structure without the
surface
network after elimination of mass coupling terms in the producing perforation
connections.
[0083] Fig. 17 illustrates an example of a pressure matrix structure.
- 9 -
CA 02873406 2016-03-09
[0084] Figs. 18, 19 are flow charts of simulation methods according to some
embodiments of the present invention.
DESCRIPTION OF SOME EMBODIMENTS
[0085] The embodiments described in this section illustrate but do not limit
the
invention. The invention is defined by the appended claims.
[0086] Some embodiments of the present invention will now be described on the
example of the multiplicative Schwarz method as in Fig. 5. The invention
includes
additive Schwarz however, and can be combined with many existing simulation
techniques. Of note, the Schwarz methods are part of a class of methods that
can be
described as follows:
[0087] ¨ set some of the variables to constant values, and solve some or all
of the
equations for the other variables;
[0088] ¨ then set some or all of those other variables to constant values, and
solve for
some or all of the variables that were constant in the previous step.
[0089] In oil and gas simulation, one type of such methods is IMPES (implicit
pressure explicit saturation): to solve the equations (1) for each time ti
given the
solution at time ti_1 (where to means the initial values), the non-pressure
variables are
set to their ti_i values, and the equations are solved to find the ti values
of the
pressures; then the pressures are set to their ti values, and the equations
are solved to
find the ti values of the non-pressure variables. Some of these and other
methods are
described in the following publications:
[0090] - Yuanlin Jiang, "Techniques for Modeling Complex Reservoirs and
Advanced Wells", Dissertation submitted to the Department of Energy Resources
Engineering and the Committee on Graduate Studies of Stanford University,
December 2007.
[0091] ¨ Wing Hong Felix Kwok, "Scalable Linear and Nonlinear Algorithms for
Multiphase Flow in Porous Media", Dissertation submitted to the program in
scientific computing and computational mathematics and the committee on
graduate
studies of Stanford University, December 2007.
100921 ¨ Denis Jose Schiozer, "Simultaneous Simulation of Reservoir and
Surface
- 10 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
Facilities", Dissertation submitted to the Department of Petroleum Engineering
and the
committee on graduate studies of Stanford University, March 1994.
[0093] ¨ Eero Vainikko, "Robust Additive Schwarz Methods ¨ Parallel
Implementations
and Applications", Thesis for the Degree of Doctor Scientiarum, Department of
Informatics, university of Bergen, March 1997.
[0094] Other such techniques include CPR (Constrained Pressure Residual).
[0095] Some embodiments are compatible with AIM (Adaptive Implicit Method),
and
others described in the documents cited immediately above. Some embodiments
are also
compatible with FIM (Fully Implicit Method), i.e. when all the variables are
treated
implicitly (none are replaced by constants) except for the Schwarz method at
the stage of
Fig. 5. However, as noted above, the invention includes pressure matrix
derivations
which are not limited to use of the Schwarz method.
[0096] Now some multiplicative Schwarz embodiments will be described with just
two
groups of variables, G1 and G2, for the models illustrated in Figs. 3 and 6
(described
below). In some such embodiments, the group G1 corresponds to the reservoir
and
perforation cells 320, the wellbore nodes 325 and segments 331 and the well
nodes 324
and segments 330 of wells 120W. The group G2 corresponds to all the wells
120W, the
wellbores, perforation cells 320 and surface network 130. The overlapping
domain is the
nodes 324 and segments 330 of wells 120W, the wellbore nodes 325 and segments
331,
and the perforation cells 320.
[0097] The global linear system (4) is, in matrix form:
Aõ xrhr
[0098] ['Emir Aww Awfl[xwpFbwl (7)
Afw Alf Xf bf
[0099] Here the subscripts r, w, f denote the reservoir, well, and surface
network 130
respectively. For subscript r, the reservoir includes the cells with
perforations 134. The
vector [br, bw, bf] on the right hand side (RHS) is the same as b in equation
(4). The
vector [xr, xw, xi] is the same as Vari .
[00100] In the matrix on the left, the rows with Arr correspond to the
reservoir
portion of equations (4), including the perforation cells 320. The rows with
Awr
-11-
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
correspond to the well equations, which include equations for flow from each
perforated
grid block 320 to the corresponding wellbore node 325. The rows with Afiv
correspond to
the facilities equations. The columns with Arr correspond to the reservoir
variables
(including the perforation grid blocks 320). The columns with Arw correspond
to the well
variables. The columns with Awf correspond to the facilities variables. The
blank entries
in the matrix are zeroes.
[00101] The Cil group consists of variables xi., xw, and the GI '
equations can be
written as:
[00102] [Aõ Arw1 xr _ [ br 1
(8)
kw, ww W[bud
[00103] These equations are obtained by setting xf to constant values
as described
above. The bw values in (8) are not necessarily the same as in (7).
[00104] The G2 group consists of variables xw, xf and the 02' equations
can be
written as:
rilw, Awl [xw] ovvi
(9)
1 A11 J Lxf J [bf
[00105] The b values in (9) are not necessarily the same as in (7).
[00106] System (9) might not exist if the surface network 130 is absent
or is not
being simulated. If the surface network exists in the simulation system, the
system (9) can
be solved using known techniques, e.g. a sparse direct solver or ILU
(Incomplete
Lower/Upper) type iterative solver. We focus on the solution of system (8),
which is
more complicated in some embodiments.
[00107] For FIM, we need to construct the CPR pressure equation on the
reservoir.
For FIM and IMPES, an approximated pressure matrix for wellbore nodes needs to
be
constructed since an implicit well modeling is used in both formulations.
[00108] Fig. 6 shows exemplary well models 680, 690, 694 which will be
used in
addition to the model of Fig. 3 to illustrate some embodiments of the present
invention.
These models are non-limiting examples. All of these models are for production
wells
(for pumping oil or gas out of the reservoir), but similar techniques apply to
injector wells
(for injecting water or other substances into the reservoir). The model 680
includes a
- 12-
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
perforation cell 320 connected to a wellbore node 325 (also called "node 1")
by
perforation 134, a "bottom hole" node 2 also marked as 324BH (the first node
of well
120W downstream of the perforation), and the surface network 130. The "q"
variables
denote mass transfer rates (mass per time) at various connections 330.
[00109] The model 690 includes a multisegment well which splits into two
segments at node 324.1. In each segment, the respective bottom-hole node 324BH
is
connected to the respective wellbore node 325 connected to perforation cells
320. Above
the common node 324.1, there are two other nodes 324.2 and 324.3. Node 324.3
is
connected to surface network 130. The models of Figs. 3 and 6 are "lumped
well"
models, i.e. the flow from the perforation cells in each wellbore is modeled
as the flow
from these cells into a single node 325.
[00110] Model 694 (gridded well model, i.e. with multiple nodes 325)
contains
three perforation nodes 325, marked as nodes 1, 2, 3 from the bottom. Node 3
is
connected to bottom-hole node 4 (marked as 324BH). Each of nodes 1, 2, 3 is
connected
to a respective perforation cells 320. Node 324BH is connected to surface
network 130.
[00111] The global matrix for implicit formulation (FIM) for model 694
is shown
in Fig. 11. In the following, we will discuss how to approximate the pressure
solution at
the well side (equations (8)) with a simple lumped well model.
[00112] The pressure drop across the connection between an upstream
node and
the next downstream node (e.g. from node 1 to node 2 in model 680) can be
determined
by a hydraulics equation as follows:
[00113] R = Pup ¨ Pdw f(q1, (12, === Pdw) (e4)
[00114] Here the subscripts up and dw represent the upstream and
downstream
nodes (e.g. 1 and 2 in model 680) at a flow connection 330 or 331;
[00115] - Pp and P dw are the pressures at these nodes;
u
[00116] - plc is the number of components in the fluid flowing through the
connection;
[00117] - for each component k, the value qk is the mass flow rate of
the
component;
[00118] - fis
the function of the hydraulics equation (e.g. based on a pressure drop
correlation, tabular input of the pressure drop, Darcy's law or some other
model).
- 13 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[00119] We define the composition of component k as
qk
Xk -= ¨ (e5)
qT
[00120] where oT is the total mass flow rate through the connection
(i.e. the sum of
qk). (The x notation for the compositions should not be confused with the x
symbol in
equations (7), (8) and (9)).
[00121] Thus,
qk = XkqT (e6)
[00122] where
nc
qT =Iqk (e7)
k=1
[00123] Combining the equations (e6) and (e7) yields
nc
1dXk = 0 (e8)
k=1
[00124] Differentiating equation (e6) gives
dqk = qTdxk + xkdoT (e9)
[00125] Plugging equations (e8) and (e9) into equation (e4) yields
dR = (1 ¨ dPõ ¨ (1 + f.19 clwdP ¨
a Pup P a Pdw
(e 10)
vnc-1 (0 f fq dx ¨ Enc x d
la o
k=1 aqk atinc T k k=1 aqk k T
[00126] So the Jacobian derivatives (partial derivatives) with respect
to the
primary unknowns are
OR f
(ell)
a Pup a Pup
- 14 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
a R Of
_________________________ = 1 _______________________________ (e12)
a Pdw a Paw
OR Off
_ _ ' (e13)
Oxk \Oq, qõ
nc
OR Of
¨ = (c14)
ach
k=1
[00127] Of note, the letter R herein corresponds to F in equations (1). The
same
symbol R is used for different scalar equations (1).
[00128] Fig. 7 is a Jacobian example. The derivatives (e 1 1)-(e14) are
shown in the
fourth row, marked "(e4)" in the right column (referring to the equations
(e4). Fig. 7
assumes that the derivatives ¨af and ¨Of are zero. This is true iff is treated
explicitly
a qup aP dw
(i.e. replaced with a value from the previous iteration), but can be true in
case of
implicitly treated f. Some derivatives (e13), (e14) are shown as h and s
symbols.
[00129] A pressure constraint may be applied to the bottom-hole node:
R = P Pspec = 0 (e15)
where P spec is a predefined constant. Alternatively, a phase flow rate
constraint can be
activate at a connection, i.e.:
np
R = anQn Qspec = 0 (e16)
n=i
[00130] Here np is the number of phases (in many applications, there
are three
phases: gas, oil, and water). For each phase n, Q71 is the corresponding
volume flow rate;
an is the corresponding predefined constant. Qspec is the maximum allowed
volume flow
rate that can flow in the connection.
[00131] Alternatively, the following constraint can be specified:
np nc np
aQn d N7, D
dR =1(4,1¨ u.qk a 1,L1 = u (e17)
n=1 k=1 aqk n aP
n=1
- 15 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[00132] Plugging equations (e8) and (e9) into equation (e17) yields
?V nc-1
dR =11 an (oar, aQn)
qTdxk
a qk a qnc)
n=1 k=1
np nc np (e18)
aQn Qn dn
+ an -- Xk a ¨ =
n
a qk n ap
n=1 k=1 n=1
[00133] Usually the rate constraint is specified at surface condition, so
the pressure
derivative term can be dropped, but it is shown herein for a more general
case. The
Jacobian derivatives with respect to the primary unknowns are obtained as
follows:
np
aR (aQr, aQn)
-= an CIT
(e19)
axk kaqk a äq7
n=1
np nc
äR 0 Q
¨ = an xk (e20)
a CIT a qk
n=1 k=1
np
aR aQT,
ap n ap
,1
[00134] The above equations are obtained by considering the flow
through a
connection 330.
[00135] For each
node 324 and 325, mass balance equations are used. In particular
(see for example nodes 2 and 3 of model 694), for each wellbore node, the
total mass
flow rate of each component into the node is equal to the sum of the mass flow
rates of
the same component from each incoming connection and each perforation. Each of
these
nodes may or may not have perforations, and let perf denote the number of
perforations at
an exemplary node "r". The node r may have multiple incoming connections. Let
"inc"
denote the number of incoming connections 330; "out" denote the number of
outgoing
connections. The fluid flows into the node through perforations and
connections, so the
total mass flow rate into the node is:
-16-
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
(inc perf nc
1CIT jr qpkr)
p k
[00136] where qiir is the total mass flow rate from connection j into
node r, and
qpkr is the mass flow rate of component k from perforationp into the node.
[00137] Let x
denote the composition of component i in the incoming flow into a
node r through connection j; and let xir be the composition of component i in
the total
incoming flow into node r. For a component i, the component's mass flow rate
out of the
node is:
[00138]
inc per! nc
Xir qTjr
j p k
[00139] The incoming mass flow rate for component i is:
inc per!
xjjqj + qpir
[00140] The mass balance equation for component i is then
inc perf inc per! Tic
R xijrqTir + qpir ¨ Xir qTir 4- qpkr = 0
(e21)
j p k
[00141] and the corresponding Jacobian derivatives are
inc perf nc
OR
¨ ¨qTi =r = + qpkr) (e22)
aXir
p k
OR
qTir (e23)
tJXijr
OR
__________________________ = x= =r ¨ x =LT (e24)
agTir y
OR =11¨ xir. if k == i
(e25)
qpkr ( ¨Xir if k #
- 17-
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[00142] The equations (1) will include the mass balance equations for
the
components i = 1, ..., nc-1.
[00143] The example of Fig. 7 shows the Jacobian for two nodes such as
in model
680. One of the nodes is node 325, also shown as node 1, whose equations
correspond to
the rows marked (e21), (e21), and (e27). The corresponding variables are shown
in the
heading row as qn (total mass flow rate from node 1 to node 2) , x11
(composition of
component 1 at node 1), x21 (composition of component 2 at node 1).
[00144] The other node is a bottom-hole node 2, corresponding to the
last three
rows. The corresponding variables are qT2 (total mass flow rate out of node
2), x21
(composition of component 1 at node 2), x22 (composition component 2 at node
2). Fig.
7 is not limited to model 680 in applicability.
[00145] The mass balance equation may be simpler if there is only one
incoming
connection and one outgoing connection and the node has no perforations; the
equation
becomes:
R = xup ¨ xd, = 0 (e26)
[00146] where xup is the composition of component fat the upstream node
of the
incoming connection into node r, and xdu, is the composition of component i at
the node
r. See the rows marked with (e26) near the bottom of Fig. 7 for example.
[00147] The total mass balance equation equates the incoming and
outgoing mass
flow rates at each node r as follows:
inc Per f nc out
R = Iqrjr = 0 (e27)
p k 1
[00148] where gni is the total mass flow rate from the node r into outgoing
connection 1. The derivatives are
aR
= 1 (e28)
u qrjr
- 18-
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
R
________________________________ = 1 (e29)
a q Tlr
R
= 1 (e30)
u qpkr
[00149] The mass transfer between a cell 320 with a perforation 134 and
a
wellbore node 325 is governed by the perforation equation as follows (here we
assume
that the perforations are producing, not injecting; the index r identifies the
reservoir's cell
320 adjacent to the perforation):
phase
qpk = CA4:13 Anr Pnrzknr (e31)
n=1
6.43 = Pr Pw + Yr (Dw ¨ Dr) (e32)
[00150] Here qk is the mass flow rate from perforation p (cell 320)
into node 325;
C is a user defined constant; ilrir is the mobility of phase n in perforation
cell 320
identified by index r; p nr is the density of phase n in the same cell; zknr
is the mass
fraction of component kin phase n in the cell; Pr is the pressure in the same
cell; /314, is
the pressure at node 325; Dr and Dw are the depths of cell 320 and node 325
respectively; and yr is the hydrostatic pressure gradient, so that the product
?r(1w-Dr)
represents the difference in pressure due to the fluid head.
[00151] For the first three rows of Fig. 7, marked (e31), the variables
in equation
(e31) are expressed in terms of primary variables: pressure in perforated grid
block 320,
and mi- the mass of component i in perforated grid block 320.
[00152] The primary variables are:
[00153] 1. In the network domain (i.e. in wells, and thus in the
overlap between G1
and G2, and possibly in the surface network and thus in group G2):
[00154] qT (total flows on various connection);
[00155] --xi (compositions of components i on connections); and
- 19 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[00156] P (pressures).
[00157] 2. In the reservoir domain (i.e. in the reservoir cells 320,
thus in group
01):
[00158] mi is the mass of component i in the perforated grid block
320;
[00159] P (pressures).
[00160] 3. At the interface between these two domains (at perforation
nodes 325,
thus in group G1):
[00161] qpk (the mass flow rate of component k through perforation p).
[00162] In equation (8):
[00163] -- The Awr rows correspond to equations (e10), (e15), (e16), (e17),
(e18).
[00164] -- The Afw rows correspond to some of the same equations, and
possibly
other equations (not shown) for the surface network.
[00165] -- The Arr rows correspond to mass balance equations for each
cell, and a
volume balance for each cell (fluid volume = pore volume); these equations can
be well
known equations, and are not shown herein.
[00166] In the following example illustrating a solution of equations
(8) for the
well model 680, we will show how to eliminate fir and qpk, and how to reduce
the
resulting linear system, in which the primary variables are xi,P in the
network domain
and inc, P in the reservoir domain, to a pressure only system. This example
includes a
single phase with three components. The Newton-Raphson Jacobian matrix at the
network side is shown in Fig. 7 (Fig. 7 is not limited to model 680 but
applies to other
Ar,
models also). In other words, Fig. 7 shows the portion of equation (8). In
equation
Aõ
(e4), the functionf is treated explicitly (its value is the value obtained in
the previous
iteration, or is an initial value if this is the first iteration) so its
derivatives are zero. The
flow rates qpk at perforation connections and total flow rates ink at all
other network
connections are eliminated as shown in Fig. 8, and only the rows initially
based on the
mass balance equations remain.
[00167] In Fig. 7, the equation number for each row is specified in the
right
column. In Figs. 8-9, the rows have been added or subtracted as described
below, but the
- 20 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
right column provides the same equation numbers as in the corresponding rows
of Fig. 7
so the rows are easily identified. The cell 320 in model 680 is assumed to
have one
perforation, and the first three rows in Fig. 7 correspond to equations (e31).
The primary
variables are given at the top of the figures, so the variables in equations
(e31) are
expressed in terms of the primary variables, i.e. the masses m1 and m2 of
components 1
and 2 in the grid block 320 and the pressure P in the grid block. P1 in Fig. 7
is the
pressure at the wellbore node 1 (Ply in equations (e31)); P2 is the pressure
at the next
node 2 up the wellbore.
[00168] The c values in the first three columns in the (e31) rows in
Fig. 7 represent
the appropriate derivatives.
[00169] The values T1 , T2, T3 are:
phase
= ¨C Xnpnzkii
n=i
[00170] where k= 1, 2, 3. See equation (e31); the r index is omitted
for simplicity.
[00171] The variable qn is the mass flow rate on the connection between
nodes 1
and 2 (i.e. nodes 320 and 324BH) in model 680. The variable q2-2 is the mass
flow rate
from node 2. The variables x1 are compositions at different nodes; the first
index, i, is the
component number, and the second index is the node number. Thus, x12 is the
composition of component 1 at node 2. The component number nc=3, and the
primary
variables include the compositions for the first two components only.
[00172] In the row for equation (e4), the values a, h11, h21 depend on
the function
fin the equation. Suitable equations are well known. Of note, f can be treated
implicitly
or explicitly.
[00173] As noted above, np = 1. The variables P1 and P2 are pressures
at nodes 1
and 2 respectively.
[00174] In the row (e17), the constraint equation can be for node 2. In
this case, xk
indicates xk2, and qk are given by equation (e6), i.e.qk ¨ . The volume
flow rates
xk2qT2
Qn are on the connection from node 2 to the surface network. As noted above,
Figs. 7-10
have a wider applicability than model 680.
-21-
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[00175] Of note, the columns x12 , x22 (node 2 compositions) have zero
entries in
the mass balance rows (e21), (e21), (e27) for node 1. However, the columns x11
, x21
(node 1 compositions) have some non-zero entries in the mass balance rows
(three last
rows (e26), (e26), (e27)) for node 2. It is desirable to eliminate these non-
zero entries to
construct a pressure matrix from the mass balance equations.
[00176] In Fig. 8, the total mass flow rates q7-2 and the perforation
connection
flow rates qpk have been eliminated. Figs. 8-10 show only the mass balance
equation
rows for nodes 1 and 2. In Fig. 9, the mass balance equations at the
perforation node 1
(the first three rows, labeled respectively (e21), (e21), (e27)) have been
added to the
respective equations of the reference node (the bottom-hole node 2 in this
example), i.e.
rows 1, 2, 3 are added to respective rows 4, 5, 6 to cancel the composition
coupling terms
between these two nodes. Before the addition, the rows 4 and 5 are multiplied
by
[00177] Err Lk qpk
f
[00178] Now the columns for node 1 compositions x11, x21 have zero
entries for
node 2 (in the last three rows). The columns for node 2 compositions x12, x22
retain zero
entries for node 1 (in the first three rows).
[00179] At the reference node 2 (324BH), the composition terms only
remain at
the diagonal block as shown below:
[00180]
pert nc
qpk
p k
peal nc
pk
Op nP nc op op nc
aQn aQn) aQn
"n an 21 an (aqi acine n qk a q 2 a q)nc a qXk
k
[00181] We can eliminate the nonzero composition terms in the last row
by left-
multiplying the last three rows by the matrix
- 22 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
1
1
[(01 (02 1 I
[00182] where
=[qnp nP nc Fpert nc 1
T2 an (¨aQn ¨ ¨aQn) / Q n Xki1 I q pk (e33)
\dqi gm! n aqk
n=1 n=1 k=1 p k
np np nc
2 = [77,2 an (a Qn oQn) / 11 a n aQn
aq2 qnc aqk
n=1 n=1 k=1
pen f nc (e34)
p k
[00183] This means that the component mass balance equations (26) for
node 2 are
added to the total mass balance equation (27) for node 2 with respective
factors 001 and
002 to eliminate the node-2 component terms xi2 and o.)2 in the total mass
balance
equation. The node-1 component terms remain zero.
[00184] The resulting matrix is shown in Fig. 10. Thus, the pressure matrix
at the
well side based on reference node 2 is obtained, as shown in the last row in
Fig. 10. In
other words, for the well nodes, only the pressure variables remain. In other
embodiments, similar results are achieved by using node 1 as the reference
node, and in
particular all the component terms are eliminated in the total mass flow
equation for node
1.
[00185] With regard to the embodiment of Fig. 10, since
Pert rtc
q T2 = qpk (e35)
p k
[00186] we have
- 23 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
np np nc
aQ aQ aQ
(01=Iaii( ¨ Xk
aq1 aq aqk (e36)
n=1 nc n=1 k=1.
np nP nc
(1)2 aqz =Ian (¨a Q ¨ ¨aQn (e37) aqn,
n=-1 n=1 k=1 k
[00187] Next we explain the physical meaning behind this pressure
matrix. For an
oil rate constraint in a black-oil system
= 1,a2 = 0,a3 = 0
[00188] The denominator in equations (e36) and (e37) can be rewritten
as
fnc
aQi =v = aQiqk Q1
(e38)
=1qk Z"'"J aqk cIT qT
k k=i
[00189] and the multipliers in equations (e36) and (e37) can be
simplified as
Qi qT
= ¨ (e39)
ach Qi
a(21 CIT
(62 = (e40)
ach
[00190] Plugging equations (e39) and (e40) into the last row in Fig.
10, the
pressure coefficient at the reference node 2 (the P1 coefficient) can be
reorganized as
follows:
Pww = ¨T1¨ T2 ¨ T3 ¨ CO 1[(x11 1)T1 X11T2 X11T3]
¨ [X21-Ti + (x21 ¨ 1)T2 X2iT3]
= ¨(1 + (01x11 + (02x21)(T1 + T2 + T3) + (4) Ti W2T2
= 601 T1 + (02 T2 (e41)
q, a Q, aQ,
=-( __________ T, ,T,")
Qi aq1 1 aq2
- 24 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[00191] where
phase
Tk = -C Anpnzkn (e42)
n=i
[00192] See equation (e31); the r index is omitted for simplicity.
[00193] It is evident that the pressure coefficient (P1 coefficient) is the
sum of
mobility of each component weighted by its contribution in the total oil phase
flow rate,
which is consistent with a simple well model. The only difference is that the
pressure
matrix is scaled by the constant - glyQi.
[00194] The coupling matrix between the reservoir and the well (the
left three
columns in the bottom row of Fig. 10) is
Pwr
qT [aaql 2 1 aQ1 a(21 aQi aQi 3Q1 (e43)
=
¨ + c21 a
C12 + c22 cip + c2p
Qi i aq2 aqi aq2 aql q2
[00195] As can be seen in this example in Figs. 7 and 8, the
composition coupling
teiiiis xii in the off-diagonal sub-matrix of downstream node 2 are the
negative of the
composition coupling terms in the diagonal sub-matrix of its upstream node 1
when.f in
the hydraulics equation (e4) is treated explicitly. ("Off diagonal" denotes
herein the
compositions at one node in the rows for another node; "diagonal" denotes
composition
at a node in the rows of the node itself) Of note, the rows (- EPp"f qpk,
0) and (0,
_ zpjerf zi7 qpk) in the (e21) lines of Fig. 7 are the same, up to a constant
multiple, as rows
(1, 0) and (0, 1) in the (e26) lines. Those composition coupling terms in the
first nc-1
rows (the component mass balance equations) for each node are the total
incoming flow
rate at the downstream node. This should be true for all nodes in a gridded
well if there is
no cross-flow or constraints applied to those nodes but the bottom-hole node.
We will
show this on the following example.
[00196] Fig. 11 shows the Jacobian matrix at the wellbore side (the
portionww
[A'}
in equation (8)) after elimination of the total flow rates qa and the
perforation
- 25 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
connection flow rates qpn (i.e. at the stage of Fig. 8) for the model 694 with
a rate
constraint applied to node 4. The first three rows correspond to the mass
balance
equations (21), (27) for node 1; the next three rows correspond to the mass
balance (MB)
equations (21), (27) for node 2; the next three rows correspond to the MB
equations (21),
(27) for node 3; the last three rows correspond to the mass balance equation
(21), (27) for
node 4. For each node, the first two MB equations are for respective
components 1 and 2;
the last MB equation is for the total flow through the node.
[00197] The columns P1, P2, P3 in the left half denote pressures in the
corresponding perforation cells 320. The same symbols in the right half refer
to pressures
at the respective nodes 325.
[00198] Figs. 12, 13, and 14 show how the pressure matrix is generated.
Figs.
12-14 show only the rows corresponding to the mass balance equations for node
4 (last
three rows in Fig. 11). For these figures, we define:
= (x1) ¨ 1.)Tij + x11T21 +
T2j= x21T21 + (x2.1 ¨ 1)T zi X21T 3j
-1; ¨ T2,1 ¨ T
[00199] where j is the perforation index, and T is calculated using
equation (e42).
[00200] In Fig. 12, the MB equations for nodes 1, 2, and 3 (first nine
rows) were
added to the equations at the bottom-hole node 4 (three last rows) to obtain
the matrix of
Fig. 12. More particularly, the first equations (for component 1) for nodes 1,
2, 3 are
added to the first equation for node 4; the second equations (for component 2)
for nodes
1, 2, 3 are added to the second equation for node 4; and the third equations
(for total
flow) for nodes 1, 2, 3 are added to the third equation for node 4. In the
composition
aR
coupling terms, corresponding to the ¨ columns, the terms
ax
--Ax remain which are similar to each other. In these terms, the values IT, E
are the
hydraulics derivatives with respect to the composition and total flow rate,
and are
calculated using equations (e13) and (e14), respectively. Ax is the difference
in
component mass composition between a node and its upstream node. As in the
lumped
- 26 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
well model, if the pressure gradient is treated explicitly, i.e. the function
f in equation (e4)
is evaluated using results from the previous iteration, all those derivatives
are zeroes, so
all the composition columns in Fig. 12 are zeroes. Then all the pressure
columns are
lumped (added) to the last column, corresponding to the pressure at the bottom-
hole
node, as shown in Fig. 13. Finally we use the same method as in the lumped
well model
680 to eliminate the diagonal composition terms, corresponding to the last
three columns
in Fig. 13, and the last row in Fig. 14 is taken as the pressure coefficients
for the well
node 4, which has the pressure as the sole unknown for node 4 at the well
side. The
pressure unknowns are decoupled from other unknowns at the well side.
1002011 In the previous two examples, if the functionf in hydraulics
equation (e4)
is treated explicitly or the pressure gradient is taken constant, we should
have
SP i = 8Pre1 (e44)
and
aR aR
ref (e45)
aPt aP
1002021 where ref denotes the reference node. We can apply the same
approach to
the coupling terms from the reservoir to the well side, and move those terms
to the
ap,
column corresponding to the reference node, as shown in Fig. 16. Note that in
those
coupling terms, aR = 0 for producing perforations. At the wellbore side, the
mass
ax
coupling terms ¨aR can be eliminated using the reduced mass balance equations
at the
am
reservoir side, and new coupling terms between network and reservoir are
introduced, as
shown in Fig. 16.
[00203] More particularly, Fig. 15 shows a global matrix structure for
implicit
formulation with implicit wells and surface network. Thick lines divide the
matrix into
nine blocks which correspond to the nine blocks of the matrix of equation (7),
i.e. the
blocks Arr, , A , etc. The "x" terms indicate possible non-zero elements. The
matrix of
- 27 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
Fig. 15 has a number of cells some of which are empty (corresponding to zero
elements)
and the others contain four x's, indicating that the cells may have any
structure.
[00204] Fig. 16 shows the matrix structure for the case of Fig. 15
after elimination
of mass coupling terms in the producing perforation connections. Only four
blocks are
shown, as in the matrix for equation (8), i.e. the blocks Arr, , Arw , Awr, ,
A . In the
Aww block, the pressures are the only non-zero elements. The "y" terms in the
A.wr block
indicate fill-ins, i.e. possible non-zero elements generated by linear
combinations of
equations.
[00205] The coefficient matrix, which is derived from mass balance
equations
(MBEs), at the reservoir side can be reduced to a pressure only matrix, shown
in Fig. 17
and having the following form:
Prr = RArrP (e46)
[00206] where the prolongation matrix P and the restriction matrix R
are defined as
follows:
Fel
P = (e47)
eN1
[00207] where
01
[00208] e ¨
1
[00209] The restriction matrix R is a diagonal matrix as follows
[LT'
R = PT (e48)
4,1-1
[00210] where L is the lower matrix of an LDU factorization of the
diagonal sub
matrix Arr, i.e.
LiDiU1
diag(Aõ) = (e49)
LNDNUN
[00211] We finally get a global pressure linear system including
reservoir and
wells, as shown in Fig. 17. More particularly, we get:
- 28 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[Prr Prwir = [dr (e50)
]
Pwr Pww U141 dw
[00212] Due to the Jacobian processing described above, the
approximation errors
in generating the approximate solutions by Newton-Raphson or other methods are
hopefully reduced, which may lead to a better pressure solution and hence a
better overall
solution of equation (7).
[00213] In the previous two examples, we assume there is a rate constraint
at the
well connection. When it is a pressure constraint, as shown in equation (e15),
or there is
no constraints at all (a constraint might be at the top of the surface
network), the method
to construct the pressure matrix is similar and simpler. The multipliers
defined in
equations (e33) and (e34) are all zeroes.
[00214] The Jacobian conditioning described above can be used in many
simulation schemes. One example is illustrated in Fig. 18. At step 1810, the
equations (1)
are selected, i.e. F(Vari = 0, to include mass balance equations for a
selected set of
wellbore nodes (such as nodes 1 and 2 in model 680, or nodes 1-4 in model
694). In some
embodiments, these nodes and the connections between them form a connected
graph. In
some embodiments, the nodes are consecutive nodes in a path of fluid. In some
embodiments, the set of nodes includes nodes of surface network 130 and/or
reservoir
cells 320. In the examples above only one set of nodes (such as nodes 1 and 2)
was
considered, but multiple sets of nodes can be provided. In each set, a
reference node is
chosen, and the mass balance equations of the other nodes are added to the
mass balance
equations of the reference node, and other processing is performed for each
set as
described above, to eliminate compositions in the reference node mass balance
equations.
[00215] The mass balance equations include, for each node in the set,
equations for
individual components and/or equations for total mass balance. Perforation
equations
may also be present.
[00216] The equations (1) are then solved for each time period ti by the
Newton-
Raphson method as follows. At step 410 (same as in Fig. 4), an initial
iterative value
Var.NT is chosen. At step 414 (same as in Fig. 4), the Jacobian is determined
at the
Var.NT value, and the RHS term "b" is determined for the equation:
[00217] JF(Var.NT)Vari = b (e50A)
-29-
CA 02873406 2019-11-12
WO 2013/180709 PCT/US2012/040023
[00218] At step 1820, the Jacobian is processed to eliminate some non-
pressure
coefficients by combining the mass balance equations in each set of nodes, and
possibly
combining different columns, as illustrated above in connection with Figs. 7-
17 and
further explained below. Let us denote the resulting "modified Jacobian" as JM
(i.e. JM
may include the matrix of Fig. 10 or 14), and denote the corresponding
variables as VarM
(the variables can change since pressure columns may have been added and some
columns may have been dropped; e.g. the columns corresponding to the total
flow rates in
the connections and the perforation flow rates may have been dropped); and
denote the
corresponding right-hand side as bM. The equation (e50A) then becomes:
[00219] JM(VarM.NT)VarMi,k = bM (e50B)
[00220] At step 418 (as in Fig. 4), the computer system solves the
equations
(e50B), possibly (but not necessarily) using a Schwarz method illustrated in
Fig. 5. Let us
denote the solution VarM*; ki.
J
[00221] At step 1824, the original variables' values Var*j,k,1 are
obtained from
VarM*.j,k J. The values Var*i jci are thus a solution of equation (e50A), and
this solution
is used as the next iterative value for Var.NT, and control returns to step
414, or
alternatively the Newton-Raphson method is stopped (steps 422, 426 are as in
Fig. 4).
[00222] Many variations are possible. For example, the modified
Jacobian
computation steps 414, 1820 can be performed in generalized form, without
plugging in
the value Var.NT, as in fact is illustrated in Figs. 7-14. Then at each
iteration, the value
Var.NT can be plugged into the generalized expression for the modified
Jacobian and
modified right-hand side to determine the equation (e50B).
[00223] Fig. 19 illustrates some embodiments of the Jacobian processing
step
1820. At step 1910, for each set of nodes, a reference node is chosen. In the
examples
discussed above in connection with models 680 and 694, the reference node was
the
bottom-hole node, but any other node in the set could be used instead. In some
embodiments, the reference node is such that the flow into the reference node
is made up
of flows from the other nodes of the set, or such that the flow out of the
reference node is
made up of flows into the other nodes of the set.
[00224] At step 1920, the total mass flow variables qT are eliminated (as
in Fig. 8)
- 30 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
for the connections within each set. Some rows can be removed from the
Jacobian.
[00225] At step 1930, within each set of nodes, the mass balance
equations for the
non-reference nodes are added to the corresponding mass balance equations for
the
reference node. In particular, the mass balance equations for each component k
are added
to the reference node's mass balance equation for the same component k; the
total mass
balance equations (for the total flow through each non-reference node) are
added to the
reference node's total mass balance equation. The equations may have to be
multiplied by
some factors before the addition. The factors depend on the equations (1) and
the choice
of primary variables. The addition is performed so as to eliminate the off-
diagonal mass
flow variables (e.g. compositions) in the reference node's mass balance
equations.
[00226] At step 1940, the reference node's MB (mass balance) equations
are
combined to eliminate the diagonal compositions in the reference node's MB
equations.
This operation can be expressed as left-multiplying the reference node's MB
equations by
the matrix L-1 where L is the lower matrix in the LU (lower-upper)
factorization of the
square matrix corresponding to the diagonal pressure and composition portion
of the
reference node's MB equations. An example of the L-1 matrix is described above
in
connection with equations (e33), (e34), i.e. the matrix:
1
[00227] L-1 ¨ 1
(01 6)2 1 I
[00228] At step 1950, the pressure columns may be combined as explained
above
in connection with equations (e44), (e45).
[00229] In some embodiments, some of these operations are omitted. For
example,
step 1950 and/or 1940 can be omitted.
[00230] As noted above, there are many ways to solve the modified
Jacobian
equation at step 418 in Fig. 18. An exemplary method, using the multiplicative
Schwarz
method of Fig. 6, is a two-stage preconditioning method performed as follows.
The
network equations (9) are solved at the pre-iteration phase (for the initial
iterative value
x.sc shown at step 620 in Fig. 6), and equations (9) are solved again at the
end of each
iteration to maintain the network residual (the residual of (9)) always at
zero in the
GMRES calculation. The preconditioning procedure consists of the following
steps:
- 31 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[00231] Pre-iteration
[00232] 1) Initialize xr = 0
-1 b
[00233] 2) Solve equation (9), i.e. [xxwil ¨ [AAwwAAwri [bw], to
determine xw
fw ff f
and xf
[00234] 3) Update residual in equation (8):
[00235]
br Ar , 0
¨ A fxf
[bw11-= ¨[Aw Awwl[xd w
[00236] Schwarz Iteration, by k, until convergence
[00237] 1) Solve equation (8) using a CPR method (Constrained Pressure
Residual) using the pressure matrix (e46) as explained below.
[00238] 2) Update residual in equation (9):
[00239] [bwl h
_wl r www Awfi rsx_
A,Sx,
bf [ k Afw Aff I- 0
b
[00240] 3) Solve equation (9), i.e. rwl =A ww Awf-I wl
Xf k+1 Afw Aff bf k+1
[00241] 4) Update residual in equation (8)
[br 1 br Arr Arw 0 A 1.
W] k+i [bwlk [A wr A wwila wl
O ¨Alwiuxf
[00242] 5) Compute the solution change Sxr= 8,cr + PSur where P is as
in (e47).
Compute the residual change a, = bir"-1 ¨ O. Feed 5x,, Sx,õ 5x1, art GMRES to
update the solutions xr, xw, xf.
[00243] Step 1) above, solving equation (8) by CPR method:
[00244] 1) Update the residual in equation (e50):
[00245] [dr br
W
k] = L vc
bj where r is a restriction matrix defined below.
[00246] 2) Solve equation (e50):{pt pPwr Situ,' idd wr
[00247] 3) Map the pressure solution to the full solution:
[00248]6x, 5x, our
[6 x [Sxd + 17 [ ouj
[00249] where )7 a prolongation matrix defined below.
- 32 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[00250] 4) Update residual of equation (8):
[00251] rbrl [hi.] rAõ Arwi r Sxr1
Lbw k+1 [bulk lAwr Avvvid LS Xyvi
[00252] 5) Solve equation (8):
[00253] r Aõ Arwl[ xr1 1'r1
LAõ Aj I4XW]k+l 1-bw-lk+1
[00254] The prolongation and restriction matrices are as follows:
[00255] tu is the restriction matrix
iRr
(e51)
[RC Rj
[00256] where Rr is the restriction matrix on the reservoir domain, and
is defined
in Equation (e48). R, is the restriction matrix at the well side, and is
defined as:
Wi
= (e52) ,
Wnwi
[00257] and wt = [col (02 1 ... coi (02 11 where (01 and (02 are
determined using equations (e33) and (e34). Here nw is the number of wells. Re
is an
operation matrix accounting for the elimination of mass coupling terms as
shown in Figs.
14 and 16.
[00258] n is the prolongation matrix:
_ [Pr
Pwl (e53)
[00259] Pr is the prolongation matrix on the reservoir domain as
defined in (e47).
Pis the prolongation matrix on the well side:
P, = (e54)
EõI
[00260] where EL for the well i is defined as:
- 33 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
0
I
1 u 1
Ill
[00261] Ei=lil
101
1.01]
[00262] Here the O's correspond to the columns with composition
unknowns, and
the l's correspond to the columns with pressure unknowns.
[00263] The invention is not limited to the embodiments described
above. The
invention includes modified Newton-Raphson and Schwarz methods. For example,
in
some Newton-Raphson modifications, the same Jacobian is used for a number of
iterations. The invention is not limited to particular variables, and while
the examples
above used the compositions of nc-1 components and the mass total flow rates
other
variables can be used instead, e.g. mass flow rates for all the components can
be used.
Further, while the embodiments described above do not use volume flow rate
equations
(e.g. volume balance equations) on the network side, other embodiments use
such
equations on the network side and use volume flow rates and volume
compositions
(volume fractions of individual components) instead of (or together with) mass
flow rates
or mass compositions. Volume flow rate terms are eliminated from volume
balance
equations using techniques described above for mass balance equations.
[00264] Some embodiments provide a computer implemented method for
determining a plurality of values (e.g. pressures, compositions, flow rates,
etc.) in a
facility comprising an oil and/or gas reservoir and a plurality of wellbores
each of which
has been made or is to be made, each wellbore providing one or more passages
for fluid
flow into or out of the reservoir. The method comprises:
[00265] obtaining, by a computer system, a data model (e.g. 680, 690,
694) which
represents the reservoir and each wellbore, the model comprising one or more
nodes at
said facility, the nodes comprising one or more sets of nodes, each set
comprising a
single node or a plurality of interconnected nodes, at least one node in each
set being in a
wellbore. In some embodiments, a set contains just one node, and diagonal
composition
terms are eliminated as in step 1940; step 1920 may be omitted. Different sets
may or
may not have nodes in respective different wellbores or wellbore segments.
- 34 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[00266] The method further comprises obtaining, by the computer system,
a
system of linear algebraic equations interrelating said values, said values
being
represented as variables in the equations, wherein said values include a
pressure for each
node in each set and include flow rate values (e.g. mass flow rates or mass
compositions
or volume flow rates or compositions) defined by flow rates into or out of the
nodes in
each set, and for each node the equations include a set of one or more
equations with
values at the node, wherein for at least one node of each set of nodes, the
set of one or
more equations comprises at least one equation with a pressure at the node and
one or
more flow rate values at the node;
[00267] performing, by the computer system, a linear transfotmation of the
system
of equations, the linear transfotmation comprising a linear transformation of
the
equations of each set of equations to eliminate at least one flow rate value
from at least
one equation in each set of equations, the linear transformation providing a
transformed
system of equations; and
[00268] solving, by the computer system, the transformed system of
equations for
said values.
[00269] Some embodiments provide a computer implemented method for
determining a plurality of values in a facility comprising an oil and/or gas
reservoir and a
plurality of wellbores each of which has been made or is to be made, each
wellbore
providing one or more passages for fluid flow into or out of the reservoir,
the method
comprising:
[00270] obtaining, by a computer system, a data model which represents
the
reservoir and each wellbore, the model comprising one or more nodes at said
facility, the
nodes comprising a set of one or more well nodes each of which is a node in at
least one
wellbore with no perforations along the wellbore between the node and a ground
surface;
[00271] obtaining, by the computer system, a system of linear algebraic
equations
(e.g. Jacobian based) interrelating said values, said values being represented
as variables
in the equations, wherein for each well node, said values include at least one
value at the
well node;
[00272] solving, by the computer system, the system of the equations for
said
values, wherein solving the system comprises:
- 35 -
CA 02873406 2016-03-09
[00273] solving, by the computer system, equations obtained from a
first
subsystem of the system of said equations (for example, the first subsystem
may be
equations (8)), to obtain values including at least one value at each well
node of said
set and at least one value in the reservoir, or for differences between such
values and
their approximations (e.g. the differences may be Sxr,dx,,Sxf as indicated
above);
[00274] solving, by the computer system, equations obtained from a
second
subsystem of the system of said equations (e.g. the second subsystem may be
equations (9)), to obtain values including at least one value at each well
node of said
set and at least one value in the reservoir, or for differences between such
values and
their approximations;
[00275] wherein the first subsystem overlaps with the second subsystem,
and at
least one value of a well node of said set, or the corresponding difference
between the
value and its approximation, is solved for in both solving the first subsystem
and
solving the second subsystem;
[00276] wherein the method further comprises constructing, by the computer
system, a solution of the system of the equations from solutions of the first
and second
subsystems.
[00277] The invention includes computer programs operable to cause a
computer system (such the system of Fig. 2) to perform any of the methods
described
above. The invention includes computer program products, e.g. computer
readable
media or computer systems, comprising a computer program operable to cause a
computer system to perform such methods. The computer programs can be stored
on
computer readable media, and can be transmitted over a computer network to a
computer system and downloaded onto such media. The invention includes methods
comprising network transmission of such computer programs. Other embodiments
and variations are within the scope of the invention, as defined by the
appended
claims.
[00278] Some information related to the above disclosure can be
obtained from
the following documents (no representation is being made as to whether or not
the
documents are prior art with respect to the present application):
[00279] - U.S. patent application publication 2012/0022841 Al (January
26,
2012), inventor Appleyard.
- 36 -
CA 02873406 2019-11-12
WO 2013/180709
PCT/US2012/040023
[00280] ¨ Martin J. Gander et al., "An optimal block iterative method
and
preconditioner for banded matrices with applications to PDEs on irregular
domains",
Report 10-05-21, pages 1-27, May 2010, available in the World Wide Web at
http://vvww.math.temple.edu/¨szyld.
[00281] ¨ PCT application WO 2011/136861 Al (November 3, 2011), applicant
EXXONMOBIL Upstream Research Company.
- 37 -
