Numerical Study for Two-Phase Flow with Gravity in Homogeneous and Piecewise-Homogeneous Porous Media

ABSTRACT

In this work we study two-phase flow with gravity either in 1-rock homogeneous media or 2-rocks composed media. These phenomena can be modeled by a non-linear scalar conservation law with continuous flux function or discontinuous flux function, respectively. Our study is essentially from a numerical point of view, we apply the new Lagrangian-Eulerian (LEH) finite difference method developed by Abreu and Pérez ^{1), (2), (3), (4), (5), (25} and the classical Lax-Friedrichs method to obtain numerical entropic solutions. Comparisons between numerical and analytical solutions show the efficiency of the methods even for discontinuous flux function. Our main contribution, is the comparison and error analysis between the new LEH and the classical Lax-Friedrichs (LF) methods, in order to show the good performance of the LEH scheme for models with discontinuous flux functions.

Keywords:
conservation laws; finite differences; Lagrangian-Eulerian approach; two-phase flow; heterogeneous porous medium

RESUMO

Neste trabalho estudamos o escoamento bifásico com gravidade em meios homogêneos de 1-rocha ou em meio composto de 2-rochas, estes fenômenos podem ser modelado por uma lei de conservação escalar não-linear com funçãao de fluxo contínua ou função de fluxo descontínua, respectivamente. Nosso estudo é essencialmente de um ponto de vista numérico, aplicamos o novo método de diferenças finitas Lagrangian-Eulerian desenvolvido por Abreu e Pérez ^{1), (2), (3), (4), (5), (25} e o método clássico Lax-Friedrichs para obter soluções numéricas entrópicas. Comparações entre soluções numéricas e analíticas mostram a eficiência dos métodos mesmo para a função de fluxo descontínua.

Palavras-chave:
leis de conservação; diferenças finitas; abordagem Lagrangiana-Euleriana; escoamento bifásico; meio poroso heterogêneo

1 INTRODUCTION

Many problems in engineering, physics and other areas of sciences lead us to the study of conservation laws. For instance, they model many physical phenomena that appear in aerodynamics, fluid mechanics, traffic flow, groundwater flow, multi-phase flow in porous media and others ^{7), (8), (9), (16), (19), (20), (26), (29), (30), (32}. In general a scalar conservation law in one dimension takes the form:

$u_t+{\left(f\left(u\right)\right)}_x=\mathrm{0,}$ (1.1)

where u is the conserved quantity and f(u) is the flux function.

A conservation law jointly with a piecewise constant initial data, having a single discontinuity, are known as Riemann problem. Riemann solutions are in general either composed by shocks, rarefaction waves or combinations of them ^{7), (17), (20), (29), (30}.

Often there is more than one weak solution to the conservation law 1.1 with the same initial data, however only one of them is physically correct, the so called entropic solution. There are distinct entropy criteria for scalar conservation laws with continuous flux function, we use the Oleinik condition which provides a geometrical construction to obtain the entropic solutions ^{17), (24), (30}.

Immiscible two-phase flow in porous media can be modeled by a non-linear scalar conservation law. Buckley and Leverett ⁷ in 1942 solved two-phase flow without gravity by means of fractional flow theory commonly used in petroleum engineering ¹². Proskurowski solved the Buckley-Leverett equation with gravity taken into account ²⁶. Others works on two-phase flow with or without gravity are ^{10), (11), (27), (28), (31}.

In the case that the medium is composed of two types of rocks, the flux function has a spatial discontinuity and the Riemann solutions contain in general an additional stationary shock ^{6), (11}. Although, for the model studied in ¹¹, the analytical entropic solutions can be obtained through an extension of Oleinik construction with appropriate connections between the fluxes ^{6), (11}, there is no an entropy condition that works in general for all discontinuous flux functions. Thus, the usage of appropriate numerical methods to obtain the Riemann solutions of hyperbolic conservation laws with discontinuous flux function gain importance.

