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 [310]. 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 [1115] 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 [1921]. 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 [2528]. 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 (CH4) 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 .

In this condition, governing equations (see (4) and (7)(10)) are solved to collect samples at different moments, using the highly accurate finite difference method (FDM) to ensure the grid-independent results for the grid number more than 60 × 60 [29]. Thus, the grid number in this paper (100 × 100) is enough. The temporal scheme is an explicit advancement for each time step () with time steps. After the total time scale of the simulation of 30 days, 2000 samples are collected. To evaluate the precision and computational speed of the POD model, the relative deviation and the acceleration ratio are defined in the following:where is the CPU time for the computation and the subscripts “POD” and “FDM” represent the results from the POD model and from the direct calculation of the governing equations via the FDM. If , the POD model can accelerate the computation. The larger represents the larger acceleration.

3.1. POD Results for Ideal Gas

First of all, the POD model is examined in the case of ideal gas (), where and in the Van der Waals equation (see (4)) are all zero. Thus, the expression of the gas compressibility (see (7)) decays to so that (31) and (33) can be simplified as

As shown in Figure 2, the first POD mode occupies the vast majority (96.86%) of the total energy. This indicates that the first mode captures the main characteristics of the whole transient process. However, the relative deviation of the reconstruction results largely fluctuates with time. The maximum deviation is as high as 210% while the minimum deviation is only % (Table 2). This phenomenon indicates that a dominant mode occupying the vast majority of energy can only obtain high-accurate reconstruction at some moments but loses the fidelity at other moments. To promote the reconstruction precision of the whole transient process, more POD modes should be included. To ensure high acceleration of the reduced-order model, the inclusion of the POD modes should be as little as possible. Thus, there is an optimal truncation number of modes. Theoretically, the number of the POD modes can be usually determined when the cumulative energy contribution achieves 100%. As shown in Figure 2, the cumulative energy contribution achieves 100% at and after the 7th POD mode. Although the maximum deviation decreases rapidly with increasing number of POD modes, the maximum deviation is still as high as 13% when the top 7 modes are used. Therefore, more modes should be considered in the reduced-order model.

Relative deviations along with time, using more POD modes, are shown in Figure 3. With the increasing number of POD modes, the deviation converges rapidly when more than top 10 modes are used. The maximum deviation decreases greatly from 13.06% at the top 7 modes to 6.39% at the top 10 modes. After the top 10 modes, the decrease of the maximum deviation becomes much slower (5.31% for the top 20 modes, 5.25% for the top 30 modes, and 5.21% for the top 40 modes). The minimum deviations are all small. Thus, it can be considered that the reconstruction results converge at the number of top 20 modes. However, the computational time of the POD model is increasing with increasing number of POD modes, causing the acceleration ratio to be largely influenced (Table 3). The largest drop of the acceleration ratio occurs between the top 10 modes and top 20 modes. More POD modes cause the dimension of the POD model to increase much faster so that the computation of a larger equation system consumes much more time. Therefore, the optimal number of POD modes needs a balance between the precision and the acceleration. Between the top 10 modes and top 20 modes, the precision is promoted a little (ε from 6.39% to 5.31%) but the acceleration ratio of the POD model decreases from 24 for the top 10 modes to 6.6 for the top 20 modes, so that the number of the top 10 modes is an appropriate optimal number of POD modes. From the comparison, it should be noted that the theoretical energy criterion is not enough to determine the optimal number of POD modes. The balance of precision and acceleration should be considered. This will include much more POD modes with very small energy contribution, but their actual contributions to precision are important.

It can also be seen from Figure 3 that the maximum deviation (6.39%) and the minimum deviation (%) occur at about 0.8 days and about 22 days when the top 10 modes are used in the POD model. The velocity field and pressure field obtained by the POD model are compared with those obtained by FDM to examine the reconstruction precision of local characteristics. From Figure 4, it is clear that the two components of Darcy velocity and pressure are all reconstructed well with tiny local deviation. From Figure 5, reconstructed velocity component and reconstructed pressure agree well with those of FDM, but reconstructed velocity component can only capture part of the features. Fortunately, this maximum deviation only exists in a very narrow range in Figure 3. The deviation at other moments decreases rapidly. We examine the local flow field at the mean deviation (3.54%) in Figure 6 and find that the local distributions of , , and are well reconstructed. Thus, the POD model using the top 10 POD modes can reconstruct the flow field in high precision for the whole transient process.

3.2. POD Results for Real Gas

For real gas, compressibility factor () and compressibility () should be computed by the equation of state (see (4)) and the equation of real gas compressibility (see (7)). POD results for the same case in Section 3.1 are obtained by solving (34). The reconstruction precision (%6.39%) is as high as that of ideal gas (Figure 3). Flow field comparisons are also very similar to those in Figures 46. Thus, these results are not redundantly shown here. The above POD model maintains high precision for real gas simulation.

