Iterative Solution of the Pressure Problem for the Multiphase Filtration

Applied problems of oil and gas recovery are studied numerically using the mathematical models of multiphase fluid flows in porous media. The basic model includes the continuity equations and the Darcy laws for each phase, as well as the algebraic expression for the sum of saturations. Primary computational algorithms are implemented for such problems using the pressure equation. In this paper, we highlight the basic properties of the pressure problem and discuss the necessity of their fulfillment at the discrete level. The resulting elliptic problem for the pressure equation is characterized by a non-selfadjoint operator. Possibilities of approximate solving the elliptic problem are considered using the iterative methods. Special attention is given to the numerical algorithms for calculating the pressure on parallel computers.


Introduction
Mathematical modeling of multiphase flows in porous media is of great importance in oil and gas recovery. Traditionally, the hydrodynamic simulators for these applications are based on three-phase black oil model [1,11]. A mathematical model of fluid dynamics in porous media includes differential equations, which express the conservation laws of mass and momentum [3,6]. First of all, the continuity equations describing the mass conservation law for each separate phase are used. The momentum equations in a porous medium are written in the form of Darcy's law, which links the velocity with the pressure. When capillary effects are omitted, the pressure is common to all phases.
Applied mathematical models of mass transfer processes in porous media are essentially nonlinear and difficult to study [20,23]. Next, it is necessary in these models to implement the closure of the system of equations via the constant sum of all saturations. Such algebraic components of the model should be taken into account in constructing computational algorithms to predict multiphase flows in porous media. [5,7].
Two classes of methods are used to solve approximately unsteady boundary value problems for coupled systems of partial differential equations. The first of them employs various implicit schemes for the initial system of equations. In this case, we face some computational problems in the transition to a new time level. The second class of methods reduces computational costs by means of using splitting schemes and solving simpler problems at the new time levelsplitting with respect to physical processes [9,15]. The above two classes of methods are presented through the fully implicit method (FIM) and the implicit pressure explicit saturation (IMPES) approach [1,5,11].
FIM is widely used in hydrodynamic modeling of oil and gas recovery [25]. The fully implicit approximation is used in FIM for all equations of the mathematical model. It allows to expect stability of the method and possibility to use large time steps. The basic drawback of the method is connected with its complexity -we have to solve a large system of nonlinear equations.
IMPES provides more efficient algorithms for solving the problem at each time level. In this approach, we formulate the problem for the pressure with implicit approximations in time. After evaluation of the pressure all other unknowns are calculated via explicit approximations. Unfortunately, the problem of stability (time step restriction) is typical for the IMPES method. Therefore, various modifications of IMPES have been developed in order to improve its stability. For example, after evaluation of the pressure we can use implicit approximations for the calculation of saturations [24] (the sequential method).
In this paper, we highlight the main features of the pressure problem which should be taken into account in constructing computational algorithms. Possibilities of obtaining the pressure equation -elliptic for incompressible media and parabolic for compressible ones are discussed here. If the condition of the constant sum of saturations is treated explicitly then the corresponding elliptic operator of the pressure problem is non-selfadjoint. This fact should be taken into account in constructing iterative algorithms.
The authors do not claim to develop algorithms for industrial applications. Our interest is academic. We have identified an important feature of applied filtration problems and discussed some possibilities to take into account peculiarities of the considered problems in their numerical solution. The analysis is based on model problems, which, of course, are very far from real-world simulations.
The paper is organized as follows. A basic system of equations is formulated in Section 2 to describe multiphase fluid flows in porous media. This mathematical model is obtained at assumptions that capillary and gravity forces are negligible. The pressure equation is derived in this section. The main features of the grid problem for the pressure are discussed in Section 3. The simplest uniform grids for the problem in a rectangle are used.
The emphasis is on iterative methods for calculating the pressure at the new time level. The two-dimensional test problem is described in Section 4. The possibility of using standard iterative methods with preconditioners is discussed in Section 5. In Section 6, we present the results of the iterative solving model pressure problems on parallel computers. Conclusions are summarized in Section 7.

