Geofluids

Volume 2018, Article ID 8482352, 15 pages

https://doi.org/10.1155/2018/8482352

## Acceleration of Gas Reservoir Simulation Using Proper Orthogonal Decomposition

^{1}National Engineering Laboratory for Pipeline Safety and MOE Key Laboratory of Petroleum Engineering and Beijing Key Laboratory of Urban Oil and Gas Distribution Technology, China University of Petroleum, Beijing 102249, China^{2}Key Laboratory of Thermo-Fluid Science and Engineering, Xi’an Jiaotong University, Ministry of Education, Xi’an 710049, China^{3}School of Mechanical Engineering, Beijing Institute of Petrochemical Technology, Beijing 102617, China^{4}Key Laboratory of Railway Vehicle Thermal Engineering, Lanzhou Jiaotong University, Ministry of Education, Lanzhou 730070, China

Correspondence should be addressed to Bo Yu; moc.361.piv@xobobuy

Received 3 July 2017; Revised 26 August 2017; Accepted 17 September 2017; Published 3 January 2018

Academic Editor: Jianchao Cai

Copyright © 2018 Yi Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

High-precision and high-speed reservoir simulation is important in engineering. Proper orthogonal decomposition (POD) is introduced to accelerate the reservoir simulation of gas flow in single-continuum porous media via establishing a reduced-order model by POD combined with Galerkin projection. Determination of the optimal mode number in the reduced-order model is discussed to ensure high-precision reconstruction with large acceleration. The typical POD model can achieve high precision for both ideal gas and real gas using only 10 POD modes. However, acceleration of computation can only be achieved for ideal gas. The obstacle of POD acceleration for real gas is that the computational time is mainly occupied by the equation of state (EOS). An approximation method is proposed to largely promote the computational speed of the POD model for real gas flow without decreasing the precision. The improved POD model shows much higher acceleration of computation with high precision for different reservoirs and different pressures. It is confirmed that the acceleration of the real gas reservoir simulation should use the approximation method instead of the computation of EOS.

#### 1. Introduction

Reservoir simulation, utilizing mathematical models to predict fluid flow in petroleum reservoirs, has been developed since the 1800s [1]. It plays more and more important roles in the development of oil exploration and production [2]. Correspondingly, mathematical theories and numerical methods have attracted extensive attention [3–10]. The computational cost of the oil reservoir simulations is huge because a very dense mesh should be used causing long-time iterations to achieve enough resolution. This kind of long-time computation is usually not endurable in engineering. Petroleum engineers attempted several approaches to speed up oil reservoir simulations, including the scale-up approach [11–15] and parallel computation [16]. Proper orthogonal decomposition (POD) can reduce computational loads efficiently via a reduced-order model. It has been applied to heated-crude-oil pipe flow [17, 18] and other fields [19–21]. Ghommem et al. [22, 23] and Efendiev et al. [24] discussed the POD and dynamic mode decomposition method for time-dependent incompressible single-phase flow in high-contrast heterogeneous porous media. They proved that the POD model can achieve good precision in an appropriate range. With fast progress of the exploration and development of gas reservoirs, gas reservoir simulation attracted more and more attention [25–28]. Additional computation should be made to calculate the compressibility of the gas, which does not exist in incompressible oil flow. Thus, the numerical simulations of gas reservoirs are more complex and time-consuming. The study on the acceleration method of gas reservoir simulation by POD is meaningful. On the other hand, the compressibility of the gas leads to a POD model with stronger nonlinearity than oil flow. This situation might cause the POD modeling process to be more difficult than the previous POD modeling for incompressible single-phase flow. To the best of the authors’ knowledge, there is no report concerning the POD model for gas reservoir simulation especially for real gas. Therefore, it is worth discussing the POD modeling method for gas reservoir simulations.

This paper is organized as follows. We revisit the governing equations of gas reservoir simulation and explain how to apply the POD to the gas reservoir simulation briefly. Then, a POD model will be established for gas reservoir simulation in single-continuum porous media. Key effects will be discussed on the acceleration ability of the POD model. Finally, conclusions and suggestions will be made.

#### 2. Numerical Methods

##### 2.1. Governing Equations of Gas Reservoir Simulation

Pure gas flow obeys mass conservation law (see (1)), Darcy’s law (see (2)), and the gas law (see (3)) as follows:where is the Darcy velocity of gas flow in reservoirs, is the diagonal permeability tensor, is the pressure, is the density of gas, is the dynamic viscosity of gas, is the porosity of porous media, is the injection or production rate, is the molecular weight of gas, is the universal gas constant, is the temperature of gas, and is the compressibility factor of gas. should be calculated by the equation of state of gas. In this paper, we use the Van der Waals equation:where is the molar volume of gas, and and are gas constants changing with different real gases. The compressibility of gas for isothermal flow can be defined as

Using pressure as the primary variable, (1) and (5) can be rewritten as follows:

