Acceleration of Gas Flow Simulations in Dual-Continuum Porous Media Based on the Mass-Conservation POD Method

Abstract

Reduced-order modeling approaches for gas flow in dual-porosity dual-permeability porous media are studied based on the proper orthogonal decomposition (POD) method combined with Galerkin projection. The typical modeling approach for non-porous-medium liquid flow problems is not appropriate for this compressible gas flow in a dual-continuum porous media. The reason is that non-zero mass transfer for the dual-continuum system can be generated artificially via the typical POD projection, violating the mass-conservation nature and causing the failure of the POD modeling. A new POD modeling approach is proposed considering the mass conservation of the whole matrix fracture system. Computation can be accelerated as much as 720 times with high precision (reconstruction errors as slow as 7.69 × 10⁻⁴%~3.87% for the matrix and 8.27 × 10⁻⁴%~2.84% for the fracture).

1. Introduction

A dual-porosity, dual-permeability (dual-continuum) model is an important conceptual model for unconventional oil/gas flow in fractured reservoirs [1,2,3,4]. Its computational consumption is very large for solving the two coupled sets of partially differential equations for fluid flows in the matrix and the fracture [5,6], which is not endurable in engineering for analyzing the long-term transient characteristics of flow in fractured porous media. Due to the high computational cost, many analytical solutions have been proposed [7,8,9,10,11]. However, analytical solutions can only be obtained under much idealized assumptions, such as slight compressibility, infinite radial flow, homogenization etc., so that their applications are very limited to simplified cases of oil flow. For gas flow, the slight compressibility assumption is no longer held, leading to nonlinear equations which cannot be analytically solved. Therefore, gas flow simulation in dual-continuum porous media has the same computational problem as has also existed with other problems [12] and has its own problems such as strong nonlinearity. It is urgent to develop an efficient numerical method for gas flow in dual-porosity, dual-permeability porous media.

Proper orthogonal decomposition (POD) is an efficient numerical method to largely accelerate the computational speed via project governing equations onto low dimensional eigen spaces to establish a reduced-order model. It has been widely utilized for many non-porous-medium flow problems [13,14,15] and also proved to be efficient for some liquid flow cases in single-continuum porous media [16,17]. In reference [16], Ghommem et al. discussed a high-precision mode decomposition method for a time-dependent incompressible single-phase flow, but did not describe the acceleration performance of the method. In reference [17], a POD Galerkin model is proposed for an incompressible single-phase flow. Only four samples and two modes used in the model can predict hundreds of cases with high-precision and fast-computation. This research is all on incompressible single-phase flow. The flow domain is a single-continuum porous media. To the best knowledge of the authors, no work has been reported for a POD model of compressible single-phase gas flow in dual-continuum porous media, which has much higher computational cost than single-continuum porous media. Thus, POD has larger potential to accelerate the numerical computation of gas flow in dual-continuum porous media. In this study, we proposed a POD acceleration method for gas flow in fractured porous media based on a dual-porosity, dual-permeability model for the first time. The principle for the POD modeling of this complex flow system is discussed to guide the improvement of the POD modeling of a dual-continuum, porous-medium compressible flow system. The precision and acceleration of the improved POD model are examined.

2. Problems Arising from the Typical POD Modeling Approach

2.1. Model Derivation via the Typical Approach

For statement convenience, we only consider ideal gas flow in two-dimensional Cartesian coordinates in this paper. Darcy’s law is as follows:

$$u_M=-K_{xxM}\frac{\partial p_M}{\partial x}, v_M=-K_{yyM}\frac{\partial p_M}{\partial y}$$

(1)

$$u_F=-K_{xxF}\frac{\partial p_F}{\partial x}, v_F=-K_{yyF}\frac{\partial p_F}{\partial y}$$

(2)

Governed by the Darcy’s law, mass conservation equations for gas flow in dual-porosity dual-permeability porous media become:

$${\varphi }_M\frac{\partial p_M}{\partial t}=\frac{\partial }{\partial x}\left(K_{xxM}{p}_M\frac{\partial p_M}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyM}{p}_M\frac{\partial p_M}{\partial y}\right)+\frac{RT}{W}{q}_M-\alpha K_M{p}_M\left(p_M-p_F\right)$$

(3)

$${\varphi }_F\frac{\partial p_F}{\partial t}=\frac{\partial }{\partial x}\left(K_{xxF}{p}_F\frac{\partial p_F}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyF}{p}_F\frac{\partial p_F}{\partial y}\right)+\frac{RT}{W}{q}_F+\alpha K_M{p}_M\left(p_M-p_F\right)$$

(4)

where $\varphi$, $p$, $q$, $K_{xx}$, $K_{yy}$ are porosity, pressure, injection/production rate and two components of conductivity (diagonal permeability tensor divided by dynamic viscosity $\mu$). The subscripts M and F represent the matrix and fracture respectively. $K_M$ is the intrinsic conductivity of the matrix for the matrix–fracture interaction term. T is temperature. $W$ is the molecular weight of gas. R is the universal gas constant (8.3147295 J/(mol·K)). $\alpha$ is the shape factor of the fracture, taking the form proposed by Kazemi [18]:

$$\alpha =4\left(\frac{1}{L_x^2}+\frac{1}{L_y^2}\right)$$

(5)

where $L_x$ and $L_y$ are the lengths of the fracture spacing in the x and y directions respectively. The last terms of Equations (3) and (4) are interaction terms between matrix and fracture.

To establish the POD model for this flow system, pressure is assumed to be the following linear combination:

$$p={\displaystyle \sum _{n=1}^N{a}_n{\phi }_n}$$

(6)

where $p$ consists of $p_M$ and $p_F$ with corresponding spatial modes for matrix pressure (${\phi }_n^M$) and fracture pressure (${\phi }_n^F$), respectively, so that $p_M={\displaystyle \sum _{n=1}^N{a}_n{\phi }_n^M}$ and $p_F={\displaystyle \sum _{n=1}^N{a}_n{\phi }_n^F}$, $a_n$ are the temporal coefficients, and N is the number of POD modes. Substitute Equation (6) to Equations (3) and (4) to obtain the variable-separated Equations (7) and (8):

