key: cord-0046007-6ygv3ap6
authors: de Paula, Filipe Fernandes; Quinelato, Thiago; Igreja, Iury; Chapiro, Grigori
title: A Numerical Algorithm to Solve the Two-Phase Flow in Porous Media Including Foam Displacement
date: 2020-05-25
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50436-6_2
sha: c7f7ae6f5de060df2a2ba113c71810cadc7fffed
doc_id: 46007
cord_uid: 6ygv3ap6

This work is dedicated to simulating the Enhanced Oil Recovery (EOR) process of foam injection in a fully saturated reservoir. The presence of foam in the gas-water mixture acts in controlling the mobility of the gas phase, contributing to reduce the effects of fingering and gravity override. A fractional flow formulation based on global pressure is used, resulting in a system of Partial Differential Equations (PDEs) that describe two coupled problems of distinct kinds: elliptic and hyperbolic. The numerical methodology is based on splitting the system of equations into two sub-systems that group equations of the same kind and on applying a hybrid finite element method to solve the elliptic problem and a high-order finite volume method to solve the hyperbolic equations. Numerical results show good efficiency of the algorithm, as well as the remarkable ability of the foam to increase reservoir sweep efficiency by reducing gravity override and fingering effects.

The enhanced oil recovery by injection of gas is a technique that is used since the 1930's [14] . The sweep efficiency of gas, however, can be affected by gravity (by the gravity override phenomenon, that occurs when the injected gas accumulates in the upper layers of the reservoir) and by the development of preferential paths (viscous fingering), due to gas lower density and viscosity. These obstacles can be surpassed by the creation of foam, that can be defined as the agglomeration of gas bubbles separated by thin liquid films (lamellae), since foam apparent viscosity is much higher than the viscosity of gas [11, 12, 23] . The usage of foam in oil recovery is mainly motivated by the reduction of the gas phase mobility [15] .

Population balance models can be used to simulate foam creation, destruction and flow through porous media. In this approach, it is common to define the foam texture (n f ), a quantity that represents the number of bubbles (or lamellae) per unit volume. Mechanisms of foam creation and coalescence play an important role in the model [15, 17] . The major hypothesis adopted for the lamella-coalescence mechanism is that bubbles collapse near the limiting water saturation (S * w ) or, equivalently, the limiting capillary pressure. In this context, a foam model based on the well-known steady-state behavior of foam in porous media was proposed by Ashoori et al. in [2] . It considers a large, nearly constant, reduction in gas mobility at high water saturation and an abrupt weakening or collapse of foam at a limiting water saturation. Foam texture in local equilibrium (n LE D ), where bubble generation and destruction reach a local equilibrium state, depends only on the water saturation (S w ):

with constant A. The dynamic foam net generation is given by a first-order approach, introduced in [24] and later related to the local-equilibrium bubble texture in [2] , with a time constant 1/K c , as follows:

where r g and r c are the foam generation and coalescence rates, respectively, n max is the maximum foam texture and n D = n f /n max is the dimensionless foam texture. This model represents a simplification of foam behavior in porous media without significantly sacrificing the physical phenomena [24] . Other models that associate bubble generation and destruction to the limiting capillary pressure or the gradient of gas pressure can be found in [15, 17] . In this model, the reduction in the mobility of gas by foam is viewed as a reduction of the gas relative permeability k rg :

where k 0 rg is the gas relative permeability when no foam is formed. Another view on the mobility reduction is based on the apparent viscosity of foam [13, 15, 17, 24] . Comprehensive reviews on the mechanisms of bubble creation, destruction, and also on the reduction of gas mobility can be found in [9, 21] .

The following system of equations describes an incompressible two-phase flow in a porous medium:

where (4) is a population balance equation for foam texture [6] and (5)-(6) account for fractional flow and hydrodynamics. We use φ to denote the porosity of the medium, and ρ α , S α , u α , and p α to denote density, saturation, superficial velocity, and pressure, respectively, of phase α. Also, K = K(x) is the intrinsic permeability tensor, and g is the gravity vector. From the fractional flow theory, λ w = k rw /μ w and λ g = k rg /μ g denote the mobility of water and gas phases, respectively, where the viscosity of water and gas are given by μ w and μ g . It is assumed that the porous medium is fully saturated, i.e., S w + S g = 1.

The numerical approach for solving this system of PDEs should be capable of handling several complexities due to discontinuity, non-linearity, stiffness, natural instabilities, among others. The numerical methods should also preserve important properties, such as local conservation of mass, shock capture, nonoscillatory solutions, accurate approximations, and reduced numerical diffusion effects.

To the extent of our knowledge, this problem is usually solved using explicitin-time finite difference schemes [1, 16, 21, 25] . Also, the use of commercial software is prevalent in the literature [7, 20 , and references within]. The most common approach in commercial software is to represent the effect of foam by a factor that reduces the mobility of the gas phase; therefore, bubble creation and destruction are not represented explicitly [9, 21] .

An effective numerical scheme to solve this kind of model and to address its inherent complexities is based on rewriting the problem in terms of global pressure, as in [4] . In this scheme, one has two distinct coupled problems: an elliptic problem and a degenerate hyperbolic problem. The next step is to decouple the system of PDEs into two subsystems of equations, each one of a different nature. In doing so, each subsystem can be solved by specialized methods, such as finite element and finite volume methods, according to their mathematical properties and the relation between precision and computational efficiency required in the resolution of each step. In this sense, for the spatial discretization of the hyperbolic problems we can employ the finite volume method, for instance; for time discretization, a common choice is a finite difference method, while a hybrid finite element method can be applied to the elliptic equations.

In this context, we develop a staggered algorithm to decouple the hydrodynamics from the hyperbolic system, resulting in a scheme that uses a locally conservative hybrid mixed finite element method to approximate the velocity and pressure fields and a high-order finite volume scheme to solve the hyperbolic equations. The two problems are solved in different time scales. Thus, the proposed staggered algorithm is employed to simulate two-phase (water and gas) flow in a heterogeneous porous medium. We compare pure gas-water injection with the gas-water-foam flow. The results show a reduction of gravity override and viscous fingering effects when the foam is present.

This work is organized as follows: in Sect. 2, we define a fractional flow formulation for Eqs. (4)-(6) using the concept of global pressure; in Sect. 3, we present an algorithm to solve the problem using a hybrid mixed finite element method for the elliptic problem and a high-order finite volume method methodology to solve the hyperbolic equations; numerical results are shown in Sect. 4. Finally, in Sect. 5, we present some concluding remarks.

To build a fractional flow model for the water-gas-foam flow in porous media we follow the global pressure approach from [4] . The global pressure is defined as

where P c (S w ) = p g − p w is the capillary pressure. From (6) and (7) the total velocity u is written as [4] : (8) where g = −9.81ĵ m/s 2 andĵ is the unitary vector in vertical direction. It follows directly from the previous definitions that

Using the hypothesis of rigid porous medium, the total fluid velocity u, the global pressure p, the water saturation S w , and the foam texture n D satisfy, in Ω × (0, T ], the following system of equations:

where i denotes a spatial direction, and

and boundary and initial conditions

with Γ N denoting the boundary region with Neumann condition (specified injection velocity), Γ D defining the boundary region with Dirichlet condition on the potential, Γ − N = x ∈ Γ N ;ū(x) < 0 and n is the unit outer normal vector to Γ . For simplicity, we assume Γ D = ∅ and Γ − N = ∅.

In this section, we introduce the sequential algorithm that combines two kinds of numerical methods to solve (9)- (11) .

The hydrodynamics (9)-(10) is approximated using a naturally stable mixed finite element method introduced in [22] . This method is locally conservative, relying on the strong imposition of the continuity of normal fluxes and on a discontinuous pressure field. The combination of a hybrid formulation with a static condensation technique reduces the number of degrees of freedom in the global problem.

