Numerical model theory

Background

There are two main types of models used to simulate the flow of fluid through porous media: analytical models and numerical models. Both types have their advantages and disadvantages.

Analytical models

Analytical models have relatively short calculation times, so you can history match with these models quickly. More than that, you can even match the model automatically using automatic parameter estimation (APE).

The main disadvantage for these models is that they can only simulate single-phase flow – either liquid or gas, but not gas and liquid flowing together. Additionally, analytical models do not fully account for changing fluid properties with pressure. In case of liquid flow, fluid properties are assumed to be constant, and in case of gas flow, the change of fluid properties with pressure is accounted for by using pseudo-pressure and pseudo-time.

Unfortunately, pseudo-time is not an exact transformation. As a result, changing gas properties are only accounted for to a certain extent. Therefore, in cases where pressure varies significantly across the reservoir (for example, pays with low permeability), long-term forecasts for gas analytical models may be inaccurate. There is no known way to fully account for the change in fluid properties when using analytical models.

Numerical models

The numerical models have been significantly improved in Harmony Enterprise. One of the key benefits is the improvement in computation time, so that the model runs almost as fast as an analytical model. These modifications include:

  • using the pseudo-pressure formulation
  • optimized "gridding" — by using the pseudo-pressure formulation, you can use larger grid cells without affecting the accuracy of calculation results; using larger grid cells also makes calculations run faster.

With numerical models, you can model multiphase flow, and account for changing properties of each phase, and the interaction between phases.

The main disadvantage for these models is their long computation time. With the numerical model, the reservoir is divided into a number of cells, and then the flow is modeled simultaneously in all cells. To be able to assume constant pressure and constant fluid properties within each cell, the size of each cell needs to be relatively small. Therefore, the number of cells required for an accurate solution becomes large, which results in a long computation time. Note that decreasing the number of cells to reduce computation time may significantly affect calculation results.

Performing a manual history match for numerical models is time consuming. In addition, an automatic history match for numerical models is complicated, and thus not widely available in commercial software.

General formulation for the single-phase model

To calculate pressure across the reservoir, we divide it into a number of grid cells. At each timestep (n+1), we formulate a material balance for each cell.

To define components of mass balance for cell i at timestep n+1:

Mass in place (MIP) for cell i at timestep n+1 can be calculated as:

Equation 1

where:

Vbi = bulk volume for cell i

pin+1 = pressure for cell i at timestep n+1

φ = porosity

ρ = density

Mass flow from cell i to the adjacent cell j during timestep n+1 can be calculated as:

Equation 2

(We used Darcy’s Law to calculate flow rate.)

In equation 2:

= flow rate from cell i to cell j at timestep n+1

Ai,j = cross-section area between cells i and j

Li,j = distance between the centers of cells i and j

μ = viscosity

k = permeability at the initial pressure

km(p) = permeability multiplier (permeability at a certain pressure is calculated as ).

Mass produced at the well during timestep n+1 (if well penetrates cell i) can be calculated as:

Equation 3

where:

WI = wellbore index (as per Peaceman 1978 or Babu and Oden, 1989)

= pressure at the wellbore at timestep n+1

Material balance for cell i at timestep n can be written as:

Equation 4

Using equations 1, 2, and 3, and bringing all the terms to the left, we can rewrite equation 4 as:

Equation 5

Assuming that the pressure distribution has been calculated for timesteps 1, 2, ... n, to calculate pressure distribution at timestep n+1, we formulate material balance (equation 5) for each cell, and solve all these equations simultaneously. Unknowns of the system are for each cell i; therefore, the number of equations is the same as the number of unknowns, and the system can be solved.

Well constraint considerations

For a well producing at a specified sandface pressure, the described system of equations can be used directly (because is known). However, for a well that produces at a specified surface q, an extra equation is required.

Equation 6

where:

i = summation index for all cells penetrated by the well

ρsc = fluid density at standard conditions

Considerations when estimating fluid properties

In the above formulation, we use rock and fluid properties estimated at certain pressures. The term is used in (equation 2) and in (equation 3).

Some logical questions are:

  • At what pressure should we estimate this term?
  • Should we use pressure in one of the cells, or some kind of average?

Fluid properties may significantly vary with pressure (especially if the modeled fluid is gas); therefore, selecting the correct way of estimating properties becomes important.

In classical numerical simulation (including numerical models in Harmony), properties are estimated at the pressure of the upstream cell (i.e., the cell with the higher pressure). For such a simulation to be accurate, you should ensure that properties in the adjacent cells are close. To achieve this, grid cells have to be small enough. Unfortunately, having smaller grid cells results in longer computation times.

