Central schemes for porous media flows

Central schemes for porous media flows

E. Abreu^I; F. Pereira^II; S. Ribeiro^III

^IInstituto Nacional de Matemática Pura e Aplicada, 22460-320 Rio de Janeiro, RJ, Brazil E-mail: [e-mail]

^IIDepartment of Mathematics and School of Energy Resources, University of Wyoming, Laramie, WY 82071-3036, USA E-mail: [e-mail]

^IIIDepartamento de Ciências Exatas, Universidade Federal Fluminense, 27255-250 Volta Redonda, RJ, Brazil E-mail: [e-mail]

ABSTRACT

We are concerned with central differencing schemes for solving scalar hyperbolic conservation laws arising in the simulation of multiphase flows in heterogeneous porous media. We compare the Kurganov-Tadmor (KT) [3] semi-discrete central scheme with the Nessyahu-Tadmor (NT) [27] central scheme. The KT scheme uses more precise information about the local speeds of propagation together with integration over nonuniform control volumes, which contain the Riemann fans. These methods can accurately resolve sharp fronts in the fluid saturations without introducing spurious oscillations or excessive numerical diffusion. We first discuss the coupling of these methods with velocity fields approximated by mixed finite elements. Then, numerical simulations are presented for two-phase, two-dimensional flow problems in multi-scale heterogeneous petroleum reservoirs. We find the KT scheme to be considerably less diffusive, particularly in the presence of high permeability flow channels, which lead to strong restrictions on the time step selection; however, the KT scheme may produce incorrect boundary behavior.

Mathematical subject classification: Primary: 35L65; Secondary: 65M06.

Key words: hyperbolic conservation laws, central differencing, two-phase flows.

1 Introduction

We are concerned with high resolution central schemes for solving scalar hyperbolic conservation laws arising in the simulation of multiphase flows in multidimensional heterogeneous petroleum reservoirs.

Many of the modern high resolution approximations for nonlinear conservation laws employ Godunov's appoach [35] or REA (reconstruct, evolve, average) algorithm, i.e., the approximate solution is represented by a piecewise polynomial which is Reconstructed from the Evolving cell Averages. The two main classes of Godunov methods are upwind and central schemes.

The Lax-Friedrichs (LxF) scheme [33] is the canonical first order central scheme, which is the forerunner of all central differencing schemes. It is based on piecewise constant approximate solutions. It also enjoys simplicity, i.e., it does not employ Riemann solvers and characteristic decomposition. Unfortunately the excessive numerical dissipation in the LxF recipe (of order ((ΔX)²/Δt)) yields poor resolution, which seems to have delayed the development of high resolution central schemes when compared with the earlier developments of the high resolution upwind methods. Only in 1990 a second order generalization to the LxF scheme was introduced by Nessyahu and Tadmor (NT) [27].They used a staggered form of the LxF scheme and replaced the first order piecewise constant solution with a van Leer's MUSCL-type piecewise linear second order approximation [8]. The numerical dissipation in this new central scheme has an amplitude of order ((ΔX)⁴/Δt). When applying these methods to multiphase flows in highly heterogeneous petroleum reservoirs or aquifers we need to use decreasing time steps as the heterogeneity increases, yielding greater numerical diffusion. Kurganov and Tadmor (KT) [3] combined ideas from the construction of the NT scheme with Rusanov's method [36] to obtain the first second order central scheme that admits a semi-discrete formulation which is then solved with an appropriate ODE solver. The resulting scheme has a much smaller numerical diffusion than the NT scheme. Due to the semi-discrete formulation, this numerical diffusion is independent of the time step used to evolve the ordinary differential equation. This property guarantees that no extra numerical diffusion will be added if the time step is forced to decrease. For this reason, the application of this central scheme results in a new numerical approach to model two-phase flows with a much lower numerical diffusion, even in the presence of a highly heterogeneous porous media.

The goals of this paper are (i) to discuss the coupling of NT and KT schemes to velocity fields approximated by Raviart-Thomas mixed finite element method (See [30]), and (ii) to compare the KT semi-discrete central scheme with the NT central scheme for numerical simulations of two-phase, incompressible, two-dimensional flows in heterogeneous formations. Both methods can accurately resolve sharp fronts in the fluid saturations without introducing spurious oscillations or excessive numerical diffusion.