The transport system (11) is solved using the KNP method, a conservative, high-order, central-upwind finite volume scheme introduced in [18] that shows reduced numerical diffusion effects. The KNP scheme is an extension of the KT method [19] that generalizes the numerical flux using more precise information about the local propagation velocities. At the same time, the KNP scheme has an upwind nature, since it respects the directions of wave propagation by measuring the one-sided local velocities. KNP is a semi-discrete method based on the REA (Reconstruct Evolve Average) algorithm of Godunov [8] . Furthermore, the KNP scheme allows for using small steps in time without requiring an excessive refinement of the spatial mesh, since the numerical diffusion does not depend on the time step. After discretization in space, the resulting system of ODEs is integrated in time using a BDF (Backward Differentiation Formula), an implicit, multi-step method that is especially indicated to solve stiff equations [5] .

The system of Eqs. (9)-(11) is strongly coupled. It is possible, however, to apply a staggered algorithm to solve an approximate problem composed of two subsystems: an elliptic one, with time step Δt u , and a hyperbolic one, with time step Δt s . Each of them is solved separately using adaptive time steps, as described in Algorithm 1, i.e., one can use smaller time steps to bound the error in the approximations under a certain tolerance. In addition, it is often recommended that Δt u > Δt s , as the time scale of the hydrodynamics is usually much slower than of the transport. In each iteration, approximations for velocity u n+1 and pressure p n+1 fields at t = t n+1 are computed from

supplemented by the boundary conditions (12) .

Then an iterative algorithm is used to find approximations for water saturation (S n+1 w ) and foam texture (n n+1 D ) by solving the following system of PDEs in Ω × (t n , t n+1 ] (for simplicity, we omit the superscript n + 1):

with boundary and initial conditions

Algorithm 1: Sequential algorithm to solve (13)- (15) .

Set initial conditions S 0 w and n 0 D ; n ← 0; t ← 0; ts ← 0; do Compute velocity (u n+1 ) and pressure (p n+1 ) fields using (13); t = t + Δtu; k ← 0; S n+1,0 ← S n ; do Compute (S n+1,k+1 ) using (14) and (15) with u n+1 ; ts = ts + Δts; k = k + 1; while ts < t; S n+1 ← S n+1,k ; n = n + 1; while t < T ;

In the following sections we comment on the methods employed to solve each problem.

When a mixed finite element formulation is used to approximate the Darcy system (13) , it is necessary to simultaneously fulfill the compatibility condition between spaces and to impose the continuity of the normal vector across interelement edges. In addition, the resulting linear system is indefinite, which can restrict the numerical solvers that could be applied. By using a hybrid formulation, the continuity requirement is imposed via Lagrange multipliers, defined on the interelement edges. Furthermore, if the local problems are solvable, it is possible to eliminate all degrees of freedom related to local problems using a static condensation technique, resulting in a considerable reduction of the computational cost, since the global system involves only the degrees of freedom of the Lagrange multiplier. Also, in this case, the global problem is positivedefinite. Once the approximation for the Lagrange multipliers is found, the original degrees of freedom (associated with velocity and pressure) can be computed in local, independent problems.

We first introduce some notations and definitions, for simplicity restricting ourselves to Ω ⊂ R 2 . The three-dimensional case follows directly. Let L 2 (Ω) denote the Hilbert space of square-integrable functions in Ω, with the usual inner product (·, ·) Ω , and let H(div, Ω) be the space of vector functions having each component and divergence in L 2 (Ω).

Assuming Ω is a polygon, we define a partition T h of Ω composed of quadrilaterals and use K to denote an arbitrary element of the partition. The set of edges of K is denoted by ∂K, the set of edges in T h is denoted by E h , and E ∂ h denotes the set of all boundary edges, i.e., those with all points in Γ . Finally, the set of interior edges is denoted by

