Numerical Simulation of Multiphase Multicomponent Flow in Porous Media: Efficiency Analysis of Newton-Based Method

Abstract

Newton’s method has been widely used in simulation multiphase, multicomponent flow in porous media. In addition, to solve systems of linear equations in such problems, the generalized minimal residual method (GMRES) is often used. This paper analyzed the one-dimensional problem of multicomponent fluid flow in a porous medium and solved the system of the algebraic equation with the Newton-GMRES method. We calculated the linear equations with the GMRES, the GMRES with restarts after every m steps—GMRES (m) and preconditioned with Incomplete Lower-Upper factorization, where the factors L and U have the same sparsity pattern as the original matrix—the ILU(0)-GMRES algorithms, respectively, and compared the computation time and convergence. In the course of the research, the influence of the preconditioner and restarts of the GMRES (m) algorithm on the computation time was revealed; in particular, they were able to speed up the program.

1. Introduction

The modeling of hydrodynamic processes that take place in oil reservoirs is one of the difficult problems of fluid mechanics. This is due to the fact that the processes occurring in reservoirs can be very complex. It is necessary to take into account phase transitions, chemical transformations, temperature effects, degassing processes when pressure drops at the bottom and in the reservoir, etc. When simulating multiphase, multicomponent flow in porous media, the phase properties can vary depending on the composition of phases, temperatures (in oil reservoirs, they can reach up to 200 °C) and pressure (in oil reservoirs, it can reach up to 150 MPa). For adequate modeling of complex flow processes that occur in oil reservoirs, it is necessary to take into account the component composition of the phases [1,2]. Accounting for convection, diffusion, phase transitions and chemical reactions makes it difficult to simulate such problems [3,4,5,6].

In [7,8], the authors used various methods and schemes (IMPES, SS, etc.) in order to simulate multiphase, multicomponent flow in porous media. Other methods and schemes can be found in [9,10,11,12]. We could find Newton’s method often in current research works in this area. Newton’s method underlies many algorithms used to solve modeling problems. The compositional modeling of oil and gas reservoirs requires the construction and solution of large, rare, non-linear systems. The solution of sparse linear systems is the most time-consuming part of many sciences and engineering applications; we usually spend over 70% of the time on solving linear systems obtained from Newton’s methods [8].

Linear solvers and preconditioners have been studied for decades, and various methods have been developed [13,14,15]. Saad et al. developed the linear solver of the generalized minimal residuals method (GMRES), which is a general-purpose solver for asymmetric linear systems [16,17]. Vinsome developed the ORTHOMIN solver (a variant of the generalized coupled residuals method with incomplete orthogonalization), which was developed initially for reservoir modeling [18]. The BICGSTAB solver (stabilized bis-conjugate gradient method) is also widely used for oil modeling. We could apply these linear solvers of the Krylov subspace methods to any kinds of linear systems. We can understand the main idea of the Krylov subspace method as the process of recursively constructing the residual vector for each iteration step. Let the linear algebraic equations of the solution be:

$$Ax=b$$

(1)

Take $K_m=K_m\left(A,r_0\right)=span\left\{r_0,A{r}_0,A^2{r}_0,\dots ,A^{m-1}{r}_0\right\}$ where $r_0=b-A{x}_0$, and $x_0$ is the initial solution, looking for an approximate solution $x_m$ from $x_0+K_m$, so that the residual vector is orthogonal to another subspace $L_m$, that is, $r_m=b-A{x}_m\perp L_m$, and then $K_m$ is the Krylov subspace, and the above method is called the Krylov subspace method.

The GMRES algorithm [17] is one of the popular Krylov subspace methods, which is used to solve the solution of asymmetric linear equations. GMRES has very small residuals, but as the number of iterations increases, the amount of calculation and storage increases linearly, which is not efficient. So, the commonly used restarted GMRES, denoted by GMRES(m) [19], performs $m$ iterations of GMRES, and then the resulting approximate solution is used as the initial guess to start another set of $m$ iterations. This process repeats until the residual norm is small enough. The group of $m$ iterations between successive restarts is referred to as a cycle, and $m$ is referred to as the restart parameter. We indicate the restart cycle number with a subscript: $x_i$ is the approximate solution after $i$ cycles, or $m\times i$ total iterations, and $r_i$ is the corresponding residual. The shortcoming of this method is that it reduces the robustness of the original GMRES algorithm and cannot guarantee the convergence of the algorithm. In addition, due to the fact that the restarted algorithm convergence may become very slow, this may affect the good convergence of the GMRES algorithm on specific issues due to the selection problem of $m$.