Our numerical experiments indicate that the KT scheme is considerably less diffusive, particularly in the presence of viscous fingers, which lead to strong restrictions on the time step selection. On the other hand the KT scheme may produce incorrect boundary behavior in a typical two-dimensional geometry used in the study of porous media flows: the quarter of a five spot.

Numerous methods have been introduced to solve two-phase flow problems in porous media. Among eulerian-lagrangian procedures we mention the Modified Method of Characteristics [25, 29], the Modified Method of Characteristics with Adjusted Advection [22], the Locally Conservative Eulerian Lagrangian Method [23] and Eulerian Lagragian Localized Adjoint Methods [26]. Additional techniques, to name just a few, include higher-order Godunov schemes [10], the front-tracking method [7], the streamline method [32, 34] the streamline upwind Petrov-Galerkin method (SUPG) [1, 4] and a second-order TVD-type finite volume scheme [31] (this procedure aims at the modeling of flow through geometrically complex geological reservoirs). Each of these procedures has advantages and disadvantages. We refer the reader to [28, 23] and references cited there for a discussion of these methods.

We remark that central schemes are particularly interesting for the numerical simulation of multiphase flow problems in porous media because they have been formulated to solve hyperbolic systems; this is not the case for several of the procedures mentioned above, which have been developed only for scalar equations.

Moreover these central schemes were also used to deal with many other applied problems: to solve Hamilton-Jacobi Equations (see [11] and [2]), to model the two-dimensional magneto hydrodynamics (MHD) equations and to study the Orszag-Tang vortex system, which describes the transition to supersonic turbulence for the equations of MHD in two space dimensions (see [9] and [20]),to mention just a couple of them.

This paper is organized as follows. In Section 2, we discuss our strategy for solving numerically the model for two-phase, immiscible and incompressible displacement in heterogeneous porous media considered here. In Section 3, we discuss the application of central differencing schemes to porous media flows. In Section 4, we present the computational solutions for the model problem considered here and our conclusions.

2 Numerical approximation of two-phase flows

We consider a model for two-phase immiscible and incompressible displacement in heterogeneous porous media. The governing equations are strongly nonlinear and lead to shock formation, and with or without diffusive terms they are of practical importance in petroleum engineering [15, 24]. See also [21] and the references therein for recent studies for the scale-up problem for such equations. The conventional theoretical description of two-phase flow in a porous medium, in the limit of vanishing capillary pressure, is via Darcy's law coupled to the Buckley-Leverett equation. The two phases will be referred to as water and oil, and indicated by the subscripts w and o, respectively. Without sources or sinks and neglecting the effects of capillarity and gravity, these equations read (See [15] for more details)

Here, v is the total seepage velocity, s is the water saturation, K(x) is the absolute permeability, and p is the pressure. The constant porosity has been scaled out by a change of the time variable. The total mobility, λ(s), and water fractional flow function, f(s), are defined in terms of the relative permeabilities k_ri(s) and phase viscosities µ_i by

2.1 Operator splitting for two-phase flow

An operator splitting technique is employed for the computational solution of the saturation equation (2) and the pressure equation (1) in which they are solved sequentially with possibly distinct time steps. This splitting scheme has proved to be computationally efficient in producing accurate numerical solutions for two-phase flows. We refer the reader to [22] and references therein for more details on the operator splitting technique; see also [16, 14, 17, 18] and [5] for applications of this strategy to three phase flows taking into account capillary pressure (diffusive effects).

Typically, for computational efficiency larger time steps are used for thepressure-velocity calculation (Equation 1) than for the convection calculation (Equation 2). Thus, we introduce two time steps: Δt_c for the solution of the hyperbolic problem for convection, and Δt_p for the pressure-velocity calculation so that Δt_p > Δt_c. We remark that in practice variable time steps are always useful, especially for the convection micro-steps subject dynamically to a CFL condition.