For the convenience of expression, we only use two-dimensional cases in this paper, but the extension to three-dimensional cases is straightforward. Thus, the above governing equations have the following expressions:where and are two components of permeability in the and directions, respectively, and and are two components of Darcy velocity in the and directions, respectively. Gas reservoir simulation can be accomplished via numerically solving (8)(10) to obtain pressure and velocity with the calculation of parameters ( and ) using (4) and (7). Here, is the implicit function of (shown in (7)) so that it can only be explicitly solved. Therefore, the whole term is treated explicitly.

##### 2.2. Basic Idea of Proper Orthogonal Decomposition

The above governing equations clearly show that the core of the whole computation is the calculation of . Once is solved, one can directly solve and . Therefore, the acceleration of the computation of is the most important in the POD modeling. The basic idea of the POD modeling for (8) is briefly stated as follows.

(1) Separate the variable in space and time using the following linear combination:where are the POD modes which are the functions of space (), are the temporal coefficients (), and are the identity basis vectors orthogonal with each other:where “” means the inner product of two vectors and is the Kronecker delta taking the value 1 for and 0 for others.

(2) can be calculated from the eigenvalue decomposition of a sample matrix, which can be obtained by aligning the samples in a time series ():where is the number of grid points (including boundary points), and is the number of samples. A kernel can be calculated as follows:where means domain. The matrix form of (14) is

The matrix can be decomposed by eigenvalue decomposition to obtain the eigenvalues and eigenvectors:where is the eigenvector matrix and are the eigenvalues of the corresponding eigenvectors. Along with the descending order of (), the importance of eigenvectors in is also descending. To quantitatively represent the importance of each eigenvector, the energy of the th eigenvector is introduced:Correspondingly, the cumulative energy of the top eigenvectors can be expressed as follows:

POD modes can be obtained as follows:where “” means the Euclidean distance.

(3) Substituting (11) into (8) and projecting the whole equation onto the Hilbert space spanned by the POD modes (), one can finally obtain a linear equation system of in the scale of , where is the number of POD modes in (11). Once the temporal coefficients are solved, the original variable can be directly reconstructed in (11), instead of the computation of (8). The computational cost of (11) is in the scale of while the computational cost of (8) is in the scale of . The number of samples is usually quite smaller than the number of grid points (). The number of effective POD modes is smaller than the number of samples (). Thus, the computational time can be largely reduced using POD. The POD model is a type of reduced-order model. The details of this step will be explained in Section 2.3.

According to the above analyses, the core computation of POD is the calculation of temporal coefficients . The equation of for gas reservoir simulation will be established in detail in the next section.

##### 2.3. Establishment of the POD Model for Gas Reservoir Simulation

Substituted by (11), (8) can be transformed toEquation (20) is projected onto :where and are the lengths of the computational domain in the and directions, respectively. As stated in (11), and are functions of time and space, respectively. Thus, the left-hand side of (21) is transformed to

Using integration by part and the property of the POD modes (see (12)), we can simplify the first term and the second term of the right-hand side of (21) as follows:

Please note that the projection of the boundary condition in (23) is dependent on boundary pressures and boundary pressure gradients. To update them, the pressure field should be firstly reconstructed using (11) in every time step once is updated by the POD model. Therefore, the projection of the boundary condition leads to additional computations of the full-order equations (see (9), (10), and (11)), consuming much more time. This is very harmful for the high acceleration of the POD model. To ensure the high-acceleration advantage of the POD model, the computations of (9), (10), and (11) should be avoided in the POD model. This can be fulfilled by further decomposing the pressures in (23). The derivations will be made in the following three steps.

(1) If all the boundaries are Dirichlet condition (known boundary pressure), boundary treatment can be expressed as follows:where the boundary pressures , , , and are constant so that they remain in the expressions, and their neighbor pressures () are changing with time so that they are decomposed into and using (11). Thus, the projection of the boundary condition can be transformed to

(2) If all the boundaries are Neumann condition (known boundary velocity), boundary treatment can be expressed as follows:The projection of the boundary condition can be transformed to

(3) For general boundary conditions, the projection of the boundary condition can be summarized via both (26) and (28):wherewhere and are 1 for Dirichlet boundary condition and 0 for Neumann boundary condition. Thus, boundary pressures only appear in the expressions of and while boundary velocities only appear in the expression of . , , and can also be calculated only once, saving computational time. LetThen, the POD model for gas reservoir simulation can be obtained in the following expression:

Equation (34) is the time evolution equation of . For each time step, (34) can be solved via an efficient linear solver (DLSARG) provided by FORTRAN coding language. Numerical results are discussed in the next section to evaluate the precision and the acceleration performances of the POD reduced-order model.

#### 3. Results and Discussions

The computational domain (100 m × 100 m) is shown in Figure 1. Permeability takes the value in the blue area (50 m × 50 m in the center of the domain) and in the red area. Natural gas with the main component of methane (CH_{4}) is the most important gas in subsurface reservoirs. Thus, it is used in the numerical cases in this paper. Initial pressure is given in the whole domain with zero injection or production. and are imposed on the left and the right boundaries. No flow boundary condition is applied to the top and bottom boundaries. All parameters are listed in Table 1, where .