Velocity Projection with Upwind Scheme Based on the Discontinuous Galerkin Methods for the Two Phase Flow Problem

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.


Introduction
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]- [3], both the pressure and saturation equations will be discritized by the PDG methods, and a process of reconstruction of the velocity needs to be done after the pressure equation is solved.In [2], an average total velocity was post processed by substituting the piecewise constants of pressure gradient and saturation gradient into the velocity-pressure expression directly.Actually such reconstructed velocity, on some level, belongs to the lowest order Raviart-Thomas finite element space.In [4], a post-processed total velocity is reconstructed in the Brezzi-Douglas-Marini (BDM) finite element spaces.But it needs that the degree of the polynomial is more than one, that is, using the linear approximation in DG method is not enough to reconstruct a velocity in BDM 1 space.A more stable and accurate reconstruction was developed in [1], in which the velocity reconstructed from the piecewise linear pressures could even belong to the first-order Raviart-Thomas finite element space.However, all the reconstructions mentioned above didn't consider the upwind scheme, which was basically used in the discretization of the equations.The property of the local mass conservation is crucial in porous media flow and transport problems.The upwind scheme has direct effect on the local mass conservation of the reconstructed velocity.We found that unless the upwind scheme and penalty terms which are used in the discretization of the two-phase flow problem are considered together into the velocity reconstruction, the error of the local mass conservation cannot reach a satisfactory level.In this paper, we present a scheme of velocity reconstruction in some H(div) spaces [5] with considering the upwind scheme totally.
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 [6] to that each coefficient in element K is evaluated as the average of the upwind value on K ∂ .This improves the stability of the numerical scheme even when an explicit scheme is used.In contrast with the explicit scheme described in [2], our explicit PDG scheme with this special approximation of coefficients can not only get rid of the extra penalties from the pressure equation but also have a robust performance in the heterogeneous media.
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.

Mathematical Formulation
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 [7].Here we use the phase formulation for which the primary variables are wetting phase pressure and saturation ( w p and w S ), and in the absence of gravity and sink/source term we have: ( ) 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; φ denotes the porosity of the me- dium; c p is the capillary pressure; n λ , w λ and t λ are wetting, nonwetting and total mobility respectively; w f is the fractional flow; ∇ and ∇ ⋅ are the gradient operator and divergence operator respectively.We will use the Brooks-Corey model [8] throughout this paper, in which some of these coefficients are non-linear functions defined below: ( ) , , , , , , on , , on , 0, on , and the other is Neumann-Drichlet boundary condition used in [9], ( ) , on , , on , , on , , on .
The whole boundary of the porous medium domain ∂Ω is divided into three mutually disjoint parts: the in-

Interface Conditions
In order to test the barrier effect phenomenon of two phase flow, the nonlinear interface condition discussed in [1] [6] [10]- [13] will be introduced here.Following [9], we assume an initially fully water saturated domain ( ) with an interface J Γ between two different sands, and the oil is injected from the inflow part of boundary in Γ , see Figure 1.In addition, we assume that I Ω stands the coarse sand and II Ω is the fine sand.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 w p is continuous on the interface, capillary pressure c p and saturation w S are discontinuous and satisfy: .
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 ( ) , the oils begin to penetrate and enter the fine sand.At this time, both w p and c p are continuous, but saturation w S is still discontinuous and satisfies: ( ) 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 we have, ( . 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,  Condition (20) is the same as that described in [1] except that the wetting phase (instead of the nonwetting phase) is used as the saturation variable.Moreover, (19) is only written for the capillary pressure and not for the wetting phase pressure, since the variable w p is always continuous in the problem discussed.Noting that if the sub-domain I Ω has a finer texture than II Ω , all the relationship above can be treated in a similar manner with superscript I and II reversed.

Approximation of Coefficients
For the approximation of coefficients,we extend the one used in [6].Let σ denote any coefficients waiting for some proper approximations.Firstly we recall the original way to approximate the coefficients, 1 , 1 .
The approach described in [6] is, (

TS
denotes the mean water saturation on the edge e of the element K, see [6] for more details.Now it is extended to the following one, where the upwind value of the side average on the interior edge is considered.The quantity σ ↑ is called the upwind flux which is done with respect to the normal component of the total velocity t u , such that for all e K K .
where the normal vector e n points from K + to K − .Throughout this paper, all the coefficients on element K and edge e are calculated by the upwind averaged constant and the integral average constant which are described in (23).

Pressure Approximation with PDG
In this section we apply the Penalty Discontinuous Galerkin (PDG) methods [14] such as Nonsymmetric Interior Penalty Galerkin (NIPG) to the pressure Equation (1).Some notations for DG methods are defined: where v ± are the restrictions of v on two adjacent elements K ± respectively.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.

Velocity Reconstruction
After solving the discrete pressure Equation (28), the total velocity will be reconstructed in the lowest-order Raviart-Thomas space ( ) 0 RT , the first-order Raviart-Thomas space ( ) RT and the first-order Brezzi-Douglas- Marini space ( ) BDM respectively, refer to [5] for more details about those spaces.The main idea of the re- construction in the current section follows the one depicted in [1], and we will extend it to the situation that the discretization of the pressure equation contains an upwind scheme.
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 w q and n q are the sink and source terms which are zeros here.Noting that if the edge e belongs to both K ∂ and D Γ on the right hand side of Equations ( 29) and (30), and the sign ± is determined by the direction of e n , for example, the sign is positive when e n is the outer normal vector with respect to K ∂ .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 0 RT , 1 RT and 1 BDM spaces, we should specify a set of degree of freedom (DOF) for these H(div) spaces and a corresponding set of basis functions.If let v be any constant in the polynomial space of degree zero 0 P , (29) will vanish and (30) will become the 0 RT space's DOF which is the integral of the normal component of velocity on each edge.Correspondingly, the set of basis functions for 0 RT on the reference element is, ( ) where K is the area of the reference element K and i a is one of its the vertices.

Let | v K and |
v e be any functions in the space of polynomial of degree one 1 P , then (29) and (30) be- come the DOFs for 2 6 2 6 12 6 ˆ, .4 6 It is noted that the choice of DOFs for the 1 BDM and 1 RT spaces is not unique, for example, the half-edge integral of the normal components of velocity is also available and applicable.

Saturation Approximation
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 ( ) is the interface condition of saturation described in (20).

DDG Methods with Some Other Projections
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, BDM ), it means the projection with considering the penalty term but without the upwind scheme, BDM ), it means the projection without considering both the penalty term and the upwind scheme, As indicated in introduction, for all DDG methods the velocity derived by the projection BDM ) preserves the local mass conservation property best, which will shown in the numerical examples.

