Numerical Solutions of the Classical and Modified Buckley-Leverett Equations Applied to Two-Phase Fluid Flow

Raphael de O. Garcia¹, Graciele P. Silveira²

¹Department of Actuarial Science, Federal University of São Paulo, Osasco, Brazil.
²Department of Physics, Chemistry and Mathematics, Federal University of São Carlos, Sorocaba, Brazil.

DOI: 10.4236/ojfd.2024.143009                

Abstract

Climate change is a reality. The burning of fossil fuels from oil, natural gas and coal is responsible for much of the pollution and the increase in the planet’s average temperature, which has raised discussions on the subject, given the emergencies related to climate. An energy transition to clean and renewable sources is necessary and urgent, but it will not be quick. In this sense, increasing the efficiency of oil extraction from existing sources is crucial, to avoid waste and the drilling of new wells. The purpose of this work was to add diffusive and dispersive terms to the Buckley-Leverett equation in order to incorporate extra phenomena in the temporal evolution between the water-oil and oil-water transitions in the pipeline. For this, the modified Buckley-Leverett equation was discretized via essentially weighted non-oscillatory schemes, coupled with a three-stage Runge-Kutta and a fourth-order centered finite difference methods. Then, computational simulations were performed and the results showed that new features emerge in the transitions, when compared to classical simulations. For instance, the dispersive term inhibits the diffusive term, adding oscillations, which indicates that the absorption of the fluid by the porous medium occurs in a non-homogeneous manner. Therefore, based on research such as this, decisions can be made regarding the replacement of the porous medium or the insertion of new components to delay the replacement.

Keywords

Computational Fluid Dynamics, Buckley-Leverett Equation, Numerical Methods, Two-phase Fluid Flow

1. Introduction

Technological advances achieved after the Industrial Revolution, followed by the development of equipment dependent on energy sources from fossil fuels, made the petroleum industry play a crucial role in the world economy. As it is a depletable natural source and a potential pollutant, efficiency in the extraction process and a detailed understanding of the phenomena involved become indispensable for the challenges of delivering a product inserted in a clean economy.

This historical role of the use of oil, natural gas and coal as the main sources of the global energy system has resulted in the current climate problem, defining as urgent the need for an energy transition to a low-carbon future via renewable sources such as wind and solar.

Although such sources are important in the energy transition, the change will not be immediate. In this way, creative and efficient uses of natural sources are essential for the transition of the energy system to a low-carbon one to be viable. Therefore, the need to maximize the efficiency of the use of these sources becomes urgent [1].

Fluid flows, including two-phase flows, represent a substantial issue in several technological applications and natural sciences. Analytical studies of mathematical models, such as those of [2], provide theoretical bases for the existence of solutions. At the same time, these topics also provide content for the areas of applied mathematics and engineering, with an emphasis on numerical analysis and scientific computing. A multidisciplinary approach, which takes into account porous media, is essential and aims to detail the dynamics of multiphase fluids, so that process efficiency results in fewer wells being drilled, both onshore and in deep water [3].

A classic modeling in this area consists of analyzing the displacement of oil (or petroleum) and/or gases through pipes filled with a porous medium, characterized by the injection of another fluid (saturated water) to help maintain the flow inside the pipe. In this sense, the mathematical model used to describe the flow of two-phase incompressible fluids is the classical Buckley-Leverett equation [4]. The nonlinearity characteristic of the partial differential equation requires attention in the use of numerical methods and computational techniques, so that the approximate solutions obtained are not harmed by spurious oscillations (numerical dispersion) or excessive numerical dissipation (numerical diffusion).

In preliminary studies, [5] investigated an essentially non-oscillatory weighted fifth-order scheme, applied to the classical Buckley-Leverett equation, and additionally, a numerical study related to the order of precision and stability was carried out. Furthermore, the authors evaluated the inclusion of a diffusive term.

Laboratory experiments on two-phase flow in pipes filled with porous material revealed complex profiles that include fluid infiltrations modeled by diffusive and dispersive terms [6]. The authors of [7] studied two-phase fluid flow between air and water (via continuity and momentum equations) in cathode channels and showed that different pore structures, which result in different infiltration profiles, significantly affect water dynamics. Findings such as these suggest that modifications to the classical Buckley-Leverett equation are necessary for a more detailed understanding of the occurrences involved in these flows and in the infiltrations.

