^{1}

^{2}

^{*}

^{1}

The upwind scheme is very important in the numerical approximation of some problems such as the convection dominated problem, the two-phase flow problem, and so on. For the fractional flow formulation of the two-phase flow problem, the Penalty Discontinuous Galerkin (PDG) methods combined with the upwind scheme are usually used to solve the phase pressure equation. In this case, unless the upwind scheme is taken into consideration in the velocity reconstruction, the local mass balance cannot hold exactly. In this paper, we present a scheme of velocity reconstruction in some H(div) spaces with considering the upwind scheme totally. Furthermore, the different ways to calculate the nonlinear coefficients may have distinct and significant effects, which have been investigated by some authors. We propose a new algorithm to obtain a more effective and stable approximation of the coefficients under the consideration of the upwind scheme.

In the context of some fields, such as modeling and simulation of fluid flows in petroleum or groundwater reservoirs, the studies of processes of the simultaneous flow of two or more fluid phases within a porous medium are of great significance. In this paper, we consider the cases of two-phase flow where the fluids are immiscible. A large number of methods, which are based on the finite difference (FD), the finite volume (FV) or the finite element (FE) methods, have been developed to deal with the two-phase flow problem. As is well known, no matter which kind of numerical methods is used, the upwind scheme is of great significance in the approximation of some problems such as the convection dominated problem, the two-phase flow problem, and so on.