Besides the restarted GMRES(m) algorithm, another way to reduce the number of iterations of the GMRES method is that we use preconditioning as an explicit or implicit modification to the system of linear equations. Therefore, similar to other iterative methods, we usually combine the GMRES algorithm with a preconditioner for accelerating convergence. We could effectively and easily realize the incomplete LU factorization for sparse matrices with no fill or ILU(0) to obtain the preconditioning matrices [16,17,20,21]. Incomplete LU factorization assumes that sparse matrix A is represented in the form A = LU + R, where the matrices L and U are lower and upper triangular matrices, respectively, and the matrix R is the residual of the decomposition. After the ILU factorization is performed, the matrix M = LU can be used as a preconditioner. ILU factorization is also the basis of many preconditions, such as domain decomposition methods.

This article considers the mathematical model describing the process of multiphase, multicomponent fluid flow linearized using Newton’s method in a simulation, and we solve the system of the algebraic equation with one of the Krylov subspace methods—the Newton-GMRES method. We especially consider a one-dimensional problem of multicomponent filtration of a liquid in a porous medium with the Newton-GMRES method. We solve and analyze the problem with the pure GMRES method, restarted GMRES(m) and preconditioned ILU(0)-GMRES method.

2. Mathematical Modeling

The considered model of multiphase, multicomponent flow in porous media makes it possible to simulate such modern methods of hydrocarbon production as the injection of coolant into an oil reservoir, combustion, injection of chemical reagents (polymer, surfactant) and others.

Let us consider a mathematical model of the multicomponent flow of three-phase (w—water, o—oil, g—gas) fluid in a porous medium. We write the mass conservation law for each component [8]:

$$\frac{\partial \left(\varphi {\xi }_w{S}_w\right)}{\partial t}+\nabla ·\left({\xi }_w{u}_{\mathit{w}}\right)=q_w,$$

(2)

$$\frac{\partial \left(\varphi \left(x_{mo}{\xi }_o{S}_o+x_{mg}{\xi }_g{S}_g\right)\right)}{\partial t}+\nabla ·\left(x_{mo}{\xi }_o{u}_o+x_{mg}{\xi }_g{u}_g\right)=q_m$$

(3)

$m=1.2,\dots N_c$, where $N_c$—the number of components; $\varphi$—the porosity of porous media; $x_{mo}, x_{mg}$—the mole fraction of component m in the oil and gas phases; ${\xi }_w$, ${\xi }_o$, ${\xi }_g$—the molar density of the water, oil and gas phases; $S_w,S_o, S_g$—the saturation of the water, oil and gas phases; $q_w, q_m$—the molar velocities of the water flow and m-th component, respectively. In Equations (2) and (3), the flow rates $u_{\mathsf{\alpha }}$ are given by the Darcy law:

$$u_{\mathsf{\alpha }}=-\frac{k_{r\alpha }}{{\mu }_{\alpha }}k\left(\nabla p_{\alpha }-{\rho }_{\alpha }\wp \nabla z\right), \alpha =w,o,g$$

(4)

where ℘—the gravitational acceleration; ${\rho }_{\alpha }$—the pressure of phase α.

In addition to the differential Equations (2)–(4), the following algebraic expressions hold:

$${\displaystyle \sum }_{m=1}^{N_c}{x}_{mo}=1, {\displaystyle \sum }_{m=1}^{N_c}{x}_{mg}=1$$

(5)

$$S_w+S_o+S_g=1$$

(6)

$$p_{cow}=p_o-p_w, p_{cgo}=p_g-p_o$$

(7)

