The MATLAB Reservoir Simulation Toolbox

B ard Skaestad

SIAM Geosciences, Long Beach, March 2124, 2011

The Matlab Reservoir Simulation Toolbox (MRST)

The toolbox has the following functionality for rapid prototyping of solvers for ow and transport: grid structure, grid factory routines, input/processing of industry-standard formats, real-life and synthetic example grids petrophysical parameters and incompressible uid models, conversion routines to/from SI and common eld units, very simplied geostatistical routines routines for setting and manipulating boundary conditions, sources/sinks, and well models reservoir state (pressure, uxes, saturations, compositions, . . . ) visualisation routines for cell and face data (scalars) Download http://www.sintef.no/MRST/ Version 2011a was released on the 22nd of February, 2011, and can be downloaded under the terms of the GNU General Public License (GPL)

2 / 21

Software goals
General framework for ow and transport in porous media Special focus on unstructured grids and multiscale methods Research: Development of new solvers and discretization schemes Application to new scientic/engineering problems Preserve know-how and promote reuse of results from previous research Contribute to accelerating production of new research (by our peers) Benchmarking compare dierent methods on standard test problems Students: Develop students intuition of porous media ow Make it easy to test, compare, and extend existing methods Textbook (with worked examples) in preparation

3 / 21

Why open-source and why in Matlab?

First of all, prototyping in a scripting language is much more eective than in traditional compiled languages (C/C++/FORTRAN) Explore alternative algorithms/implementations close to mathematics Gradually replace individual (or bottleneck) operations with accelerated editions callable from matlab Direct access to matlab environment and prototype whilst developing replacement components More experienced in Matlab/Octave than in e.g., Python, Why free and open-source software: Research funded by Norwegian Research Council should be freely available Combined with publications: a means to collaborate and disseminate results, while protecting IP rights Support reproducible research and promote replicability

4 / 21

How is MRST designed?

The fundamental object in MRST is the grid: Data structure for geometry and topology Several grid factory routines Input of industry-standard (proprietary) format(s) Physical quantities dened as dynamic objects in matlab Properties of medium (, K, net-to-gross, . . . ) Reservoir uids (, , kr , PVT, . . . ) Driving forces (wells, boundary conditions, sources) Reservoir state (pressure, uxes, saturations, etc) All MRST operations accept, manipulate and produce objects of these types. Physical quantities are assumed to be in SI units.

5 / 21

How is MRST designed?

MRST core routines for creating and manipulating grids and physical properties basic ow and transport solvers (sequential splitting) for incompressible and immiscible ow Functionality is stable and not expected to change in future releases Modules Similar to matlabs toolboxes. Implements more advanced solvers and tools: adjoint methods, experimental multiscale, fractures, MPFA, upscaling black-oil models, three-phase ow, vertically integrated models, . . . streamlines, (ow-based) coarsening, . . . Octave support, C-acceleration, . . . Some are stable. Some are constantly changing to support ongoing research. New modules initiated by others are much welcome
5 / 21

Grids in MRST
Complex reservoir geometries
Corner point: Tetrahedral: PEBI:

Accurate modelling of features such as faults, fractures, or erosion requires grids that are exible with respect to geometry. Industry-standard grids are often nonconforming and contain complex grid-cell connectivities, as well as skewed and degenerate cells There is a trend towards unstructured grids with general polyhedral cells Standard discretization methods produce wrong results on skewed and rough cells condition numbers The representation of (un)structured grids aects the eciency of numerical methods

6 / 21

Grids in MRST
Cell geometries are challenging from a discretization point-of-view
Skewed and deformed blocks: Non-matching cells:

Many faces:

Small interfaces:

Dicult geometries: (Very) high aspect ratios:

800 800 0.25 m

7 / 21

Grids in MRST
All grids are assumed to be unstructured

Basic representation: cells

2
6 5

nodes
cells.faces = 1 10 1 8 1 7 2 1 2 2 2 5 3 3 3 7 3 2 4 8 4 12 4 9 5 3 5 4 5 11 6 9 6 6 6 5 7 13 7 14 : :

faces
faces.nodes = 1 1 1 2 2 1 2 3 3 1 3 4 4 1 4 7 5 2 5 3 6 2 6 6 7 3 7 4 8 3 8 5 9 3 9 6 10 4 10 5 : : faces.neighbors = 0 2 2 3 3 5 5 0 6 2 0 6 1 3 4 1 6 4 1 8 8 5 4 7 7 8 0 7

6
9

4
8 7

12

2
1

1
10

7 14