where h e is the diameter of the edge e ∈ ∂K and h, the mesh parameter, is the element diameter. For each edge of an element K we associate a unit outward normal vector n K .

The (discontinuous) RT spaces of index k [22] are here denoted by U k h × P k h . We define the following sets of functions on the mesh skeleton:

where p k (e) denotes the set of polynomial functions of degree up to k on e. From these definitions, we can write the following hybrid mixed formulation for the hydrodynamics problem (13) 

where A = A(S n w , n n D ) = (Kλ(S n w , n n D )) −1 .

To solve the hybrid formulation (18)- (19) we apply the static condensation technique that consists in a set of algebraic operations, done at the element level, to eliminate all degrees of freedom corresponding to the variables u h and p h , leading to a global system with the degrees of freedom associated with the multipliers only.

We can observe that static condensation causes a major reduction in the size of the global problem, which is now rewritten in terms of the multiplier only. Also, the new system of equations is positive-definite, allowing for using simpler and more robust solvers. In the end, a hybrid formulation associated with static condensation leads to a great reduction of the computational cost required to solve the global problem. In this work, the deal.II library [3] is used to solve this hydrodynamics problem.

The numerical methodology used to approximate the water saturation and bubble texture Eqs. (14) and (15) is a high-order non-oscillatory central-upwind finite volume method proposed in [18] and here referred to as KNP. Like many other finite volume methods, the KNP scheme is based on a grid of control volumes (or cells).

The upwind nature of KNP is because it respects the directions of wave propagation by measuring the one-sided local speeds, given by

on direction i and a cell of index l, where l+1/2 is the right (resp. top) face and j− 1/2 is the left (resp. bottom) face of a cell, S − l±1/2,i is the local reconstruction of S at the left (resp. bottom) side of a face, and S + l±1/2,i is the local reconstruction of S at the right (resp. top) side of a face; Λ max X and Λ min X are the maximum and minimum eigenvalues, respectively, of the Jacobian ∂f i /∂S at S = X. The result of spatial discretization using KNP is the system of ODEs in conservative form:

where Φ l = Φ (S l ), h i is the cell size in the i-th direction, with the convective numerical fluxes given by

and diffusive numerical fluxes given by

whereD l±1/2,i is defined as the harmonic mean of (DB) l and (DB) l±1,i . The scheme (21)-(23), combined with minmod reconstruction of the typẽ

, whereS l is a piecewise linear approximation to the solution at time t n , i.e.,S l (x) ≈ S n l (x). Then we can use the fact that x l±1/2,i = x l,i ± h i /2 to find S ± l±1/2,i . Various numerical methods can be used to solve the system of ODEs (21) . In this work, a variable order, adaptive step Backward Differentiation Formula (BDF) was chosen. This stable, implicit scheme allows for taking larger time steps than an explicit method would require, which reduces computational cost. In our numerical simulations we used the implementation of the BDF scheme from the CVode package, available in the SUNDIALS library [10] .

Applying the numerical methods described in Sect. 3, we now present results of numerical experiments that aim to assess the influence of foam and gravity effects in two-phase flow. Two scenarios are simulated: flow without and with foam. In the first scenario, we consider that only a mixture of water and gas is flowing through the porous medium, setting n D = 0 in the hyperbolic problem and solving only for S w in (11) . The hydrodynamics and the mobility of the gas phase remain unchanged (k rg = k 0 rg ). In the second scenario, we assume that surfactant is readily available in the water phase, allowing for foam creation and changes in the mobility of gas phase. This scenario is simulated using the full problem (9)- (11) . In both scenarios, the capillary pressure and relative permeabilities are