where $p_{cow}$—the capillary pressure between the oil and water phases; $p_{cgo}$—the capillary pressure between the gas and oil phases.

For the closure of the system (1)–(7), we set the conditions for the equilibrium relation below, which is described by Gibbs law for the composition system:

$$f_{mo}\left(p_o,x_{1o},x_{2o,}\dots ,x_{N_{c}o}\right)=f_{mg}\left(p_g,x_{1g},x_{2g},\dots ,x_{N_{c}g}\right),$$

(8)

$$f_{m\alpha }=p_{\alpha }{x}_{m\alpha }{\phi }_{m\alpha }, m=1.2,\dots N_c,=,o,g$$

(9)

where ${\phi }_{m\alpha }$—the fugacity coefficient of component m in phase α which is oil or gas, and $f_{mo}$ and $f_{mg}$ are fugacity functions of m-th component in the gas and liquid phases, respectively.

We described the distribution of chemical components in the hydrocarbon phase by the K-value method (instantaneous evaporation coefficient) and determined each component by the following expression:

$$K_m=x_{mg}/x_{mo }, m=1,2,\dots ,N_c$$

We described the thermodynamic behavior of fluids under reservoir conditions using Peng–Robinson equations of state [22], in which the relationship between pressure, temperature and phase volume is determined.

Initial and boundary conditions of system (2)–(9) in a one-dimensional field are as follows [0, 1]:

$$S_w\left(x,0\right)=S_w^0\left(x\right) S_o\left(x,0\right)=S_o^0\left(x\right) p_o\left(x,0\right)=p_o^0\left(x\right)$$

$$p_o\left(0,t\right)=p^{inj} p_o\left(1,t\right)=p^{prod}$$

$$S_w\left(0,t\right)=S_w{}^{inj} \frac{\partial S_w}{\partial x}\left(1,t\right)=0$$

$$S_o\left(0,t\right)=0 \frac{\partial S_w}{\partial x}\left(1,t\right)=0$$

where inj—injection; prod—production. $S_g, p_w,$ $p_g$ are calculated using (5) and (6).

3. Numerical Method

In the considered system of differential equations, several components in the hydrocarbon system are provided. From the computational point of view, the computation time and the number of operations for solving this system increase with increasing the number of components. Therefore, it is necessary to group several components with the same chemical compositions into pseudo-components. We wrote the system of equations for three phases and two components as below (neglecting gravitational forces):

$$\varphi \frac{\partial \left({\xi }_w{S}_w\right)}{\partial t}+div\left({\xi }_w{u}_{\mathrm{w}}\right)=q_w$$

$$\varphi \frac{\partial \left(x_{1o}{\xi }_o{S}_o+x_{1g}{\xi }_g{S}_g\right)}{\partial t}+div\left(x_{1o}{\xi }_o{u}_{\mathrm{o}}+x_{1g}{\xi }_g{u}_{\mathrm{g}}\right)=q_1$$

$$\varphi \frac{\partial \left(x_{1o}{\xi }_o{S}_o+x_{2g}{\xi }_g{S}_g\right)}{\partial t}+div\left(x_{2o}{\xi }_o{u}_{\mathrm{o}}+x_{2g}{\xi }_g{u}_{\mathrm{g}}\right)=q_2$$

$$u_{\alpha }=-\frac{k_{r\alpha }}{{\mu }_{\alpha }}k\nabla p_{\alpha }, \alpha =o, g, w$$

(10)

$$x_{1o}+x_{2o}=1, x_{1g}+x_{2g}=1, S_o+S_w+S_g=1$$

$$p_{cow}=p_0-p_w, p_{cgo}=p_g-p_0$$

$$f_{1o}\left(p_o, x_{1o}, x_{2o}\right)=f_{1g}\left(p_g, x_{1g}, x_{2g}\right), f_{2o}\left(p_o, x_{1o}, x_{2o}\right)=f_{2g}\left(p_g, x_{1g}, x_{2g}\right)$$

