Finite Element Simulation of Multiphase Flow in Oil Reservoirs – Comsol Multiphysics as Fast Prototyping Tool in Reservoir Simulation

Reservoir simulation is a powerful tool to mimic the formation behaviour during primary production and later on for planning enhanced oil recovery (EOR) pattern. However, all available commercial and developed scientific/academic software for this purpose is based on either finite difference method (FDM) or finite volume method (FVM). Recently finite element method started to gain more attention in the scientific and commercial practices due to its robust results and the ability to deal with complex boundaries. COMSOL Multiphysics is a finite element method (FEM)-based software, having very special features, which are different from standard reservoir engineering software packages like Eclipse or CMG, which are black box-type software. The most important feature of the COMSOL is that user can see equation and modify it – customize for specific conditions and objectives, as well as couple different physics together and apply different solvers, which are under user’s disposal. In this paper, short background of FEM will be illustrated and then the mathematical models of two-phase immiscible flow of water and heavy oil will be reviewed and simulated using COMSOL Multiphysics on the famous inverted five-spot model. The comparison between the results of Comsol Multiphysics and Eclipse shows good agreement. This study is the first step in applying Comsol Multiphysics to reservoir simulation. Further steps will involve simulating thermal enhanced oil recovery using steam flooding technique and coupling Comsol Multiphysics with CMG software package to enhance simulation inputs and outputs.


Introduction
Reservoir simulation combines mathematics, physics, reservoir engineering, and computational science to develop a tool for predicting oil and gas reservoirs performance at different operating patterns [1]. To describe the fluid flow inside an oil or gas reservoir a set of partial differential equations (PDE's) must be solved with the use of consistent set of initial and boundary conditions. Different numerical methods are applied to solve these equations, but the most applicable methods in commercial and scientific software are finite difference method (FDM) followed by finite volume method (FVM). The application of finite element method (FEM) is still very limited in this branch of science.
In this paper, short background of FEM will be illustrated and then the mathematical models of two-phase immiscible flow of water and heavy oil will be reviewed and simulated using COMSOL Multiphysics on the famous inverted five-spot model. The comparison between the results of Comsol Multiphysics and Eclipse shows good agreement.

Finite element method
While FDM and FVM can be considered as easier methods in terms of programming language, FEM is a little bit difficult to program due to its speciality in equation discretization. There are different procedures to solve PDE's in FEM, but the most famous method is the Galerkin method. To solve PDE using Galerkin approach the following steps should be followed [2]:  multiply the original equation with a test function;  integrate the equation and apply boundary conditions to produce weak form of the original equation;  write the finite element solution as linear sum of a set of basis function (in the Galerkin method, test functions and basis functions are identical);  apply the finite element solution into the weak form equation;  solve the set of algebraic equations produced at the previous step;  perform the error analysis.

Computational methods 3.1. Two-phase immiscible flow model
Two-phase flow in any porous medium can be described by continuity and momentum equations for each phase [2,4,5].

Mass conservation (continuity equation)
Let the porous medium fills a domain . Conservation equation for each phase α is written as: Darcy's law (momentum equation). It is also defined for each phase : Define the mobility for each phase too: Fractional flow for each phase In the case of two phase flow for oil and water, the total velocity is: To close the equations' system, two customary equations for saturation and capillary pressure are introduced as follows:

Fractional flow formulation of immiscible two-phase model
Equation (1) can be simplified by the following assumptions [3]:  temperature is constant in the domain;  there are two phases: water (w) and oil (o);  there are two components: water (w), only in the water phase, and oil (o), only in the oil phase;  the rock (porous matrix) and the fluids are incompressible;  the solid matrix is not poroelastic; this means that the available pore space (porosity) is constant.
Applying equation (2) into (1) results in: When applying the above assumptions and adding oil and water continuity equations, equation (1) can be written for the two-phase system as: Subtracting the continuity equation for the two phases produces: This will lead to the total velocity and phase velocities: Neglecting capillary forces leads to our model that may be applied in Comsol Multiphysics:

Comsol implementation of immiscible two-phase model
Simulation steps using COMSOL Multiphysics ® software are different from the traditional reservoir simulation software packages owing to the speciality in discretization procedure used in FEM that is significant difference compared to FDM and FVM. In general, the simulation procedure is as flows:  setting model environment;  creating geometric objects;  specifying materials properties;  defining physics boundary conditions;  creating the mesh;  running simulation;  post-processing of the results. Mathematical module is used due to its merits in stabilizing the solution, and we use general form and coefficient form PDE to solve pressure (p) and saturation (s) equations, respectively. The numerical study is applied on the famous SPE case study quarter five-spot inverted model shown in Fig. (1) and meshed in Fig. (2).

Pressure equation
General form time dependent PDE from mathematical module in Comsol is applied to pressure equation (4) with p as an independent variable. The components of the total velocity are the components of conservative flux vector. There is no source term, equation (11).  

Initial conditions
Initial pressure in the whole reservoir

Simulation findings
Selected data from [1] are chosen to be applied in this study. Figs. (3) to (6) illustrate the saturation profile in the domain at different periods through the life of the reservoir where the oil saturation decreases significantly toward the production well to reach minimum at 2350 days. Fig. 7 shows the saturations (So, Sw) as the function of time. Fractional flow (fo, fw) is shown in Fig. 8. The comparison between the daily production from Comsol and Eclipse-100 shows almost identical results for the first period before water breakthrough, Fig. 9.

Heavy oil model
Heavy oil model can be derived from classical continuity equation (1) and Darcy's law [4] as follows: Saturation equation: Equations (14) and (15) can be applied in Comsol Multiphysics in the same way as in twophase flow model.

Conclusions and future work
Principle of FEM has been reviewed and mathematical models of two-phase immiscible flow model and heavy oil model have been formulated to be easily applied in Comsol Multiphysics to illustrate saturation profile in the reservoir. The results show good agreement with commercial software Eclipse 100.
This study is the first step in applying Comsol Multiphysics ® to reservoir simulation. Further steps will involve simulating thermal enhanced oil recovery using hot water and steam flooding techniques and coupling Comsol Multiphysics with CMG software package to enhance simulation inputs and outputs.