DDG Method without Explicit Projections
In [2] the velocity is used directly as the combination of the gradient of the solutions and coefficients, as follows, 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 0 RT space.The velocity derived from (41) and ( 42) is close to the velocity projection in

RT
space which is constructed by (39) and (40).But their value in each element is different.Furthermore, the DDG method with using the velocity reconstruction presented in this subsection has certain differences in contrast to what proposed in [2], which are reflected in two aspects below: 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.

Numerical Examples
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 [2].In test 3 we consider the displacement of the wetting phase by the non-wetting phase which is used to simulate the barrier effect in [9].The domains used in the experiments are the square ( ) 2 0, 2 with two corners be cut off, and for the mesh used in the discontinuous problem a small square with different rock property is fixed inside the domain, see Figure 2(a The errors of the local mass conservation at selected times are listed in Table 2 and Table 3.

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

RT
in a homogeneous media.For the results with using projection ( )

BDM and
( ) RT , they are similar with using ( ) RT so are omitted.The mesh used in test 2 is the same as the previous test.The initial and boundary conditions (43)-( 48) are used in this test.To make sure that the water front stays inside the domain, the final time is set to T = 180 s.A constant time step is used, and the ratio of the time step to the space step's square is about  1, Test 2. The contours of wetting phase saturation in the homogeneous medium at selected times are presented in Figure 3.

Test 3
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      1, Test 3. When the oil flows from coarse sand to fine sand with the injection of oil from the inflow boundary in Γ , more and more oil approaches and accumulates at the front of the interface of the fine sand.When the accumulation reaches a critical point, that is, when the capillary pressure at the coarse side of the interface is greater than at the fine side, the accumulated oil will penetrate the interface and enter the fine sand area.By contrast, in

Figure 1 .
Figure 1.The interface (dashed line) between two subdomains with different rock properties.
) and Figure2(b).In each test, we use the Nonsymmetric Interior Penalty Galerkin (NIPG) method with the penalty parameters 1 In order to prevent the oscillations, a slope limiter procedure described in[15] is used.If considering the mixed-Neumann type boundary (8)-(9) for the saturation equation, the following initial and boundary conditions are used: property and Brooks-Corey model are listed in Table

Figure 4 .Γ see Figure 4 .
IIΩ is the fine sand and I Ω is the coarse sand, so the oil-trapped phenomenon will appear on the interfaces J + The critical point in (19) is * 0.44 w S ≈ , and the oil will penetrate the interface J

Figure 3 .
Figure 3.The contours of wetting phase saturation in the homogeneous medium at selected times in the Test 2.

Figure 4 .
Figure 4. Discontinuous quarter-five spot problem.used in this test is Figure 2(b).The initial and boundary conditions (49)-(55) are used in this test.The initial and boundary conditions (49)-(55) are used in this test.To make sure that the water front stays inside the domain, the final time is set to T = 200 s.A constant time step is used, and the ratio of the time step to the space step's square is about 2 4.5 dt h  .The parameters of rock property and Brooks-Corey model are listed in Table1, Test 3. When the oil flows from coarse sand to fine sand with the injection of oil from the inflow boundary in Γ , more and more oil approaches and accumulates at the front of the interface of the fine sand.When the accumulation reaches a critical point, that is, when the capillary pressure at the coarse side of the interface is greater than at the fine side, the accumulated oil will penetrate the interface and enter the fine sand area.By contrast, in be any functions in 1 P polynomial space, thus the basis functions of

Table 1 .
Parameters used in the numerical simulations.
Since there is no sink and source terms the exact local mass is zero on each elements, that is,

Table 2 .
Numerical errors of the local mass conservation in test 1 at times t = 40 s and t = 80 s.

Table 3 .
Numerical errors of the local mass balance in test 1 at times t = 120 s and t = 160 s.