To achieve stable numerical computations in the simulation of two-phase flow problem, an accurate approximation of the flux is one of the most important and desirable ingredients. If we use Penalty Discontinuous Galerkin (PDG) methods to discretize the pressure equation, like in [_{1} space. A more stable and accurate reconstruction was developed in [

The different ways to calculate the nonlinear coefficients may have distinct and significant effects, which have been investigated by some authors. For the approximation of the coefficients, we extend the one used in [

The rest of the article is organized as follows: In addition to the introduction and conclusion, we divide the text of this document into four parts. Section 2 is the first part and consists of two subsections, in which we introduce governing equations of two-phase flow problem and the corresponding interface conditions in Subsections 2.1 and 2.2 respectively. The second part, Section 3, comes in four subsections. In Subsection 3.1, the upwind average approximations of coefficients are introduced. In Subsections 3.2 and 3.3, the PDG methods are used for the pressure equation and the velocity reconstruction is presented respectively. In Subsection 3.4, the PDG methods are used for the saturation equation. The third part, Section 4, consists of two subsections. In Subsections 4.1 and 4.2, we introduce all the possible projection schemes with respect to the velocity reconstruction and the scheme without any explicit projections. In the last part, Section 5, several numerical examples in two dimensions are provided.

We consider two immiscible incompressible fluids in porous media and there is no mass transfer between the phases. Various and alternative model equations for two-phase flow problem can be found in reference [

where the denotations and meanings of each coefficient and their relationships are defined as follows: D is the absolute permeability tensor and is discontinuous in heterogeneous media;

where

Some notations for the mesh are given below:

and the other is Neumann-Drichlet boundary condition used in [

The whole boundary of the porous medium domain

In order to test the barrier effect phenomenon of two phase flow, the nonlinear interface condition discussed in [

The process of the phenomenon is described briefly below. First, oil approaches the material interface but cannot penetrate it and begin to accumulate. In this case, only water pressure

Then, when more and more oils accumulate at the interface and the capillary pressure on the coarse side exceeds the entry pressure of the other side

We note that a critical point of saturation can be found when the capillary pressure on coarse side increases to the value equivalent to the threshold pressure on fine side. That is, deducing from

This point will be used to judge whether the nonwetting phase can or cannot penetrate the material interface. So the interface conditions can be rewritten in the form below. For capillary pressure,

and for wetting phase saturation,

Condition (20) is the same as that described in [

For the approximation of coefficients,we extend the one used in [

The approach described in [

where

where the upwind value of the side average on the interior edge is considered. The quantity

where the normal vector

In this section we apply the Penalty Discontinuous Galerkin (PDG) methods [

where

The pressure Equation (1) discretized by PDG reads as follows.

Indeed, the PDG methods are only applied to the wetting-phase pressure term, for the capillary pressure term a traditional DG method with the upwind scheme is used.

After solving the discrete pressure Equation (28), the total velocity will be reconstructed in the lowest-order Raviart-Thomas space

A proper reconstruction of velocity stems from the local mass conservation law as shown in the following description. Firstly, we recast the variational Equation (28) on element K into two parts as follows,

Combing (29) and (30), it is easily seen that the local mass is conserved,

where

Secondly, using (29) and (30) as the degree of freedom for some H(div) spaces, the total velocity will be obtained as some appropriate projections or interpolations in these spaces. In order to have a proper interpolation in

where

Let

Let

It is noted that the choice of DOFs for the

The spatial discretization of the saturation equation is similar to that of the pressure equation given in (28). The diffusion term of the saturation equation is discretized by the PDG methods, and the advective term is discretized by a traditional DG method with using the upwind scheme. An Euler scheme in time is used. The saturation Equation (2) equipped with Mixed-Neumann boundary conditions (8)-(11) could be written as:

The variational form in terms of Neumann-Dirichlet boundary conditions (12)-(15) reads:

where

The abbreviation DDG means that DG methods are used for both pressure and saturation equations. For a clear comparison in the numerical experiments, we list all the possible and feasible projections below. Firstly, we denote RT^{(1)} (or BDM^{(1)}) as the the velocity space projected into RT (or BDM) space by (29) and (30). Secondly, RT^{(2)} (or BDM^{(2)}) means the projection into RT (or BDM) space with considering the upwind scheme but without the penalty term,

Thirdly, for

At last, for

As indicated in introduction, for all DDG methods the velocity derived by the projection

In [

where the average of the total velocity is used in the interior edges. Although it doesn’t use any projections explicitly, the velocities constructed from (41) and (42) are some kind of implicit projections into

1) The variational form of the saturation equation doesn’t incorporate any additional penalties from the pressure equation.

2) The approximations of the coefficients are totally different.

In this section, we present some computer experiments to examine the proposed methods on two dimensional spaces. Both two boundary conditions with different types are used in the examination of all the methods. In tests 1 and 2 we consider the displacement of the non-wetting phase by the wetting phase, which is similar to the so called quarter-five spot problem introduced in [

If considering the mixed-Neumann type boundary (8)-(9) for the saturation equation, the following initial and boundary conditions are used:

When the Neumann-Dirichlet type boundary (12)-(13) is used for the saturation equation, the following initial and boundary conditions are considered:

The parameters including rock and fluid properties used in the simulation are summarized in

In test 1, we examine the property of the local mass conservation law. For this purpose, we solve the so-called quarter-five spot problem on a homogeneous medium and check the numerical local mass of the reconstructed velocity. All the projection methods discussed above will be used and compared. The domain used in the experiment is the square

Test 1 | |
---|---|

porosity | |

permeability [m^{2}] | |

viscosity[kg/(ms)] | |

residual saturation | |

Brooks-Corey | |

Test 2 | |

porosity | |

permeability [m^{2}] | |

viscosity[kg/(ms)] | |

residual saturation | |

Brooks-Corey | |

Test 3 | |

porosity | |

permeability [m^{2}] | |

viscosity[kg/(ms)] | |

residual saturation | |

Brooks-Corey |

explicit projections. Since there is no sink and source terms the exact local mass is zero on each elements, that is,

errors of the local mass conservation under the vector norms

The errors of the local mass conservation at selected times are listed in

In this test, we show the numerical solutions solved by our scheme with using projection

In the last test, we examine our scheme in a discontinuous media. We assume that the domain used here is initially fully water saturated and with the interfaces between two different sands, see

t = 40 s | t = 80 s | |||
---|---|---|---|---|

2.2913e−04 | 5.7110e−04 | 2.3767e−04 | 5.8819e−04 | |

2.8284e−09 | 1.6994e−08 | 2.2655e−09 | 1.6564e−08 | |

2.2117e−04 | 5.2564e−04 | 2.2860e−04 | 5.3194e−04 | |

2.0947e−05 | 7.6105e−05 | 1.8192e−05 | 6.8195e−05 | |

2.2087e−04 | 5.2553e−04 | 2.2808e−04 | 5.3135e−04 | |

2.6577e−09 | 1.6747e−08 | 2.1773e−09 | 1.5875e−08 | |

2.2117e−04 | 5.2564e−04 | 2.2860e−04 | 5.3194e−04 | |

2.1408e−05 | 7.8417e−05 | 1.7849e−05 | 6.9403e−05 | |

2.2086e−04 | 5.2552e−04 | 2.2808e−04 | 5.3135e−04 | |

2.0632e−09 | 1.6436e−08 | 2.5179e−09 | 1.7404e−06 | |

2.2933e−04 | 5.7035e−04 | 2.3798e−04 | 5.8939e−04 | |

2.1631e−05 | 7.9100e−05 | 1.7994e−05 | 7.0000e−05 | |

2.2913e−04 | 5.7110e−04 | 2.3766e−04 | 5.8819e−04 |

t = 120 s | t = 160 s | |||
---|---|---|---|---|

2.4309e−04 | 5.9636e−04 | 2.4690e−04 | 6.0363e−04 | |

2.4551e−09 | 1.7281e−08 | 2.3734e−09 | 1.8916e−08 | |

2.3367e−04 | 5.4152e−04 | 2.3644e−04 | 5.4898e−04 | |

2.1099e−05 | 7.1351e−05 | 1.4933e−05 | 6.7183e−05 | |

2.3305e−04 | 5.3781e−04 | 2.3576e−04 | 5.4316e−04 | |

2.4264e−09 | 1.7214e−08 | 2.4283e−09 | 1.7782e−08 | |

2.3368e−04 | 5.4152e−04 | 2.3644e−04 | 5.4898e−04 | |

2.1464e−05 | 7.2709e−05 | 1.5027e−05 | 6.7311e−05 | |

2.3305e−04 | 5.3780e−04 | 2.3576e−04 | 5.4316e−04 | |

2.0203e−09 | 1.6556e−08 | 2.4811e−09 | 1.8123e−08 | |

1.5926e−04 | 5.9819e−04 | 2.4740e−04 | 6.0948e−04 | |

2.1777e−05 | 7.3186e−05 | 1.5423e−05 | 6.6493e−05 | |

2.4309e−04 | 5.9636e−04 | 2.4691e−04 | 6.0363e−04 |

used in this test is

the reversed direction the oil immediately penetrate the interface, that is, the oil-trapped phenomenon will not happen if the oil flows from fine sand to coarse sand. The contours of wetting phase saturation in the discontinuous media at selected times are presented in

The velocities reconstructed from projections

The work is supported by the National Natural Science Foundation of China (No.11371288 and No.11401467).