The Pressure Problem
In this section we formulate the basic mathematical model for fluid flows in porous media. The system of governing equations for multiphase flows includes the continuity equation for each phase, where α = 1, 2, . . . , m -the phase index. The mass conservation law for each particular phase is expressed by the following equation Here φ stands for the porosity, b α is the phase density, S α -the phase saturation, u α -the velocity, and q α -the volumetric mass source. For simplicity, we neglect the capillary and gravitational forces. In this simplest case the equation of fluid motion in porous media has the form of Darcy's law, where the velocity is directly determined by the common pressure: In (2.2), k is the absolute permeability (in general, symmetric second-rank tensor), k α is the relative permeability, µ α is the phase viscosity and p is the pressure. The unknowns in the system of equations (2.1), (2.2) are the phase saturations S α , α = 1, 2, . . . , m and the pressure (m + 1 unknowns in all). In the simplest case, the coefficients in equations (2.1), (2.2) are defined as some relations Let us consider more convenient forms of system (2.1)-(2.3), which lead to the typical problems of mathematical physics for the pressure. It should be noted that such equivalent formulations do exist only at the differential level. At the discrete level such equivalence of formulations is broken even for linear problems. So, the choice of the initial form of the equations is essential for calculations.
The most natural way to derive the equation for the pressure is the following. Divide each equation (2.1) by φ b α > 0 and add them together, which gives With the natural assumption for compressible fluids d(φ b α ) dp > 0, α = 1, 2, . . . , m equation (2.4) for the pressure is the standard parabolic equation of second order. In particular, the maximum principle holds for its solutions [8]. When using equation (2.4), the basic system of equations for flows in porous media can include m equations Note that the above forms of equations for multiphase flows in porous media are algebraically equivalent only at the differential level. Discussions of problems of preserving the equivalence at the discrete level are presented in [22]. In the case of variable coefficients φ b α , the elliptic operator for the pressure equation (2.4) is non-selfadjoint. This fact leads to some problems in using implicit schemes for equation (2.4). That is why some modifications are employed for the pressure equation. For instance, we can obtain the pressure equation via the direct summation of equations (2.1) taking into account equation In this case we have the selfadjoint elliptic operator for the pressure. However, this approach has some drawbacks. In particular, this system of equations (2.1), (2.2), (2.6) is not closed, because the basic algebraic relation (2.3) is not involved. The most powerful argument in favor of using the model (2.4) in compare with models (2.6) or (2.5) is related to the correctness of the corresponding parabolic problem. In contrast to the model (2.4), we cannot guarantee for models (2.6) or (2.5) the positivity of the coefficient at the pressure derivative with respect to time throughout the whole computational domain. We can have in some subdomains a parabolic equation with the inverse time. In this case we obtain, in general, the ill-posed problems, which require special and more sophisticated computational algorithms for their numerical solution [17].