For the pressure solution we use a (locally conservative) hybridized mixed finite element discretization equivalent to cell-centered finite differences [21, 22], which effectively treats the rapidly changing permeabilities that arise from stochastic geology and produces accurate velocity fields. The pressure and Darcy velocity are approximated at times t^m = mΔt_p, m = 0, 1, 2, ....The linear system resulting from the discretized equations is solved by a preconditioned conjugate gradient procedure (PCG) (See [22] and the references therein). The saturation equation is approximated at times = t^m + κ Δt_c for t^m < < t^m+1. We remark that we must specify the water saturation at t = 0.

3 Central differencing schemes for porous media flows

In this section, we shall study the family of high resolution, non-oscillatory, conservative central differencing schemes introduced by Nessyahu and Tadmor (NT) and Kurganov and Tadmor (KT). They will be applied to the numerical approximation of the scalar hyperbolic conservation law modeling the convective transport of fluid phases in two-phase flows. For the associated elliptic problem (Eq. (1)), we use the lowest order Raviart-Thomas [30] locally conservative mixed finite elements. These central schemes enjoy the main advantage of Godunov-type central schemes: simplicity, i.e., they employ neither characteristic decomposition nor approximate Riemann solvers. This makes them universal methods that can be applied to a wide variety of physical problems, including hyperbolic systems. In the following sections we will discuss the main ideas of the NT and KT central schemes coupled to the mixed finite element discretization mentioned above. We will not repeat here all the details involved in the development of the NT and KT schemes; instead, we refer the reader to [27] and [3] for this material.

3.1 The Nessyahu-Tadmor scheme for two-phase flows

Consider the following scalar hyperbolic conservation law,

subject to prescribed initial data, s(x, y, 0) = S₀(x, y). Here ^xv = ^xv(x, y, t) and ^yv = ^yv(x, y, t) denote the x- and y-components of the velocity field v (see Eq. 1). To approximate (3) by the NT scheme, we begin with a piecewise constant solution of the form

is the approximate cell average at t = associated with the cell C_j,k = I_j × I_k = [x_j-1/2, x_j+1/2] × [y_k-1/2, y_k+1/2] and c_j,k(x, y) is a characteristic function of the cell C_j,k.

We first reconstruct a piecewise linear approximation of the form

In Eq. (4), the discrete slopes along the x and y directions satisfy

to guarantee second-order accuracy.

The reconstruction (4) retains conservation, i.e.:

Let {s(x, y, t), t > } be the exact solution of the conservation law (3), subject to the reconstructed piecewise-linear data (4) at time t = . The evolution step in the NT scheme consists of approximating this exact solution at the next time step t = + Δt_c, by its averages over staggered cells, C_j+1/2,k+1/2 : = I_j+1/2 × I_k+1/2. See dashed grid in Figure 1 (denote κ + 1 : = + Δt_c).

Let

These new staggered cell averages are obtained by integrating the conservation law (3) over the control volumes C_j+1/2,k+1/2 × [, + Δt_c] following the same manipulations as described in [19] (denote α_x≡ and α_y≡ ):

The cell average

has contributions from the four cells C_j,k, C_j+1,k, C_j+1,k+1, and C_j,k+1:

Computing these integrals exactly yields

To approximate the four flux integrals on the right hand side of (8), we use the second-order rectangular quadrature rule for the spatial integration and the midpoint quadrature rule for second-order approximation of the temporal integrals. For instance, letting κ + 1/2 be + Δt_c/2,

Since these midvalues are computed at the center of the cells, C_j,k, where the solution is smooth, provided an appropriate CFL condition is observed, we can use Taylor expansion together with the conservation law (3) to get

Here, and are one-dimensional discrete slopes in the x and y directions, respectively. They satisfy the conditions

in order to produce a second order scheme for the approximation of (3). To avoid spurious oscillations, it is essential to reconstruct the discrete derivatives given by Equations (5) and (13) with built-in nonlinear limiters. In this work we use the following MinMod limiter

where Δ is the centered difference, We refer the reader to [27] and [3] and the references therein for the various options for the form of such discrete derivatives.

In our sequential scheme, when solving for the saturation in time, the total velocity v is given by the solution of the velocity-pressure equation. Recall that the solution of Eq. (1) is approximated the lowest order Raviart-Thomas mixed finite element method. Thus, the computed total velocity v is discontinuous at the vertices of the original non-staggered grid cells. This constitutes a difficulty for the staggered scheme (8), which requires the values of the total velocity v at these vertices at every other time step. To avoid this difficulty we use the non-staggered version of the NT scheme.