$$\begin{array}{c}{\varphi }_M{\displaystyle \sum _{n=1}^N{\phi }_n^M\frac{d{a}_n}{dt}}={\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l\left[\frac{\partial }{\partial x}\left(K_{xxM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial y}\right)\right]}}\\ +\frac{RT}{W}{q}_M-\alpha K_M{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l{\phi }_l^M\left({\phi }_n^M-{\phi }_n^F\right)}}\end{array}$$

(7)

$$\begin{array}{c}{\varphi }_F{\displaystyle \sum _{n=1}^N{\phi }_n^F\frac{d{a}_n}{dt}}={\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l\left[\frac{\partial }{\partial x}\left(K_{xxF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial y}\right)\right]}}\\ +\frac{RT}{W}{q}_F+\alpha K_M{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l{\phi }_l^M\left({\phi }_n^M-{\phi }_n^F\right)}}\end{array}$$

(8)

According to the typical procedure of POD modeling, one should project the matrix equation, Equation (7), onto ${\phi }_m^M$ and project the fracture equation, Equation (8), onto ${\phi }_m^F$, respectively (m = 1~N), and add these two projection equations together to solve the unique temporal coefficients $a_n$ (n = 1~N). The procedure will be explained as follows. First of all, the projections of Equations (7) and (8) are expressed as:

$$\begin{array}{l}{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^M{\varphi }_M{\displaystyle \sum _{n=1}^N{\phi }_n^M\frac{d{a}_n}{dt}}dxdy}}\\ ={\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^M{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l\left[\frac{\partial }{\partial x}\left(K_{xxM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial y}\right)\right]}}dxdy}}\\ +{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^M\frac{RT}{W}{q}_{M}dxdy}}-{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^M\alpha K_M{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l{\phi }_l^M\left({\phi }_n^M-{\phi }_n^F\right)}}dxdy}}\end{array}$$

(9)

$$\begin{array}{l}{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^F{\varphi }_F{\displaystyle \sum _{n=1}^N{\phi }_n^F\frac{d{a}_n}{dt}}dxdy}}\\ ={\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^F{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l\left[\frac{\partial }{\partial x}\left(K_{xxF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial y}\right)\right]}}dxdy}}\\ +{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^F\frac{RT}{W}{q}_{F}dxdy}}+{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^F\alpha K_M{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l{\phi }_l^M\left({\phi }_n^M-{\phi }_n^F\right)}}dxdy}}\end{array}$$

(10)

where l_x and l_y are the domain length in the x and y directions, respectively. Note that $\phi$ is only the spatial function and $a$ is only the temporal function. Using this property, the above two equations can be transformed to:

$$\begin{array}{cc}{\displaystyle \sum _{n=1}^N\frac{d{a}_n}{dt}{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\varphi }_M{\phi }_m^M{\phi }_n^{M}dxdy}}}& ={\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^M\left[\frac{\partial }{\partial x}\left(K_{xxM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial y}\right)\right]dxdy}}}}\\ & +\frac{RT}{W}{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^M{q}_{M}dxdy}}-\alpha {\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{K}_M{\phi }_m^M{\phi }_l^M\left({\phi }_n^M-{\phi }_n^F\right)dxdy}}}}\end{array}$$

(11)

$$\begin{array}{cc}{\displaystyle \sum _{n=1}^N\frac{d{a}_n}{dt}{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\varphi }_F{\phi }_m^F{\phi }_n^{F}dxdy}}}& ={\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^F\left[\frac{\partial }{\partial x}\left(K_{xxF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial y}\right)\right]dxdy}}}}\\ & +\frac{RT}{W}{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^F{q}_{F}dxdy}}+\alpha {\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l}}{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{K}_M{\phi }_m^F{\phi }_l^M\left({\phi }_n^M-{\phi }_n^F\right)dxdy}}\end{array}$$

(12)

The first integrations on the right hand side of Equations (11) and (12) can be transformed using integration by part:

$$\begin{array}{l}{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^M\left[\frac{\partial }{\partial x}\left(K_{xxM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyF}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial y}\right)\right]dxdy}}\\ ={\displaystyle {\int }_0^{ly}\left[{\left({\phi }_m^M{K}_{xxM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial x}\right)}_{lx}-{\left({\phi }_m^M{K}_{xxM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial x}\right)}_0\right]dy}\\ +{\displaystyle {\int }_0^{lx}\left[{\left({\phi }_m^M{K}_{yyM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial y}\right)}_{ly}-{\left({\phi }_m^M{K}_{yyM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial y}\right)}_0\right]dx}\\ -{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_l^M\left(K_{xxM}\frac{\partial {\phi }_n^M}{\partial x}\frac{\partial {\phi }_m^M}{\partial x}+K_{yyM}\frac{\partial {\phi }_n^M}{\partial y}\frac{\partial {\phi }_m^M}{\partial y}\right)dxdy}}\end{array}$$

(13)

$$\begin{array}{l}{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^F\left[\frac{\partial }{\partial x}\left(K_{xxF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial y}\right)\right]dxdy}}\\ ={\displaystyle {\int }_0^{ly}\left[{\left({\phi }_m^F{K}_{xxF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial x}\right)}_{lx}-{\left({\phi }_m^F{K}_{xxF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial x}\right)}_0\right]dy}\\ +{\displaystyle {\int }_0^{lx}\left[{\left({\phi }_m^F{K}_{yyF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial y}\right)}_{ly}-{\left({\phi }_m^F{K}_{yyF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial y}\right)}_0\right]dx}\\ -{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_l^F\left(K_{xxF}\frac{\partial {\phi }_n^F}{\partial x}\frac{\partial {\phi }_m^F}{\partial x}+K_{yyF}\frac{\partial {\phi }_n^F}{\partial y}\frac{\partial {\phi }_m^F}{\partial y}\right)dxdy}}\end{array}$$

(14)

By substituting these two expressions back to the first term of right hand side of Equations (11) and (12), one can obtain the following expressions:

$$\begin{array}{l}{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l}}{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^M\left[\frac{\partial }{\partial x}\left(K_{xxM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial y}\right)\right]dxdy}}\\ ={\displaystyle {\int }_0^{ly}\left[{\left({\phi }_m^M{K}_{xxM}{p}_M\frac{\partial p_M}{\partial x}\right)}_{lx}-{\left({\phi }_m^M{K}_{xxM}{p}_M\frac{\partial p_M}{\partial x}\right)}_0\right]dy}\\ +{\displaystyle {\int }_0^{lx}\left[{\left({\phi }_m^M{K}_{yyM}{p}_M\frac{\partial p_M}{\partial y}\right)}_{ly}-{\left({\phi }_m^M{K}_{yyM}{p}_M\frac{\partial p_M}{\partial y}\right)}_0\right]dx}\\ -{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_l^M\left(K_{xxM}\frac{\partial {\phi }_n^M}{\partial x}\frac{\partial {\phi }_m^M}{\partial x}+K_{yyM}\frac{\partial {\phi }_n^M}{\partial y}\frac{\partial {\phi }_m^M}{\partial y}\right)dxdy}}}}\end{array}$$

(15)

$$\begin{array}{l}{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_m^F\left[\frac{\partial }{\partial x}\left(K_{xxF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yyF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial y}\right)\right]dxdy}}}}\\ ={\displaystyle {\int }_0^{ly}\left[{\left({\phi }_m^F{K}_{xxF}{p}_F\frac{\partial p_F}{\partial x}\right)}_{lx}-{\left({\phi }_m^F{K}_{xxF}{p}_F\frac{\partial p_F}{\partial x}\right)}_0\right]dy}\\ +{\displaystyle {\int }_0^{lx}\left[{\left({\phi }_m^F{K}_{yyF}{p}_F\frac{\partial p_F}{\partial y}\right)}_{ly}-{\left({\phi }_m^F{K}_{yyF}{p}_F\frac{\partial p_F}{\partial y}\right)}_0\right]dx}\\ -{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{\phi }_l^F\left(K_{xxF}\frac{\partial {\phi }_n^F}{\partial x}\frac{\partial {\phi }_m^F}{\partial x}+K_{yyF}\frac{\partial {\phi }_n^F}{\partial y}\frac{\partial {\phi }_m^F}{\partial y}\right)dxdy}}}}\end{array}$$

(16)

For the Dirichlet boundary condition, the boundary pressure gradients in Equations (15) and (16) can be expressed as:

$$\begin{array}{l}{\frac{\partial p_k}{\partial x}|}_{lx}=\frac{{\left(p_k\right)}_{nx+1,j}-{\left(p_k\right)}_{nx,j}}{\Delta x/2}=\frac{{\left(p_k\right)}_{nx+1,j}-{\displaystyle \sum _{n=1}^N{a}_n}{\left({\phi }_n^k\right)}_{nx,j}}{\Delta x/2}\\ {\frac{\partial p_k}{\partial x}|}_0=\frac{{\left(p_k\right)}_{1,j}-{\left(p_k\right)}_{0,j}}{\Delta x/2}=\frac{{\displaystyle \sum _{n=1}^N{a}_n}{\left({\phi }_n^k\right)}_{1,j}-{\left(p_k\right)}_{0,j}}{\Delta x/2}\\ {\frac{\partial p_k}{\partial y}|}_{ly}=\frac{{\left(p_k\right)}_{i,ny+1}-{\left(p_k\right)}_{i,ny}}{\Delta y/2}=\frac{{\left(p_k\right)}_{i,ny+1}-{\displaystyle \sum _{n=1}^N{a}_n}{\left({\phi }_n^k\right)}_{i,ny}}{\Delta y/2}\\ {\frac{\partial p_k}{\partial y}|}_0=\frac{{\left(p_k\right)}_{i,1}-{\left(p_k\right)}_{i,0}}{\Delta y/2}=\frac{{\displaystyle \sum _{n=1}^N{a}_n}{\left({\phi }_n^k\right)}_{i,1}-{\left(p_k\right)}_{i,0}}{\Delta y/2}\end{array}$$

(17)

where i and j are the index of the grid points in the x and y directions, respectively; nx and ny are grid numbers; and $\Delta x$ and $\Delta y$ are corresponding grid spaces. K = M for the matrix and F for the fracture.

For the Neumann boundary condition, boundary pressure gradients can be expressed via Darcy’s law (Equations (1) and (2)). Substitute Equations (15) and (16) and the boundary conditions to Equations (11) and (12), and the projection equations can be expressed as follows:

Projection equation of the matrix:

$${\displaystyle \sum _{n=1}^N\frac{d{a}_n}{dt}}H{U}_{m,n}^M=H{S}_m^M+H{B}_m^M+{\displaystyle \sum _{n=1}^N{a}_n}\left(HB{V}_{m,n}^M-HB{P}_{m,n}^M\right)-{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l}}\left(H{D}_{m,n,l}^M+\alpha H{I}_{m,n,l}^M\right)$$

(18)

Projection equation of the fracture:

$${\displaystyle \sum _{n=1}^N\frac{d{a}_n}{dt}}H{U}_{m,n}^F=H{S}_m^F+H{B}_m^F+{\displaystyle \sum _{n=1}^N{a}_n}\left(HB{V}_{m,n}^F-HB{P}_{m,n}^F\right)-{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l}}\left(H{D}_{m,n,l}^F-\alpha H{I}_{m,n,l}^F\right)$$

(19)

where $H{U}_{m,n}^k={\displaystyle \sum _{j=1}^{ny}{\displaystyle \sum _{i=1}^{nx}{\left({\varphi }_k{\phi }_m^k{\phi }_n^k\right)}_{i,j}}}$, $H{S}_m^k=\frac{RT}{W}{\displaystyle \sum _{j=1}^{ny}{\displaystyle \sum _{i=1}^{nx}{\left({\phi }_m^k{q}_k\right)}_{i,j}}}$,

