8.
NUMERICAL CONTROLS
In order the simulation to be performed the flow equations have to be solved. These equations are partial non
linear (terms existed in the equation depend on the unknown value) differential equations and are solved
numerically by the simulator as the direct solution of these equations is too difficult and requires huge
computational times. The numerical method that the simulator uses to approximate the exact solution is the
Newton-Raphson method.
The first step for the solution of the flow equations is their discretization in time (eg. implicit) and space (to
each block) by using a finite difference method. Since their discretization has been done then the system of
equations that is required to be solved by the simulator is of the form: where A is a matrix in which
the coefficients (eg. density, permeability) are stored, x is the vector of the unknowns (pressure and
saturation) and b is a vector updated at each time step and includes the total compressibility. Finally, i t has to
be noted that the numerical section is not recommended to be changed by the user without first consulting
the technical support department of CMG.
The numerical control section is constituted by the following sub-sections:
Newton-Raphson method
In the Newton-Raphson method sub-section the required tolerance values for the convergence of the
Newtonian iterations are defined.
Explicit, implicit and IMPES methods
In the explicit, implicit and IMPES methods sub-section which method will be used for the simulation is
defined.
Adaptive implicit method
In the adaptive implicit method sub-section the criterion based on which the change of IMPES blocks to
implicit and the opposite is defined.
Linear solver
In the linear solver sub-section the way the procedure that it is followed in order the system of equations to
be solved is defined.
Parallel processing
In the parallel processing sub-section the parameters that have to be defined in order the system of equations
to be solved in parallel.
Numerical dispersion control
In the numerical dispersion control sub-section the way by which the numerical dispersion will be minimized
is defined.
Geomechanics
In the geomechanics sub-section the required numerical controls for the deformation modeling are defined.
8.1 NEWTON-RAPHSON METHOD
The Newton-Raphson method at each timestep starts making an initial guess which has to be close to the
actual solution and uses the tangent line of the function that the solution is wanted to find a point that is
closer to the actual solution (interception point of the tangent line with the x axis) than the initial guess and if
the iterations continue then the Newton-Raphson method gets closer and closer to the actual solution. If this
solution deviates from the actual solution in an accepted range (convergence tolerance) then the simulator
assumes that the equations have been converged and the provided solution is accepted. On the other hand if
the provided solution is not in the accepted deviation range then the simulator continue iterating until the
provided solution to be accepted or until to reach the maximum allowed number of iterations. In the case that
the Newton-Raphson method iterates until the maximum allowed number of iterations and is not able to
provide an accepted solution then the timestep has to be cut in smaller timesteps and the simulator has to
repeat the calculation procedure.
The Newton-Raphson algorithm which provides the solution of the system of equations for each iteration is
given by the following formula:
, where x() and x(+1) are the unknowns solved for the and +1 iterations of
the Newton-Raphson method performed at each time step, J is the matrix constituted by the derivatives of the
F(x) and the F matrix is defined as follow:
8.1.1 Convergence Tolerance
In order the simulator to proceed to the next time step of the simulation the Newtonian iteration performed
for the previous time step has to be converged. The convergence of the Newtonian iterations is required so as
the simulator to be able to provide the pressure and saturation values that will be used for the solution of the
next time step. Although, the convergence of the iterations is constraint by some specific limits in order the
calculations to be performed as accurately as possible and provide more realistic results. In order the
simulator to provide realistic results takes into consideration the normal variation and the maximum allowed
variation of properties (eg. pressure, saturation, saturation pressure, composition) per timestep.
Moreover, the material balance of the flow equations is not solved to zero and as a result an error occurs and
propagates from time step to time step (Computer Modelling Group Ltd, 2010). The magnitude of this error
is called residual value (R) and it is calculated by the following equation:
, as sink and source term are referred the outflow
and inflow to a block accordingly whereas as accumulation term is referred the mass that remains in the
block. The simulator in order to limit this error propagation sets limiting constraints for the maximum
allowed changes in pressure and saturations during the time step or the maximum allowed residual values by
the solution of the material balance of the flow equations.
8.1.2 Relaxation Of The Accumulation Term
For large simulation problems with many blocks direct methods are difficult to be applied so as the
accumulation term to be calculated and as a result an iterative method namely relaxation should be used
(Peaceman, 1977). The relaxation option is applied to the accumulation term of the mass balance.
Relaxation is nothing more than the convergence tolerance for this term. The simulator for the relaxation of
the accumulation term offers the following two options:
Over relaxation
This keyword enables the over relaxation in accumulation term and may be helpful in cases of gas
appearance and disappearance problems.
Under relaxation
The under relaxation of the accumulation term prevents the disappearance of gas from a block for three
Newtonian iterations that means that if the simulator predicts the disappearance of gas for more than three
iterations only then the gas will be disappeared. This may lead to faster convergence in cases that the gas is
being liberated in the reservoir but it can cause difficulties in cases where free gas is present and the reservoir
pressure is rising (Computer Modelling Group Ltd, 2010).
8.2 EXPLICIT, IMPLICIT AND IMPES METHODS
Explicit and implicit methods are different ways used in order to discretize the partial differential equations
in time. The calculation of the value of a property (pressure, saturation) for a specific block, using the
explicit method depends on the value of the property of this block and its neighboring blocks calculated by
the previous time step. On the other hand the calculation of a property value for a specific block, using the
implicit method depends on the property value of this block calculated by the previous time step and the
values of this property of the neighboring blocks at this time step. Moreover, except the implicit and explicit
methods for the calculation of the pressures and saturation a widely used method in reservoir simulation is
the IMPES (implicit in pressure-explicit in saturation) method. By using the IMPES method the pressures are
calculated by using the saturations of the old time step and then these pressure values are used for the explicit
calculation of the saturations (Fanchi, 2006).
Due to the fact that implicit methods calculate the properties at the current time step are more stable (do not
maximize the errors occurred during the simulation from time step to time step) than explicit methods
(Research Gate, 2017). As a result, implicit methods can use greater time steps than explicit methods.
Although, implicit methods require much more computer resources and time than the explicit methods. On
the other hand, the IMPES method requires smaller memory and time than the implicit method but small
time steps have to be used (Computer Modelling Group Ltd, 2010).
8.3 ADAPTIVE IMPLICIT METHOD
The two methods that are commonly used in reservoir simulation for the treatment of the blocks are the fully
implicit and the IMPES methods. In a field the majority of the blocks can be solved using the IMPES method
whereas only a small amount of blocks is required to be solved using the implicit method.
The blocks that need to be treated implicitly are blocks where sharp changes in pressures and saturation
values occur. This situation is apparent at the well's blocks, the blocks near the wells and in cases of flooding
fronts propagating inside the reservoir. As a result during a simulation some blocks which were initially set
to be solved using the IMPES method may have to be solved using the implicit method for some time period
and the opposite. The method that allows this switching between the two solution methods during the
simulation is the adaptive implicit method (AIM). By using the adaptive implicit method the simulator uses
some criteria so as to check if blocks have to switch from IMPES method to implicit (forward switching) or
from implicit to IMPES (backward switching) (Wig, Codln, Nghlem, & Group, 1989). The available criteria
of the simulator based on which will be determined if a block will switch from IMPES to implicit and vice
versa are the following:
Stability switching criterion
By the stability switching criterion the IMPES blocks are checked if they have to change to implicit mode by
using the stability criterion. The stability criterion defines if the used method provides the required stability.
The stability check can be done either for the IMPES neighbors of implicit blocks or all the blocks of the
field which is more reliable but increases a lot the CPU time in large simulation problems with many blocks.
Threshold switch criterion
When the threshold switch criterion is used the switch is based on pressure and saturation or polymer
concentration. When this method is used a fraction for pressure and saturation or polymer concentration has
to be introduced. This fraction is by turn multiplied with the value entered in normal variation per time step
and their product determines the value of the threshold criterion. The switching of blocks to implicit mode is
done when the calculated value by the simulator changes over the threshold criterion value.
Stability switching criterion and threshold switch criterion
In the case that the stability switching criterion and the threshold switch criterion the simulator is used the
simulator checks the IMPES blocks which are in contact with implicit blocks for stability switching and at
the same time checks all the blocks of the field for threshold switching. When this option is selected then
also the threshold pressure and threshold saturation or polymer concentration fractions have to be defined.
Stability switching criterion using backward switching
By using this option the switch of a block from implicit to IMPES is based on the stability criterion. When
this option is selected then the frequency of checking for backward switching has also to be introduced into
the simulator. This check is done periodically and if a block is found that can provide th e required stability
by using IMPES method then the block switches from implicit to IMPES.
8.6 LINEAR SOLVER
The linear solver is used during each Newtonian step in order to solve the partial differential equations which
describe mass balances of individual components (Petrowiki, 2017f).
8.6.1 Matrix ordering
In order the simulation to be performed a system constituted by all the equations for all the field blocks has
to be solved. This system of equations is stored in a matrix and each cell of the matrix represents a block of
the field. In that matrix many elements are equal to zero (sparce matrix) indicating that there is no connection
between these blocks or in other words that they are not neighboring blocks. These zero elements can be
eliminated from the matrix in order to minimize the store requirements (Golub & Van Loan, 1996). This
elimination of the zero elements can be done by re-ordering the cells of the matrix trying to concentrate the
non zero elements close to each other. This re-ordering of the matrix is done by using specific techniques
also known as ordering algorithms.
The ordering techniques provided by the simulator are the red black ordering, the natural ordering and the
reverse Cuthill-McKee ordering. By using the red-black ordering method the even cells of the matrix are the
red points whereas the odd cells of the matrix are the black points (Young, 2014). This elimination method is
directed towards the minimization of the black points (University of Utah, 2017). By using the natural
ordering the cells of the matrix are ordered in increasing row (MathWorks, 2017b). The reverse Cuthill-
McKee ordering technique is used when the matrix presents a symmetric pattern of the zero elements
(Cuthill & McKee, 1969). The reverse Cuthill-McKee ordering technique aims to the creation of a matrix
with the smallest possible bandwidth (concentrate all the non zero elements close to the diagonal)
(Mathworks, 2017)).
8.6.2 Matrix Preconditioning
The matrix preconditioning is called the procedure during which the previously re-ordered matrix is
decomposed to a product of matrices which is equal to the original matrix. This preconditioning of the matrix
is done in order to result in a system that provides the same solution as the original system but it can be
solved much more easily by the iterative solver.
This decomposition can be done by using an orthogonal matrix (Q) and an upper triangular matrix (R) whose
their product is equal to the initial matrix A ( ). As a result, after this decomposition the original
system of equations that has to be solved is transformed as following:
which can be solved easier as the R is a triangular matrix.
Moreover, GEM simulator also offers and an alternative way to decompose the matrix A by using incomplete
LU factorization. In this case the matrix A is decomposed to a lower and an upper triangular matrices by
using incomplete Gaussian elimination. Gaussian elimination is an algorithm by which the system of
equations is reduced by performing row and column operations in the matrix so as to create as many zero
entries as possible and for computational reasons the elimination is stopped before the matrix become totally
reduced (Golub & Van Loan, 1996). By using the complete LU factorization the original system of
equations is transformed as follow: . By setting the product U and x equal to y the problem can
be solved in two steps, firstly by calculating and then by calculating . Although, this cannot
be done in the case that incomplete elimination is done as the L and U product of the triangular matrices is
not exactly equal to the matrix A ( ) and as a result the system cannot be solved exactly
but approximately. Therefore, the product of the lower and upper triangular matrices is set equal to a matrix
M( ).
8.6.3 Generalized Minimal Residual Method (GMRES)
As soon as the decomposition has take place the system is solved by using the generalized minimal residual
method (GMRES). The generalized minimal residual method is an iterative solver which makes a first
estimate for the solution of the system (a first estimate for the unknown vector x) and by using iterations it is
improving this estimation approximating the exact solution of the system. In order the simulator to assume
that the system is solved the provided solution is required to minimize the residual values (r) of the system
which are defined as follow: (Chen, 2015).
8.7 PARALLEL PROCESSING
In cases of large simulation problems in order to avoid extremely high computational times, clusters of
computers are used so as to solve in parallel the simulation problem. When the simulation is to be performed
in parallel then the field has to be divided in different domains (classes) that are not communicate directly so
as the simulator can solve them simultaneously. Moreover, the number of the computer processors (threads)
that will be used has to be defined.
8.8 NUMERICAL DISPERSION CONTROL
Numerical dispersion depends on the discretization in time and space of the partial differential equations and
it is caused by the fact that numerical methods do not provide an exact solution but they provide an
approximation of this solution (Fanchi, 2006). In general, the coarser is the grid (big disretization in space)
and the bigger is the timestep (eg. using the implicit method) the more pronounced is the numerical
dispersion as the occurred errors from the solution of the equations are greater. These errors propagate during
the simulation time inside the reservoir and look like physical dispersion, as a result the simulator cannot
provide an accurate description of a front movement (Fanchi, 2006). In order the GEM simulator to
minimize the numerical dispersion offers the option to use a flux limmiter which imposes a limit to the
spatial derivatives in order to be maintained in a realistic range of values.
The two flux limmiters that the simulator offers are the Fromm limiter and the minimally extended scheme
(MES) limiter have to be selected. The Fromm limiter dispersion control is preferred to be used in cases
where the grid is in general uniform with almost equal blocks to any direction whereas the minimally
extended scheme limiter is preferred in cases of non uniform grids (Sammon, Masanori, & Jialing, 2001).
Finally, both of methods increase the accuracy of the recovery predictions and the movement of the fluid
fronts.
8.9 GEOMECHANICS
For the solution of the geomechanic deformations the simulator creates a matrix with the force balance
equations and solves it with the same manner as the flow equations. In the case that deformations are
modeled a pressure tolerance multiplier has to be introduced so as to be multiplied with the pressure
convergence tolerance of the Newtonian iterations. This multiplier has to be applied as deformations may
cause large pressure changes and cause convergence problems. Moreover, an allowed tolerance regarding the
displacement of the block nodes due to deformation has to be determined. Finally, if the calculated node
displacement is less than the tolerance value then the equations is thought to have converge.