The permeability is assumed isotropic K = κ(x)I, where κ(x) is the permeability field of layers 1 (case A, Fig. 1(a) ) and 36 (case B, Fig. 1(b) ) of the SPE10 project 1 , rotated to the xy plane. The right boundary is chosen to be the Dirichlet type (Γ D ) withp = 0, while left, top and bottom boundaries are set to Neumann condition (Γ N ) withū < 0 for the left boundary andū = 0 for the top and bottom boundaries. Coefficients and numerical parameters used in the simulations are shown in Table 1 . 

The water saturation profiles for case A at t = 2 000 s and t = 10 000 s are shown in Figs. 2 and 3 , respectively. The gravity effects are much more pronounced in the no-foam simulation. Also, as expected, the water phase displacement occurs more slowly in the foam presence due to the gas mobility reduction caused by foam. Note that, without foam, the gas breakthrough has already taken place at t = 10 000 s (Fig. 3) , which does not occur when foam is present. Moreover, for the foam model adopted [2] , viscous fingering and gravity override are reduced with foam as time advances. As a result, a better sweep efficiency of the medium is observed when foam is present, as can be seen in Fig. 6(a) . In our experiments, foam injection increased total water recovery by approximately 100%. In case B, the permeability field has a more evident preferential channel in the lower region (see Fig. 1(b) ). The results for this channelized porous formation (Figs. 4 and 5) reinforce the foam's ability to reduce the effects of gravity override and viscous fingering, according to the model used, even though this case presents a more pronounced preferential path. The water cumulative production curves for case B (Fig. 6(b) ) show that gravity effects on the production are much more pronounced when no foam is present; this is due to the influence of gravity on diverting the flow from the high permeability channel, resulting in higher sweeping efficiency. Moreover, foam injection in this case increases total water production by approximately 50%, when gravity is considered, and by about 100% when gravity is neglected. 

In this work, we presented a locally conservative numerical algorithm to solve the gas-water flow including foam injection. The system of PDEs that models this phenomenon was derived considering a fractional flow formulation based on the global pressure. The numerical staggered approach proposed combines a high-order central upwind finite volume method for the hyperbolic equations adopting BDF time integrations with a hybrid finite element method to solve the Darcy's problem employing Raviart-Thomas spaces. The proposed methodology was applied to simulate regimes with pure gaswater injection and gas-water-foam flow. In this context, we have established a comparison between these two regimes considering two layers of SPE10 project with different heterogeneous permeability fields. The results, based on the model proposed in [2] , point to the foam's ability to reduce the gravity override and viscous fingering even in cases of porous media with rather pronounced preferential channels.

Mechanistic foam modeling and simulations: gas injection during surfactant-alternating-gas processes using foam-catastrophe theory

Roles of transient and local equilibrium foam behavior in porous media: traveling wave

II -a general-purpose objectoriented finite element library

Mathematical Models and Finite Elements for Reservoir Simulation: Single Phase, Multiphase and Multicomponent Flows Through Porous Media

Integration of stiff equations

Development of a mechanistic foam simulator: the population balance and generation by snap-off

Effect of permeability on implicit-texture foam model parameters and the limiting capillary pressure

A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics

Foam flow in porous media: concepts, models and challenges

SUNDIALS: suite of nonlinear and differential/algebraic equation solvers

A review of the steam foam process mechanisms

The steam-foam process

Mechanisms of foam flow in porous media: apparent viscosity in smooth capillaries

Foam flow in a model porous medium: I. The effect of foam coarsening

Improved mechanistic foam simulation with foam catastrophe theory

Dynamic simulations with an improved model for foam generation

A mechanistic population balance model for transient and steady-state foam flow in boise sandstone

Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi equations

New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations

Non-uniqueness, numerical artifacts, and parameter sensitivity in simulating steady-state and transient foam flow through porous media

Modeling techniques for foam flow in porous media

A mixed finite element method for 2-nd order elliptic problems

Surfactant-Based Mobility Control

A new stochastic bubble population model for foam in porous media

Numerical analysis of a new stochastic bubble population foam model

The authors are thankful to Professor Pacelli P. L. Zitha for fruitful preliminary discussions.