• $\begin{array}{ll}H{B}_m^k& =\frac{2}{{\left(\Delta x\right)}^2}{\displaystyle \sum _{j=1}^{ny}\left[{\left(K_{xxk}{\phi }_m^k{p}_k^2\right)}_{nx+1,j}Diri{X}_{nx+1,j}^k+{\left(K_{xxk}{\phi }_m^k{p}_k^2\right)}_{0,j}Diri{X}_{0,j}^k\right]}\\ & +\frac{2}{{\left(\Delta y\right)}^2}{\displaystyle \sum _{i=1}^{nx}\left[{\left(K_{yyk}{\phi }_m^k{p}_k^2\right)}_{i,ny+1}Diri{Y}_{i,ny+1}^k+{\left(K_{yyk}{\phi }_m^k{p}_k^2\right)}_{i,0}Diri{Y}_{i,0}^k\right]}\end{array}$,
• $\begin{array}{ll}HB{V}_{m,n}^k& =\frac{1}{\Delta x}{\displaystyle \sum _{j=1}^{ny}\left[u_{0,j}^k{\left({\phi }_m^k{\phi }_n^k\right)}_{0,j}\left(1-Diri{X}_{0,j}^k\right)-u_{nx,j}^k{\left({\phi }_m^k{\phi }_n^k\right)}_{nx+1,j}\left(1-Diri{X}_{nx+1,j}^k\right)\right]}\\ & +\frac{1}{\Delta y}{\displaystyle \sum _{i=1}^{nx}\left[v_{i,0}^k{\left({\phi }_m^k{\phi }_n^k\right)}_{i,0}\left(1-Diri{Y}_{i,0}^k\right)-v_{i,ny}^k{\left({\phi }_m^k{\phi }_n^k\right)}_{i,ny+1}\left(1-Diri{Y}_{i,ny+1}^k\right)\right]}\end{array}$,
• $\begin{array}{ll}HB{P}_{m,n}^k& =\frac{2}{{\left(\Delta x\right)}^2}{\displaystyle \sum _{j=1}^{ny}\left[{\left(K_{xxk}{\phi }_m^k{p}_k\right)}_{nx+1,j}{\left({\phi }_n^k\right)}_{nx,j}Diri{X}_{nx+1,j}^k+{\left(K_{xxk}{\phi }_m^k{p}_k\right)}_{0,j}{\left({\phi }_n^k\right)}_{1,j}Diri{X}_{0,j}^k\right]}\\ & +\frac{2}{{\left(\Delta y\right)}^2}{\displaystyle \sum _{i=1}^{nx}\left[{\left(K_{yyk}{\phi }_m^k{p}_k\right)}_{i,ny+1}{\left({\phi }_n^k\right)}_{i,ny}Diri{Y}_{i,ny+1}^k+{\left(K_{yyk}{\phi }_m^k{p}_k\right)}_{i,0}{\left({\phi }_n^k\right)}_{i,1}Diri{Y}_{i,0}^k\right]}\end{array}$,
• $H{D}_{m,n,l}^k={\displaystyle \sum _{j=1}^{ny}{\displaystyle \sum _{i=1}^{nx}{\left[{\phi }_l^k\left(K_{xxk}\frac{\partial {\phi }_n^k}{\partial x}\frac{\partial {\phi }_m^k}{\partial x}+K_{yyk}\frac{\partial {\phi }_n^k}{\partial y}\frac{\partial {\phi }_m^k}{\partial y}\right)\right]}_{i,j}}}$,
• $H{I}_{m,n,l}^M={\displaystyle \sum _{j=1}^{ny}{\displaystyle \sum _{i=1}^{nx}{\left[K_M{\phi }_m^M{\phi }_l^M\left({\phi }_n^M-{\phi }_n^F\right)\right]}_{i,j}}}$, $H{I}_{m,n,l}^F={\displaystyle \sum _{j=1}^{ny}{\displaystyle \sum _{i=1}^{nx}{\left[K_M{\phi }_m^F{\phi }_l^M\left({\phi }_n^M-{\phi }_n^F\right)\right]}_{i,j}}}$.

k = M for the matrix and F for the fracture. Diri is 1 for the Dirichlet boundary and 0 for the Neumann boundary. To uniquely solve the value of the POD coefficients, Equations (18) and (19) should be added to obtain the final POD projection equation of the dual-porosity, dual-permeability system:

$${\displaystyle \sum _{n=1}^N\frac{d{a}_n}{dt}}H{U}_{m,n}=H{B}_m+H{S}_m+{\displaystyle \sum _{n=1}^N{a}_n}\left(HB{V}_{m,n}-HB{P}_{m,n}\right)+{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n{a}_l\left(\alpha H{I}_{m,n,l}-H{D}_{m,n,l}\right)}}$$

(20)

where $H{U}_{m,n}=H{U}_{m,n}^M+H{U}_{m,n}^F$, $H{B}_m=H{B}_m^M+H{B}_m^F$, $H{S}_m=H{S}_m^M+H{S}_m^F$, $HB{V}_{m,n}=HB{V}_{m,n}^M+HB{V}_{m,n}^F$, $HB{P}_{m,n}=HB{P}_{m,n}^M+HB{P}_{m,n}^F$, $H{D}_{m,n,l}=H{D}_{m,n,l}^M+H{D}_{m,n,l}^F$, $H{I}_{m,n,l}=H{I}_{m,n,l}^F-H{I}_{m,n,l}^M$.

2.2. Numerical Methods and Parameters

Equation (20) is the POD model for the gas flow in dual-porosity, dual-permeability porous media using the typical modeling approach. The temporal advancement of coefficients $a_n$ can be expressed by the following linear algebraic equations:

$${\displaystyle \sum _{n=1}^N{A}_{m,n}{a}_n\left(t+\Delta t\right)}=B_m$$

(21)

where $B_m=H{B}_m+H{S}_m+{\displaystyle \sum _{n=1}^N{a}_n\left(t\right)}\left(A_{m,n}+HB{V}_{m,n}-HB{P}_{m,n}\right)+{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{l=1}^N{a}_n\left(t\right)a_l\left(t\right)\left(\alpha H{I}_{m,n,l}-H{D}_{m,n,l}\right)}}$, $A_{m,n}=H{U}_{m,n}/\Delta t$. Starting with initial coefficients $a_n\left(t_0\right)$, the coefficients for every time step $a_n\left(t+\Delta t\right)$ can be calculated via Equation (21). $a_n\left(t_0\right)$ comes from the inner product of initial pressure onto POD modes:

$$a_n\left(t_0\right)=\left[\begin{array}{l}{p}_M\\ p_F\end{array}\right]\cdot \left[\begin{array}{l}{\phi }_n^M\\ {\phi }_n^F\end{array}\right]$$

(22)

The symbol “·” means inner product. Equation (22) exists because of an important property of POD modes, namely that they are unitary vectors orthogonal to each other, i.e., $\left[\begin{array}{l}{\phi }_n^M\\ {\phi }_n^F\end{array}\right]\cdot \left[\begin{array}{l}{\phi }_m^M\\ {\phi }_m^F\end{array}\right]=\{\begin{array}{l}1n=m\\ 0n\ne m\end{array}$. To calculate the temporal coefficients of the POD model in Equation (21), $A_{m,n}$ and $B_m$ should be firstly calculated by POD modes. The modes are generated from the following sampling and eigenvalue decomposition process:

(1)

Through the numerical computation of Equations (3) and (4), a sample matrix of pressure can be collected at different moments as:

$$S=\left[\begin{array}{ccc}{p}^M\left(t_1\right)& \cdots & p^M\left(t_{Ns}\right)\\ p^F\left(t_1\right)& \cdots & p^F\left(t_{Ns}\right)\end{array}\right]$$

(23)

where Ns is the number of samples. To ensure stability and accuracy at a large time step, we used the semi-implicit finite difference method (FDM) to compute Equations (3) and (4). With the samples, a kernel for the eigenvalue decomposition can be calculated as follows:

$$C=\frac{1}{Ns}{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}{S}^{T}Sdxdy}}$$

(24)

(2)

Take the eigenvalue decomposition for the kernel to obtain eigenvalues and eigenvectors:

$$CV=\left[\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_{Ns}\end{array}\right]V$$

(25)

where ${\lambda }_1>{\lambda }_2>\dots >{\lambda }_{Ns}$ are the eigenvalues; and $V$ is the corresponding eigenvector matrix.

(3)

Calculate the POD modes using the eigenvectors and samples:

$$\left[\begin{array}{l}{\phi }_1^M{\phi }_2^M\cdots {\phi }_{Ns}^M\\ {\phi }_1^F{\phi }_2^F\cdots {\phi }_{Ns}^F\end{array}\right]=SV/\Vert SV\Vert$$

(26)

where “$\Vert \Vert$” means L₂ norm.

A numerical case is designed for the gas flow in the dual-porosity, dual-permeability porous media. The computational domain and boundary conditions are shown in Figure 1. Permeability is assumed to be isotropic with spatial distribution (Figure 2). The gas property was chosen to be methane as it is the main composition of natural gas. Uniform mesh was used. Sampling was made every 10 time steps. Other parameters are listed in

Table 1. Additional parameters for computation.

Parameter	Value	Unit
${\varphi }_M$	0.5	/
${\varphi }_F$	0.02	/
$p_M\left(t_0\right)$	1,013,250	Pa
$p_F\left(t_0\right)$	1,013,250	Pa
$p_1$	2,026,500	Pa
$p_2$	101,325	Pa
$q_M$	0	Kg/(m3·s)
$q_F$	0	Kg/(m3·s)
$K_M$	8.9177127 × 10−11	m2/(Pa·s)
W	16 × 10−3	Kg/mol
R	8.3147295	J/(mol·K)
T	298	K
$\mu$	11.067 × 10−6	Pa·s
nx	100	/
ny	100	/
Ns	2433	/
lx	100	m
ly	100	m
Lx	0.2	m
Ly	0.2	m
$\Delta x$	1	m
$\Delta y$	1	m
$\Delta t$	1296	s
Simulation time scope	365	days

.

2.3. Problem Analyses

The typical POD model for gas flow in dual-porosity, dual-permeability porous media (Equations (6) and (21)) is computed using the above parameters. To examine the deviation quantitatively, the relative error for pressures in the matrix and the fracture are defined as follows:

$${\epsilon }_M=\Vert p_{POD}^M-p_{FDM}^M\Vert /\Vert p_{FDM}^M\Vert \times 100\%$$

(27)

$${\epsilon }_F=\Vert p_{POD}^F-p_{FDM}^F\Vert /\Vert p_{FDM}^F\Vert \times 100\%$$

(28)

where $p_{FDM}^M$ and $p_{FDM}^F$ are matrix pressure and fracture pressure computed by FDM in Equations (3) and (4), $p_{POD}^M$ and $p_{POD}^F$ are corresponding pressures computed by the POD model. The details of the FDM algorithm are as follows. The temporal advancement for the whole equation system is a semi-implicit method. The diffusion terms in Equations (3) and (4) are discretized using the second-order cell-centered finite difference scheme. A staggered grid, where velocity components are placed on the cell surfaces of the pressure grid, is used to ensure the computational stability. The discrete equations via these FDM algorithms can be easily solved using the linear solver of algebraic equations.

It is unusual that the typical POD model can only be computed when the top one POD modes of matrix and fracture (${\phi }_1^M$ and ${\phi }_1^F$) are used. If more POD modes are used in the POD model, the computation will be broke up. For the POD results using the top one POD modes, the relative errors are very large all over the whole transient process (18.89~38.67% for the fracture pressure and 16.26~38.98% for the matrix pressure), as shown in Figure 3. It is further confirmed in Figure 4 and Figure 5 that the reconstructed pressure fields for both the fracture and matrix cannot reproduce pressure fields accurately at the two moments: t₁ and t₂. Thus, the typical POD model has quite low robustness and low precision, and it should be improved.