The Properties of the Grid Operators
Let us consider stationary and unsteady model problems, which are linear prototypes for the pressure problem in modeling multiphase flows. Consider the two-dimensional problem in the rectangle In accordance with (2.4) we solve in Ω the boundary problem for the equa- . . , m, and elliptic operators L α are defined by under the standard assumptions 0 < κ α ≤ k α ≤ κ α . This equation is supplemented with homogeneous Dirichlet boundary conditions In addition, the initial condition is given in the following form In some cases (incompressible media) it is reasonable to consider the stationary problem. The boundary value problem is formulated for the equation which is supplemented by the boundary conditions (3.3). Peculiarities of the considered problems are connected with the differential formulation of the problem and can be taken into consideration when applying various approximations in space. For clarity, we restrict ourselves to the simplest finite-difference approximations on regular grids. The approximate solution is given at the nodes of the uniform rectangular grid in Ω: and let ω be the set of internal nodes (ω = ω ∪∂ω). For grid functions y(x) = 0, x ∈ ∂ω we define the Hilbert space H = L 2 (ω) with the inner product and norm Approximation in space for problem (3.1)-(3.4) will be performed at the assumption that the coefficients and solution are sufficient smooth. For the elliptic operator L α we put into the correspondence the grid operator Λ α : for all α = 1, 2, . . . , m. In H the operators Λ α , α = 1, 2, . . . , m are selfadjoint and positive definite [13,14]: where E is the identity operator, and The grid operator for the pressure problem can be represented as In general (non-constant coefficients a α (x), α = 1, 2, . . . , m) the operator A is non-selfadjoint. It approximates the corresponding differential operator with the error of O(|h| 2 ), where |h| 2 = h 2 1 + h 2 2 . After discretization in space we go from (3.1)-(3.4) to the differentialoperator equation considered on the set of grid functions y(t) ∈ H. The initial condition is taken in the form For the stationary problem (3.3), (3.5) the grid analog has the form To solve approximately problem (3.9), (3.10), we use the standard two-level schemes. Let τ be the fixed time step and y n = y(t n ), t n = nτ , n = 0, 1, . . . , N , N τ = T . Equation (3.9) is approximated by the two-level scheme with weights y n+1 − y n τ + A σy n+1 + (1 − σ)y n = ϕ n , n = 0, 1, . . . , N − 1, (3.12) where, for example, ϕ n = f (σt n+1 + (1 − σ)t n ). It is supplemented by the initial condition The difference scheme (3.12), (3.13) has the approximation error in time O(τ 2 + (σ − 0.5)τ ). If we employ the fully implicit scheme (σ = 1), then the transition to the new time level is performed through solving the grid problem (3.14) The main subject of our consideration is the methods of solving grid problems (3.11) and (3.14), which are linear prototypes for stationary and unsteady problem for the pressure. The primary question here is the non-selfadjoint property of the grid operator A.
For the grid problem (3.11) the maximum principle holds [13]. With regard to considered approximations on the five-point stencil, we formulate it as follows [16]. Consider the difference equation which is supplemented by boundary conditions We assume, that the coefficients of the difference scheme (3.15) satisfy the conditions Let in the difference scheme (3.15)-(3.16) we have ϕ(x) ≥ 0 for all x ∈ ω (or ϕ(x) ≤ 0 for x ∈ ω). Then for we have (the grid maximum principle) y(x) ≥ 0 for all x ∈ ω (y(x) ≤ 0 for x ∈ ω). In our case (see (3.6), (3.8)) fulfillment of the sufficient conditions (3.18) can be verified directly. Because of this, for the grid operator at the new time level (3.14) we have the strict diagonal dominance.
To study properties of operators A and A in Hilbert spaces H = L 2 (Ω) and H = L 2 (ω), it is convenient to treat A and A as the corresponding convectiondiffusion operators. In this case it is possible to employ in our research the results from [10,16]. Based on these studies, we can, in particular, to construct approximations in time, to prove the convergence of an approximate solution to the exact one as well as to construct iterative algorithms for solving grid problems. In our work we have restricted ourselves to using the standard iterative methods without employing advanced results on iterative methods for solving convection-diffusion problems.
Taking into consideration (3.2) and (3.5), we have the representation where The effective diffusion coefficient and the convection velocity for the separate phase α are d α = k α a α , w α = k α grad a α .
Then the pressure operator takes the form of convection-diffusion operator with the convective term in the non-divergent form. Note that application of equation (2.6) to evaluate the pressure corresponds to using only the diffusion part (3.20) of the operator (3.19). Operators of diffusion in the above assumptions about the coefficients are self-adjoint and positive definite in H = L 2 (Ω). Next, we present some facts about the properties of convective transport operators. A detailed discussion of these issues is given in the book [16]. We have the following representation where C α is the operator of convective transport in the symmetric form: The operator C α is skew-symmetric in H: for any w α (x), x ∈ Ω.
From (3.21) and (3.22) we directly obtain the estimate for the energy of the convective transport operator C α : It is interesting to consider the subordination estimate for the operator of convective transport with respect to the diffusion operator. In our model twodimensional problem the corresponding estimate has the following form These properties of differential operators of diffusion and convection (3.21), (3.22), (3.23), (3.25) are inherited not only for the difference operators on rectangular grids [16], but also for difference operators on irregular grids with the Delaunay triangulation [21]. We consider this issue here for the grid operator (3.6), (3.8). First of all, we are interested in the grid analog of (3.19)- (3.20).

similarly (3.19) we obtain
The grid diffusion operator has the form Similarly (3.7) in H = L 2 (ω) we have Approximation of the convective part of the grid operator A is conducted via setting the coefficients (effective velocity w α ) on the grids shifted in the corresponding direction on the half-step. Let us define with accuracy of O(|h| 2 ) the components of the grid analog of w α using the following relations The convective transport operator in representation (3.27) has the form The grid analogue of (3.11) can be written as For the skew-symmetric part C α = −C * α we have The following grid analog of (3.23) takes place: where now (see (3.24)) The subordination inequality (see (3.25) and (3.26)) has the form The fundamental issue here is that for these approximations the constants M α and M α are the complete grid analogues of the corresponding constants M α and M α for the differential problem.

