POD-Galerkin Model for Incompressible Single-Phase Flow in Porous Media

2.1 Description of the Original Governing Equations

Incompressible single-phase flow in porous media are described by the mass balance equation (Eq. (1)) and the Darcy’s law (Eq. (2)) as follows:

where u=ui→+vj→ is Darcy velocity, k=kxxkyy is permeability tensor, p is pressure, g is the gravitational acceleration, q is injection or production rate, ρ and µ are density and dynamic viscosity of fluid respectively, z is the depth. The components of Darcy velocity and permeability are u, k_xx in the x direction and v, k_yy in the y direction. Gravitational acceleration can be included into the piezometric head (Φ = p -ρgz) so that the treatment is the same as that without gravity. It is neglected here for convenience. Substituting Eq. (2) without gravity to Eq. (1), we can obtain:

Eq. (3) shows that computational time is primarily used in the computation of pressure, usually requiring long-time iterations for large systems. POD-Galerkin model is a good choice for greatly reducing the computational time of Eq. (3). The main procedure to establish the model will be stated below.

2.2 POD-Galerkin Model

POD method is data dependent so that the first step is collection of samples through computation of governing equations (Eq. (2) and Eq. (3)). Finite difference method (FDM) is used here to obtain discrete equations of the governing equations, which are solved via the Gauss-Seidel iteration method combined with the successive over relaxation method. The samples can be expressed as follows:

where S₁, S₂, S₃, · · · , S_N represent the serials of samples, N is the number of samples, L is the number of grid points (including boundary). All the samples compose an L × N sample matrix:

The number of grid points is usually much larger than the number of samples, i.e. L>>N. Thus, the snapshot POD [19, 20] should be used to obtain POD modes:

where ϕ_n are referred to the POD modes (n = 1~ N), σ_n and V_n are the singular values and eigenvectors obtained from singular value decomposition (SVD) of matrix S^TS.

With the known singular values, we can define:

where e_n is the energy contribution of the nth POD mode to the whole energy spectrum of POD, E_n is the cumulative energy contribution of the first n POD modes (ϕ₁ ~ ϕ_n) to the whole energy spectrum. The higher energy contributions indicate more features of samples captured by POD modes so as to more importance of these modes. This property will be used to determine the selection of POD modes in Section 3.

With the known modes ϕ_n, pressure can be reconstructed by the linear combination:

where M(≤ N) is the number of used POD modes, c_n are unknown coefficients, independent of space, to be calculated from the POD model. To derive this model, we substitute Eq. (8) to Eq. (3) and obtain:

Project Eq. (9) onto the POD mode ϕ_m (m = 1~ M):

Simplify Eq. (10) via the integration by part:

The undefined symbols in Eqs. (10) and (11) are illustrated in Fig. 1. Uniform mesh and staggered grid method are utilized with grid number of 100 × 100 and domain size of 100 m × 100 m in this paper.

Figure 1

Domain and staggered grid

For a Dirichlet boundary condition (known boundary pressure), any one of the following equations could be used:

For a Neumann boundary condition (known boundary velocity), any one of the following equations could be used:

Combining the Dirichlet and Neumann boundary conditions in Eq. (11), the final expression for POD-Galerkin model is:

where

Dirix and Diriy are 1 for Dirichlet boundaries and 0 for Neumann boundaries.

For a prediction case that is different from samples, the coefficients can be calculated from Eq. (14) so that the pressure under the prediction condition can be calculated directly from Eq. (14) instead of Eq. (3). From the above derivation, we know the dimension of Eq. (14) is much smaller than that of Eq. (3), i.e. M ≤ N ≪ L. Thus, computational time can be reduced largely via POD-Galerkin model.

3.1 A Homogeneous Isotropic Porous Medium

Firstly, a homogeneous isotropic porous medium is considered. Permeability components are the same over the whole domain with the value of 100 md (1 md = 9.869233 × 10⁻¹⁶ m²). The parameters and boundary conditions are shown in Fig. 2. Each boundary pressure takes two different values: p₁ = {50, 60} mH₂O, p₂ = {20, 40} mH₂O (1 mH₂O = 9800 Pa). Thus, the number of the samples is 2 × 2 = 4.

Figure 2

A homogeneous isotropic porous medium