Before the improvement of the above POD model, the reason for the low robustness and low precision should be revealed firstly to clarify the direction of improvement. We revisited the governing equations (Equations (3) and (4)) of the gas flow in dual-porosity, dual-permeability porous media and realized that these two equations actually reflect the mass transfer of gas between the matrix and fracture. The term $-\alpha K_M{p}_M\left(p_M-p_F\right)$ in Equation (3) represents the mass leaving the matrix. The term $\alpha K_M{p}_M\left(p_M-p_F\right)$ in Equation (4) represents the mass entering the fracture. The mass of gas in the matrix and fracture changes every second, but the total mass for the whole system does not change. Mass is only transferred from the matrix part to the fracture part so that the net mass transfer should be zero. Therefore, the projection of the net mass transfer should also be zero. However, in the typical POD modeling approach, the transfer term in Equation (3) is projected onto ${\phi }^M$ while the transfer term in Equation (4) is projected onto ${\phi }^F$ so that the projection of the net mass transfer $I=-{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}\alpha K_M{p}_M\left(p_M-p_F\right){\phi }^{M}dxdy}}+{\displaystyle {\int }_0^{ly}{\displaystyle {\int }_0^{lx}\alpha K_M{p}_M\left(p_M-p_F\right){\phi }^{F}dxdy}}$ could not be zero. This non-zero term is generated from the projection process but not from its physical nature. Thus, the typical POD modeling is correct in mathematics but incorrect in physics because it induces artificial mass transfer of the whole system, which does not exist in nature.

To verify the above theoretical analyses, we examine the value of $I$ and compare its magnitude with other terms. It is obvious that $I=\alpha H{I}_{m,n,l}\Delta x\Delta y$ according to the definition of integration. Therefore, $\alpha H{I}_{m,n,l}$ in Equation (20) can represent the projection of the net mass transfer $I$. For convenience of comparison, we define a ratio of the net mass transfer to the diffusion term:

$$r=\alpha {\displaystyle \sum _{l=1}^N{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{m=1}^{N}H{I}_{m,n,l}}}}/{\displaystyle \sum _{l=1}^N{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{m=1}^{N}H{D}_{m,n,l}}}}$$

(29)

The ratio is −3.62 × 10², −1.98 × 10³ and 4.85 × 10³ for N = 1, 2, 3 respectively in

Table 2. Magnitude comparison between the net mass transfer and the diffusion term.

N	1	2	3
$r$	−3.62 × 102	−1.98 × 103	4.85 × 103

. This means that the effect of the matrix–fracture interaction is hundreds or thousands of times larger than the effect of the diffusion, so that the behavior is indeed controlled by the matrix–fracture interaction in the typical POD modeling. However, this interaction effect is artificially generated from the projection of the transfer term $-\alpha K_M{p}_M\left(p_M-p_F\right)$ and $\alpha K_M{p}_M\left(p_M-p_F\right)$ separately in the matrix equation and fracture equation. Thus, the large unphysical mass transfer term causes the computation of the typical POD model to have very low precision (N = 1) or even to be broken up (N > 1). It is confirmed by the above comparison that the previous theoretical analyses are correct quantitatively. To obtain a reliable and accurate POD model for the dual-porosity, dual-permeability system, the artificial mass transfer should be avoided.

3. A New POD Modeling Approach Based on System Mass Conservation

3.1. Establishment of the New POD Model

According to the theoretical analysis and numerical comparison in Section 2.3, we realized that the separate projections of the interaction terms in matrix Equation (3) and fracture Equation (4) in the typical POD modeling approach introduces artificial mass transfer between matrix and fracture which does not exist in nature. Although it is correct in mathematics, the typical approach violates the mass conservation of the whole dual-continuum system in physics. A new POD modeling approach should be developed obeying both mathematics and physics.

From the analyses in Section 2.3, we recall that the mass is only exchanged between the matrix part and the fracture part, but the net mass transfer for the whole system is zero. This kind of mass conservation is automatically ensured in FDM but not automatically ensured in the typical POD modeling approach. To guarantee the mass conservation in the POD modeling, we propose a new approach in which the summation of Equations (3) and (4) is made before their projection. After the summation, the mass conservation equation of the whole system is:

$$\begin{array}{ll}\frac{\partial \left({\varphi }_M{p}_M+{\varphi }_F{p}_F\right)}{\partial t}& =\frac{\partial }{\partial x}\left(K_{xxM}{p}_M\frac{\partial p_M}{\partial x}+K_{xxF}{p}_F\frac{\partial p_F}{\partial x}\right)\\ & +\frac{\partial }{\partial y}\left(K_{yyM}{p}_M\frac{\partial p_M}{\partial y}+K_{yyF}{p}_F\frac{\partial p_F}{\partial y}\right)+\frac{RT}{W}\left(q_M+q_F\right)\end{array}$$

(30)

This equation can be projected onto either ${\phi }_m^M$ or ${\phi }_m^F$. Using a general symbol ${\phi }_m^{*}$ (* = M or F representing the matrix or fracture respectively) in a similar projection process to that in Section 2.1, the projection equation of Equation (30) can be expressed as follow:

$${\displaystyle \sum _{n=1}^M\frac{d{a}_n}{dt}}H{U}_{m,n}^{new}=H{S}_m^{new}+H{B}_m^{new}+{\displaystyle \sum _{n=1}^M{a}_n}\left(HB{V}_{m,n}^{new}-HB{P}_{m,n}^{new}\right)-{\displaystyle \sum _{n=1}^M{\displaystyle \sum _{l=1}^M{a}_n{a}_{l}H{D}_{m,n,l}^{new}}}$$

(31)

where $H{U}_{m,n}^{new}={\displaystyle \sum _{j=1}^{ny}{\displaystyle \sum _{i=1}^{nx}{\left({\varphi }_M{\phi }_m^{*}{\phi }_n^M+{\varphi }_F{\phi }_m^{*}{\phi }_n^F\right)}_{i,j}}}$, $H{S}_m^{new}=\frac{RT}{W}{\displaystyle \sum _{j=1}^{ny}{\displaystyle \sum _{i=1}^{nx}{\left[{\phi }_m^{*}\left(q_M+q_F\right)\right]}_{i,j}}}$,