System of Equation (10) has dependent variables (x_1o, x_2o, x_1g, x_2g, u_α, p_α, S_α). By changing variables and using several arithmetic operations (described in detail in the article [23]), you can obtain a nonlinear system of 6 equations with 6 unknowns (x_1o, L, z, F, S, p), where F—the total mass variable; L and V—the mass fractions, z_m—the total mole fraction:

$$F={\xi }_o{S}_o+{\xi }_g{S}_g$$

$$L=\frac{{\xi }_o{S}_o}{F}, V=\frac{{\xi }_g{S}_g}{F}$$

$$z_m=L{x}_{mo}+\left(1-L\right)x_{mg}, m=1,2,\dots ,N_c$$

This system is linearized by Newton’s method. For function $v$, we set:

$$v^{n+1,l+1}=v^{n+1,l}+\Delta v$$

where l refers to the number of Newton’s iteration, and $\Delta v$ represents the increment in this iteration step; n—the time approximation level.

In the (l + 1)-th layer of the Newton iterative process, the value of the unknowns is updated by the following law:

$$x_{mo}^{n+1,l+1}=x_{mo}^{n+1,l}+\Delta x_{mo}^{n+1,l+1}, m=1,2, \dots ,N_c-1,$$

$${ɀ}_{mo}^{n+1,l+1}={ɀ}_{mo}^{n+1,l}+\Delta {ɀ}_{mo}^{n+1,l+1}, m=1,2, \dots ,N_c-1,$$

$$L^{n+1,l+1}=L^{n+1,l}+\Delta L^{n+1,l+1},$$

(11)

$$F^{n+1,l+1}=F^{n+1,l}+\Delta F^{n+1,l+1},$$

$$S^{n+1,l+1}=S^{n+1,l}+\Delta S^{n+1,l+1},$$

$$p^{n+1,l+1}=p^{n+1,l}+\Delta p^{n+1,l+1}$$

The calculation algorithm is as follows [23]:

• Equation (11) is substituted into the finite difference form of a system of nonlinear equations of 6 variables.
• We obtain a system of linear equations with unknown $\Delta p_i$ and solve it using the GMRES method.
• Further, using the known $\Delta p_i$, we find the unknown $\Delta x_{1oi}$ ($\Delta x_{2oi}$ is found by relation (5)), $\Delta L_i$, $\Delta z_i$, $\Delta F_i$, $\Delta S_i$.
• Next, $\Delta x_{1gi}$ and $\Delta x_{2gi}$ are calculated according to (5).
• Next, the unknowns $x_{1oi}$, $L_i$, $z_i$, $F_i$, $S_i$, are found using (10).

As mentioned earlier, the resulting system of linear equations was also solved by the ILU(0)-GMRES method. There is a limitation in the incomplete LU decomposition; namely, problems may arise when the diagonal element of the matrix is equal to zero. This problem can be solved by selecting the reference element or, in other words, by pivoting [24]. In our problem, the diagonal elements of the matrix had a nonzero value, and finding the preconditioner did not require additional operations.

In the Figure 1, you can see a diagram of the computational algorithm.

4. Results and Discussion

To find Δp_i, the algorithms GMRES, restarted GMRES(m) and preconditioned ILU(0)-GMRES were used. The matrix of coefficients of this system of linear equations is shown in

Table 1. Matrices information for numerical experiments.

Matrices Number	#1	#2	#3	#4	#5	#6
Number of nodes	3,240,000	1,2960,000	29,160,000	81,000,000	129,960,000	190,440,000
Nonzero elements	8338	16,738	25,138	41,938	53,138	64,338

for various mesh sizes. The characteristics of the computer are the processor with 4-core Intel Core i7 3770 and RAM with 8 GB. In order to perform the calculation, we created six different dimensional matrices as shown in Table 1, which shows the number of points in each matrix and the number of nonzero elements. We have obtained a comparative graph of the time for calculating one Newton-GMRES iterative calculation (Figure 2) and the number of the GMRES inner iteration (Figure 3). In

Table 2. Acceleration obtained by preconditioning.