Thus, the objective of this work was to add diffusive terms and dispersive terms to carry out a study on the impact that these have on numerical solutions, when compared with standard solutions of the classical Buckley-Leverett equation. For this purpose, an essentially non-oscillatory weighted scheme was used, coupled with a three-stage Runge-Kutta method and a finite difference scheme centered on the discretization of the modified Buckley-Leverett equation. Using proprietary codes developed in Octave, numerical solutions capable of representing the temporal evolution of different initial scenarios were performed, focusing on transitions between oil-water and water-oil, where the fluid dynamics undergo changes depending on the diffusive and dispersive parameters adopted. Finally, the implemented simulations showed the impacts generated by the variation of the dispersive term, producing more detailed results from the point of view of the solutions of the classical equation.

2. Mathematical Modeling

For centuries, a problem that has been of global interest is the extraction of petroleum underground through a tube filled with a porous medium. After drilling the soil to the underground oil reservoir, a certain amount is drained due to the high pressure that the oil is found, but as the extraction progresses, there is a decrease in pressure with consequent interruption of the flow, still leaving a lot of petroleum in the subsoil. A standard method subsequent to the initial extraction is to pump water into rest oil reservoir to force the continuation of extraction. In this case, the fluid is two-phase, oil and water, and the flow is restarted in the porous medium consisting of rock or sand.

The mathematical modeling that will be described consists of representing the oil one-dimensional flow through the tube filled with porous material, by water pumping [8]. Such a model was initially proposed by Buckley and Leverett (1942) [4], in researches on the flow of two-phase incompressible fluids in porous media.

Let $0\le q\left(x,t\right)\le 1$ be the fraction of saturated water and $1-q\left(x,t\right)$ the fraction of oil contained in a pipe filled with a porous material. Such fluids are essentially incompressible, which ensures that the total flow between the pipe ends is equal to any smaller portion of the pipe.

In this way, in regions of the tube where $q=0$ (pure oil) or $q=1$ (pure water), the velocities are constant and distinct, but when $0<q<1$ , the difference between the surface tensions of fluids causes them to move and mix. Buckley and Leverett [4] proposed a model in which the rate of change of q over time ($q_t$ ) is described by the following conservation law,

$q_t+f{\left(q\right)}_x=0,$ (1)

in which $f\left(q\right)=\frac{q^2}{q^2+a^2{\left(1-q\right)}^2}$ is the water flux, $0<a<1$ represents the porosity of the medium and $1-f\left(q\right)$ is the oil flux, with $q=q\left(x,t\right)$ . Equation (1) models a flow from left to right, in which the tube thickness does not influence the dynamics in question.

The reestablishment of the oil flow, from left to right, can be done by filling part of the pipe on the left with saturated water, allowing the resumption of oil extraction. An initial condition for modeling that procedure is