To turn the staggered scheme (8) into a non-staggered scheme, we re-average the reconstructed values of the underlying staggered scheme, thus recovering the cell averages of the central scheme over the original non-staggered grid cells. First we reconstruct a piecewise bilinear interpolant at the time step κ + 1 := + Δt_c

as in Equation (4), through the staggered cell averages given by (8), and re-average it over the original grid cells, giving the following non-staggered scheme:

3.2 The Kurganov-Tadmor scheme for two-phase flows

The first multidimensional extension of the KT scheme was presented in [3]. This extension used the dimension by dimension approach, that is, the numerical fluxes computed along the x and y directions are viewed as a generalization of the one-spatial-dimension numerical fluxes. This approach consists of the following steps: at each time step and at each cell I_j,k,

(i) Compute the difference of the one-dimensional numerical flux in one spatial dimension in the x direction keeping y constant and equal to y_k. Denote this difference by

The numerical flux is

where

are the corresponding right and left intermediate values of (x, t^κ) at (x_j+1/2, y_k).

The local speed of wave propagation (t) is estimated at the cell boundaries (x_j+1/2, y_k) as the upper bound

where ω is a value between (t) and (t). The velocity field used in the KT scheme is obtained directly from the Raviart-Thomas space on the cell edges:

^xv_j+1/2,k(t) := (v_r)_jk(t),

^xv_j-1/2,k(t) := (v_l)_jk(t),

where v_r and v_l stand for the velocity on the one-dimensional "right" and "left" faces of the cells.

(ii) Analogously, compute the difference of the one-dimensional numerical flux in the y direction keeping x constant and equal to x_j. This difference is denoted by

The one dimensional numerical flux in the y direction is

In a similar way, the correspondent "up" and "down" intermediate values of (x, t^κ) at (x_j, y_k+1/2) are

The local speed of wave propagation (t) in the y direction is estimated at the cell boundaries (x_j, y_k+1/2) as the upper bound

where ω is a value between (t) and (t). Analogously the velocity field in the y direction is obtained directly from the Raviart-Thomas space on the cell edges:

^xv_j,k+1/2,k(t) := (v_u)_jk(t),

^xv_j-1/2,k(t) := (v_l)_jk(t),

where v_u and v_d stand for the velocity on the "upper" and "lower" faces of the cells.

(iii) The cell average in the next time step + Δt_c is then the solution of the following differential equation

The numerical derivatives are computed using the MinMod limiter given by Equation (14). In our numerical experiments, the parameter θ assumes values 1 < θ < 1.8.

The two-dimensional semi-discrete formulation (21) comprises a system of nonlinear ordinary differential equations for the discrete unknows {S_j,k(t)}. To solve it, we integrate in time introducing a variable time step Δtⁿ. Although the forward Euler scheme can be used, it may be advantageous to use higher order discretizations in numerical simulations. The numerical examples presented below use third-order Runge-Kutta ODE solvers based on convex combinations of forward Euler steps. See [12] and [13] for more details on a whole family of such schemes.

4 Two-dimensional numerical experiments

We present and compare the results of numerical simulations of two-dimensional, two-phase flows associated with two distinct flooding problems using the KT and NT schemes.

In all simulations, the reservoir contains initially 79% of oil and 21% of water. Water is injected at a constant rate of 0.2 pore volumes every year. The viscosity of oil and water used are µ_o = 10.0 cP and µ_w = 0.05 cP. The relative permeabilities are assumed to be:

K_ro(s) = (1 – (1 – s_ro)^-1s)²,

k_rw(s) = (1 – s_rw)^-2 (s – s_rw)²,

where s_ro = 0.15 and s_rw = 0.2 are the residual oil and water saturations, respectively.

For the heterogeneous reservoir studies we consider a scalar absolute permeability field K(x) taken to be log-normal (a fractal field, see [6] and [21] for more details) with moderately large heterogeneity strength. The spatially variable permeability field is defined on a 256 × 64 grid with three different values of the coefficient of variation CV (CV = 0.5, 1.2, 2.2) given by the ratio between the standard deviation and the mean value of the permeability field.