• $\begin{array}{cc}H{B}_m^{new}& =\frac{2}{{\left(\Delta x\right)}^2}{\displaystyle \sum _{j=1}^{ny}\left\{\begin{array}{l}{\left({\phi }_m^{*}\right)}_{nx+1,j}\left[{\left(K_{xxM}{p}_M^2\right)}_{nx+1,j}Diri{X}_{nx+1,j}^M+{\left(K_{xxF}{p}_F^2\right)}_{nx+1,j}Diri{X}_{nx+1,j}^F\right]\\ +{\left({\phi }_m^{*}\right)}_{0,j}\left[{\left(K_{xxM}{p}_M^2\right)}_{0,j}Diri{X}_{0,j}^M+{\left(K_{xxF}{p}_F^2\right)}_{0,j}Diri{X}_{0,j}^F\right]\end{array}\right\}}\\ & +\frac{2}{{\left(\Delta y\right)}^2}{\displaystyle \sum _{i=1}^{nx}\left\{\begin{array}{l}{\left({\phi }_m^{*}\right)}_{i,ny+1}\left[{\left(K_{yyM}{p}_M^2\right)}_{i,ny+1}Diri{Y}_{i,ny+1}^M+{\left(K_{yyF}{p}_F^2\right)}_{i,ny+1}Diri{Y}_{i,ny+1}^F\right]\\ +{\left({\phi }_m^{*}\right)}_{i,0}\left[{\left(K_{yyM}{p}_M^2\right)}_{i,0}Diri{Y}_{i,0}^M+{\left(K_{yyF}{p}_F^2\right)}_{i,0}Diri{Y}_{i,0}^F\right]\end{array}\right\}}\end{array}$,
• $\begin{array}{cc}HB{V}_{m,n}^{new}& =\frac{1}{\Delta x}{\displaystyle \sum _{j=1}^{ny}\left\{\begin{array}{l}{\left({\phi }_m^{*}\right)}_{0,j}\left[u_{0,j}^M{\left({\phi }_n^M\right)}_{0,j}\left(1-Diri{X}_{0,j}^M\right)+u_{0,j}^F{\left({\phi }_n^F\right)}_{0,j}\left(1-Diri{X}_{0,j}^F\right)\right]\\ -{\left({\phi }_m^{*}\right)}_{nx+1,j}\left[u_{nx,j}^M{\left({\phi }_n^M\right)}_{nx+1,j}\left(1-Diri{X}_{nx+1,j}^M\right)+u_{nx,j}^F{\left({\phi }_n^F\right)}_{nx+1,j}\left(1-Diri{X}_{nx+1,j}^F\right)\right]\end{array}\right\}}\\ & +\frac{1}{\Delta y}{\displaystyle \sum _{i=1}^{nx}\left\{\begin{array}{l}{\left({\phi }_m^{*}\right)}_{i,0}\left[v_{i,0}^M{\left({\phi }_n^M\right)}_{i,0}\left(1-Diri{Y}_{i,0}^M\right)+v_{i,0}^F{\left({\phi }_n^F\right)}_{i,0}\left(1-Diri{Y}_{i,0}^F\right)\right]\\ -{\left({\phi }_m^{*}\right)}_{i,ny+1}\left[v_{i,ny}^M{\left({\phi }_n^M\right)}_{i,ny+1}\left(1-Diri{Y}_{i,ny+1}^M\right)+v_{i,ny}^F{\left({\phi }_n^F\right)}_{i,ny+1}\left(1-Diri{Y}_{i,ny+1}^F\right)\right]\end{array}\right\}}\end{array}$,
• $\begin{array}{cc}H{B}_m^{new}& =\frac{2}{({\Delta x}^2}{\displaystyle \sum _{j=1}^{ny}\left[\begin{array}{l}{\left({\phi }_m^{*}\right)}_{nx+1,j}\left[{\left(K_{xxM}{p}_M\right)}_{nx+1,j}{\left({\phi }_n^M\right)}_{nx,j}Diri{X}_{nx+1,j}^M+{\left(K_{xxF}{p}_F\right)}_{nx+1,j}{\left({\phi }_n^F\right)}_{nx,j}Diri{X}_{nx+1,j}^F\right]\\ +{\left({\phi }_m^{*}\right)}_{0,j}\left[{\left(K_{xxM}{p}_M\right)}_{0,j}{\left({\phi }_n^M\right)}_{1,j}Diri{X}_{0,j}^M+{\left(K_{xxF}pF\right)}_{0,j}{\left({\phi }_n^F\right)}_{1,j}Diri{X}_{0,j}^F\right]\end{array}\right]}\\ & +\frac{2}{{\left(\Delta y\right)}^2}{\displaystyle \sum _{i=1}^{nx}\left\{\begin{array}{l}{\left({\phi }_m^{*}\right)}_{i,ny+1}\left[{\left(K_{yyM}{p}_M\right)}_{i,ny+1}{\left({\phi }_n^M\right)}_{i,ny}Diri{Y}_{i,ny+1}^M+{\left(K_{yyF}{p}_F\right)}_{i,ny+1}{\left({\phi }_n^F\right)}_{i,1}Diri{Y}_{i,ny+1}^F\right]\\ +{\left({\phi }_m^{*}\right)}_{i,0}\left[{\left(K_{yyM}{p}_M\right)}_{i,0}{\left({\phi }_n^M\right)}_{i,1}Diri{Y}_{i,0}^M+{\left(K_{yyF}{p}_F\right)}_{i,0}{\left({\phi }_n^F\right)}_{i,1}Diri{Y}_{i,0}^F\right]\end{array}\right\}}\end{array}$,
• $H{D}_{m,n,l}^{new}={\displaystyle \sum _{j=1}^{ny}{\displaystyle \sum _{i=1}^{nx}{\left[\frac{\partial {\phi }_m^{*}}{\partial x}\left(K_{xxM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial x}+K_{xxF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial x}\right)+\frac{\partial {\phi }_m^{*}}{\partial y}\left(K_{yyM}{\phi }_l^M\frac{\partial {\phi }_n^M}{\partial y}+K_{yyF}{\phi }_l^F\frac{\partial {\phi }_n^F}{\partial y}\right)\right]}_{i,j}}}$,

Equation (31) does not contain the projection of interaction terms, such as $H{I}_{m,n,l}$ in Equation (20), so that no artificial mass transfer exists. Thus, it represents the mass conservation of the whole system in physical nature. The derivation also obeys mathematical laws strictly.

3.2. Model Verification

A new POD model based on the mass conservation of the whole system is established in Section 3.1 theoretically. Its precision and computational speed should be verified by the same numerical case in Section 2.2. First of all, the projection operator ${\phi }_m^M$ or ${\phi }_m^F$ needs to be discussed to determine the exact model expression. For convenience of comparison, the relative errors in Equations (27) and (28) are averaged in time. The time-averaged errors of the matrix pressure and the fracture pressure are listed in

Table 3. Time-averaged error of matrix pressure for different projection operators.

$\overline{{\mathit{\epsilon }}_{\mathit{M}}}$ (%)	N = 1	N = 2	N = 3	N = 4	N = 5	N = 6
Project onto ${\phi }_m^M$	2.7527	1.3319	2.5884	0.9123	0.8826	0.9110
Project onto ${\phi }_m^F$	2.7469	1.8863	1.9840	1.1248	N/A	N/A

and

Table 4. Time-averaged error of fracture pressure for different projection operators.

$\overline{{\mathit{\epsilon }}_{\mathit{F}}}$ (%)	N = 1	N = 2	N = 3	N = 4	N = 5	N = 6
Project onto ${\phi }_m^M$	1.9049	0.9131	1.8079	0.6702	0.7182	0.7062
Project onto ${\phi }_m^F$	1.8990	1.3435	1.3795	0.7912	N/A	N/A

. The time-averaged errors have slight differences between the projections onto ${\phi }_m^M$ and ${\phi }_m^F$ when the numbers of the POD modes are from 1 to 4. However, the POD model based on the projection of ${\phi }_m^F$ cannot obtain results for $N\ge 5$ due to the break-up of the computation. This is because the gas in the matrix flows much more slowly than the gas in the fracture, i.e., $p_F$ changes much faster than $p_M$ in the same time scope. Thus, the values of $p_F$ in samples are in a very small range, leading the base ${\phi }^F$ to contain more numerical errors than ${\phi }^M$ (the numerical errors are generated from the POD decomposition). Thus, the blow up using more modes ($N\ge 5$) indicates that only the first four modes of the base ${\phi }^F$ contain information on fracture pressure while the following modes ($N\ge 5$) are all numerical errors with no physical meaning. Thus, the new POD modeling approach can only use the projection of matrix mode, i.e., the superscript “*” in Equation (31), and is “M”.

After the determination of the projection onto ${\phi }_m^M$, the number of POD modes should also be determined because it affects the precision of the POD model. The importance of the mode number can be represented by the energy spectrum:

$${\eta }_n={\displaystyle \sum _{j=1}^n{\lambda }_j}/{\displaystyle \sum _{i=1}^N{\lambda }_i}\times 100\%$$

(32)

where ${\eta }_n$ is the cumulative energy contribution of the top n POD modes (${\varphi }_1~{\varphi }_n$) to total energy. Once ${\eta }_n$ achieves 100%, it indicates that all information has been captured by these top n POD modes. As shown in Figure 6, the cumulative energy contribution achieves 100% from the 5th or 6th mode so that the POD modes at least after the 6th mode should be neglected. The exact mode number needs an examination of the relative error for the whole transient process.

The temporal evolutions of relative errors using different POD modes are shown in Figure 7. Relative errors for the matrix and fracture pressures fluctuate at the initial stage and then rapidly decrease to a stable value at any number of POD modes. They decrease with increasing number of POD modes and become stable after five modes. The addition of the 6th mode causes the error to increase slightly (Figure 7f), indicating that the 6th mode has no positive contribution to precision promotion. Thus, the top five modes should be adopted in the new POD model.

To further examine the local precision of the new POD model, it is necessary to compare the local flow fields at two typical moments (0.15 days and 365 days). At 0.15 days, the reconstruction errors of the POD are maximum. At 365 days, the reconstruction errors of the POD stay low and stable. As shown in Figure 8 and Figure 9, the reconstructed pressure fields and velocity fields in the fracture and matrix from the new POD model always agree well with those from FDM. The local distributions of the reconstructed flow fields by POD only have small deviations with the flow fields computed by FDM, regardless of whether the reconstruction errors are large (${\epsilon }^F=2.84\%$, ${\epsilon }^M=3.87\%$ at 0.15 days) or small (${\epsilon }^F=0.83\%$, ${\epsilon }^M=1.02\%$ at 365 days). Thus, it is confirmed that the new POD model can reproduce all the transient characteristics of pressure and velocity for gas flow in the dual-porosity, dual-permeability system with high precision.

The aim of POD modeling is not only high-precision reconstruction, but also a large-acceleration of computation. Hence, the verification of computational speed is important. The comparison of code running time, i.e., CPU time, is listed in

Table 5. Acceleration ability of the new POD model.

	FDM	New POD Model
CPU time	3600 s	5 s

. For the same period of simulation time (365 days in this case), simulation via FDM needs one hour (i.e., 3600 s), but the simulation via the new POD model only needs 5 s. The acceleration of computational speed is as high as 720 times, indicating that the computational time is largely saved. The new POD modeling approach is a good way to decrease CPU time for numerical computation for gas flow in dual-porosity, dual-permeability porous media.

4. Conclusions

POD modeling for gas flow in dual-porosity, dual-permeability porous media was studied for the purpose of the acceleration of computation with high-precision for potential engineering applications. Some conclusions related to the principle of POD modeling for this type of flow can be made as follows:

(1)

For dual-porosity, dual-permeability porous media, the typical method should be avoided to project the matrix equation and fracture equation separately. Otherwise, an artificial mass transfer term, which is 10³~10² times larger than the diffusion term, will be generated to cause the failure of the POD modeling, because it violates the mass conservation of the whole system.

(2)

A mass conservation POD modeling method is proposed to ensure that no artificial mass transfer is generated by the POD projection process. Original governing equations should be projected onto the POD modes of matrix pressure to maintain a robust POD model.

(3)

The new POD model obeying the mass-conservation nature of the whole system can promote computational speed as much as 720 times under high precision: ${\epsilon }_{\max }^F=2.84\%$, ${\epsilon }_{\max }^M=3.87\%$; ${\epsilon }_{\min }^F=8.27\times {10}^{-4}\%$, ${\epsilon }_{\min }^M=7.69\times {10}^{-4}\%$; ${\epsilon }_{\mathrm{stable}}^F=0.83\%$, ${\epsilon }_{\mathrm{stable}}^M=1.02\%$.