Matrices Number	GMRES Runtime(s)	ILU(0)-GMRES Runtime(s)	Acceleration
#1	7.03	1.23	5.715447154
#2	53.54	11.07	4.836495032
#3	169.28	36.86	4.592512208
#4	746.1	164.54	4.534459706
#5	1641.1	349.72	4.692611232
#6	3009.45	643.71	4.675164282

, we can see the acceleration obtained by preconditioning.

As one can see, the GMRES algorithm with ILU(0) preconditioning works 4.5–5 times faster than the algorithm without preconditioning for all matrices. This result allows us to conclude that this speedup will persist as the size of the task increases

As one can see from Figure 2, the GMRES(m) algorithm takes less time than the other two computing algorithms. However, it was selected experimentally by repeatedly giving values ranging from 100 to 1500 to obtain the appropriate value of six different matrices. For example, experimental calculations for matrix #3 are shown in

Table 3. Numerical experiments for matrix 3 with different m values.

Matrix #3	CPU time(s)	Residual	Inner iteration
GMRES(100)	23.374	0.00982457	5286
GMRES(200)	56.605	0.00997898	5502
GMRES(300)	49.395	0.0099962	2942
GMRES(400)	50.373	0.00999233	2184
GMRES(500)	16.83	0.00994755	512
GMRES(600)	24.017	0.00939816	604

; as we have seen from this table, when m = 500, algorithms obtain the minimum computational time and number of iterations, and for matrix #4, as can be seen from

Table 4. Numerical experiments for matrix 4 with different m values.

Matrix #4	CPU time(s)	Residual	Inner iteration
GMRES(100)	35.4249	0.00993622	5788
GMRES(200)	184.447	0.00993313	10513
GMRES(300)	281.269	0.00999677	10751
GMRES(400)	200.089	0.00999497	5471
GMRES(500)	218.66	0.00999747	4473
GMRES(600)	260.9249	0.00999968	4435

, the algorithm is executed faster when m = 100.

As can be seen from Table 4, in some cases, the total number of iterations is fewer than when m = 100, but the algorithm is slower. This is due to the fact that at m = 100, at each restart, the program builds a smaller orthogonalization basis.

The appropriate value of m for the matrix of Table 2 is given in

Table 5. Appropriate m value for different size of matrices.

Matrices Number	#1	#2	#3	#4	#5	#6
Appropriate m value	200	300	500	100	100	100

. As we have noted, in different matrix sizes, the value of m for the GMRES(m) algorithm is different. However, the optimal m value of matrix #4, #5 and #6 is 100, but with substantial inner iterations (Figure 3). Although the ILU(0)-GMRES algorithm is slower than the GMRES(m) algorithm, there is no need to carry out repeated experiments to obtain the optimal value of m. Moreover, with the ILU(0)-GMRES algorithm, we had better convergence then GMRES(100) for matrix #4, #5 and #6. Moreover, the GMRES algorithm is the slowest and most assured the number of iterations.

5. Conclusions

This article discusses the Newton-ILU0-GMRES method for solving the problem of multiphase, multicomponent flow in a porous medium. The equations describing the motion of a multicomponent, multiphase flow in a porous medium were linearized using Newton’s method. The resulting system of linear equations was solved using the GMRES algorithm. In addition, a preconditioner and a modification of the method with restarts (GMRES(m)) were used to speed up the method. The analysis and comparison (calculation time, convergence) of the Newton-GMRES, GMRES(m) and ILU0-GMRES methods were carried out. The numerical results demonstrate that for the one-dimensional problem of multicomponent flow in a porous medium, the GMRES algorithm is the most time-consuming. The GMRES algorithm with ILU(0) preconditioning works 4.5–5 times faster than the algorithm without preconditioning. For the problem size equal or bigger than matrix #4 for the GMRES(m) algorithm, the suitable m value is 100, but convergence is slow. The ILU(0)-GMRES algorithm calculated slower than the GMRES(m) at the larger size of the matrix but with better convergence in solving the given problem.

Similar results can be obtained for two-dimensional and three-dimensional equations by applying to them an analogous method of solving, but with the condition that the matrices of the resulting systems of linear equations satisfy the constraints of Incomplete LU factorization.