Spontaneous Imbibition Simulation with Rayleigh-Ritz Finite Element Method
Authors: Kaul Sandeep P. , Erwinsyah Putra , Schechter David S. 
Published: December 2007

Proposal
In highly heterogeneous reservoirs, where fractures and vugs are dominant, the role of spontaneous imbibition in displacement of oil is important. Apart from the complex phenomenon of fluid exchange in them, varied reservoir geometries make simulation in these types of reservoirs a challenging task. Except for small or medium sized reservoirs, numerical simulation of these comes under the purview of high speed or multi processor computers. Finite Element (FEM) numerical simulation has an edge over Finite Difference (FDM) schemes as it can depict complex geometries of these reservoirs with accuracy without resorting to fine scale gridding. Many researchers have already presented various approaches of fluid flow using finite elements method but nearly all of them are iterative. This is because the diffusion saturation equation is non-linear partial differential equation and resulting iterative FEM schemes puts a burden on computer memory and renders large-scale simulation unsuitable. This paper attempts to present a new mathematical model that makes the finite element method non-iterative by applying classical Rayleigh-Ritz method. The classical method cannot be applied directly. For this we have to truncate the non-linear terms, decouple and solve analytically the dependent variables from the primary variable - saturation, which reduces the non-linear FEM problem to a simple weighted integral weak form that can be solved with Rayleigh-Ritz method. We compared our numerical model with the analytical solution of this diffusion equation. A FDM numerical model was validated using X-Ray Tomography (CT) experimental data of a single-phase spontaneous imbibition experiment, where two simultaneously varying parameters of weight gain and CT water saturation were used. The results of FEM model were then compared to that of FDM model. A two-phase field size synthetic example was developed and results from a commercial simulator were compared to the FEM model.

Introduction
The search for hydrocarbons has entered into an era where conventional reservoirs are all depleted or at the end of their economical life. There arises a need now for exploiting highly heterogeneous and unconventional reservoirs, which have huge economic potential, but are too complex when, it comes to understanding the physics of fluid flow. Naturally fractured reservoirs belong to such category and the major process, which determines the amount of hydrocarbons that can be extracted from them, is spontaneous imbibition.
At the heart of spontaneous imbibition process is suction generated by capillary forces. But fluid flow in porous media, determined primarily by capillary forces, is relatively difficult phenomenon to quantify, the details of which can be had from reference 1, and much of the research effort directed in this direction appears in the works of Handy2, Garg et al.3, Babadagli and Ershaghi4, Li and Horne5, Akin and Kovscek6, Reis and Cil7, Zhou et al.8, to name a few. Spontaneous imbibition process affects the fluid flow in fractured porous media as well, and has also received much attention by other researchers. In case of naturally fractured reservoirs, ever since the model proposed by Barenblatt et al9 and extended by Warren and Root10, various researchers such Mattax and Kyte11, deSwaan12, Kazemi et al13, Civan14 to name a few, have helped to understand recovery from such reservoirs.
Finite element method and its more popular part, the Galerkin method, is an iterative numerical method. As compared to it, the finite difference method is a direct method because the coefficients of the finite difference equation are all linear. As a result, the matrix that is used for inversion. single most time consuming operation in any numerical simulation. is strictly in banded form. The Galerkin method is an iterative method where the coefficients of the weak form have to be ascertained prior to the actual inversion of the matrix. This severely taxes the computer memory. In comparison to the Galerkin FEM, the Rayleigh-Ritz FEM is simple as its formulation incorporates the boundary condition along with a form of differential equation. Thus the solution is not iterative. Also, in general, higher-order interpolation functions are to be used in Galerkin FEM since its formulation does not include any specified boundary conditions. The main advantage of this is that we can incorporate all the geological details into an efficient mathematical model which is less mathematically intense than existing finite element methods. This is achieved without having to resort to large-scale refinements, to capture reservoir heterogeneities and details, as has to be done by finite difference method.