We now discuss the simulations in the slab geometry. We consider two-dimensional flows in a rectangular, heterogeneous reservoir (slab geometry) having 256m × 64m with three different sizes of computational grid: 256 × 64, 512 × 128 and 1024 × 256 cells. The boundary conditions and injection and production specifications for the two-phase flow equations (1)-(2) are as follows. The injection is made uniformly along the left edge of the reservoir and the production is taken along the right edge; no flow is allowed along the edges appearing at the top and bottom of the reservoir.

Figures 3, 4, and ⁵ refer to a comparative study of the NT and KT schemes, showing the water saturation surface plots after 275, 250 and 225 days of simulation for the three different CV values (CV = 0.5, 1.2, and 2.2). The results obtained with the NT scheme were computed using three computational grid: the coarsest grid with 256 × 64 cells, and two levels of refinement denoted by NTr and NTrr with 512 × 128 and 1024 × 256 cells, respectively (See the first three pictures from top to bottom of Figures 3, 4, and ⁵). At the same time, the bottom pictures in those figures are the results presented by the KT dimension by dimension scheme on the coarsest computational grid of 256 × 64 cells.

For each heterogeneity, we computed the difference between the results produced by the NT scheme in the three computational grids with the corresponding result produced by the KT scheme in the coarsest grid. We consider the solution of the NT scheme in the finer grid as the reference solution. The differences are then computed using the L² norm relative to this reference solution as follows

Here F stands for the KT solution and G for a NT solution. The graph in Figure 2 shows these differences. Note that as we refine the solution of the NT scheme, the differences become smaller indicating a comparable accuracy for the simulations produced by the KT scheme in the coarsest grid and the result produced by the NT scheme in the finest grid. These differences indicate that one has to refine twice the grid used in the NT scheme to produce an equivalent solution to the one produced by the KT scheme using the coarsest grid.

We now turn to the discussion of the set of simulations performed in a five-spot pattern. In case of a five-spot flood discretized by a diagonal grid (Figure 6), injection takes place at one corner and production at the diametrically opposite corner; no flow is allowed across the entirety of the boundary. In case of a five-spot flood discretized by a parallel grid (Figure 7), injection takes place at two opposite corners (say, bottom left and top right), and production is through the remaining two corners (say, bottom right and top left). Figures 6 (diagonal grid) and 7 (parallel grid) show the saturation level curves after 260 days of simulation obtained with the NT and KT schemes for two levels of spatial discretization.

In both Figures 6 and 7, the pictures on the left column are the results obtained with the NT scheme and the ones on the right were computed with the KT scheme. In these Figures, the grids are refined from top to bottom and have 64 × 64 and 128 × 128 cells in the diagonal pattern and 90 × 90 and 180 × 180 cells in the parallel grid.

It is clear that the KT scheme (right column pictures in Figure 6 and in Figure 7) is producing incorrect boundary behavior. Moreover, as the computational grid is refined (right column and bottom picture in Figures 6 and 7) this problem seems to be emphasized.

5 Conclusions

In one spatial dimension the KT scheme is a small modification of the NT scheme which uses more precise information about the local speed of propagation. This approach leads to a very simple numerical recipe producing numerical solutions more accurate then those provided by the NT scheme. On the one hand, in two spatial dimensions the KT scheme uses numerical fluxes in the x and y directions that can be viewed as generalizations of one-dimensional numerical fluxes. This is called the dimension by dimension approach. The NT, in the other hand, uses a genuinely two-dimensional configuration. In the case of a slab geometry, the fluid flows mostly in one direction. For this reason, this flow may be viewed as a one-dimensional flow and the KT scheme is expected to produce very accurate solutions like those presented in Figures 3, 4 and ⁵. In the five-spot problem, the fluid flows in both x- and y-directions, causing a genuinely two-dimensional displacement. The KT scheme produces incorrect boundary behaviors in the five-spot numerical examples. We remark that this incorrect behavior is not present in the results produced by the NT scheme. The dimension by dimension approach of the KT scheme might be a source of numerical errors for a class of problems with an intrinsic two-dimensional geometry. These numerical errors may lead to incorrect behavior like those in the five-spot problem. The authors are currently working on an improvement of these schemes in order to compute more precisely a genuinely two-dimensional numerical flux.