With the numerical model, we use the pseudo-pressure formulation to estimate fluid properties and calculate mass transfers.

Numerical model pseudo-pressure definition

The numerical model pseudo-pressure is defined as:

Equation 7

Note:    For gas, the definition of numerical model pseudo-pressure is very similar to a traditional pseudo-pressure. To see this, express density as (z = gas compressibility factor, T = temperature, R = gas constant, M = molar mass). This results in the following equation:

Equation 8

Therefore, the numerical model pseudo-pressure and traditional pseudo-pressure are different by a constant factor.

Numerical model pseudo-pressure formulation

As was mentioned in considerations when estimating fluid properties, the calculation of mass flow from one cell to another given in equation 2 has a deficiency: it does not account for a variation of fluid properties with pressure. The numerical model formulation modifies equation 2 to get a relationship between the mass flow rate and pressure drop when properties are changing with pressure.

At each cross-section at any given distance (x), the pressure gradient across the cross-section can be expressed using the differential form of Darcy’s Law:

Equation 9

where is the mass flow rate across the cross-section. This can be re-arranged as:

Equation 10

Integrating both parts from 0 to L with respect to x, and using the left-side of the equation is a constant with respect to x, we get the following equation:

Equation 11

Note:    We used the definition given in equation 7 for the last transformation.

To summarize, when fluid properties are changing with pressure, the mass flow rate through the cross-section can be calculated as:

Equation 12

Therefore, equation 2 can be modified to account for changing properties. The mass flow from cell i to cell j during timestep n+1 is calculated as:

Equation 13

Equation 3 can be modified in a similar fashion. Mass produced at the well during timestep n+1 is calculated as:

As a result, the material balance equation for each cell (equation 5) is modified to:

Equation 14

Assuming that the pressure distribution has been calculated for timesteps 1, 2, ... n, to calculate pressure distribution at timestep n+1, we formulate a modified material balance (equation 14) for each cell, and solve all these equations simultaneously. Unknowns of the system are for each cell i; therefore, the number of equations is the same as the number of unknowns, and the system can be solved.

Note:    Well constraints are treated in a similar way to general numerical model formulations.

Modifications for the multiphase model

In the case where there are two or three phases flowing together, the multiphase solution follows similar logic as single-phase flow, with a few modifications.

  • The mass balance equation (equation 4) has to be modified to account for the fact that fluid can change from one phase to another (e.g., oil can evaporate and become gas, or gas can condense to liquid).
  • In the single-phase case, the mass balance equation is written for each cell; therefore, we get n equations (where n is the number of cells). In the multiphase case, the mass balance equation is written for each cell for each phase; therefore, we get n•"number of phases" as equations.
  • In the single-phase case, the unknowns of the system are pressures in each cell; therefore, we get n unknowns. In the multiphase case, the unknowns of the system are pressures in each cell (n), and saturations for each cell; therefore, we get total n•"number of phases" as unknowns. (We use the fact that the sum of all saturations is equal to 1 for any of the grid cells.)
  • When calculating a mass transfer from one cell to another (equation 2) and a mass transfer from the cell to the well (equation 3), we have to account for relative permeability.
  • For the single-phase case, we re-write the mass transfer term in terms of pseudo-pressure (equation 13). In the multiphase case, we can do the same thing, but there will be different pseudo-pressure functions for different phases. (The pseudo-pressure formulation for each phase follows the definition given in equation 7, but ρ, μ, and k are different for different phases.)
  • For the single-phase case, the system of equations can be solved in terms of pseudo-pressures (the unknowns are pseudo-pressures in each cell). For the multiphase case, equations have to be solved in terms of pressure (the unknowns are pressures and saturations in each of the cells). The pseudo-pressure formulation is only used to calculate mass transfer between cells.

Advantages of using the pseudo-pressure formulation

Pseudo-pressure formulation makes simulation significantly faster than classical numerical simulation modeling for the following reasons:

1. Smaller number of grid cells: In classical numerical simulation, grid cells have to be small enough to consider constant fluid properties within each cell. Pseudo-pressure formulation accounts for a variation of fluid properties between the centers of two adjacent cells; therefore, it is possible to have larger grid cells. By having a smaller number of grid cells, calculation speed increases.

2. Faster solution for a non-linear system: While performing numerical modeling, equation 5 is solved at each timestep. This system of equations is non-linear; therefore, it is solved iteratively, using the Newton-Raphson method. This method involves calculating derivatives of each matrix element with respect to each unknown. Calculating these derivatives is faster for pseudo-pressure formulation (equation 14), because we have to deal with one fluid property function () as opposed to a combination of three functions () for traditional formulation. More than that, due to the integral nature of , its derivative is calculated easily.