According to Section 2, the maximum number of POD modes is equivalent to the number of samples, i.e. four in this case. The selection of the mode number depends on the actual precision of sample reconstruction. The precisions using different number of POD modes are shown in Table. 1. It is obvious that the optimal number of POD modes is two with the relative deviation as low as 0%~ 2 × 10⁻⁴%. The inclusion of ϕ₃ and ϕ₄ generates very large deviations. The reason is analyzed in Fig. 3 and Table. 2. In Fig. 3, the smooth distributions of modes ϕ₁ and ϕ₂ refect the sample distribution correctly (Fig. 3(a) & 3(b)) while the fluctuated distributions of modes ϕ₃ and ϕ₄ indicate the unphysical numerical errors because flows in porous media are low-speed laminar flow without fluctuations. This is further verified in Table. 2 that ϕ₃ and ϕ₄ do not contain any energy indicating that they do not contain any information about the samples. These two modes are only numerical errors produced by the SVD. The sole inclusion of ϕ₁ could not generate accurate enough reconstructed results (e = 2.82% ~ 10.11% in Table. 1) although its energy contribution is very high (e_n = 99.35%). The inclusion of ϕ₂ supplements a small energy contribution (e_n = 0.65%) to make the POD model capture the whole important information (E_n = 100%) and promote model precision largely (e = 0% ~ 2 × 10⁻⁴% in Table. 1). Thus, ϕ₂ cannot be neglected even if its energy contribution is quite lower than ϕ₁. The same situations also occur in Sections 3.2 ~ 3.5 according to the computational results so that the discussion on the selection of POD modes will not be repeated in the following sections and ϕ₁ and ϕ₂ will be used directly.

Figure 3

Distribution of POD modes

To validate the precision of POD model for predicting cases different from the samples, prediction conditions are designed as p₁ = {10, 20, 30, 40, 50, 60, 70, 80, 90, 100} mH₂O, p₂ = {5, 15, 25, 35, 45, 55, 65, 75, 85, 95} mH₂O so that the number of prediction cases is 10 × 10 = 100. The relative deviations of 100 cases are shown in Fig. 4 with the maximum value of 4.0 × 10⁻⁴% and the minimum value of 0%, indicating good performance of solutions globally. Three typical cases (a, b, c) are selected to show the local solutions of the flow field. The comparisons are shown in Fig. 5. Local solutions of p, u, v using POD coincide very well with those using FDM, no matter whether the boundary conditions fall within (a) or out of (b, c) the sample scope. Even if the boundary condition of case c (p₁ < p₂) is opposite to that of samples (p₁ > p₂), the deviation is also small enough. Therefore, the POD model is proved to be very accurate for the homogeneous isotropic porous medium.

Figure 4

Relative deviation of the prediction cases for a homogeneous isotropic porous medium

Figure 5

Comparison of flow fields for 3 typical prediction cases for a homogeneous isotropic porous medium: from left to right are p, u, v respectively; black solid line-FDM, red dotted line-POD

To validate the computational speed of the POD model, computational time is compared in Table. 3. The computational time of the training stage (both the 4 sample cases and the SVD) only costs once for the POD model. Once the training stage is completed, we can use the POD model to predict any number of cases without further sampling and decomposition. Thus, the one-time cost of training time should not be considered into the acceleration ratio of POD model, even if it is small. For the 100 prediction cases, the computational times of FDM and POD are 678 s and 0.2s respectively. The acceleration ratio is high (3390), indicating excellent acceleration ability of POD model.

3.2 A Homogeneous Anisotropic Porous Medium

In this section, a homogeneous anisotropic porous medium is considered as shown in Fig. 6. The permeability tensor is unique all over the domain, but the component k_xx is 100 times larger than component k_yy. Different boundary conditions are adopted to avoid particularity. Four sample cases are taken with the same boundary pressures of Section 3.1. The 100 prediction cases are designed with the boundary pressures:

Figure 6

A homogeneous anisotropic porous medium

p₁ = {10, 20, 30, 40, 50, 60, 70, 80, 90, 100} mH₂O,

p₂ = {15, 25, 35, 45, 55, 65, 75, 85, 95, 105} mH₂O.

Fig. 7 shows the high precision of POD model for 100 cases. The relative deviations are low (1.0 × 10⁻⁴% ~ 5.2 × 10⁻³%). The precision is further examined using two cases (a, b) with the maximum and minimum deviations. For these cases, POD results agree very well with FDM results (as shown in Fig. 8). Thus, POD model is very accurate for homogeneous anisotropic porous medium.

Figure 7

Relative deviation of the prediction cases for a homogeneous anisotropic porous medium

Figure 8

Comparison of flow fields for 2 typical prediction cases for a homogeneous anisotropic porous medium: from left to right are p, u, v respectively; black solid line-FDM, red dotted line-POD

For the 100 prediction cases, the computational times are compared in Table. 4 for FDM and POD. The acceleration ratio is still very high (r = 2300), showing the good acceleration ability of the POD model.

3.3 An Inhomogeneous Isotropic Porous Medium

In this section, an inhomogeneous isotropic porous medium is designed as shown in Fig. 9. The distribution of permeability is not uniform with a small value in the “H”-shape zone and a large value in other zone, but the two components k_xx and k_yy are always equal to each other. Sample cases are designed as p₁ = {50, 55} mH₂O, p₂ = {60, 65} mH₂O. Prediction cases are designed as p₁ = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90} mH₂O, p₂ = {10, 20, 30, 40, 50, 60, 70, 80, 90, 100} mH₂O. It should be noticed that there are no flows for the 9 cases when p₁ = p₂ so those are not discussed here. Therefore, the total effective prediction cases must be 10 × 10 - 9 = 91.