$q\left(x,0\right)=\{\begin{array}{ll}1,\hfill & \text{if}x<x_0\hfill \\ 0,\hfill & \text{if}x>x_0\hfill \end{array}\left(\text{scenario}1\right).$ (2)

The initial condition, Equation (2), can be approximated by the following function:

$q\left(x,0\right)=1-\left[\frac{1+\tanh \left(\alpha \left(x-x_0\right)\right)}{2}\right],$ (3)

in which $x_0$ is a parameter associated with the position of the function and $\alpha$ corresponds to how fast the function changes from zero to one.

Another situation that can occur is, after an injection of water, the flow is restarted and a second oil portion enters the pipeline, soon after the interruption in the supply of saturated water. This is the desired circumstance, when only an amount of water injected is enough for the continuity of the extraction. The condition to represent that scenario is the conservation law

$q\left(x,0\right)=\{\begin{array}{ll}0,\hfill & \text{if}{x}_1<x\hfill \\ 1,\hfill & \text{if}{x}_1<x<x_2\hfill \\ 0,\hfill & \text{if}x>x_2\hfill \end{array}\left(\text{scenario}2\right).$ (4)

Equation (4) can also be approximated by an expression involving hyperbolic tangent functions:

$q_b\left(x,0\right)=\left[\frac{1+\tanh \left({\alpha }_1\left(x+x_1\right)\right)}{2}\right]+\left[\frac{1-\tanh \left({\alpha }_2\left(x-x_2\right)\right)}{2}\right]-1,$ (5)

where ${\alpha }_1$ represents how fast the transition between oil and water is, in the neighborhood of $x=x_1$ and ${\alpha }_2$ describes how fast the water-oil transition is in the neighborhood of $x=x_2$ .

The third situation considered arises when following the injection of water, a small amount of oil is extracted and then there is a new interruption of flow, making it necessary to inject more saturated water after a small amount of oil is extracted. So,

$q\left(x,0\right)=\{\begin{array}{ll}1,\hfill & \text{if}{x}_1<x\hfill \\ 0,\hfill & \text{if}{x}_1<x<x_2\hfill \\ 1,\hfill & \text{if}x>x_2\hfill \end{array}\left(\text{scenario}3\right).$ (6)

Equation (6) is rewritten as:

$q_p\left(x,0\right)=1-q_b\left(x,0\right),$ (7)

therefore, there is a water-oil transition in $x=x_1$ and the oil-water transition in $x=x_2$ .

Considering the spatial domain $\left[x_i,x_f\right]$ , the boundary condition on the left was of the Dirichlet type and on the right of the radiation type, that is, when the dynamics reach the boundary, the flow simply goes through the boundary, and is not being affected by the edge.

Taking into account the effects of infiltration in porous media, van Duijin, Peletier and Pop (2007) [9] proposed the following modification to Equation (1):

$q_t+f{\left(q\right)}_x=\epsilon q_{xx}+{\epsilon }^2\kappa q_{xxt},$ (8)

in which $\epsilon$ is a diffusibility coefficient and $\kappa$ is a dispersive coefficient.

As one of the objectives of this work is to observe the effect of dispersion and diffusion in the Buckley-Leverett equation, we could not consider that $\epsilon \gg {\epsilon }^2\kappa$ , as in [5].

Thus, from the point of view of numerical methods, it is interesting to rewrite Equation (8) in the form,

${\left(q-{\epsilon }^2\kappa q_{xx}\right)}_t+f{\left(q\right)}_x=\epsilon q_{xx}$ (9)

and then,

$\{\begin{array}{l}{p}_t+f{\left(q\right)}_x=\epsilon q_{xx}\hfill \\ p=\left(q-{\epsilon }^2\kappa q_{xx}\right)\hfill \end{array}.$ (10)

The choice of numerical methods to solve Equation (10) is essential, as a numerical scheme that has a high degree of numerical dissipation can stand out and mask the real effect of the dispersive and diffusive terms of Equation (8), impairing the interpretation of the results.

3. Numerical Methods

Choosing a numerical method requires care that depends on the differential equation and the initial conditions imposed by the model under study. In problems containing hyperbolic partial differential equations, even initial conditions with smooth and sufficiently differentiable functions can evolve to discontinuous solutions. Furthermore, an inappropriate choice of numerical methods may result in excessive numerical dissipation and dispersion (oscillation), harming the quality of the numerical solution. Thus, for this work, a weighted essentially non-oscillatory scheme (WENO-5 method) was chosen for spatial discretization and the third-order Runge-Kutta TVD method (Total Variation Diminishing) was chosen for temporal discretization.

3.1. Weighted Non-Oscillatory Schemes (WENO-5)

The essentially non-oscillatory (ENO) schemes proposed by Harten & Osher (1987) [10] and Harten et al. (1987) [11] proved to be efficient in terms of decreasing numerical dissipation, in addition to avoiding oscillations in the discontinuity regions of the solutions. The main feature of these schemes is to choose the smoothest stencil among the candidates, in order to preserve high order of acurrate.

In general, a rth-order ENO scheme chooses the smoothest stencil among r possibilities and uses only the chosen one to approximate the flow [12]. The idea of weighted essentially non-oscillatory (WENO) schemes is to find a convex combination of all candidate stencils, for the numerical flow approximation. A weight is assigned to each stencil, representing its contribution to the process.

In smooth regions, the weights can be defined by certain optimal weights, while maintaining the high order of accuracy. In non-smooth regions, weights close to zero are assigned for stencils that contain discontinuities [12].

The WENO schemes technique is based on the flow version of the ENO schemes, considering a one-dimensional conservation law $u_t+f{\left(u\right)}_x=0$ . The spatial operator that approximates $-f{\left(u\right)}_x$ in $x_j$ is given by

$L=-\frac{1}{\text{Δ}x}\left({\widehat{f}}_{j+1/2}-{\widehat{f}}_{j-1/2}\right),$(11)

in which $\text{Δ}x$ is the size of the spatial discretization and ${\widehat{f}}_l$ is the numerical flux [12].

The r candidate stencils are denoted by $S_k$ , $k=0,1,\cdots ,r-1$ , of the form $S_k=\left\{x_{j+k-r+1},x_{j+k-r+2},\cdots ,x_{j+k}\right\}$ , which defines the amount of points of $S_k$ , which will be used to approximate the value of ${\widehat{f}}_{j+1/2}$ . Thus, a 3-order ENO scheme ($r=3$ ) has three candidate stencils, $S_k$ , $k=0,1,2$ , given by $S_k=\left\{x_{j+k-2},x_{j+k-1},x_j\right\}$ , which results in $S_0=\left(x_{j-2},x_{j-1},x_j\right)$ , $S_1=\left(x_{j-1},x_j,x_{j+1}\right)$ and $S_2=\left(x_j,x_{j+1},x_{j+2}\right)$ .

ENO schemes approximate ${\widehat{f}}_{j+1/2}$ through a polynomial interpolation at the points of each stencil [12] [13]. This approximation is given by

${\widehat{f}}_{j+1/2}=q_k^r\left(f_{j+k-r+1},\cdots ,f_{j+k}\right)\text{with}{q}_k^r\left(g_0,\cdots ,g_{r-1}\right)={\displaystyle \sum }_{l=0}^{r-1}{a}_{k,l}^r{g}_l.$(12)

Let $g\left(x\right)$ be a smooth function. The average approximation of $g\left(x\right)$ in cell $I_j$ is defined by

${\overline{g}}_j=\frac{1}{\text{Δ}{x}_j}{\displaystyle {\int }_{x_{j-1/2}}^{x_{j+1/2}}g\left(\xi \right)\text{d}\xi }\text{where}\xi \in I_j=\left(x_{j-1/2},x_{j+1/2}\right).$(13)

To obtain the constants $a_{k,l}^r$ in (12), consider the primitive function of $g\left(x\right)$ defined by $G\left(x\right)={\displaystyle {\int }_{-\infty }^{x}g\left(\xi \right)\text{d}\xi }$ . The value of $G\left(x_{j+1/2}\right)$ can be written as

$G\left(x_{j+1/2}\right)={\displaystyle \sum _{i=-\infty }^j{\displaystyle {\int }_{x_{i-1/2}}^{x_{i+1/2}}g\left(\xi \right)\text{d}\xi }}={\displaystyle \sum _{i=-\infty }^j{\overline{g}}_i\text{Δ}{x}_i}.$(14)

Equation (14) means that, once the average approximations of the cells ${\overline{g}}_i$ are known, the values of $G\left(x\right)$ at the boundary of cell $I_i$ are also known. So, the constants $a_{k,l}^r$ are determined by interpolating $G\left(x_{j+1/2}\right)$ by a r degree polynomial $P\left(x\right)$ , at most. Therefore,

${\widehat{f}}_{j+1/2}=p\left(x_{j+1/2}\right)=P^{\prime }\left(x_{j+1/2}\right)={\displaystyle \sum }_{l=0}^{r-1}{a}_{k,l}^r{g}_l,$(15)

in which $a_{k,l}^r$ are obtained from the Lagrange interpolating polynomial [13], with the data in Table 1.

Table 1. $a_{k,l}^r$ coefficients.

r	k	$l=0$	$l=1$	$l=2$
3	0	1/3	−7/6	11/6
	1	−1/6	5/6	1/3
	2	1/3	5/6	−1/6

A rth-order ENO scheme implies to a (2r − 1)th-order WENO scheme, so a 3rd-order ENO scheme leads to a 5th-order WENO scheme. In WENO schemes, for each candidate stencil $S_k$ , $k=0,1,\cdots ,r-1$ , a weight ${\omega }_k$ is assigned and these are used to calculate the numerical flux

${\widehat{f}}_{j+1/2}={\displaystyle \sum }_{k=0}^{r-1}{\omega }_k{q}_k^r\left(f_{j+k-r+1},\cdots ,f_{j+k}\right).$(16)

The weight ${\omega }_k$ for the stencil $S_k$ is defined by

${\omega }_k=\frac{{\alpha }_k}{{\alpha }_0+\cdots +{\alpha }_{r-1}}\text{with}{\alpha }_k=\frac{C_k^r}{{\left(\epsilon +I{S}_k\right)}^p},k=0,1,\cdots ,r-1.$ (17)

Taking $p=r$ , the coefficients $C_k^r$ are optimal values to determine ${\omega }_k$ [13]. The term $I{S}_k$ is an indicator of smoothness and for $r=3$ we have

$\begin{array}{l}I{S}_0=\frac{13}{12}{\left(f_{j-2}-2f_{j-1}+f_j\right)}^2+\frac{1}{4}{\left(f_{j-2}-4f_{j-1}+3f_j\right)}^2\\ I{S}_1=\frac{13}{12}{\left(f_{j-1}-2f_j+f_{j+1}\right)}^2+\frac{1}{4}{\left(f_{j-1}-f_{j+1}\right)}^2\\ I{S}_2=\frac{13}{12}{\left(f_j-2f_{j+1}+f_{j+2}\right)}^2+\frac{1}{4}{\left(3f_j-4f_{j+1}+f_{j+2}\right)}^2\end{array}$

This measure was introduced by Jiang and Shu (1996) [12], with the aim of achieving high accurate for the case where $r=3$ . Note that as $I{S}_k$ increases, the smoothness decreases and, consequently, ${\alpha }_k$ becomes close to zero as does ${\omega }_k$ , meaning that a weight close to of zero will be assigned to non-smooth solutions.

3.2. Third-Order Runge-Kutta (RK3-TVD)

Once the spatial discretization is concluded, a method for temporal discretization that maintains the non-oscillatory characteristics achieved is necessary.

Numerical methods belonging to the TVD class (Total Variation Diminishing) have the property of avoiding oscillations that are not typical of the phenomenon under study [8]. A good alternative is the high-order Runge-Kutta TVD methods, which were developed by Shu and Osher (1988) [14] in research related to efficient implementations for ENO’s schemes.

A method is called Total Variation Diminishing (TVD) if, for any data set $U^n$ , the values $U^{n+1}$ computed by the method satisfy $TV\left(U^{n+1}\right)\le TV\left(U^n\right)$ , where

$TV\left(U^n\right)={\displaystyle \sum }_{i=1}^N\left|U_i^n-U_{i-1}^n\right|$ (18)

is the total variation. In this work, the third-order Runge-Kutta TVD (RK3-TVD) method was chosen, whose expressions are given by

$\begin{array}{l}{u}^{\left(1\right)}=u^n+\text{Δ}tL\left(u^n\right)\\ u^{\left(2\right)}=\frac{3}{4}{u}^n+\frac{1}{4}{u}^{\left(1\right)}+\frac{1}{4}\text{Δ}tL\left(u^{\left(1\right)}\right)\\ u^{n+1}=\frac{1}{3}{u}^n+\frac{2}{3}{u}^{\left(2\right)}+\frac{2}{3}\text{Δ}tL\left(u^{\left(2\right)}\right)\end{array}$ (19)

in which L is the spatial operador of differential equation.

The WENO-5 and RK3-TVD methods numerically solve the classical Buckley-Leverett equation. To discretize the diffusive and dispersive terms incorporated in the temporal part of the modified Buckley-Leverett equation, Equation (10), we use the fourth-order centered finite difference scheme.

3.3. Fourth-Order Central Finite Difference Scheme (CFDS-4)

There are at least two ways to numerically add the U term in the conservation law: one is to discretize the diffusive term in the flux context [15],

$q_{xx}\left(x_{i+1/2},t\right)\cong \frac{{\overline{q}}_{i+3/2}-2{\overline{q}}_{i+1/2}+{\overline{q}}_{i-1/2}}{\text{Δ}{x}^2},$

and the other in the finite difference context [16],

$q_{xx}\left(x_i,t\right)\cong \frac{-{\overline{q}}_{i-2}+16{\overline{q}}_{i-1}-30{\overline{q}}_i+16{\overline{q}}_{i+1}-{\overline{q}}_{i+2}}{12\left(\text{Δ}{x}^2\right)},$(20)

In this work, both discretizations were performed, as for the interpretation of the results there were no significant differences between them, we chose to keep the fourth-order central finite difference scheme (CFDS-4) in our codes.

This method was also used to discretize the second expression in Equation (10), differently of Wang’s work [16] that use a staggered mesh and numerical integration to incorporate such equation into WENO scheme. In this way, the dispersive term enters the temporal evolution by replacing the term $u^n$ of the RK3-TVD method by $v^n=u^n-{\epsilon }^2\kappa q_{xx}\left(x_i,t^n\right)$ , in which $q_{xx}\left(x_i,t^n\right)$ is defined by Equation (20). Details about stability of those methods can be founded in a study of stability analysis presented in [5].

Figure 1 contains a schematic structure of the adopted methods.

4. Simulations and Results

The simulations were performed by initial conditions defined in Section 2. For all scenarios, we have $x\in \left[-1,1\right]$ with 128 subintervals, $\text{Δ}x=1/64$ , and $\text{Δ}t=0.1\text{Δ}x$ . Furthermore, $a=0.5$ was assigned to the constant that characterizes the porous medium, Equation (1). All codes were written by the authors themselves in Octave, and the discretizations used satisfy the stability criteria described in [5].

Figure 1. Flowchart. Source: The authors.

In all simulations, the fluids mix is representative when q assume values between $q=1$ (pure water) and $q=0$ (pure petroleum). While the fluid dynamics develops to on the right, the time evolutions reveal the emergence of a region where the fluids mix, between water and oil.

4.1. Initial Conditions

In the simulations of this section the following initial conditions were adopted

$q_b\left(x,0\right)=\left[\frac{1+\tanh \left(50\left(x+0.5\right)\right)}{2}\right]+\left[\frac{1-\tanh \left(50\left(x-0.0\right)\right)}{2}\right]-1,$ (21)

for barrier (Figure 2(a)) and

$q_p\left(x,0\right)=1-q_b\left(x,0\right),$ (22)

(Figure 2(b)).

In Figure 2(a) we have oil in the region $\left[-1,-0.5\right)$ , water in the region $\left(-0.5,0.0\right)$ and oil in $\left(0.0,1\right]$ . Thus, from left to right, we have an oil-water transition in the vicinity of $x=-0.5$ and a water-oil transition close to $x=0.0$ . And, in Figure 2(b) the situation is reversed, we have a water-oil transition close to $x=-0.5$ and an oil-water transition at $x=0.0$ .

(a)

(b)

Figure 2. Initial Conditions. (a) scenario 2; (b) scenario 3. Source: The authors.

The first two simulations consider the temporal evolution of the initial conditions (Figure 2(a) and Figure 2(b)) via the classical Buckley-Leverett equation, that is, $\epsilon =\kappa =0$ , which results in Equation (1).

4.2. Standard Simulation 1

In this scenario, we started with the initial condition being a barrier function, Equation (21) and Figure 2(a).

After $n=64$ iterations, the emergence of two mixing profiles between the fluids is noted, one with linear behavior starting at $x=-0.5$ and the other with non-linear behavior at $x=0.0$ (Figure 3(a)).

(a)

(b)

Figure 3. Numerical solution after 64 and 128 iterations. (a) n = 64; (a) n = 64. Source: The authors.

Another interesting point to highlight is that these characteristics will remain in simulations, that is, when we have a transition from left to right between oil and water, the mixture evolves with linear behavior and, when it is a transition between water and oil, the evolution of the mixture will have non-linear behavior.

In Figure 3(b), simulation for $n=128$ , the mixing initiated at $x=-0.5$ reaches the vicinity of $x=0.0$ , leaving only a small region with only water.

(a)

(b)

Figure 4. Numerical solution after 180 and 250 iterations. (a) n = 180; (b) n = 250. Source: The authors.

With $n=180$ , Figure 4(a), there is no longer a region with $q=1$ (only water) and the transition that evolved from $x=0.0$ reaches the neighborhood of $x=1.0$ . In Figure 4(b), the fluid dynamics crosses the boundary without adding errors to the solution. Note that this fact indicates that the computational program and the method are adequate for the chosen boundary condition.

4.3. Standard Simulation 2

In this scenario, we started with the initial condition being a well function, Equation (22) and Figure 2(b).

(a)

(b)

Figure 5. Numerical solution after 64 and 128 iterations. (a) n = 180; (b) n = 250. Source: The authors.

Similar to what occurred in Subsection 4.3, we obtained two mixing profiles between the fluids (Figure 5). However, in this scenario, the nonlinear profile is at $x=-0.5$ and the linear one at $x=0.0$ . Also, the only region initially containing oil is $\left(-0.5,0.0\right)$ .

With $n=64$ , Figure 5(a), it is still possible to see a small region with only oil ($q=0$ ). However, in Figure 5(b), it can be seen that the non-linear transition reaches the linear transition and, then, we no longer have a region with only oil.

The simulation continues for the values of $n=180$ and $n=250$ (Figure 6(a) and Figure 6(b)), respectively.

(a)

(b)

Figure 6. Numerical solution after 180 and 250 iterations. (a) n = 180; (b) n = 250. Source: The authors.

These standard simulations are similar to those presented in [5] and have served as references for the discussions in Subsections 4.4 and 4.5.

4.4. New Simulations 1

In this subsection, we consider again the initial condition, Equation (21), represented by Figure 2(a), but the numerical solutions obtained are for the modified Buckley-Leverett equation, Equation (8).

(a)

(b)

Figure 7. Numerical solution after 64 and 128 iterations. (a) n = 64; (b) n = 128. Source: The authors.

The simulations in this section considered $\epsilon =0.04$ for the diffusive term, adopting the values of $\kappa =0.00$ , $\kappa =0.50$ , $\kappa =0.70$ , $\kappa =0.90$ , $\kappa =0.95$ , $\kappa =0.97$ and $\kappa =0.99$ for the dispersive term.

With $\epsilon =0.04$ and $\kappa =0.00$ , numerical solutions in blue in Figure 7 and Figure 8, the action of the diffusive term can be seen. The transitions between oil-water (proximities of $x=-0.5$ ) and water-oil (proximities of $x=0.0$ ) are smoothed, representing that the mixtures between the fluids are more spread out when compared to the standard simulations (Figure 3 and Figure 4).

(a)

(b)

Figure 8. Numerical solution after 180 and 250 iterations. (a) n = 180; (b) n = 250. Source: The authors.

In Figure 7(a), by increasing the $\kappa$ value to 0.50 (red color) and 0.70 (yellow color), the dispersive term starts to inhibit the diffusive effect, leaving the transitions less spread out.

When the values of $k=0.90$ (green color), $k=0.95$ (magenta color), $k=0.97$ (cyan color) and $k=0.99$ (black color) are adopted (Figure 7(a)), dispersion inhibits diffusion even more, making the solution closer to the classical solution (Figure 3(a)). However, new phenomena arise in the transitions between fluids, the dispersive term starts to add oscillatory effects in the transitions between oil-water and water-oil.

The continuation of these simulations, Figure 7(b), Figure 8(a) and Figure 8(b), allows us to visualize the increase in the emergence of oscillations.

For $k=0.99$ the oscillations are more present, allowing us to state that the higher the value of k, the more intense the oscillatory behavior.

Therefore, in these simulations we identified that the dispersive term is what actually aggregates oscillatory phenomena, in two-phase fluid dynamics via Buckley-Leverett equation, Equation (8). Such oscillations are in some ways similar to those presented in [7], whose work had a different focus, but which took into account the issue of infiltration. On the other hand, [17] studied two-phase fluids incorporating several physical parameters (with the exception of infiltration) and the results did not present oscillations.

4.5. New Simulations 2

In this Subsection we consider again the initial condition, Equation (22), represented by Figure 2(b), but the numerical solutions obtained are for the modified Buckley-Leverett equation, Equation (8).

(a)

(b)

Figure 9. Numerical solution after 64 and 128 iterations. (a) n = 64; (b) n = 128. Source: The authors.

The value for the diffusive term and the values for the dispersive term are the same as those defined in Subsection 4.4, that is, $\epsilon =0.04$ for the diffusive term, and $\kappa =0.00$ , $\kappa =0.50$ , $\kappa =0.70$ , $\kappa =0.90$ , $\kappa =0.95$ , $\kappa =0.97$ and $\kappa =0.99$ for the dispersive term.

The simulations for $\epsilon =0.04$ and $\kappa =0.00$ highlight the diffusive character, smoothing the transitions between the fluids, Figure 9 and Figure 10, graphs in blue color.

(a)

(b)

Figure 10. Numerical solution after 180 and 250 iterations. (a) n = 180; (b) n = 250. Source: The authors.

In this temporal evolution, the inhibition of the diffusive term by the dispersive term can also be seen, by increasing the value of $\kappa$ , Figure 9 and Figure 10.

In Figure 10(a) and Figure 10(b), two distinct oscillations can be seen, one with a larger amplitude appearing in the vicinity of $x=-0.5$ and one with a smaller amplitude close to $x=0.0$ . Both oscillations propagate and spread along the pipe.

It is important to highlight the performance of the combination between the WENO-5 and RK3-TVD methods, which during the temporal evolution of the transitions between oil-water and water-oil, remained stable, providing solutions sufficiently appropriate to represent the diffusive and dispersive (oscillatory) phenomena.

5. Conclusions

In this study, the WENO-5 and RK3-TVD methods were addressed, together with a central finite difference scheme. The methods were applied to the classical and modified Buckley-Leverett Equations, in order to investigate the temporal evolution of scenarios that present transitions between oil-water and water-oil under two different initial conditions.

First at all, standard scenarios with the classical Buckley-Leverett Equations were developed, and then, compared with the solutions obtained via modified Buckley-Leverett Equations.

By adopting the values of $\epsilon =0.04$ and $\kappa =0.00$ , the diffusion phenomenon becomes evident. As the value of k increases, the diffusive term is inhibited and the dispersive characteristics are visually evidenced by the oscillations.

Thus, the insertion of the diffusive term ($\epsilon$ parameter) to the classical equation adds the characteristic of smoothing the solutions, spreading the transitions between oil-water and water-oil. Finally, the dispersive term (k parameter) has the property of incorporating oscillations, representing the inhomogeneity in the absorption of fluids by the porous medium.

In future work, we intend to carry out simulations in three-dimensional pipelines, inserted in an environment that influences the dynamics of two-phase fluid flow.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Nikolova, C. and Gutierrez, T. (2020) Use of Microorganisms in the Recovery of Oil from Recalcitrant Oil Reservoirs: Current State of Knowledge, Technological Advances and Future Perspectives. Frontiers in Microbiology, 10, Article 2996. https://doi.org/10.3389/fmicb.2019.02996
[2]	Kou, W. (2022) Global Existence of Solutions for Baer-Nunziato Two-Phase Flow Model in a Bounded Domain. Open Journal of Applied Sciences, 12, 631-649. https://doi.org/10.4236/ojapps.2022.124043
[3]	Khan, M.Y. and Mandal, A. (2021) Improvement of Buckley-Leverett Equation and Its Solution for Gas Displacement with Viscous Fingering and Gravity Effects at Constant Pressure for Inclined Stratified Heterogeneous Reservoir. Fuel, 285, Article ID: 119172. https://doi.org/10.1016/j.fuel.2020.119172
[4]	Buckley, S.E. and Leverett, M.C. (1942) Mechanism of Fluid Displacement in Sands. Transactions of the AIME, 146, 107-116. https://doi.org/10.2118/942107-g
[5]	Garcia, R. O. and Silveira, G. P. (2024) Essentially Non-Oscillatory Schemes Applied to Buckley-Leverett Equation with Diffusive Term. Latin-American Journal of Computing, 11, 44-54. https://zenodo.org/records/10402169
[6]	DiCarlo, D.A. (2004) Experimental Measurements of Saturation Overshoot on Infiltration. Water Resources Research, 40, 4215.1-4215.9. https://doi.org/10.1029/2003wr002670
[7]	Wang, K., Liu, Z. and Yang, J. (2021) Numerical Study of Gas-Liquid Two-Phase Flow in PEMFC Cathode Flow Channel with Different Diffusion Layer Surface Structure. Journal of Power and Energy Engineering, 9, 106-117. https://doi.org/10.4236/jpee.2021.911006
[8]	LeVeque, R.J. (2002) Finite Volume Methods for Hyperbolic Problems. Cambridge University Press. https://doi.org/10.1017/cbo9780511791253
[9]	van Duijn, C.J., Peletier, L.A. and Pop, I.S. (2007) A New Class of Entropy Solutions of the Buckley-Leverett Equation. SIAM Journal on Mathematical Analysis, 39, 507-536. https://doi.org/10.1137/05064518x
[10]	Harten, A. and Osher, S. (1987) Uniformly High-Order Accurate Nonoscillatory Schemes. I. SIAM Journal on Numerical Analysis, 24, 279-309. https://doi.org/10.1137/0724022
[11]	Harten, A., Engquist, B., Osher, S. and Chakravarthy, S.R. (1987) Uniformly High Order Accurate Essentially Non-Oscillatory Schemes, III. Journal of Computational Physics, 71, 231-303. https://doi.org/10.1016/0021-9991(87)90031-3
[12]	Jiang, G. and Shu, C. (1996) Efficient Implementation of Weighted ENO Schemes. Journal of Computational Physics, 126, 202-228. https://doi.org/10.1006/jcph.1996.0130
[13]	Gerolymos, G.A., Sénéchal, D. and Vallet, I. (2009) Very-High-Order Weno Schemes. Journal of Computational Physics, 228, 8481-8524. https://doi.org/10.1016/j.jcp.2009.07.039
[14]	Shu, C. and Osher, S. (1988) Efficient Implementation of Essentially Non-Oscillatory Shock-Capturing Schemes. Journal of Computational Physics, 77, 439-471. https://doi.org/10.1016/0021-9991(88)90177-5
[15]	Silveira, G.P. and de Barros, L.C. (2013) Numerical Methods Integrated with Fuzzy Logic and Stochastic Method for Solving PDEs: An Application to Dengue. Fuzzy Sets and Systems, 225, 39-57. https://doi.org/10.1016/j.fss.2013.04.003
[16]	Wang, Y. and Kao, C. (2013) Central Schemes for the Modified Buckley-Leverett Equation. Journal of Computational Science, 4, 12-23. https://doi.org/10.1016/j.jocs.2012.02.001
[17]	Nyariki, E.M., Kinyanjui, M.N. and Abonyo, J.O. (2023) Heat and Mass Transfers in a Two-Phase Stratified Turbulent Fluid Flow in a Geothermal Pipe with Chemical Reaction. Journal of Applied Mathematics and Physics, 11, 484-513. https://doi.org/10.4236/jamp.2023.112030