However, the acceleration ratio is as low as 0.9 (i.e., 6240 s/7256 s in Table 4), which means the POD model causes the simulation to be even slower than FDM. This is quite different from the high acceleration ratio for ideal gas ( in Table 3). The main difference between the computations of real gas and ideal gas is that the equation of state (EOS) should be computed for real gas. We analyze the different contributions of the total computational time for FDM and POD in Table 4 and find two important phenomena: (1) the largest part of the total computational time is mainly consumed on EOS for both FDM and POD (92% for FDM and 78% for POD); the EOS cannot be expressed as an explicit function of pressure so that it cannot be decomposed by the POD; it can only be calculated locally for both FDM and POD and thus the main computational time of real gas simulation is not reduced in the POD model; (2) the computational time for flow equations in POD (1596 s) is much longer than that in FDM (493 s) with much higher contribution (22%) than that of FDM (8%). This means the POD model is still slower than FDM even not considering the computational time of EOS. This can be explained using the expression of the POD model. In the model (see (34)), all terms are constant for time advancement except the projection of the unsteady term (see (31)) and the projection of the source term (see (33)). They should be calculated in every time step because they contain the variables changing with time (, , and ). This kind of calculation largely increases the total computations.

The above analyses indicate that the improvement of the acceleration ability of the POD model should reduce the time contribution of EOS and avoid the frequent computations of the two projection terms in the POD computation. According to these two points, we propose a new method: EOS is only calculated once at the initial time. Then, the initial compressibility, the initial compressibility factor, and the initial pressure are used to approximate the two projection terms in (31) and (33), which are treated as constants in the following transient computation of the POD model. Through this treatment, the calculations of EOS and the two projection terms in every time step are all avoided for the POD computation so that the computational time of the POD model is expected to be reduced greatly. The approximations of the two projection terms are shown as follows:where , , and are the compressibility, the compressibility factor, and the pressure in the initial condition (), respectively.

Figure 7 shows that the relative deviation of the above treatment is almost the same as that updating the compressibility every time step, indicating that this treatment does not affect the precision of the POD model. Table 5 shows that the total computational time of the POD model has been largely reduced because the computational time of EOS and flow equations are all reduced greatly. The time decrease of EOS is because the EOS is only calculated once. The time decrease of flow equations is because the two projection terms are not calculated with time advancement of the POD model. These results demonstrate that large acceleration ratio () of POD for real gas simulation can only be achieved when the time contribution of EOS is greatly reduced.

3.3. Verification of the Improved POD Model in a More Realistic Case

The improved POD model in the above section is applied to a more complex case of real gas flow simulation to confirm its capability. Permeability field is shown in Figure 8, where the red area has a permeability 100 md and the blue area has a permeability 1 md. The large permeability represents the main flow path of the gas in the reservoir such as soil or sand. The small permeability represents the main obstacles such as rock. Higher boundary pressure on the left border (100 atm) is used because gas pressure is usually high in engineering. Other parameters are the same as the previous case.

As shown in Figure 9, the relative deviations are also as low as the previous case. The maximum and minimum deviations are 6.65% and 0.13%, respectively, indicating high precision of the improved POD model. It is further confirmed in Figures 10 and 11 that the local flow fields are reconstructed well at the maximum and minimum deviations. A little larger deviation of pressure in Figure 10(c) may be caused by the approximation of gas compressibility. However, the main features of the pressure field are correct while the velocity field has much smaller deviation. Thus, the precision is acceptable.

Table 6 shows that the computational time is reduced from 5823 s for FDM to 25 s for POD. The acceleration ratio is as high as 233. Therefore, the improved POD model can still greatly save computational time for more complex reservoir and higher pressure.

4. Conclusions

Proper orthogonal decomposition is utilized in gas reservoir simulation to accelerate the simulation speed of gas flow in single-continuum porous media. High-precision reconstruction can be achieved using only 10 POD modes with the deviations as low as %6.39% for ideal gas and %6.55% for real gas. The acceleration ratio is as high as 24 for the typical POD model of the ideal gas flow. However, the computational speed of the typical POD model of the real gas flow is even slower than FDM. Two key points for improving the computational speed of the POD model are discussed and verified:(1)The computation of EOS should be avoided in the solving process of the POD model because the total computational time is dominated by EOS.(2)POD projection terms containing compressibility of gas should not be updated in every time step. Otherwise, the computational time of flow equations will also be longer than FDM.

According to the two points, we proposed a new method to approximate the projection terms in all time steps using the initial compressibility so that EOS only needs to be calculated once at the initial condition. After this treatment, it is verified in two different cases that the computational speed of the POD model is largely promoted while high precision is retained (0.13%6.65%). The computational time of the real gas reservoir simulation is reduced from 6240 s and 5823 s of FDM to 23 s and 25 s of POD for the two cases. The acceleration ratios are 271 and 233, respectively.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

Acknowledgments

The work presented in this paper has been supported by the National Natural Science Foundation of China (NSFC) (nos. 51576210, 51325603, and 51476073) and the Science Foundation of China University of Petroleum, Beijing (nos. 2462015BJB03, 2462015YQ0409, and C201602). This work is also supported by the Project of Construction of Innovative Teams and Teacher Career Development for Universities and Colleges Under Beijing Municipality (no. IDHT20170507) and the Foundation of Key Laboratory of Thermo-Fluid Science and Engineering (Xi’an Jiaotong University), Ministry of Education, Xi’an 710049, China (KLTFSE2015KF01).