Project topics and materials
Get Complete Project Material File(s) Now! »
Governing equations
In this chapter, we review some of the main concepts and formulate the driving equations of the flow and transport in porous media, to provide an understanding of these equations in single-phase and multiphase flows in the reservoir.
Porosity and permeability are the main properties of a porous medium. The ratio of the volume of void spaces to the total volume is defined as the porosity, . With this definition, varies between zero and one. The absolute permeability is an intrinsic property of the porous medium that measures its ability to allow a single fluid to flow through it. The permeability is generally correlated to the porosity. The general form of permeability is a tensor, meaning that the permeability in each direction depends on the permeabilities in other directions. In 2D, the permeability tensor is a matrix of the form [ ], in which the off-diagonal terms represent this dependence. When the off-diagonal terms are all zero and diagonal terms are equal, the permeability becomes a scalar property that does not depend on the direction and is called isotropic. In other cases, the permeability is anisotropic.
Natural porous media are generally heterogeneous at all length scales from the pore scale to the reservoir scale in the order of several kilometres. Full-field numerical simulation of flow problems at the pore-scale is not computationally feasible. Moreover, obtaining data at the microscale is too challenging and too precise for most practical situations. A common scale is the representative elementary volume (REV). REV, usually in the order of several hundred to thousands of pores, is the smallest volume over which a property average does not change with a small change in its size.
Single-phase flow
One of the most fundamental laws in reservoir simulations is Darcy’s law, by the French engineer Henry Darcy (Darcy, 1856). In one dimensional single-phase flow, Darcy’s law can be represented as:
– negligible gravitational effects,
– constant fluid viscosity,
– isothermal reservoir condition,
– incompressible rock and fluid.
With these assumptions, and using Darcy’s law in Eq. 2.4, the single-phase flow equation or the “pressure” equation simplifies to: ∇. = −∇.( ) = . Eq. 2.5
To close the model, boundary conditions are imposed: = Eq. 2.6 { . = 0 which correspond to Dirichlet and homogeneous Neumann boundary conditions on the boundaries of the domain.
Multiphase flow
In most applications of reservoir simulations, two or more phases are generally present, e.g., water or gas displacing oil or CO2 displacing the brine. In multiphase flows, each phase occupies a certain fraction of the pore space and new properties are introduced in the flow model (Marle, 1981). The pore volume fraction occupied by each phase, , is defined as its saturation, . We assume that all the fluid phases fill the void space, i.e. ∑ =1. Eq. 2.7
Another property in a multiphase flow model is the relative permeability which measures how the presence of one phase hinders the ability of the flow of another phase. In the presence of other phases, each phase experiences a reduced permeability called the effective permeability, defined as: = Eq. 2.8
where the relative permeability, , generally is a function of phase saturation, , but their dependence cannot de defined exactly. It is a popular choice to use some empirical laws to relate the saturation and the relative permeability.
At the pore scale, due to the interfacial tension between two immiscible fluids, the pressure in these two phases is generally different. This pressure difference across the fluids interface is called the capillary pressure, : = − Eq. 2.9 where and denote the pressure in the and phases, usually being the nonwetting and wetting phases, respectively.
In this work, we consider an immiscible two-phase flow model, for example, water and oil in a waterflood problem in a porous medium. We neglect capillary effects, along with the other assumptions mentioned in single-phase flow model. We use the power law, one of the most common empirical relations to define the relative permeability that reads in terms of reduced water saturation, ∗ = 1− − (Brooks and Corey, 1964): = ( ∗ ) Eq. 2.10 = (1− ∗) and are irreducible water saturation and residual oil saturation, respectively. , or the maximum relative permeability of each phase , and , the exponent of each phase, should be generally found by fitting to the measured data. In the numerical examples throughout this manuscript, we set both exponents, and to 2, a popular choice in the industry, and the maximum relative permeabilities to 1, = 1. The corresponding relative permeability curves are shown in Figure 2.1.
Like in single-phase flows, when modelling the behaviour of multiple fluids flowing in porous media, mass conservation and Darcy’s law are used to derive the system of equations. Considering the fluids to be incompressible, the conservation equation for each phase is: +∇.( )= , = , Eq. 2.11 where is the phase velocity, and is the flow rate of the source and sink in each phase. We write the two mass conservation equations into a more manageable system of equations, composed of a pressure or flow equation and a saturation or transport equation.
Saturation equation
Writing the conservation equation of water in terms of total velocity and water fractional flow leads to the ‘saturation equation’: + ∇. ( ( ) ) = Eq. 2.17 to be solved with the initial and boundary conditions: { = 0 , Eq. 2.18 . <0 with 0 the initial saturation in the domain and the Dirichlet condition on the inlet boundary of the domain. The fractional flow function = ⁄ measures the water fraction of the total flow.
The system of equations described through Eq. 2.16 to Eq. 2.18 is nonlinearly coupled. The mobility term in the pressure equation is saturation dependent and the flux function term in the saturation equation is pressure dependent. There are several approaches to numerically solve this coupled system of equations. In a fully implicit approach, all the variables are solved simultaneously. This approach is unconditionally stable, but its implementation and solution are computationally expensive, especially for large problems like full-field reservoir simulations. Another class of methods, sequential approaches, aims at solving the pressure and saturation equations separately and sequentially. The main advantage of these approaches lies in the reduction of the size of the linear systems to be solved, and therefore superior computational efficiency. Another important advantage is that it allows mixing different discretization methods in the same system, which can be very beneficial as the flow and transport equations have different mathematical characters, the elliptic character of the flow versus the hyperbolic character of the fluid transport (Bell et al., 1986). Sequential approaches have different forms, like the so-called IMPES (implicit pressure – explicit saturation) and the IMPIMS (implicit pressure – implicit saturation) schemes (Coats, 2000; Fagin and Stewart, 1966; Stone and Garder, 1961). Some of the main disadvantages of sequential approach are the additional memory requirements and the splitting error related to the decoupling of flow and transport equations. One proposed remedy to decoupling error is to add an extra step after solving the transport equation to revisit the flow problem (Hajibeygi and Tchelepi, 2014; Lee et al., 2008). Adaptive implicit method (AIM) is another method which tries to combine advantages of sequential and fully implicit methods by implementing each of them adaptively in the domain (Cao and Aziz, 2002; Spillette et al., 1973; Thomas and Thurnau, 1982). For example, a fully implicit method is most suitable for local areas with large variations in the saturation, while computationally efficient IMPES scheme can be applied elsewhere. The main challenge in AIM is the difficulty of finding a suitable switching criterion (Marcondes et al., 2009).
Literature review
This chapter presents a review of literature, starting with different approaches to describe single-phase and multiphase flows in porous media, from stochastic to numerical methods and fast diagnostic tools. Next, a review and analysis of upscaling methods, in single-phase and multiphase flow modelling, known as one of the classical approaches to accelerate flow simulations, is detailed followed by a brief review of some alternative techniques. In the last section, we review some of the methods used to locate the saturation discontinuity, one of the most important local features of an immiscible multiphase flow model.
Stochastic approaches
In the stochastic approach, input parameters such as permeability and porosity are considered as being random functions of the position. Geostatistical tools provide methods allowing to generate as many random maps honouring imposed statistical properties as desired (Guadagnini et al., 2018). In this thesis, FFTMA (Le Ravalec et al., 2000) will be used to generate spatially correlated Gaussian random maps.
The idea of stochastic methods is to be able to estimate the ensemble average of the properties of interest, including dynamic data such as pressure, saturation, and cumulated oil production that are themselves random. Indeed, the randomness of the input parameters propagates through the solution of the partial differential equations that govern the flow. In practice, this ensemble average could also be obtained by averaging many independent simulations: the so-called Monte Carlo approach. Although easy to parallelize, it remains an expensive approach in the case where many simulations need to be performed, especially if each simulation is computationally costly and long. Moreover, the number of realizations that are required is itself a part of the issue.
The core idea of stochastic methods is to look for the output parameters average, such as the pressure 〈 ( , )〉 and the saturation 〈 ( , )〉, and to find if an effective set of equations driving these averages may be found by some methods. The associated variance may be studied by analogous methods, to provide some information about the related uncertainties.
Although conceptually simple, some difficulties may be anticipated. First, the form of the equations may be changed, even in the simplest single-phase context by considering ensemble-averaged Darcy’s law. It may be shown under a quite general hypothesis (Nœtinger and Gautier, 1998) that the linear relation between averaged flux and pressure gradient becomes an integrodifferential operator that only keeps the overall linearity of the original problem, a strong requirement. The case of a tracer advected in a random flow-field leads to the so-called macro-dispersion (Gelhar and Axness, 1983). In two-phase flows, few works followed this approach (Artus et al., 2004; Langlo and Espedal, 1994; Spesivtsev and Teodorovich, 2007; Zhang and Tchelepi, 1999).
Another difficulty is that even if the form of the ensemble-averaged equation is known, the question on how to relate the associated coefficients to the input geostatistical parameters remains. In the context of single-phase problems, considering the ensemble averaging of Darcy and advection/dispersion equations, the main approach is the small perturbation expansion (Dagan, 1989). Considering that the local permeability may be decomposed as: ( ) = 〈 〉 + ( ) Eq. 3.1 With the perturbation value ( ) ≥ 0, the ensemble-average 〈 〉, and the tuning parameter , one may expand the solution of Laplace and tracer equations in powers of . It may be anticipated that the second-order term provides contributions involving averages of 〈 ( ) ( ′)〉 corresponding to permeability correlations. The quantity 〈 ( ) ( )〉 = 2 is the permeability covariance function. It allows expressing the effective permeability up to the second order as: = 〈 〉(1 − 1⁄ 2 ⁄〈 〉2 + ⋯) Eq. 3.2 in which D=1,2,3 is the number of dimensions available to flow. Further developments were performed to get higher-order terms, generalizing the approach to higher variances using field theoretical approaches, Feynman diagrams, and so on (Dagan, 1993; Indelman and Abramovich, 1994; King, 1987; Noetinger, 1994).
In the case of the passive tracer transport equation that reads: Eq. 3.3 + ∇. ( ) = ∇. ( . ∇ ) in which is the concentration, is the local Darcy velocity that depends on the underlying disorder, and is the diffusivity tensor that represents the local molecular diffusion/dispersion tensor. The stochastic approach leads to an effective equation that reads asymptotically: 〈 〉 + ∇. (〈 〉〈 〉) = ∇. ( . ∇ ). Eq. 3.4
The tensor , proportional to 2 ⁄〈 〉2 with the correlation length , may have a non-zero value even if the input D vanishes. This corresponds to the so-called macro-dispersion phenomena (Dagan, 1989; Gelhar and Axness, 1983) which represents the mixing due to the fluctuations of streamlines of local velocity . Once again, methods such as renormalization group techniques were proposed to increase the validity range of the results (Jaekel and Vereecken, 1997). Good agreement was obtained using field data.
For completeness, let us mention another class of stochastic techniques that were developed for modelling flow in fracture networks. The framework is mainly percolation theory and involves different methods. The interested reader can look at the textbooks of Berkowitz and Balberg (1993) and Hunt et al. (2014).
Previously listed methods were developed for problems driven by linear equations, such as Darcy equation with random coefficients, and passive tracer advection/dispersion. The situation is far more complex when two-phase flow, such as water/oil in porous media, is considered. The generalization of stochastic methods to two-phase flows is difficult due to the coupling between the pressure equation and the saturation conservation equation which leads to a coupling between viscous effects and channelling effects. Partial decoupling approximations (Langlo and Espedal, 1994; Zhang and Tchelepi, 1999) can lead to erroneous results (Nœtinger et al., 2005). The technical difficulty is that the existence of the Buckley Leverett discontinuity leads to technicalities for setting up the perturbation theory. A change of variable proposed by King and Dunayevski (1989) allowed to get some results up to the second order in the permeability variance for a limited range of validity (Artus et al., 2004; Noetinger et al., 2004; Spesivtsev and Teodorovich, 2007). The main result is that finding an effective transport equation sharing some similarities with non-linear advection/dispersion equation driving the ensemble-averaged saturation remains to be done.
Numerical methods
Numerical modelling of multiphase flows in porous media aims at solving the governing equations through their discretization. The main challenges regarding the flow equation are the accuracy of the flux that enters the transport equation and dealing with the heterogeneous permeability that can vary by several orders of magnitude over the domain. On the other hand, the main difficulties for numerically solving the hyperbolic problems of conservation laws lie in the presence of discontinuities in the exact solution and the complicated coupling of different origins near these discontinuities. In the following, we review the main classical and most widely used methods as well as some advanced discretization methods. To review their mathematical base idea and the strengths and weaknesses, we consider the elliptic flow equation, described in Eq. 2.16. For a more thorough review and comparison of discretization methods in geosciences, see Di Pietro and Vohralík (2014).
Finite difference and finite volume methods
Perhaps the simplest and oldest approach is the finite difference (FD) method. Finite difference methods work based on replacing the derivatives in the partial differential equation (PDE) with an approximation based on Taylor series. Even though this approach is very easy to implement and can be very efficient, the need for an orthogonal grid structure limits their application in more general and complex domain geometries.
In many engineering processes, including flows in porous media, mass conservation is an important property. In these contexts, finite volume (FV) methods are very popular. The reason is that finite volume methods, unlike finite difference methods, have physical principles behind them and can be derived only based on the conservation of the desired quantity over control volumes, independently from the continuous equation. In FV methods the unknowns are approximated by constant values for each cell or control volume, that represent the volume-averaged values. This another difference with finite difference methods, in which the physical quantities are calculated at specific nodes, gives FV methods more geometric flexibility. Despite these differences in derivations, FV and FD methods are usually used as identical in the literature. A common approach to derive a finite volume representation of the flow equation is to integrate Eq. 2.16 over a grid cell or control volume Ω in the domain Ω, and use the divergence theorem to obtain: ∫ .= ∫ , Eq. 3.5 Ω Ω where is the outward unit normal vector on the boundaries of the grid cell Ω . Different finite volume schemes approximate the flux through interfaces from a set of neighbouring cell-averaged pressure values. The two-point flux-approximation (TPFA) scheme is the simplest form of a finite volume discretization. It works based on approximating the flux across the face, = Ω ∩ Ω , using the values in the two neighbouring cells, Ω and Ω , on the sides of the face : = | | − Eq. 3.6
where | | denotes the area of the face, denotes the mobility defined at the face, and is the distance between the neighbouring cell centres. This scheme is simple to implement and computationally efficient, and that is why it is widely used in industrial software. However, the main constraint in this scheme is its conditional consistency. If the grid does not satisfy some orthogonality condition, called K-orthogonality, which is usually the case in unstructured grids in porous media with complex geology, the solution of the TPFA method will have errors. The severity degree of the error depends on the direction of the interfaces with respect to the direction of the permeability tensor. These artefacts or errors, called “grid orientation effects”, can lead to the convergence to a completely wrong solution if the mentioned orthogonality condition failure is dominant (Aavatsmark, 2007; Faille, 1992; Wu and Parashkevov, 2009).
Table of contents :
Abstract
Acknowledgement
List of Figures
List of Tables
Nomenclature
Thesis outline
Governing equations
Single-phase flow
Multiphase flow
Pressure equation
Saturation equation
Literature review
Stochastic approaches
Numerical methods
Finite difference and finite volume methods
Finite element and discontinuous Galerkin methods
Flow diagnostic tools
Upscaling
Single-phase upscaling
Multiphase upscaling
Alternative approaches
Dual grid techniques
Multiscale methods
Adaptive gridding
Localization of the saturation discontinuity
Front tracking methods
Level set method
Fast marching method
Discretization methods
Finite volume method
Single-phase flow:
Multiphase flow
Numerical errors
Discontinuous Galerkin method
Discretizing the flow equation
Adaptation of FESTUNG for reservoir simulations
Reconstruction of conservative fluxes
Using the underlying permeability values
Numerical errors
Fast front tracking method
Buckley Leverett problem
FFrT formulation
Test case validation
Buckley Leverett problem
Layered medium
Random generated isotropic medium
SPE10 benchmark test
Combining FFrT with adaptive coarsening: A proposed approach for immiscible multiphase flows
Formulation
Complexity analysis
Numerical examples
Numerical example one: Pressure solution on the base coarse grid
Numerical example two: Adaptive coarsening criterion
Numerical example three: Domain decomposition
Numerical example four: Non-uniform coarsening
Numerical example five: Intermediate resolution level
Numerical example six: Favourable displacement
Numerical example seven: Unfavourable displacement
Conclusion
Concluding remarks
Future work
References