Although there are several good numerical methods to deal with discontinuous flux function in conservation laws theory (see for instance ^{13), (14), (21), (22} among others), up to now there are only two finite difference schemes whose convergence to entropic solutions have been rigorously proved for general space-discontinuous flux functions ^{1), (13), (25}: the classical LaxFriedrichs scheme ¹³ and the new Lagrangian-Eulerian (LEH) scheme developed by Abreu and Pérez ^{1), (2), (3), (4), (5), (25} which is a shock-capturing and high-resolution method for first order hyperbolic problems (maybe with forcing terms, because LEH method has also the “well-balance property when adapted for balance laws). The new scheme provides high-order resolution in smooth regions and prevents the creation of spurious oscillations near discontinuities. The new Lagrangian-Eulerian scheme satisfies Harten’s requirements ²⁶ of being a second-order accurate TVD scheme. All these features turn the LEH method a very interesting option to deal with conservation or balance laws with discontinuous flux functions. Our main contribution in this paper, is the comparison and error analysis between the new LEH and the classical LaxFriedrichs (LF) methods, in order to show the good performance of the LEH scheme for models with discontinuous flux functions.

Motivated by this, in this work we solve some Riemann problems by both the LagrangianEulerian and the Lax-Friedrichs methods for two-phase flow with gravity in two cases; first we test our codes for the homogeneous rock case and then we obtain the numerical solutions for a piecewise-homogeneous rock (two-rocks composed media). In the homogeneous test-case we compare the numerical solutions with the analytical solutions obtained via the classic Oleinik construction. In the case of two-rocks composed media, we use the analytical entropic solutions appearing in ¹¹ and briefly explained in Section 3 of this paper, to calculate the errors of the methods.

2 IMMISCIBLE TWO-PHASE FLOW WITH GRAVITY

We study the one-dimensional flow in porous media of two immiscible fluids: water and oil, under the effect of gravity force. Initially there is an impermeable interface, in $z=0$, separating the mixture of the fluids and we study the flow near the interface, once the interface disappears. This is equivalent to solve a class of Riemann problems for a scalar conservation law.

We assume also that the flow is isothermal, there are no sources or sinks, the porosity is constant, compressibility and capillary pressure effects are neglected. Equation (2.1) expresses the conservation of mass:

$\frac{\partial }{\partial t}\left(\varphi s_i\right)+\frac{\partial }{\partial z}{v}_i=\mathrm{0,}\begin{array}{c}\end{array}i=w,o,$ (2.1)

$v_i=-K\frac{k_i}{{\mu }_i}\left(\frac{\partial p_i}{\partial z}-{\rho }_{i}g{e}_z\right),\begin{array}{c}\end{array}i=w,o.$ (2.2)

Here, ϕ denotes porosity, s _i , v _i , respectively denote saturation and seepage velocity of the phase i. The seepage velocity is given by Darcy’s law (2.2), in which K is the absolute permeability of the rock, k _i , µ _i , p _i and ρ _i are the fluid relative permeability, viscosity, pressure and density for each phase i. Also in (2.2) g denotes the gravitational constant. The viscosities µ _i are assumed constants, while the permeabilities k _i are assumed to be functions of the saturations.

Under the assumption that the porous medium is totally saturated $\left(s_w+s_o=1\right)$ the equation (2.1)-(2.2) reduce to the following scalar conservation law.

$\frac{\partial s_w}{\partial t}+\frac{\partial }{\partial z}\left(f\left(s_w\right)\right)=\mathrm{0,}$ (2.3)

with flux function

$f\left(s_w\right)=\frac{s_w^2}{s_w^2+\mu {\left(1-s_w\right)}^2}\left[v+{\left(1-s_w\right)}^2\mu \left(1-\rho \right)\right].$ (2.4)

here we denote $\mu ={\mu }_w/{\mu }_o$ and $\rho ={\rho }_o/{\rho }_w,s_w$ is the saturation of the phase water. The dimensionless parameter v is related to the pressure gradients.

For two-rock composed media, the porous medium is piecewise-homogeneous, so we define the permeability of the rock by the following equation