3
4
3

8
11

13

5
4

Choices in grid representation guided by utility and convenience in low-order nite-volume methods. Available geometric information typically limited to centroids, normals, areas, and volumes
8 / 21

Grids in MRST
Examples of basic grid types

% Rectillinear grid dx = 10.5cos((1:0.1:1)pi); x = 1.15+0.1cumsum(dx); G = tensorGrid(x, sqrt(0:0.05:1)); plotGrid(G);

% Extrude a standard MATLAB dataset load seamount g = triangleGrid([x(:) y (:)]); P = pebi(g); V = makeLayeredGrid(P, 5); plotGrid(V), view(40, 60), axis off

% Make and read a simple Eclipse input le grdecl = simpleGrdecl([20, 10, 5], 0.12) G = processGRDECL(grdecl); plotGrid(G,'FaceAlpha',0.8); plotFaces(G,find(G.faces.tag>0), 'FaceColor','red'); view (40,40), axis off

9 / 21

Grids in MRST
Examples of exible gridding strategies

Hybrid grid areal grid consisting of the following components radially rened grid at wells Cartesian along boundary hexahedral in interior polyhedral transition cells extruded to 3D along vertical lines

Hierarchical grid triangular cells adapted to a curved line extruded to 3D, with throw along curved surface coarse grid constructed by collecting cells from ne grid

9 / 21

Discretization of ow equation
General family of conservative methods for: v = q, v = Kp Basic formulation u i = T i ( ei p i i ) , ei = (1, . . . , 1)T

pi the pressure at the center of cell i ui the vector of outward face uxes i the vector of face pressures Ti the one-sided transmissibilities
pi

c ik Ak nk k

Special cases: The standard two-point method: T ii = ni K ci /|ci |2 Multipoint ux-approximation methods (MPFA) Mixed nite-element methods Mimetic methods
10 / 21

Discretization of ow equation
Linear system: mixed hybrid form
ij

B 4C T DT

C 0 0

2 3 D v 0 0 5 4p5 = 4g 5 , 0 0

32

pi

vij

pj

B denes an inner product. The matrix blocks read, Z Z Z bij = i T 1 j dx, cik = k i dx, dik =

|i nk | dx

11 / 21

Discretization of ow equation
Linear system: mixed hybrid form
ij

B 4C T DT

C 0 0

2 3 D v 0 0 5 4p5 = 4g 5 , 0 0

32

pi

vij

pj

B denes an inner product. The matrix blocks read, Z Z Z bij = i T 1 j dx, cik = k i dx, dik =

|i nk | dx