The Test Problem
Capabilities of iterative methods for approximate solving the pressure equation in modeling multiphase flows in porous media are illustrated here using the test grid problem. We consider equation (3.14), which corresponds to the calculation of one time step in the numerical solution of problem (3.1)-(3.4). Numerical experiments are conducted for problem (3.14) with f = 1 in the unit square (l β = 1, β = 1, 2) on the grid h = h 1 = h 2 (N 1 = N 2 ). Particular attention should be given to the coefficients of equations (3.1), (3.2) in order to take into account peculiarities of these problems, namely, inhomogeneity of a α (x), α = 1, 2, . . . , m. Taking into account (2.4), we set , α = 1, 2, . . . , m.
Compressibility of the second phase is defined as follows: The diffusion part of operator (3.20) is The properties of the considered problems are defined (see (3.20)) by the vectors w α , α = 1, 2, . . . , m. For the test problem we have In this case div w 2 = −4ηξ.
For the constants in the estimates (3.23) and (3.25) we obtain Thus, the governing numerical parameters for this problem are η and ξ. The sign of ξ can be any, moreover, it defines the fundamental difference in the behavior of the solution (the pressure) in the vicinity of the production or injection well.

Iterative Solution of the Problem
For numerical solving the test problem we use iterative methods. In the corresponding grid equation (3.14) the operator A is non-selfadjoint. Therefore, we use iterative methods for grid problems with unsymmetric matrices [12,14]. The standard Generalized Minimal Residual Method (GMRES) with different preconditioners has been employed.
To solve the test problem, the PETSc library [2] has been used. The PETSc library, developed in the Argonne National Laboratory, is a powerful set of freely available multi-platform compatible tools for the solution of large-scale problems governed by partial differential equations. Experiments were carried out with the following preconditioners: none -without preconditioning; jacobi -the Jacobi method; sor -the successive overrelaxation method; ilu -the incomplete LU factorization; mg -the multigrid method.
We start the numerical experiments with the steady-state problem, and the influence of the time step is investigated separately. Of course, the stationary problems are of little interest from the practical point of view for the discussed filtration problems. The stationary problem is nothing but the starting point for our methodological investigations of boundary value problems for a model parabolic equation (3.1).  Table 1 shows the dependence of the computational cost (the number of iterations) on the physical parameters of the problem. Features of the problem are clearly defined by the parameters η and ξ. The calculations were performed using the unpreconditioned GMRES method on the grid with 256 × 256 unknowns. We see that with increasing of η and/or |ξ| the number of iterations decreases. The condition of termination for the iterative process is formulated as achieving the relative residual value less than 10 −5 . The same is true for negative and positive values of ξ.
Effect of preconditioning on different grids is shown in Table 2. Calculations were performed at η = 1. It is easy to see that the multigrid preconditioner is the best. In addition, it is interesting to look at the effect of the time step τ . The unpreconditioned GMRES method was used with the grid of 256 × 256 unknowns. From Table 3 we see that the number of iterations decreases with τ .
At large time steps (for τ ≥ 0.1) the spatial discrete operator A is dominant in the operator of the problem 1 τ E + A. Reducing the time step we increase the positive diagonal part and decrease the condition number.

Parallel Implementation
The parallel formulation is based on the domain decomposition methods [18,19]. The main idea of these methods is to divide the original computational domain Ω into subdomains where p is the number of processes. Review of preconditioners that are used for parallel solution of large sparse systems of equations given in [4].
In our case, we use standard 2D decomposition algorithm. Due to the stencil of discretization, the computational domains of processes are taken from the overlap. For distributed memory computers, Message Passing Interface (MPI) is used for sending messages (overlapping elements) between neighbouring processes.
The systems of linear equations are solved by the parallel version of the preconditioned GMRES algorithm. In our computations the none, bjacobi (doing the ILU-factorization of a local part of the matrix at each processor) and multigrid preconditioners were used. The calculations were performed at η = 1.
The parallel code was run on a cluster of North-Eastern Federal University. The cluster consists of four computing nodes, each node has two quad-core processors Intel Xeon E5450 (3.00 GHz) with 16 Gb RAM.
In Table 4, the results of calculations according to the preconditioner are presented. The presented results showed good efficiency of the parallel algorithm for none and bjacobi preconditioners, the number of iterations is almost independent of the number of running processes.
Calculated grid 512 × 512 is the too small for what we have observed the effectiveness of parallelization multigrid iterative method. In this case, overhead costs is sufficiently large. To see an acceleration of the computation time with increasing number of processes we give Table 5. It presents the computation time (number of iterations is always equal to 5) for a sequence of grids. The effectiveness of parallelization on large grids (for example, on the grid 2048 × 2048) is obvious, and the acceleration effect is significant.

Conclusions
1. The basic features of the pressure problem associated with the nonselfadjoint operator are considered for multiphase flows in porous media.
2. It was found that the computational cost of solving the model pressure problem does not depend strongly on ξ, more pronounced dependence is on the physical parameter η. This means that the number of iterations depends basically on the various properties of the phases than on the value of external sources.
3. Parallel computations have been performed using standard techniques with various preconditioners.