Figure 9

An inhomogeneous isotropic porous medium

With the 2 POD modes extracted from the 4 samples, predictions for 91 cases are made. The relative deviations are very low (7.0 × 10⁻⁴%~ 1.4 × 10⁻²%) as shown in Fig. 10. Pressure and velocity coincide very well with each other for both POD and FDM in the two cases with the maximum and minimum deviations (Fig. 11). Thus, POD model is again proven to be very accurate for inhomogeneous isotropic porous medium.

Figure 10

Relative deviation of the prediction cases for an inhomogeneous isotropic porous medium

Figure 11

Comparison of flow fields for 2 typical prediction cases for an inhomogeneous isotropic porous medium: from left to right are p, u, v respectively; black solid line-FDM, red dotted line-POD

Table. 5 shows the acceleration ratio of the POD model decreases to 880. The decrease is due to the shorter computational time for FDM caused by different permeability field and boundary conditions. Nevertheless, the acceleration ratio is still attractive.

3.4 An Inhomogeneous Anisotropic Porous Medium

To examine the POD model in a more complex case, we let k_xx > k_yy in the “H”-shape zone and k_xx < k_yy in the other zone so that an inhomogeneous anisotropic porous medium is designed as shown in Fig. 12. All other sampling and prediction conditions are the same as Section 3.3.

Figure 12

An inhomogeneous anisotropic porous medium

Figure 13

Relative deviation of the prediction cases for an inhomogeneous anisotropic porous medium

The relative deviations of the prediction cases are apparently larger than the previous situations, indicating that the precision of POD model may decrease along with the increasing complexity of the problems. However, the precision is still high (ɛ = 1.2 × 10⁻²% ~ 2.2 × 10⁻¹%) in Fig. 13. The flow fields still coincide with each other very well for both POD and FDM results, in the two cases with maximum and minimum deviations (Fig. 14). Thus, the POD model has also proven to be very accurate for inhomogeneous isotropic porous medium.

Figure 14

Comparison of flow fields for 2 typical prediction cases for an inhomogeneous anisotropic porous medium: from left to right are p, u, v respectively; black solid line-FDM, red dotted line-POD

Figure 15

A random porous medium

Table. 6 shows the acceleration ratio of the prediction cases. It is 98454, which exceeds all other cases in the above sections. The extremely high acceleration ratio is due to longer computational time for FDM due to much higher complexity of the problem.

3.5 A Random Porous Medium

From the comparisons in Sections 3.1 ~ 3.4, we confirmed that POD model has high fidelity in precision and large acceleration in computational time for either homogeneous or inhomogeneous and either isotropic or anisotropic porous media. However, the computational settings are basically ideal, e.g. the “H”-shaped permeability distribution. Permeability fields in real engineering are more complex than these structures. In order to achieve stronger conclusions, we use randomly distributed, inhomogeneous and anisotropic permeability field and recall other settings in the case in Section 3.1. The domain and parameters are shown in Fig. 15 while the permeability fields are shown in Fig. 16.Two components of permeability are even in different ranges: k_xx for 1 md ~ 10 md and k_yy for 10 md ~ 100 md. The sampling conditions are: p₁ = {30, 60} mH₂O, p₂ = {20, 80} mH₂O. The prediction conditions are:

p₁ = {10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110} mH₂O,

p₂ = (0, 15, 30, 45, 60, 75, 90, 105, 120, 135) mH₂O.

Figure 16

Random distribution of the inhomogeneous-anisotropic permeability field

As shown in Fig. 17, the relative deviations are 4.5 × 10⁻³% ~ 2.3 × 10⁻¹% even if 107 prediction cases are predicted by only 2 POD modes for the random permeability field. This overall precision is high enough to be accepted in engineering. To identify local precision, 4 typical cases (a, b, c, d) among these 107 cases are selected to compare the pressure field and the velocity field obtained by POD and FDM. If the same style of the previous comparison (e.g. style in Fig. 5, Fig. 8, Fig. 11, Fig. 14) is adopted, the local details cannot be recognized clearly. Therefore, we change to a new comparison style in Fig. 18, where the figures on the first row represent results by FDM while the figures on the second row represent the results by POD. The comparison shows that POD results agree very well with the FDM results, even though the prediction conditions of the 4 cases are quite different from the sampling conditions. Acceleration ratio on the computational time is also surprisingly high (r = 14200).

Figure 17

Relative deviation of the prediction cases for a random porous medium

Figure 18

Comparison of flow fields for 4 typical prediction cases for a random porous medium: from left to right are p, u, v respectively; first row-FDM, second row-POD