Acknowledgements. The authors wish to thank F. Furtado (University of Wyoming, USA) for several enlightening discussions and many suggestions during the preparation of this work. S. Ribeiro wishes to thank CAPES/Brazil for a Doctoral fellowship and also F. Pereira and the University of Wyoming for supporting a Pos-Doctoral study. E. Abreu thanks the financial support from CNPq through grants 150114/2007-9 and 473213/2007-9.

Received: 10/VI/08. Accepted: 24/XII/08.

#CAM-11/08.

• [1] A.F.D. Loula, A.L.G.A. Coutinho, E.L.M. Garcia and J.L.D. Alves, Solution of miscible and immiscible flows employing element-by-element iterative strategies. 3rd SPE Latin American/Caribean Engin. Conference, (1994), 431-444.
• [2] A. Kurganov and E. Tadmor, New high-resolution semi-discrete central scheme for Hamilton-Jacobi Equations. Journal of Computational Physics, 160 (2000), 720-742.
• [3] A. Kurganov and E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations. Journal of Computational Physics, 160 (2000), 241-282.
• [4] A.L.G.A. Coutinho and I.D Parsons, Finite element multigrid methods for two-phase immiscible flow in heterogeneous media. Communications in Numerical Methods in Engineering, 15 (1999), 1-7.
• [5] A.S. Francisco, F. Pereira, H.P. Amaral Souto and J. Aquino, A two-stage operator splitting algorithm for the numerical simulation of contaminant transport in porous media. International Journal of Computational Science, 2, 3 (2008) 422-436.
• [6] B. Lindquist, F. Pereira, J. Glimm and Q. Zhang, A theory of macrodispersion for the scale up problem. Transport in Porous Media, 13 (1993), 97-122.
• [7] B. Lindquist, J. Glimm, L. Padmanhaban and O. McBryan, A front tracking reservoir simulator: 5-spot validation studies and the water coning problem. Frontiers in Applied Mathematics, H.T. Banks, 1 (1983).
• [8] B. Van Leer, Towards the ultimate conservative difference scheme. v. a second order sequel to Godunov's method. Journal of Computational Physics, 32 (1979), 101-136.
• [9] C.C. Wu, E. Tadmor and J. Balbas, Non-oscillatory central schemes for one- and two-dimensional MHD Equations. Journal of Computational Physics, 201 (2004), 261-285.
• [10] C.N. Dawson, G.B. Shubin and J.B. Bell, An unsplit high-order Godunov scheme for scalar conservation laws in two dimensions. Journal of Computational Physics, 75 (1988), 1-24.
• [11] C.T. Lin and E. Tadmor, High-resolution non-oscillatory central scheme for Hamilton-Jacobi Equations. SIAM Journal on Scientific Computations, 21 (2000), 2163-2186.
• [12] C.W. Shu, Total variation diminishing time discretizations. SISSC, 6 (1988), 1073.
• [13] C.W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes. Journal of Computational Physics, 77 (1988), 439.
• [14] E. Abreu, D. Marchesin, F. Furtado and F. Pereira, Transitional Waves in Three-Phase Flows in Heterogeneous Formations. Computational Methods for Water Resources, I.5 Unsaturated and Multiphase Flow, Edited by C.T. Miller, M.W. Farthing, W.G. Gray and G.F. Pinder, Series: Developments in Water Science, ELSEVIER, I (2004), 609-620.
• [15] D.W. Peaceman, Fundamentals of Numerical Reservoir Simulation. Elsevier (1977).
• [16] E. Abreu and F. Pereira, Desenvolvimento de Novas estratégias Para a Recuperação Avançada de Hidrocarbonetos em Reservatórios de Petróleo. Boletim Técnico da Produção de Petróleo (PETROBRAS-CENPES), 2 (2007), 83-106.
• [17] E. Abreu, F. Furtado and F. Pereira, Three-phase immiscible displacement in heterogeneous petroleum reservoirs. Mathematics and computers in simulation, 73 (2006), 2-20.
• [18] E. Abreu, J. Douglas, F. Furtado and F. Pereira, Operator splitting based on physics for flow in porous media. International Journal of Computational Science, 2 (2008), 315-335.
• [19] E. Tadmor and G.-S. Jiang, Non-oscillatory central schemes for multidimensional hyperbolic conservation laws. SIAM Journal on Scientific Computing, 19 (1998), 1892-1917.
• [20] E. Tadmor and J. Balbas, A central differencing simulation of the Orszag-Tang vortex system. IEEE Transactions on Plasmas Science. The 4^th Triennal Issue on Images in Plasma Science, 33(2) (2005), 470-471.
• [21] F. Furtado and F. Pereira, Crossover from nonlinearity controlled to heterogeneity controlled mixing in two-phase porous media flows. Comput. Geosci., 7(2) (2003), 115-135.
• [22] F. Furtado, F. Pereira and J. Douglas Jr., On the numerical simulation of waterflooding of heterogeneous petroleum reservoirs. Comput. Geosci., 1(2) (1997), 155-190.
• [23] F. Pereira, J. Douglas Jr. and Li-Ming Yeh, A locally conservative Eulerian-Lagrangian method for flow in a porous medium of a mixture of two components having different densities. In: "Numerical treatment of multiphase flows in porous media (Beijing, 1999)", volume 552 of "Lecture Notes in Phys.", pages 138-155. Springer, Berlin (2000).
• [24] G. Chavent and J. Jaffré, Mathematical Methods and Finite Elements for Reservoir Simulation., volume 17. Studies in Mathematics and its Applications, North-Holland, Amsterdam, studies in mathematics and its applications edition (1986).
• [25] H.K. Dahle, M.S. Espedal, O.Sævareid and R.E. Ewing, Characteristic adaptive subdomain methods for reservoir flow problems. Numer. Methods Partial Differential Equations, 6(4) (1990), 279-309.
• [26] H.K. Dahle, R.E. Ewing and T.F. Russell, Eulerian-Lagrangian localized adjoint methods for a nonlinear advection-diffusion equation. Comput. Methods Appl. Mech. Engrg., 122(3-4) (1995), 223-250.
• [27] H. Nessyahu and E. Tadmor, Nonoscillatory central differencing for hyperbolic conservation laws. J. Comput. Phys., 87(2) (1990), 408-463.
• [28] H. Wang and R.E. Ewing, A summary of numerical methods for time-dependent advection-dominated partial differential equations. J. Comput. Appl. Math., 128(1-2) (2001), 423-445. Numerical analysis 2000, Vol. VII, Partial differential equations.
• [29] J. Douglas Jr. and Thomas F. Russell, Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures. SIAM J. Numer. Anal., 19(5) (1982), 871-885.
• [30] J.M. Thomas and P.A. Raviart, A mixed finite element method for 2^nd order elliptic problems. In: "Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975)", pages 292-315. Lecture Notes in Math., Vol. 606. Springer, Berlin (1977).
• [31] L.J. Durlofsky, A triangle based mixed finite element-finete volume technique for modeling two phase flow through porous media. Journal of Computational Physics, 105 (1993),252-266.
• [32] M.J. Blunt, M.R. Thiele and R.P. Batycky, A 3d field scale streamline-based reservoir simulator. SPE Reservoir Engineering, November (1997), 246-254.
• [33] P.D. Lax, Weak solutions of non-linear hyperbolic equations and their numerical computation. Comm. Pure Appl. Math., 7 (1954), 159-193.
• [34] R.P. Batycky, A Three-Dimensional Two-Phase Field Scale Streamline Simulator. PhD thesis, Stanford University, Dept. of Petroleum Engineering, Stanford, CA (1997).
• [35] S.K. Godunov, A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics. Mat. Sb. (N.S.), 47(89) (1959), 271-306.
• [36] V.V. Rusanov, The calculation of the interaction of non-stationary shock waves with barriers. Z. Vycisl. Mat. i Mat. Fiz., 1 (1961), 267-279.

Publication Dates

• Publication in this collection
  30 Mar 2009
• Date of issue
  2009

History

• Received
  10 June 2008
• Accepted
  24 Dec 2008