Positive-denite system obtained by a Schur-complement reduction ` T 1 D B D F T L1 F = F T L1 g, F = C T B 1 D, L = C T B 1 C.

Reconstruct cell pressures and uxes by back-substition, Lp = q + F T , Bv = Cp D.

11 / 21

Discretization of ow equation
Herein: a mimetic method, Brezzi et al., 2005

M u = ( ep )

u = T ( ep )
T

Requiring exact solution of linear ow (p = x a + k): MNK = C NK = T C

C vectors from cell to face centroids. N area-weighted normal vectors Family of schemes (given by explicit formulas): M=
T 1 CK 1 C T + Q N S M QN |i | T 1 T = N KN T + Q C SQC |i |

T Q N is an orthonormal basis for the null space of N , and S M is any positive denite matrix. Herein, we use null-space projection T P N = QN S M QN = I QN QN

12 / 21

Discretization of ow equation
Example: grid orientation eects Homogeneous domain with Dirichlet boundary conditions (left,right) and no-ow conditions (top, bottom) computed with three dierent pressure solvers in MRST.
TPFA MFD MPFA

100 100

Homogeneous permeability with anisotropy ratio 1 : 1000 aligned with the grid.
13 / 21

Discretization of ow equation
Example: lack of monotonicity Homogeneous domain with Dirichlet boundary conditions (left,right) and no-ow conditions (top, bottom) computed with three dierent pressure solvers in MRST.
TPFA MFD MPFA

10 10

Homogeneous permeability with anisotropy ratio 1 : 1000 rotated by /6.

13 / 21

Application examples
A (very) simple ow solver
K p + g z

v = q,

v=

Vertical well and Dirichlet boundary

% Grid and rock parameters nx = 20; ny = 20; nz = 10; G = computeGeometry(cartGrid([nx, ny, nz])); rock.perm = repmat(100 millidarcy, [G.cells.num, 1]); fluid = initSingleFluid('mu', 1centipoise, ... ' rho ' , 1014kilogram/meter3); gravity reset on % Fluid sources and boundary conditions c = (nx/2ny+nx/2 : nxny : nxnynz) .'; src = addSource([], c , ones(size(c)) ./ day ()); bc = pside ([], G , ' LEFT', 10barsa()); % Construct components for mimetic system S = computeMimeticIP(G, rock, 'Verbose', true); % Solve the system and convert to bars rSol = initResSol(G, 0); rSol = solveIncompFlow(rSol,G,S,fluid,'src',src,'bc' , bc ); p = convertTo(rSol.pressure(1:G.cells.num), barsa() ); From tutorial: simpleSRCandBC.m

Source term and boundary condition

Pressure distribution
14 / 21

Application examples
Modelling faults in consistent schemes
Water cuts from producer #1

Faults modelled as internal boundaries, with internal jump conditions u f =

Tf (f

0.5 0.4

without multiplier

f )

0.3 0.2 0.1 0 1

Gives an extended hybrid system. In addition, method to convert TPFA multipliers to fault transmissibility Tf

with multiplier

+ = TPFA solution

= mimetic

1000

100

10

0.1

15 / 21

Application examples
Multiscale module: bypassing the need for upscaling?

Key idea of multiscale methods: Local decomposition P p(x) = i pi i (x) P v (x) = ij vij ij (x) i coarse grid block. ij = i j Multiscale basis functions:

homogeneous (RT0)

heterogeneous

16 / 21

Application examples
Multiscale module / black-oil module: primary production
Rate in well perforation (m3 /day)
560

Shallow-marine reservoir (realization from SAIGUP) Model size: 40 120 20 Initially lled with gas, 200 bar Single producer, bhp=150 bar Multiscale solution for dierent tolerences compared with ne-scale reference solution.

555

Reference 5102 5104 5106

550

5107

545

540

535 100

200

400

600

800

1000

17 / 21

Application examples
Coarsening module: (ow-based) coarsening by amalgamation
Amalgamation of cells ow-adapted grids simple and exible coarsening adaptive gridding schemes ecient model reduction Algorithmic primitives working on a partition vector segment (ow) indicator into bins merge small cells rene large cells intersect partitions sanity checks to ensure connected blocks, etc great exibility Triangular grids: Cartesian grid:

18 / 21

Application examples
Coarsening module: (ow-based) coarsening by amalgamation

Dynamic adaption

Dierent partitioning:

Adapting to geology

18 / 21

Application examples
Coarsening module: (ow-based) coarsening by amalgamation

Flow-adapted vs METIS:

Model reduction of real-eld model:

ne grid 11 864 cells

Watercut curves
1 0.9 0.8

Reference solution 1581 blocks 854 blocks 450 blocks 239 blocks 119 blocks

ow-based 127 blocks

0.7 0.6 0.5 0.4 0.3 0.2 0.1

METIS 175 blocks

0 0 0.1 0.2 0.3

Pore volume injected

0.4

0.5

0.6

0.7

0.8

0.9

18 / 21

Application examples
Upcoming adjoint module: production optimization
Specialized simulator: using dierent grids for pressure and transport Multiscale pressure solver:

Water-ood optimization:

Reservoir geometry from a Norwegian Sea eld

14 12 Cum.Prod. [m3] 10 8 6 4 2 0 0 2 4 6 Time [years] 8 10 x 10
6

Oil initial Oil optimized Water initial Water optimized

Transport on ow-adapted grid:

Forward simulations: 44 927 cells, 20 time steps, < 5 sec in Matlab

19 / 21

Application examples
Vertically integrated module: educational and research tool

Vertically integrated module: part of a forthcoming numerical CO2 laboratory

20 / 21

Application examples
Vertically integrated module: educational and research tool

Vertically integrated module: part of a forthcoming numerical CO2 laboratory

20 / 21

Application examples
Vertically integrated module: educational and research tool

Vertically integrated module: part of a forthcoming numerical CO2 laboratory

20 / 21

Application examples
Vertically integrated module: CO2 migration at Sleipner
3D model made by Statoil to study CO2 in the top layer of Utsira. MRST used to study physical assumptions and numerical methods: Simple to modify the code Flexible and simple upscaling (homogeneous model) Fast response time, simulation 210 min in Matlab Seismic data (2006) / 3D simulation (tough2) VE simulation

Fra: Chadwick, Noy, Arts & Eiken: Latest time-lapse seismic data from Sleipner yield new insights into CO2 plume development, Energy Procedia (2009), 21032110.

21 / 21