$K\left(z\right)=\{\begin{array}{c}{K}_l,\text{if}z<0\\ K_r,\text{if}z>0\end{array}$

which leads to a discontinuous flux function of the form

$f\left(s_w,z\right)=\{\begin{array}{c}{f}_l,\text{if}z<\mathrm{0,}\\ f_r,\text{if}z>\mathrm{0,}\end{array}$ (2.5)

where

$f_l=\frac{s_w^2}{s_w^2+\mu {\left(1-s_w\right)}^2}\left[v+{\left(1-s_w\right)}^2\mu \left(1-\rho \right)\right]$ (2.6)

and

$f_r=\frac{s_w^2}{s_w^2+\mu {\left(1-s_w\right)}^2}\left[v+{\left(1-s_w\right)}^2\mu \left(1-\rho \right)\frac{K_r}{K_l}\right].$ (2.7)

The initial Riemann data for equations (2.3)-(2.5) are denoted by

$s_w\left(z\mathrm{,0}\right)=\{\begin{array}{c}{s}_l\begin{array}{c}\end{array}z\le \mathrm{0,}\\ s_r\begin{array}{c}\end{array}z>0.\end{array}$ (2.8)

3 ANALYTICAL SOLUTIONS FOR SCALAR CONSERVATION LAWS

For non-linear scalar conservation laws, discontinuous or non-smooth solutions commonly appear even for smooth initial data. Thus a weak formulation of the conservation law and the corresponding “weak” concept of solutions are needed to deal with this kind of equations. Specially for Riemann problems the weak solutions are combinations of shocks and rarefaction waves. The shocks (even for smooth initial data) happen when characteristics cross, at the time where the characteristics first cross, the function s(z, t) has an infinite slope the wave “breaks” and a shock forms ¹⁷. The shock waves joining states s _l and s _r travel with constant speeds σ which satisfy the Rankine-Hugoniot condition:

$\sigma =\frac{f\left(s_r\right)-f\left(s_l\right)}{s_r-s_l}.$ (3.1)

Unfortunately weak solutions are not unique and additional conditions are needed in order to select the physically correct one (entropic solution). For instance we consider the Oleinik entropy condition which apply to general non convex flux f: the solution s(z, t) is the entropy solution if all discontinuities have the property that ^{17), (23}

$\frac{f\left(s\right)-f\left(s_l\right)}{s-s_l}\ge \sigma \ge \frac{f\left(s\right)-f\left(s_r\right)}{s-s_r},$ (3.2)

for all s between s _l and s _r and σ is gives by equation (3.1).

The above exposed, is the basement of the Oleinik’s geometric construction which allows to obtain the entropic analytic solutions for any Riemann problem in the one-dimensional scalar case for continuous flux functions (the case of two-phase flow in homogeneous porous medium). Thus the entropic solution can be constructed by using the convex or concave hulls of the flux function, depending on the relation between the initial Riemann left and right states, see the Figure 1. If part of this hull is a tangent or secant line to the flux function (for example the dotted and dashed lines in the Figure 1) this represents a shock wave, with shock speed coinciding with the inclination of the line. Everywhere the hull coincides with the flux function (for example, if $s>s_1$ when $s_r=0$ or if $s<s_0$ when $s_r=1$ in Figure 1), the solution is a rarefaction wave, with spreading speeds varying along with the derivative of the flux function. For the case of combinations of shocks and rarefactions, the junction points between them can be obtained by calculating the zero of the equation (3.3)

$\frac{f\left(s^{+}\right)-f\left(s^{-}\right)}{s^{+}-s^{-}}-f\text{'}\left(s^{+}\right)=\mathrm{0,}\text{with}{s}^{-}\text{fixed}$ (3.3)

for a right characteristic shock wave or the zero of the equation

$\frac{f\left(s^{+}\right)-f\left(s^{-}\right)}{s^{+}-s^{-}}-f\text{'}\left(s^{-}\right)=\mathrm{0,}\text{with}{s}^{+}\text{fixed}$ (3.4)

for a left characteristic shock wave.

Figure 1:
Example of flux function f and entropic shock waves as limits of travelling waves.

In ¹¹ Kaasschieter considered a generalization for discontinuous flux functions in space variable z

$\begin{array}{l}\frac{\partial s}{\partial t}+\frac{\partial }{\partial z}{f}_l\left(s\right)=\mathrm{0,}z<\mathrm{0,}\\ \frac{\partial s}{\partial t}+\frac{\partial }{\partial z}{f}_r\left(s\right)=\mathrm{0,}z>\mathrm{0,}\\ f_l\left(S_L\right)=f_r\left(S_R\right).\end{array}$ (3.5)

Here $f_l,f_r:\left[\mathrm{0,1}\right]\to ℝ$ are assumed to be twice differentiable such that $f_l\left(0\right)=f_r\left(0\right)$ and $f_l\left(1\right)=f_r\left(1\right)$. The states

$S_L\left(t\right)=\underset{z\uparrow 0}{\lim }s\left(z,t\right),\begin{array}{c}\end{array}{S}_R\left(t\right)=\underset{z\downarrow 0}{\lim }s\left(z,t\right),t>0$ (3.6)

are the left and right states of the stationary shock always appearing in the solutions. They are called fluxes connections in literature ^{6), (11}. The last expression in (3.5) represents the Rankine-Huguniot condition for stationary shocks.

Kaasschieter obtained entropy conditions for a model of two-phase flow with gravity in two-rocks of the form (3.5) with flux functions in (2.5)-(2.7). We will use the same entropy conditions to construct our analytical solutions in Section 5.2.

To obtain the entropy conditions, Kaasschieter ¹¹ regularised the problem by adding the capillary diffusion term, he considered similar solutions $s\left(z,t\right)=s(\xi )$, with $\xi =z/t$ for the regularized problem, he assumed some hypothesis like $S_L\ne S_R$ and $p_{c,L}\left(S_L\right)\ne p_{c,R}\left(S_R\right)$ (distinct capillary pressures at the connection states), between others, and after some calculations he proved that:

• (i) All non-stationary waves must satisfy the classical Oleinik entropy condition for one of the flux functions f _l or f _r . Of course, waves with negative speeds are obtained through the convex or concave hulls of f _l , while waves with positive speeds are obtained through the convex or concave hulls of f _r .

• (ii) There is a stationary shock separating the non-stationary waves groups. The stationary shock connecting states S _L and S _R can be obtained geometrically by using a horizontal jump line connecting the fluxes.

• (iii) The appropriate (unique) selection of these connection states, requires an entropy condition. For this problem the states S _L and S _R must satisfy the entropy inequality

$f_l^{\text{'}}\left(S_L\right)\underset{\_}{\succ }0\text{or}{f}_r^{\text{'}}\left(S_R\right)\underset{\_}{\prec }0$ (3.7)

• where the notation $f_l^{\text{'}}\left(S_L\right)\underset{\_}{\succ }0$ is used if $f_l^{\text{'}}\left(S_L\right)\ge 0$ and ${\lim }_{s\uparrow S_L}sign\left(f_l^{\text{'}}\left(s\right)\right)=1$ or ${\lim }_{s\downarrow S_L}sign\left(f_l^{\text{'}}\left(s\right)\right)=1$. Analogously the notation $f_r^{\text{'}}\left(S_R\right)\underset{\_}{\prec }0$ means that $f_r^{\text{'}}\left(S_R\right)\le 0$ and ${\lim }_{s\uparrow S_R}sign\left(f_r^{\text{'}}\left(s\right)\right)=-1$ or ${\lim }_{s\downarrow S_R}sign\left(f_r^{\text{'}}\left(s\right)\right)=-1$.

Thus, in order to obtain the entropic solutions for our model with discontinuous fluxes, the key fact is the appropriate selection of the connection states S _L and S _R of the stationary shock by using a horizontal jump connecting the fluxes while satisfying the entropy inequality in (3.7), the other waves can be obtained through the Oleinik construction for f _l and f _r , for illustration see Figures 8(a), 9(a),10(a), 11(a), 12(a), 13(a) in Subsection 5.2.

4 FINITE DIFFERENCE METHODS

Finite difference methods are numerical methods for solving differential equations, in particular very useful for conservation laws. This kind of method consists on replacing the derivatives in the differential equations by finite difference approximations obtained through Taylor expansions of the functions.

In order to preserve the conservation property of the solutions, finite differences schemes in conservative forms ¹⁵ are used:

$U_j^{n+1}=U_j^n-\frac{k}{h}\left(G\left(U_j^n,U_{j+1}^n\right)-G\left(U_{j-1}^n,U_j^n\right)\right).$ (4.1)

The function G is called the numerical flux.

Lax-Friedrichs Scheme

The classical Lax-Friedrichs (LF) scheme for conservation laws with continuous flux functions in conservative form (4.1) has the following numerical flux, see ¹⁸:

$G\left(U_{j-1}^n,U_j^n\right)=\frac{1}{2}\left(f\left(U_{j-1}\right)+f\left(U_j\right)\right)-\frac{1}{2v}\left(U_j^n-U_{j-1}^n\right),$

where

$v=k/h.$

For this method to be convergent, it must satisfy the CFL condition ¹⁸

${\max }_j\left|f\text{'}\left(U_j^n\right)\right|\frac{k}{h}\le 1.$

For conservation laws with discontinuous flux functions the numerical flux of the method can be adapted to

$G\left(U_{j-1}^n,U_j^n\right)=\frac{1}{2}\left(f\left(U_{j-1},x_{j-1}\right)+f\left(U_j,x_j\right)\right)-\frac{1}{2v}\left(U_j^n-U_{j-1}^n\right),$

where

$f\left(U_j,x_j\right)=\{\begin{array}{c}{f}_l\left(U_j\right),\text{if}{x}_j<\mathrm{0,}\\ f_r\left(U_j\right),\text{if}{x}_j>\mathrm{0,}\end{array}$ (4.2)

with CFL condition

${\max }_j\left(\left|f_l^{\text{'}}\left(U_j^n\right)\right|,\left|f_r^{\text{'}}\left(U_j^n\right)\right|\right)\frac{k}{n}\le 1.$

Lagrangian-Eulerian Scheme

The Lagrangian-Eulerian (LEH) method developed by Abreu and Pérez ^{1), (2), (3), (4), (5), (25} is a scheme based on a Eulerian central scheme finite volume formulation, but in a space-time Lagrangian-Eulerian framework. These scheme is derived from the divergence forms of the equations ²⁵, and it is an order high-resolution scheme for numerically solving nonlinear conservation law problems. This scheme provides high-order resolution in smooth regions and prevents the creation of spurious oscillations near discontinuities. The convergence of the numerical solutions to the entropic solutions for models with discontinuous flux functions was proved in ^{25), (1}.

This method for hyperbolic conservation laws with continuous flux functions is given by the following numerical formulation, see ^{1), (25} for the details:

$U_j^{n+1}=\frac{1}{4}\left(U_{j-1}^n+2U_j^n+U_{j+1}^n\right)-\frac{k}{2h}\left(f\left(U_{j+1}^n\right)-f\left(U_{j-1}^n\right)\right),$ (4.3)

and the CFL condition for this scheme is given by

${\max }_j\left|f\text{'}\left(U_j^n\right)\right|\frac{k}{h}\le \frac{\sqrt{2}}{2}.$

The scheme (4.3) can be written in the conservative form (4.1), where the numerical flux is, see ^{1), (25}

$G\left(U_j^n,U_{j+1}^n\right)=\frac{1}{4}\left[\frac{h}{k}\left(U_j^n-U_{j+1}^n\right)+2\left(f\left(U_{j+1}^n\right)+f\left(U_j^n\right)\right)\right].$

For conservation laws with discontinuous flux functions the method use the following numerical flux

$G\left(U_j^n,U_{j+1}^n\right)=\frac{1}{4}\left[\frac{h}{k}\left(U_j^n-U_{j+1}^n\right)+2\left(f\left(U_{j+1}^n,x_{j+1}\right)+f\left(U_j^n,x_j\right)\right)\right],$

with f in the form (4.2) and CFL (see ^{1), (25} for details):

${\max }_j\left(\left|f_l^{\text{'}}\left(U_j^n\right)\right|,\left|f_r^{\text{'}}\left(U_j^n\right)\right|\right)\frac{k}{h}\le \frac{\sqrt{2}}{2}.$

5 NUMERICAL SOLUTIONS

5.1 Two-Phase Flow in One-Rock Media

First we consider a test-case of homogeneous media, i.e., constant absolute permeability of the rock K. We proceed to compare the numerical solutions in the following two cases:

• (i) Purely gravitational flow (i.e., $v=0$);

• (ii) Mixed flow with dominant gravity (v small and non-zero, for instance $v=0.01$). For both cases we set the following parameters $\mu =0.25$ and $\rho =0.8$.

In order to illustrate the solution for the pure gravitational problem (case i), we consider initial Riemann data $s_l=1$ and $s_r=0$, see Figure 2. The solution consists in a shock wave with negative speed joining $s_l=1$ with $s_1=0.4732$ traveling upwards in the reservoir. There is also a rarefaction wave with positive and negative speeds joining the states $s_1=0.4732$ with $s_2=0.274$, preceding a last shock wave with positive speed joining $s_2=0.274$ with $s_r=0$. This last shock moves to the right (bottom in the reservoir).

Figure 2:
Case (i) with the Riemann data $s_l=\mathrm{1,}{s}_r=0$.

A numerical refinement study for this case is shown in Figures 2b-2d to illustrate that accuracy increases when the mesh is refined. In the Table 1, we show the errors between the approximate solution obtained by Lagrangian-Eulerian and Lax-Friedrichs methods when compared to the analytical Oleinik solution.

Table 1:
Errors between the approximate solutions U and exact solutions u of case (i) with Riemann data $s_l=\mathrm{1,}{s}_r=0$.

Cells	$LEH{\Vert u-U\Vert }_{l_h^1}$	$LEH{\Vert u-U\Vert }_{l_h^2}$	$LF{\Vert u-U\Vert }_{l_h^1}$	$LF{\Vert u-U\Vert }_{l_h^2}$
128	1.57 × 10−2	5.04 × 10−2	2.52 × 10−2	6.79 × 10−2
256	9.80 × 10−3	4.00 × 10−2	1.60 × 10−2	5.38 × 10−2
512	6.10 × 10−3	3.23 × 10−2	1.03 × 10−2	4.32 × 10−2
1024	3.50 × 10−3	2.37 × 10−2	6.20 × 10−3	3.27 × 10−2
2048	2.20 × 10−3	1.97 × 10−2	3.90 × 10−3	2.65 × 10−2
4096	1.30 × 10−3	1.54 × 10−2	2.30 × 10−3	2.07 × 10−2

The variation of the error with respect to the time is shown in Figure 3, notice that as time increases the error tends to increase, and this is smaller in refined meshes.

Figure 3:
Errors Lax-Friedrichs and Lagrangian Eulerian of case (i) with Riemann data $s_l=\mathrm{1,}{s}_r=0$.

In Figures 4 and 5 the solutions for other Riemann data in the case (i) are illustrated.

Figure 4:
Case (i) with the Riemann data $s_l=\mathrm{0.4,}{s}_r=0$.

Figure 5:
Case (i) with the Riemann data $s_l=\mathrm{1,}{s}_r=0:3$.

In the case (ii), for mixed flow with dominant gravity, the Riemann solutions for initial data $s_l=1$ and $s_r=0$, are shown in the Figure 6 and in Figure 7 for initial data $s_l=0$ and $s_r=0.8$.

Figure 6:
Case (ii) with the Riemann data $s_l=\mathrm{1,}{s}_r=0$.

Figure 7:
Case (ii) with the Riemann data $s_l=\mathrm{0,}{s}_r=0.8$.

5.2 Two-Phase Flow in Two-Rocks Media

For two-rocks composed media the conservation law has a discontinuous flux function (2.5)-(2.7). We consider some distinct combination of values for µ, ρ, K _l , K _r and we only study the pure gravitational flow, i.e. $v=0$.

• Case 1: $\rho =3/2$ (heavy oil), permeabilities $K_l=0.5$ and $K_r=1$, so $f_l>f_r$;

• Case 2: $\rho =3/2$ (heavy oil), permeabilities $K_l=1$ and $K_r=0.5$, so $f_l<f_r$;

• Case 3: $\rho =0.8$ (light oil), permeabilities $K_l=0.5$ and $K_r=1$, so $f_l<f_r$;

• Case 4: $\rho =0.8$ (light oil), permeabilities $K_l=1$ and $K_r=0.5$, so $f_l>f_r$.

In the case 1, three simulations were performed for distinct values of µ and initial Riemann data, see the Figures 8, 9, 10. To illustrate we explain the solution in Figure 8: it consists of two groups of waves, the first group has negative speeds and travels upward in the reservoir, it is composed by a shock wave joining the states $s_l=0$ and $s_1=0.298$ followed by a rarefaction wave up to state $S_L=0.4095$. The second group consists in a solely shock wave with positive speed travelling downwards in the reservoir connecting the states $S_R=0.6659$ and $s_l=1$. There is a stationary shock joining the states $S_L=0.4095$ with $S_R=0.6659$ separating the two groups of waves.

Figure 8:
Case 1 with $t=\mathrm{0.5,}\mu =1/3$ and Riemann data $s_l=\mathrm{0,}{s}_r=1$.

Figure 9:
Case 1 with $t=\mathrm{1,}\mu =0.1$ and Riemann data $s_l=\mathrm{0,}{s}_r=0.25$.

Figure 10:
Case 1 with $t=\mathrm{1,}\mu =0.1$ and Riemann data $s_l=\mathrm{0.4,}{s}_r=1$.

The errors of numerical solutions are shown in the Table 2, notice that Lagrangian-Eulerian scheme had better performance (more accurate) than the Lax-Friedrichs method and also that the error decreases when the mesh is refined, see the Figures 8b-8d.

Table 2:
Errors between the numerical solutions U and exact solutions u of case 1 with $t=\mathrm{0.5,}\mu =1/3$ and Riemann data $s_l=\mathrm{0,}{s}_r=1$.

Cells	$LF{\Vert u-U\Vert }_{l_h^1}$	$LF{\Vert u-U\Vert }_{l_h^2}$	$LEH{\Vert u-U\Vert }_{l_h^1}$	$LEH{\Vert u-U\Vert }_{l_h^2}$
256	1.48 × 10−2	5.59 × 10−2	1.00 × 10−2	4.76 × 10−2
512	8.90 × 10−3	4.28 × 10−2	5.80 × 10−3	3.55 × 10−2
1024	5.30 × 10−3	3.27 × 10−2	3.30 × 10−3	2.66 × 10−2
2048	3.00 × 10−3	2.41 × 10−2	1.80 × 10−3	1.91 × 10−2

The solution in Figure 9 has a structure similar to the solution in Figure 8. The last simulation corresponding to the case 1, see Figure 10, shows a solution without the negative-speed group of waves and only a stationary shock followed by a positive speed shock.

For the cases 2, 3 and 4 the simulations shown in the Figures 11, 12 and 13 respectively, illustrate again the well performance of the methods in particular the high accuracy of the LagrangianEulerian scheme when compared with Lax-Friedrichs method. For all the cases the analytical solutions were obtained using the entropic condition (3.7) developed in ¹¹ to choose the appropriate connection between the fluxes (stationary shock).

Figure 11:
Case 2 with $t=\mathrm{0.5,}\mu =1/3$ and Riemann data $s_l=\mathrm{0,}{s}_r=0.4$.

Figure 12:
Case 3 with $t=\mathrm{1,}\mu =1/3$ and Riemann data $s_l=\mathrm{1,}{s}_r=0$.

Figure 13:
Case 4 with $t=\mathrm{2,}\mu =1/3$ and Riemann data $s_l=\mathrm{0.2,}{s}_r=0.6$.

6 CONCLUSIONS

The Lax-Friedrichs and the Lagrangian-Eulerian methods were used to obtain the approximate Riemann solutions for two-phase flow with gravity in both homogeneous and piecewisehomogeneous porous media (two-rocks media). For two-rocks media the flux function has a discontinuity in space variable z so find analytical solutions is not simple. We used the entropy criterion (3.7) developed in ¹¹ to compare with our numerical results. Both methods LF and LEH present a diffusive behavior, however the Lagrangian-Eulerian scheme showed to be less diffusive and more accurate as reflected in the table and figures shown in Section 5. Notice also that the error decreases when the mesh is refined.

REFERENCES

• ¹
  E. Abreu, V. Matos, J. Perez & P. Rodriguez-Bermudez. A class of Lagrangian-Eulerian shock-capturing schemes for first-order hyperbolic problems with forcing terms. sent, (2019).
• ²
  E.C. Abreu, W. Lambert, J. Perez & A. Santo. A new finite volume approach for transport models and related applications with balancing source terms. Mathematics and Computers in Simulation, 137 (2017), 2-28.
• ³
  E.C. Abreu & J.A. Perez. A Lagrangian-Eulerian algorithm scheme for hyperbolic conservation laws and balance laws. HYP2014XV International Conference on Hyperbolic Problems, (2014). Acesso em 31 de agosto de 2018.
• ⁴
  E.C. Abreu & J.A. Perez. A new locally Conservative Lagrangian Eulerian method for hyperbolic and Balance laws. VIII Pan-American Workshop Applied and Computational Mathematics, (2014).
• ⁵
  E.C. Abreu & J.A. Perez. A fast, robust, and simple Lagrangian-Eulerian solver for balance laws and applications. Computers & Mathematics with Applications, 77(9) (2019), 2310-2336.
• ⁶
  B. Andreianov & C. Cance`s. Vanishing capillarity solutions of Buckley-Leverett equation with gravity in two-rocks medium. Comput. Geosci., 17(3) (2013), 551-572.
• ⁷
  S. Buckley & M. Leverett. Mechanims of fluid displacement in sands. Trans. AIME, 146 (1942), 187-196.
• ⁸
  P. Castaneda, E. Abreu, F. Furtado & D. Marchesin. On a universal structure for immiscible threephase flow in virgin reservoirs. Computational Geosciences, 20 (2016), 171-185.
• ⁹
  C. Dafermos. “Hyperbolic Conservation Laws in Continuum Physics”. Springer (2000).
• ¹⁰
  M. Hayek, E. Mouche & C. Mugler. Modeling Vertical Stratification of CO2 injected into deeep layered aquifier. In “Advances in Water Resours”, volume 32 (2009), pp. 450-462.
• ¹¹
  E. Kaasschieter. Solving the Buckley Leverett Equation with gravity in a Heterogeneous Porous Medium. Comp. Geosci., 3 (1999).
• ¹²
  A. Kantzas, J. Bryan & S. Taheri. Fundamentals of Fluid Flow in Porous Media (S.I.). URL http://perminc.com/resources/fundamentals-of-fluid-flow-in-porous-media/
  » http://perminc.com/resources/fundamentals-of-fluid-flow-in-porous-media/
• ¹³
  K.H. Karlsen & J.D. Towers. Convergence of the Lax-Friedrichs Scheme and Stability for Conservation Laws with a Discontinuous Space-Time Dependent Flux. Chinese Annals of Mathematics Ser. B, 25(3) (2004), 287-318.
• ¹⁴
  K.H. Karlsen & J.D. Towers. Convergence of a Godunov scheme for conservation laws with a discontinuous flux lacking the crossing condition. J. Hyper. Differential Equations, 14(4) (2017), 671-701.
• ¹⁵
  D. Kröner. “Numerical Schemes for Conservation Laws”. John Wiley Sons, Stuttgart (1997).
• ¹⁶
  P.D. Lax. Hyperbolic systems of conservation laws II. Comm. Pure Appl. Math., 10 (1957), 537-566.
• ¹⁷
  R.J. LeVeque. “Numerical Methods for Conservation Laws”. Birkhauser, Basel; Boston; Berlin, 2 ed. (1992).
• ¹⁸
  R.J. LeVeque. “Finite Difference Methods for Ordinary and Partial Differential Equations”. SIAM (2007).
• ¹⁹
  H.W. Liepmann & A. Roshko. “Elements of gasdynamics”. John Wiley & Sons (1959).
• ²⁰
  D. Marchesin, E. Isaacson, B. Plohr & J.B. Temple. Multiphase Flow Models with Singular Riemann Problems. Computation and Applied Mathematics, 11(2) (1992), 147-166.
• ²¹
  S. Mishra. Numerical Methods for Conservation Laws With Discontinuous Coefficients. Handbook of Numerical Analysis, 18 (2017), 479-506.
• ²²
  S. Mishra. On the convergence of numerical schemes for hyperbolic systems of conservation laws. Proc. Int. Cong. of Math., Rio de Janeiro, 3 (2018), 3625-3652.
• ²³
  O. Oleinik. Discontinuous solutions of nonlinear differential equations. Amer. Math. Soc. Trans. Ser., 26 (1957), 95-172.
• ²⁴
  O.A. Oleinik. On the uniqueness of generalized solution of Cauchy problem for non linear system of equations occurring in mechanics. Uspekhi Mat. Nauk, 12 (1957), 3-73.
• ²⁵
  J.A. Perez. “Lagrangian-Eulerian approximation methods for balance laws and hyperbolic conservation laws”. Ph.d. thesis, UNICAMP (2015).
• ²⁶
  W. Proskurowski. A note on solving the Buckley-Leverett equation in the presence of gravity. J. Comput. Phys., 41 (1981), 136-141.
• ²⁷
  A. Riaz & H. Tchelepi. Dynamics of vertical displacement in porous media associated with CO2 sequestration. SPE Journal, (2008).
• ²⁸
  A. Riaz & H.A. Tchelepi. Stability of two-phase vertical flow in homogeneous porous media. Phys. Fluids, 19(7) (2007). DOI:10.1063/1.2742975.
  » https://doi.org/10.1063/1.2742975
• ²⁹
  P. Rodríguez-Bermúdez & D. Marchesin. Riemann Solutions for vertical flow of three phases in porous media: simple cases. Journal of Hyperbolic Differential Equations, 10 (2013), 335-370.
• ³⁰
  C.J. van Duijn. “An Introduction to Conservations Laws: Theory and Applications to Multi-Phase Flow”. Lecture Notes, Delft University of Technology (2002).
• ³¹
  M. Wangen. Vertical migration of hydrocarbons modeled with fractional flow theory. Geophysical Journal International, 115 (1993), 109-131.
• ³²
  G.B. Whitham. “Linear and Nonlinear Waves”. Wiley-Interscience (1999).

Publication Dates

• Publication in this collection
  30 Apr 2020
• Date of issue
  Jan-Apr 2020

History

• Received
  10 Dec 2018
• Accepted
  04 Sept 2019