A numerical study of two-phase miscible flow through porous media with a Lagrangian model

Abstract

The multiphase flow mechanism in miscible displacement through porous media is an important topic in various applications, such as petroleum engineering, low Reynolds number suspension flows, dusty gas dynamics, and fluidized beds. To simulate such flows, volume averaging spatial operators are considered to incorporate pressure drag and skin friction experienced by a porous medium. In this work, a streamline-based Lagrangian methodology is extended for an efficient numerical approach to handle dispersion and diffusion of solvent saturation during a miscible flow. Overall pressure drag on the diffusion and dispersion of solvent saturation is investigated. Numerical results show excellent agreement with the results obtained from asymptotic analysis. The present numerical simulations indicate that the nonlinear effects due to skin friction and pressure drag cannot be accurately captured by Darcy’s method if the contribution of the skin friction dominates over that of the pressure drag. Moreover, mass conservation law is investigated, which is an important feature for enhanced oil recovery, and the results help to guide a good agreement with theory. This investigation examines how the flow regime may be optimized for enhanced oil recovery methods.

Introduction

The study of multiphase flow mechanism in porous media is an important topic in both academia and industry.^1–3 Examples include enhanced oil recovery techniques in petroleum engineering, which may be used at the tertiary stage of oil production to extract the remaining oil that is 50%–80% of the original oil in the reservoir. ⁴ In such techniques, if carbon dioxide is injected into deep oil reservoirs, it becomes supercritical and miscible with oil, leading to a reduction in the oil viscosity and in the tension force between the oil and the surrounding rock. Clearly, a numerical study of multiphase flow is useful to complement costly experimental approaches for enhanced oil recovery techniques.

In computational studies of miscible two-phase flow, e.g. Peaceman and Rachford, ⁵ the momentum equation may be expressed in the form of Darcy’s equation, ⁶ or in the form of a volume-averaged Navier–Stokes equation.^7,8 Derived empirically (see Whitaker ⁶ ), the former relates the flow rate linearly with the pressure gradient. The latter is robust and nonlinear; it incorporates Darcy’s results and provides a complete description of the nonlinear interaction between the fluid and the porous medium.^7–10 In this article, we report a numerical simulation of the miscible displacement in a confined reservoir using the volume-averaged Navier–Stokes equation. One objective of this work is to demonstrate a numerical model for dispersion of the displacing fluid (i.e., solvent), particularly when the inertial force strengthens because of miscibility. We also investigate how dispersion is influenced by the nonlinear drag (skin friction and pressure drag) of a porous medium, and by the viscous stress adjacent to impermeable regions. We present a methodology for modeling the nonlinear interaction between the porous medium, fluid motion, and reservoir boundary. Our numerical model introduces an algorithm for solving the governing equations, which is an extension of the classical IMPES method (IMplicit Pressure and Explicit Saturation). ¹¹ To improve numerical modeling of the nonlinear effects of dispersion, we have studied a Lagrangian method to compute the saturation of the displacing fluid efficiently. Based on our simulations, the Lagrangian method appears robust with respect to an equivalent Eulerian method.

Following is a brief review of related work that highlights the contribution of the current research and motivates the numerical modeling techniques presented in this article.

Motivation and related work

Peaceman and Rachford ⁵ used the theory of Fick’s law to model miscible two-phase flow of carbon dioxide and oil in a porous medium. In their model, ⁵ carbon dioxide and oil remain distinct and occupy their own volumes, where miscibility diminishes the tension force at the interface. The numerical implementation considers the two-phase flow in terms of the saturation of one of the phases, and the interface condition is derived from the assumption that the surface tension force vanishes. The comparison between the experimental results of Blackwell et al. ¹² and the corresponding numerical results of Peaceman and Rachford ⁵ supports the use of Fick’s law to model miscible two-phase flow accurately in porous media. A similar approach was also adopted by other researchers to study pore-scale viscous fingering.^13,14

In contrast, Koval ¹⁵ considered the theory of immiscible two-phase flow based on the Buckley–Leverett equations because for a miscible two-phase flow the time to sweep out oil by solvent is much smaller than the characteristic time scale of molecular diffusion.^16,17 An excellent review of the use of Fick’s law versus immiscible theory for simulating miscible flow in porous media is given by Udey and Spanos.¹⁷ In the present work, we have considered the model used by Peaceman and Rachford ⁵ ; however, we have also considered the volume-averaged Navier–Stokes equation because the Reynolds number of the flow may be increased during the enhanced oil recovery process. A detailed derivation of the volume-averaged Navier–Stokes equations in porous media is given by Breugem and Rees ⁹ (see the methodology section).

To solve the volume-averaged Navier–Stokes equation, there are two commonly used computational fluid dynamics approaches. One is the fully coupled solver, which resolves the dependence between the pressure and the saturation accurately at each iterative step by simultaneously discretizing all variables into one nonlinear system.^11,18–20 However, the computational cost and memory requirement may restrict the application of the fully coupled approach, particularly when the number of degrees of freedom is large. The coupled solver is often robust if an appropriate pre-conditioner is available. ²⁰ The other approach is the segregated solver, which treats the entire problem as a collection of sub-problems, each of which may be solved independently by individual schemes. The segregated solver is often flexible with respect to the fully coupled solver. In the petroleum industry, a very popular and powerful segregated algorithm for solving a two-phase flow is the IMPES method. ²¹ The IMPES algorithm is simpler and more efficient than some fully coupled solvers. ¹¹ However, for it to be stable, the classical IMPES requires extremely small time steps for the saturation, which is prohibitive, particularly for long time integration problems with small grid-blocks. The IMPES method can be improved by optimizing the pressure and saturation calculation because this is where most of the computational time is spent by the classical IMPES method. This article aims at advancing our knowledge in this direction.

Present work

The objectives of the present study are twofold: capturing the physics of the flow (i.e. flow in the reservoir scale), and optimizing the computational cost. We use the Darcy–Forchheimer equation to model the pressure drag and skin friction. We have combined the benefits of the multigrid method with that of the Lagrangian method to optimize the computational cost, which is in contrast with the model of Peaceman and Rachford. ⁵ We have demonstrated the performance of a Lagrangian method for transport in porous media. The present results indicate that the influence of skin friction on the overall transport mechanism is better resolved with the Darcy–Forchheimer equation when the flow rate is enhanced by miscibility. It is worth mentioning that the present idealized simulations help understand the development of optimal computational techniques for multiphase flows using the Lagrangian and the multigrid methods.

The article is organized as follows. The next section outlines the volume-averaged Navier–Stokes equation, and demonstrates the calculation of the drag force and its incorporation into the volume-averaged Navier–Stokes equation. After this, the segregated solver based on multigrid and Lagrangian methodologies is summarized. Then the results of numerical simulation are presented. Finally, the findings of the research, as well as further extensions of this development, are discussed.

Methodology

The use of multiphase flow in enhanced oil recovery is illustrated schematically in Figure 1(a). A region of homogeneous isotropic porous medium and a representative elementary volume of that medium are illustrated in Figure 1(b) and 1(c), respectively. A similar idealized model was also considered by Bottero et al. ²² (see Figure 2 therein). Consider now a porous matrix with pores completely filled with oil that is to be displaced by a solvent. Throughout the displacement process, the interface between the solvent and the oil is represented by a continuous phase indicator function c( x , t). However, the pores are not resolved by the mesh, where the overall pressure drag and skin friction exerted by the porous matrix is calculated through the volume-averaged method.

Figure 1.

(a) Schematic description of the oil recovery method by injecting a second fluid such as carbon dioxide or water. The flow direction in a porous medium is shown by red arrows. (The drawing is obtained from Google image repository; http://www.kgs.ku.edu/CO2/evaluation/slide03.html.) (b) The porous medium is represented by a collection of uniformly distributed solid bodies. (c) Representative elementary volume that contains both the solid and the fluid.

Figure 2.

A validation of the procedure, where the numerical solution u(x, y, 1), v(x, y, 1) and P(x, y, 1) (left column) is compared with the exact solution (right column). Here, (a) and (b) are the u velocity component, (c) and (d) are the v velocity component, and (e) and (f) are the pressure.

The volume averaging technique

Defining the porosity as $\phi =\Delta V\text{f}/\Delta V$ , where ΔV is the volume of the representative elementary volume and $\Delta V\text{f}$ is the volume of fluid in the representative elementary volume, the flux of fluid per unit area (the “Darcy velocity”) is given by the volume average

$$u_i^{\text{D}}=\frac{1}{\Delta V}\int \Delta V\tilde{u}i(xi,t)\text{d}V$$

which is related to the intrinsic average through the pore spaces

$$ui\equiv 〈ui〉=\frac{1}{\Delta V\text{f}}\int \Delta V\text{f}\tilde{u}i(xi,t)\text{d}V$$

such that $u_i^{\text{D}}=\phi 〈ui〉$ .^7,8,10,23 We may drop the notation 〈·〉 for simplicity. It should be noted that ΔV may be taken as the same as a computational cell

Let us assume the following. ⁹

All pores in the medium are connected and the porosity is constant in space and time.

The volume-averaged quantities are well-behaved; i.e. if $\tilde{u}i=〈ui〉+u\text{'}i$ , then $〈u\text{'}i〉=0$ .

The system is in local thermal equilibrium and the geothermal effects are compensated by the geo-pressure gradient.

The production and dissipation of $〈(1/2)u\text{'}iu\text{'}i〉$ are in local balance.

The dispersion and transport of each of the fluid phases are balanced locally by the rate of change of mass within a given volume.

Calculation of the drag force

Based on assumptions (1) to (4), an average drag per unit volume is^9,10

$$fi=\stackrel{\text{pressuredrag}}{\overbrace{-\frac{1}{\Delta V}\int \int \partial S[\tilde{P}ni]\text{d}S}}+\underset{\text{skinfriction}}{\underbrace{\frac{\mu }{\Delta V}\int \int \partial S[\frac{\partial \tilde{u}i}{\partial n}]\text{d}S}}$$

(1)

where the momentum equations are averaged over a representative elementary volume. In the literature, the drag in equation (1) is approximated by an empirical linear formulation ⁹

$$fi=-\frac{\mu \phi }{K}〈ui〉$$

which compensates the overall drag force derived from Darcy’s experiment; i.e.

$$\frac{\mu \phi }{K}〈ui〉=-\nabla P+\rho g$$

Using a nonlinear component, known as Forchheimer drag, ¹⁰ the total drag may be expressed by

$$fi=-\frac{\mu }{K}{u}_i^{\text{D}}+\frac{c\phi \rho \left|u_i^{\text{D}}\right|u_i^{\text{D}}}{\sqrt{K}}$$

where $c\phi$ is referred to in the literature as the nonlinear Forchheimer coefficient.

In a confined reservoir, the Forchheimer drag term may be modified to account for the rate of shear and dispersion in a flow. For sediment transport applications, the shear velocity (u_τ) is typically evaluated at the lower boundary by $u\tau =\sqrt{\tau 0/\rho }$ . For a laminar channel flow, the shear stress at wall (τ₀) vanishes linearly as the center of the channel is approached. In a confined reservoir, the shear stress would act against the drag of the porous medium in addition to that of the impermeable boundary, which in terms of the shear velocity is given by equation (1) of Nield et al. ²⁴

$$\mu \nabla 2u\tau -\frac{\mu }{K}u\tau =0$$

(2)

The fluid above an impermeable surface tends to push the fluid underneath forward, whereas the surface tends to drag the upper fluid backward. From equation (2), it can be shown that shear stress is confined to a boundary layer of thickness $\mathcal{O}(\mathcal{D}a-1/2)$ as $\mathcal{D}a\to 0$ . ²⁴ As a result, the overall effect of the pressure drag and skin friction can be evaluated by

$$fi=-\frac{\mu }{K}{u}_i^{\text{D}}+\frac{c\phi \rho \left|u\tau \right|u_i^{\text{D}}}{\sqrt{K}}$$

(3)

In equation (3), the shear velocity, $u\tau$ , computed from equation (2), is used in the Forchheimer drag term.

Computational modeling of miscible displacement

As derived by Peaceman and Rachford, ⁵ the mass balance equation for a two-phase flow in porous media is

$$\frac{\partial (\phi \rho kSk)}{\partial t}+\frac{\partial }{\partial xj}[\rho kSk{u}_j^{\text{D}}-D\frac{\partial (\rho kSk)}{\partial xj}]=0\sum _{k=0}^{1}Sk=1$$

(4)

where subscripts k = 0 and k = 1 denote oil and solvent, respectively, and D is the dispersion coefficient.^5,11,17 As can be seen from equation (4), the individual saturation S_k is not conserved. However, the total mass with respect to a given volume containing oil, solvent, and rock is conserved (e.g. ρ = ρ₀ + ρ₁)

$$\frac{\partial (\phi \rho )}{\partial t}+\frac{\partial \rho u_j^{\text{D}}}{\partial xj}=0$$

(5)

It is useful to define the saturation $c(x,y,t)=S1=1-S0$ as a phase indicator, which leads to the following transport equation^5,10,11

$$\frac{\partial (\phi \rho c)}{\partial t}+\frac{\partial }{\partial xj}[c\rho u_j^{\text{D}}-\rho D\frac{\partial c}{\partial xj}]=0$$

(6)

Combining equation (5) with equation (6) and using the volume-averaged Navier–Stokes equation, one obtains a set of governing equations for miscible displacement in porous media.

The Darcy–Forchheimer Navier–Stokes model of miscible flow

Following the work of Peaceman and Rachford, ⁵ we have combined equations (5) and (6) into a single equation that models the dynamics of solvent and oil. We have used the convective terms in the rotational form in the volume-averaged Navier–Stokes equation to calculate the velocity pressure coupling effectively, ²⁵ where the drag force is obtained according to equations (2) and (3). With respect to a length scale H and a velocity scale U, it is convenient to express all variables in the dimensionless form. Since φ is assumed constant in this work, the equations governing a miscible two-phase flow in porous media can be simplified in terms of the intrinsic velocity (using u_i for 〈u_i〉) and the saturation

$$\frac{\partial ui}{\partial xi}=0$$

(7)

$$\frac{\partial ui}{\partial t}+uj(\frac{\partial ui}{\partial xj}-\frac{\partial uj}{\partial xi})=-\frac{\partial }{\partial xi}(P+\frac{1}{2}uiui)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+\frac{1}{\Re e}\frac{\partial }{\partial xj}\frac{\partial ui}{\partial xj}+fi$$

(8)

and

$$\frac{\partial c}{\partial t}+uj\frac{\partial c}{\partial xi}=\frac{1}{\Re e\mathcal{S}c}\frac{\partial }{\partial xj}\frac{\partial c}{\partial xj}$$

(9)

Equations (7) to (9) contain three dimensionless parameters: the Reynolds number $\Re e=\mathrm{UH}/\nu$ , the Darcy number $\mathcal{D}a=K/H2$ , and the Schmidt number $\mathcal{S}c=\nu /D$ . The dimensionless form of the drag term f_i in equation (8) is given by

$$fi=\underset{\text{linear}}{\underbrace{-\frac{\phi ui}{\Re e\mathcal{D}a}}}-\stackrel{\text{quadratic}}{\overbrace{\frac{C\varphi \phi \left|u\tau \right|ui}{\sqrt{\mathcal{D}a}}}}$$

(10)

In a study of an oil reservoir using the volume-averaged Navier–Stokes equation, the computational challenge of solving the coupled system of equations (7) to (9) requires fast numerical algorithms.

The segregated algorithm

Calculation of pressure and velocity

The velocity pressure coupling can be described by a pressure Poisson equation, which is derived by taking the divergence of equation (8) and requiring that equation (7) is satisfied at each time, step such that

$$\frac{\partial 2}{\partial xj\partial xj}(\underset{pn+1}{\underbrace{P+\frac{uiui}{2}}})=\frac{\partial }{\partial xi}(fi-u_j^n(\frac{\partial u_i^n}{\partial u_j^n}-\frac{\partial u_j^n}{\partial u_i^n}))$$

(11)

For a two-phase flow in porous media, pressure changes less rapidly in time than saturation. In equation (11), it is thus appropriate to take a larger time step for the pressure than for the saturation. Solving equation (11) with a multigrid method is asymptotically optimal, which means that the CPU time does not increases drastically when the number of grid-blocks increases (see Alam and Bowman ²⁵ ). When the Reynolds number is not too large, explicit treatment of advection and drag force, and implicit treatment of pressure and viscous force, uses a near-optimal algorithm for solving equation (8), where the elliptic system

$$\frac{2u_i^{n+1}}{\Delta t}-\frac{1}{\Re e}\frac{\partial 2u_i^{n+1}}{\partial xj\partial xj}=\frac{2u_i^n}{\Delta t}+\frac{1}{\Re e}\frac{\partial 2u_i^n}{\partial xj\partial xj}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}-2[\frac{\partial }{\partial xi}(\underset{pn+1}{\underbrace{P+\frac{u_i^n{u}_i^n}{2}}})+u_j^n(\frac{\partial u_i^n}{\partial xj}-\frac{\partial u_j^n}{\partial xi})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}-\frac{\phi u_i^n}{\Re e\mathcal{D}a}+\frac{C\varphi \phi \left|u\tau \right|u_i^n}{\sqrt{\mathcal{D}}a}]$$

(12)

is solved by a multigrid method.

Calculation of the saturation

Let us now outline the streamline-based Lagrangian model. Like the treatment of the momentum equation (12), the linear dispersion term in equation (9) is treated implicitly (convenient for large values of $\mathcal{S}c$ ), and the nonlinear advection term is treated explicitly, where a fractional time stepping method is considered (Chapter 7 in Chen et al. ¹¹ ).

Given the saturation cⁿ and the velocity $u_i^n$ from a recent time step, equation (9) can be solved as a pure advection problem without the dispersion term; let us denote this solution c^*. It is convenient to leave the presentation of the advection problem for c^* to the next section. Now, the effect of dispersion can be incorporated into c^* by solving

$$\frac{2cn+1}{\Delta t}-\frac{1}{\Re e\mathcal{S}c}\frac{\partial }{\partial xj}(\frac{\partial cn+1}{\partial xj})=\frac{2c*}{\Delta t}+\frac{1}{\Re e\mathcal{S}c}\frac{\partial }{\partial xj}(\frac{\partial c*}{\partial xj})$$

(13)

using the multigrid method with the same time step Δt considered for the pressure equation.

In such a splitting method, as discussed in detail by Chen et al., ¹¹ the saturation c^* could be calculated using a time step smaller than Δt in order to satisfy the stability condition. The approach is still good because the costly pressure equation is not solved for every smaller time step. However, the convective derivatives of saturation must be discretized with an upwind biased method. An artifact of such discretization is the damping of high-frequency modes of the numerical solution. In contrast, we have calculated c^* with a Lagrangian method that overcomes the artificial damping of high-frequency modes by modeling the convective transport using the streamlines. As a result, the only restriction on the time step comes from the consideration of accuracy in contrast with stability.

The Lagrangian method for saturation

The streamlines associated with a given velocity field are curves at which the velocity field is tangential. Such streamlines can be parametrized, and traced by solving the differential equation

$$\frac{\text{d}si}{\text{d}\xi }=ui$$

(14)

For an incompressible flow, streamlines do not intersect, begin, or end inside the fluid. If k  = (k₁, k₂) is the index of a grid point, only one streamline passes through this point. Knowing a reference $\mathbf{s}k(0)$ of the streamline that passes through the k th grid point, the entire streamline can be traced by solving equation (14) using an appropriate ordinary differential equation solver routine. Since the velocity field is tangential to the streamlines, the velocity vector $v=[u1,u2]\text{T}$ is given by the directional derivative of a streamline, s , along the direction of s . Thus, from equation (14), the directional derivative along a streamline can also be parametrized, such that

$$\frac{\text{d}}{\text{d}\xi }=uj\frac{\partial }{\partial xj}$$

and hence, the directional derivative of saturation can be written as

$$\frac{\partial c}{\partial \xi }=uj\frac{\partial c}{\partial \xi }$$

A technical detail of the streamline-based Lagrangian method is also given by Alam and Penney. ²⁶

The directional derivative is beneficial because the advective evolution of saturation in equation (9) takes the form

$$\frac{\partial c}{\partial t}+\frac{\partial c}{\partial \xi }=0$$

(15)

In equation (14), ξ is a time-like parameter and, for each value of ξ. we have a corresponding value for the streamline s (ξ). However, in equation (15), ξ is a space-like parameter describing the advection of c along a streamline. In the ξ−t plane, the solution of equation (15) is unique because the characteristic curves of equation (15) never intersect, and the solution of equation (15) is obtained analytically such that $c*(\mathbf{s}(\xi ),\Delta t)=cn(\mathbf{s}(\xi -\Delta t)).$

It is clear from this discussion that the costly computational task for the multidimensional advection of saturation is replaced with the analytical solution of a simple one-dimensional problem. However, we only need to trace a portion of the streamline passing through each grid point at each time step. Starting a streamline from a grid point, its terminal position may be traced within a grid-block that is not a grid point. As a result, the evaluation of the right hand side of equation (13) requires interpolation of c^* between grid points and streamlines. Alam and Penney ²⁶ developed an interpolation technique that conserves mass. Since only the saturation is interpolated, satisfying conservation of mass by the interpolation method is important to simulate a two-phase flow problem.

Implementation of the methodology

The proposed segregated algorithm has been implemented with the object-oriented C++ programming paradigm. We have extended an existing C++ code that implements the multigrid solver presented by Alam and Bowman ²⁵ to solve equations (11) to (13). The implementation of the proposed Lagrangian solver (e.g. equations (14) and (15)) takes advantage of the $\mathcal{O}(\log \mathcal{N})$ complexity of some data structures (e.g. std::map) from the C++ standard template library. The entire solver used for this research has been implemented as a module of the computational fluid dynamics code, AWCM++, which solves the Navier–Stokes equations and was fully verified for several other applications.^19,20 The following algorithm outlines the implementation of the methodology proposed in this article:

For each time step,

‐ solve equation (11) for P;

‐ solve equation (12) for u ;

‐ solve equation (13) for C;

‐ solve equation (14) for the streamlines;

‐ solve equation (15) for advection;

• goto next time step.

Numerical simulation, results, and discussion

In this section, we analyze the current methodology for simulating miscible displacement in a fluid saturated porous channel.^5,22 The first set of simulations aims to understand the computation of forces between the fluid and the porous medium. The second set of simulations aims to understand how accurately the proposed Lagrangian method captures an interface between two nearly immiscible fluids; how this interface moves and deforms in porous media. More specifically, conserving the mass of an individual phase of miscible fluids is investigated by the Lagrangian method and the results are compared with numerical experiments for the fluid flows in porous media. Finally, we discuss a potential application of the methodology in reservoir engineering, where flow paths in heterogeneous reservoirs are analyzed using methods called pressure tests and tracer tests.

Parameters and numerical verification

The accuracy of our method for modeling the velocity pressure coupling is validated using a common test case: an incompressible flow in a cavity for which the exact solutions of equations (7) and (8) are given by

$$u(x,y,t)=-\cos (2\pi x)\sin (2\pi y)\text{e}-8\pi 2t/\Re ev(x,y,t)=\sin (2\pi x)\cos (2\pi y)\text{e}-8\pi 2t/\Re eP(x,y,t)=-\frac{1}{4}(\cos (4\pi x)+\cos (4\pi y))\text{e}-16\pi 2t/\Re e$$

In the limit $\mathcal{D}a\to \infty$ , the numerical solutions of equations (7) to (8) have been compared with the exact solution in Figure 2. One notes the excellent agreement between the exact solution and the numerical solution. For a more detailed understanding, the exact solution of the horizontal velocity along the line y = 0.5 is compared with the numerical solution u(x, 0.5, 1). Similarly, the vertical velocity is considered along the line x = 0.5. Both u(x, 0.5, 1) and v(0.5, y, 1) are presented in Figure 3. This quick verification demonstrates the accuracy of the methodology in simulating a single phase flow in a clear region without a porous matrix. Let us now discuss the parameters for simulating the interaction between a two-phase flow and a porous medium.

Figure 3.

(a) Numerical solution u(0.5, y, 1), compared with exact solution. (b) Numerical solution u(x, 0.5, 1), compared with exact solution.

As seen from Table 1, the physical properties of fluid flows in petroleum reservoirs vary widely. We have derived a set of characteristic scales (see Table 2) from the data listed in Table 1. For clarity, these parameters are listed in Tables 3 and 4. Note that the results presented in the following sections are discussed using the dimensionless numbers listed in Table 2. A porous channel of depth H is decomposed into a mesh of 513 × 129 nodes with Δx = 2.9 m and Δy =3.9 m, unless stated otherwise. On this mesh, a V-cycle multigrid iteration with 9 × 3 grid points on the coarsest level and 513 × 129 grid points on the finest level is used. On each level, a red–black ordered Gauss–Seidel method is employed to damp high-frequency modes of the error. The divergence-free velocity is obtained at each time step through the staggered arrangement of velocity and pressure (see Harlow and Welch) ²⁷ . The rate of convergence is tested by changing the pre- and post-smoothing sweeps, where the number of maximum V-cycle iterations is set to 10. We have found that about 3–5 pre- and post-smoothing sweeps are necessary to reduce the residual for both the pressure and the velocity below 10⁻⁶. In the presentation of the results, yellow is chosen to represent oil, and red is used for solvent.

Table 1.

Characteristic parameters collected from relevant references.

Reference	Velocity, U, m/s	Reservoir thickness, H, ft	Oil viscosity, μ, mPa s	Density, ρ, g/cm3	Permeability, κ, mD	Porosity, φ, %
28	1.52 × 10−6 − 1.35 × 10−5	–	4.3–53.6	0.831–0.895	14, 71, 136	20.8–21.9
29	–	60.96	–	–	0.0001	6
30	–	30.48–48.77	70–1500	–	500–7500	20–36
31	3.53 × 10−6	15.24–51.82	–	–	20–150	–
32	–	152.4–182.88	400–1500	0.95– 0.98	–	29–35
(Lloydminster oil reservoir, CANADA)

Table 2.

Parameters and dimensionless numbers.

Parameter	Symbol
Velocity scale	U (ms−1)
Length scale	H (m)
Time scale	T (s)
Permeability	K (m2)
Kinematic viscosity	ν (m2s−1)

Dimensionless parameters	Symbol
Reynolds number	$\Re e=\frac{\mathrm{UH}}{\nu }$
Darcy number	$\mathcal{D}a=\frac{K}{H2}$
Schmidt number	$\mathcal{S}=\frac{\nu }{D}$

Table 3.

Parameters (available in the literature) used in the simulations for an idealized reservoir model.

Parameter	Value used in the simulations	Reference
U	3.5 × 10−6m/s	31
H	500 m	32
ν	1.75 × 10−3 m2/s	32
φ	18%	33

Table 4.

Parameters for simulations with Reynolds number, $\Re e\ge 1$ and Darcy number $10-4<\mathcal{D}a<106$ .

Parameter	Value
L1 × L2	3 × 1
n1 × n2	512 × 128
Δt	10−2
∇ P 0	2
α	1
$\Re \mathrm{eSc}$	20, 000
φ	18%

Effect of shear stress and permeability on miscible displacement

Peaceman and Rachford ³⁴ carried out experiments to investigate the overall flow rate in a porous channel. Accordingly, in a confined reservoir or aquifer, the flow rate is dependent on the shear stress and permeability of the medium, which may also be characterized by the Darcy friction factor, $f\text{D}$ , a dimensionless number that measures the wall shear stress in a porous channel. $f\text{D}$ is inversely proportional to the Reynolds number, and is also related to the pressure drop in the channel, i.e.

$$\frac{\Delta P}{L}=\frac{f\text{D}\rho U2}{2H}$$

As observed in experiments by Beavers and Joseph, ³⁴ the shear stress is confined within a boundary layer. Zhang and Prosperetti ³⁵ have further discussed the experimental findings of Beavers and Joseph. ³⁴ For large $\mathcal{D}a$ and small $\Re e$ , the mass flow rate in such a porous channel of width H is given by^34,35

$$M=-\frac{\rho H3}{12\mu }\frac{\text{d}P}{\text{d}x}$$

which implies that the dimensionless mean streamwise velocity in the channel is about $\Re e/12$ . ³⁵

To demonstrate these phenomena, two sets of simulations have been considered here: one is for $\Re e=1$ and the other is for $\Re e=102$ (see Table 4). Figure 4 presents the result of eight numerical experiments, where the dimensionless value of the pressure dropwas kept fixed, the Darcy number was varied, $10-4<\mathcal{D}a<106$ , and the Forchheimer effect was turned off, i.e. C_ϕ = 0. The streamwise and spanwise velocities are computed at a distance x = 1.5H from the inlet when the solvent phase has traveled a distance of approximately 2H along the channel. These data demonstrate how the flow rate is affected by the Darcy number ( $\mathcal{D}a$ ) and the Darcy friction factor ( $f\text{D}$ ). In Figure 4, the simulated velocities are compared with respect to the variation of the Darcy number and the Reynolds number. Clearly, as the Darcy friction factor weakens at higher $\Re e$ , the drag of porous media weakens with respect to a given $\mathcal{D}a$ , and vice versa. At a small Darcy number, e.g. $\mathcal{D}a\le \phi \times 10-2$ , the pore-scale flow and transport are dominated by the capillary pressure drag and skin friction. We see that a reduction of $f\text{D}$ by a factor of 10² has increased the flow rate by a factor 10² (e.g. Figure 4(b) and (c)). For $\mathcal{D}a>\phi \times 10-2$ , however, the influence of the pressure drag and skin friction on flow rate has been diminished. Clearly, at $\mathcal{D}a\ge \phi \times 106$ , the flow rate is primarily dependent on the hydrodynamic pressure drop and viscous stress. ³⁵ In Figure 4(a), we see that the dimensionless mean streamwise velocity approaches $\Re e/12$ as expected. One also notices from Figure 4 that the parabolic nature of the horizontal velocity profile turns into a relatively flat shape when the Darcy friction factor is reduced. The vertical velocity profiles indicate clearly that the hydrodynamic transverse dispersion is influenced by the dominant inertial effect when the Darcy friction factor is reduced.

Figure 4.

Velocity profiles u₁(1.5, y, T) and u₂(1.5, y) for $\Re e=1$ and $\Re e=102$ and for $10-2\le \mathcal{D}a/\phi \le 106$ . Plots with increased $\Re e$ correspond to those with decreased $f\text{D}$ . Here, the parameter values are: (a) u₁(1.5, y) at $\Re e=1$ , (b) u₁(1.5, y) at $\Re e=102$ , (c) u₂(1.5, y) at $\Re e=1$ , and (d) u₂(1.5, x₂) at $\Re e=102$ .

Parabolic shape of the dispersed solvent phase

Corresponding to the velocity field presented in Figure 4, let us now present the simulated moving front of the solvent at $\Re e\mathcal{S}c=2\times 104$ . The initial solvent phase is displayed in Figure 5(a). The time evolution of this solvent front is presented in Figure 5(b) to (e). To investigate the effects of the pressure drag and skin friction on the transport of the solvent phase, we have simulated the flow with various values of $\Re e$ and $\mathcal{D}a$ . Note that the dimensionless porous term in the equation without the Forchheimer drag is $\phi ui/\Re e\mathcal{D}a$ . We have compared the evolution of the solvent phase for $\phi /\Re e\mathcal{D}a=1$ (Figure 5(b) and (c)) and $\phi /\Re e\mathcal{D}a=10-6$ (Figure 5(d) and (e)) when both parameters ( $\Re e$ and $\mathcal{D}a$ ) are varied. The values of the Darcy number for the plots in Figure 5 are: (b) $\mathcal{D}a=\phi \times 100$ ; (c) $\mathcal{D}a=\phi \times 10-2$ ; (d) $\mathcal{D}a=\phi \times 106$ ; and (e) $\mathcal{D}a=\phi \times 104$ . Between Figure 5(b) and 5(c), the impact of the reduced value of the Darcy friction factor is clear. However, between the second and third rows in Figure 5, the change in the Darcy number by a factor of 10⁶ does not exhibit a notable change in the flow pattern.

Figure 5.

Miscible flow of carbon dioxide and oil; red and yellow represent nonzero and zero concentrations of carbon dioxide, respectively. (a) Initial concentration, where carbon dioxide has been emplaced in the region of area A. (b) Displacement of oil by carbon dioxide at $\Re e=1$ , $\mathcal{D}a=\phi$ , and t = 24.4H/U, where the rectangle represents the region of area A, which is expected to be filled with carbon dioxide. (c) $\Re e=102$ , $\mathcal{D}a=\phi$ , and t = 2.3H/U. (d) $\Re e=1$ , $\mathcal{D}a=\phi \times 106$ , and t = 22.5H/U. (e) $\Re e=102$ , $\mathcal{D}a=\phi \times 106$ , and t = 2.28H/U.

Note that the time for Figure 5(b) and (d) is t = 22.5H/U, while that for Figure 5(c) and (e) is t = 2.3H/U. It is the individual time when the solvent phase has traveled a distance of approximately 2H. At $\Re e=1$ , the Darcy friction factor is strong and the front moves approximately 10 times slower than the case at $\Re e=102$ where the Darcy friction factor is weak. One observes from these data that the distortion of the solvent phase toward the plug flow regime (i.e. a flat shape of the front) from the pipe flow (i.e. parabolic) is strongly influenced by the Darcy friction factor, in contrast to the Darcy number. For a qualitative measure of the distortion of the initial shape of the sample, a rectangle has been drawn in each of the plots in Figure 5. We see that the proportion of solvent within the marked region has increased at $\Re e=102$ , compared with $\Re e=1$ , indicating that field-scale sweep efficiency would improve if the Darcy friction factor weakens. It is worth mentioning that these results need to be interpreted carefully because experimental measurements are not available at this stage of the study.

To provide further insight into the effect of the Darcy friction factor, two sets of results are displayed in Figure 6 for three values of the Darcy number, $\mathcal{D}a=\phi \times 10-2$ , $\mathcal{D}a=\phi \times 10-1$ , and $\mathcal{D}a=\phi \times 100$ , and two values of the Reynolds number, $\Re e=1$ and $\Re e-102$ . The plots in Figure 6 indicate that the solvent moves faster if the permeability increases, as expected. However, this speed-up is also influenced strongly by the Darcy friction factor, and the distortion of the solvent sample into a parabolic shape is not primarily dependent on the permeability, but on the wall shear stress.

Figure 6.

Temporal evolution of flow with respect to variations of $\Re e$ and $\mathcal{D}a$ with $\Re e\mathcal{S}c=20000$ and φ = 18%. Left column (a, b, c) $\Re e=1$ , t = 30H/U. Right column (d, e, f): $\Re e=102$ , t = 2H/U. Variation of Darcy number: (a, d) $\mathcal{D}a=\phi \times 10-2$ ; (b, e) $\mathcal{D}a=\phi \times 10-1$ ; (c, f) $\mathcal{D}a=\phi$ .

Dispersion and diffusion during miscible displacement

As a miscible flow progresses, we would like to investigate the nature of the dispersion experienced by the flow. We compare this dispersion and diffusion with the theoretical prediction of convection–diffusion theory. ¹⁷ A key feature of miscible flow is that the dispersion can be described in analogy with Fick’s law of diffusion. ⁵ Such an analogy between theoretical prediction and experimental data is given in Figure 1 of Udey and Spanos. ¹⁷ We have solved equations (7) to (9) to predict the dispersion of saturation for several values of $\Re e$ and $\mathcal{S}c$ . We have compared the numerical solution with the exact solution based on the convection–diffusion theory. ¹⁷ We have also compared the current Lagrangian solution with an equivalent Eulerian solution. A technical detail of the Eulerian method is given by Ahammad. ³⁶

In Figure 7(a), saturation c(x, 0.5, 24.5) is predicted by the Eulerian method and the Lagrangian method, and compared with the exact solution for $\Re e=1$ . Similar data c(x, 0.5, 2.55) for $\Re e=102$ are presented in Figure 7(b). The Lagrangian prediction has excellent agreement with the convection–diffusion theory. The convection–diffusion theory for dispersion in porous media was extensively analyzed and compared with experimental data by Udey and Spanos. ¹⁷ Although we do not have experimental data corresponding to our present simulation, it is clear from Figure 7 that the Lagrangian simulation is extremely accurate with respect to the theory. Figure 8 shows domain-averaged maximum saturation at t = 24.5U/H for $\Re e=1$ and at t = 2.55U/H for $\Re e=102$ , and the results have been compared for various values of $\Re \mathrm{eSc}$ . The Schmidt number, $\mathcal{S}c$ , represents the ratio of the viscous diffusion rate to the molecular diffusion rate. As $\mathcal{S}c$ increases, the rate of molecular diffusion decreases. The effect of molecular diffusion is further illustrated in Figure 9, in which Lagrangian simulations are compared with Eulerian simulations. Figures 7 to 9 indicate that the Lagrangian solution agrees closely with the theoretical solution. As $\mathcal{S}c$ increases, the impact of artificial numerical diffusion for the Eulerian solution is also clear. For technical details on the impact of artificial numerical diffusion of Eulerian methods, readers are referred to Pletcher. ³⁷

Figure 7.

Comparison of the prediction of saturation between the Lagrangian method, the Eulerian method, and the theory for $\Re e\mathcal{S}c=2\times 104$ : (a) $\Re e=1$ ; (b) $\Re e=102$ .

Figure 8.

Domain-averaged maximum saturation as a function of $\Re \mathrm{eSc}$ : (a) $\Re e=1$ ; (b) $\Re e=102$ .

Figure 9.

Comparison of saturation between Eulerian and the Lagrangian calculations: (a) $\Re e=1$ , Eulerian; (b) $\Re e=1$ , Lagrangian; (c) $\Re e=102$ , Eulerian; (d) $\Re e=102$ , Lagrangian.

Forchheimer effect on dispersion and diffusion

Let us now demonstrate the Forchheimer effect on the dispersion of saturation. The simulation is similar except that we now switch on the Forchheimer term by setting C_ϕ = 1. For the plots in Figure 10, a small Darcy number is used ( $\mathcal{D}a=\phi \times 10-8$ ). To compare the impact of the Darcy friction factor, we have chosen $\Re e=1$ and $\Re e=102$ . Note that the effect of the linear component of the drag is neglected for this test. As a result, we can see a competition between the Darcy friction factor and the Forchheimer drag. As can be seen from Figure 10(a) and (b), the flow rate is faster at $\Re e=1$ with respect to $\Re e=102$ . Note that we have used ΔP = 0.5 for the dimensionless pressure drop in the case of $\Re e=1$ and ΔP = 0.15 in the case of $\Re e=102$ . This is the consequence of an important feature of reservoir performance: the stronger the viscous effect, the higher the pressure gradient. ³⁸

Figure 10.

Comparison between the Darcy friction factor and nonlinear Forchheimer drag. (a) $\Re e=1$ , t = 16.5U/H; (b) $\Re e=102$ , t = 50U/H; (c) center line of the concentration field, as (a), for different values of $\Re \mathrm{eSc}$ ; (d) center line of the concentration field, as (b), for different values of $\Re \mathrm{eSc}$ .

The simulations have been repeated for various values of $\Re e$ and $\mathcal{S}c$ , and representative plots are demonstrated in Figure 10(c) and (d). To illustrate further, we have computed c(x, 0.5, t), where t is the time when the saturation has traveled an approximate distance 1.5H. We have found that t = 24.5U/H for $\Re e=1$ and t = 2.55U/H for $\Re e=102$ . The result is compared between the Lagrangian and the Eulerian calculations in Figure 9. Clearly, the Eulerian calculation exhibits a large amount of diffusion as the saturation front is dispersed in the porous channel.

In Figure 11(a) to (c), the saturation is compared for a few values of $\Re \mathrm{eSc}$ , where the Forchheimer drag is used. We observe that a flat profile of saturation exists except for a thin boundary layer. This plug flow regime results from the reduction of shear stress by the nonlinear drag term. In addition, the maximum value of saturation in Figure 11(a) to (c) has been compared with the asymptotic convection–diffusion theory and the results are exhibited in Figure 11(d). An excellent agreement between the Lagrangian and theoretical prediction is observed. This test also ensures the conservation of mass fraction for miscible displacement in a multiphase flow system.

Figure 11.

(a) to (c) Migration of solvent as a flat band with reducing Darcy friction factor near the reservoir boundaries (Forchheimer drag is applied). (d) Maximum values of the saturation show agreement with asymptotic convection–diffusion theory. For clarity, the maximum value of the saturation predicted by the Eulerian method is also compared. The Eulerian simulation does not converge to the present Lagrangian simulation. (a) $\Re \mathrm{eSc}=104$ ; (b) $\Re \mathrm{eSc}=2\times 104$ ; (c) $\Re \mathrm{eSc}=105$ . (d) Mass conservation analysis for the different values of $\mathcal{S}c$ .

Conclusions and future direction

There has been great interest in miscible displacement of oil by a solvent in the field of petroleum engineering, and the topic of carbon dioxide injection for the enhanced oil recovery process has received increased attention. Recently, there is a growing propensity for employing high-performance computational fluid dynamics reservoir simulation techniques with the hope that this may lead to the discovery of substantial currently unrecognized opportunities for increasing the economic recovery of hydrocarbons. ³⁹ Intrinsic volume-averaged techniques are used for the development of a mathematical model to capture the flow at the field scale. A detailed computational fluid dynamics investigation helps to understand how injected solvent displaces oil during a multiphase flow system through porous materials in an oil reservoir, and provides further insight when field results fall short of laboratory performance. Diffusion and dispersion of the interface during miscible displacement have been studied extensively with a streamline-based Lagrangian technique. The numerical model was examined by two test cases to investigate its efficiency, and numerical predictions are in good agreement with theoretical data. ³⁶ Many researchers have investigated non-Darcy flow in oil reservoir.^40,41 We extensively discuss the Darcy effect and the Forchheimer effect (non-Darcy) on the flow, i.e. non-Darcy drag force on the flow. During the miscible displacement, it is found that flow rate is dependent on shear stress along the boundary layer and the permeability of the medium. Phase behavior phenomena of the solvent dissolution during miscible displacement are identified as well. Moreover, the mass conservative law is studied and explained, and a good agreement is seen with the analytical solution. The field-scale sweep efficiency is one of the most important factors affecting the economic recovery of hydrocarbons, for which an improved upscaling model may be beneficial to industries instead of studying the details of pore-scale flow. This article reiterates such a growing trend, and outlines some significant developments in this direction.

To analyze the sweep efficiency for a miscible displacement, overall flow characteristics can be studied through an upscaling model. In particular, the idealized influence of an isolated impermeable region partially embedded in a porous matrix has been studied. Note that the averaging process employed in this article is commonly used to model turbulent flow in porous media. We have outlined how to extend and redesign an intrinsic space–time average operator to model small-scale transient variability in a laminar miscible flow, which is a distinct fundamental development of this work. This research has clearly explained the regime of the viscous shearing stress, which plays a negligible role on large-scale reservoir flow but is important for accurate upscaling of small-scale physics.

Some important future developments include the following. An extension to a three-dimensional multilevel mesh employing a massively parallel algorithm would add potential benefits toward the needs of petroleum industries. We are currently developing data structures using the object-oriented C++ programming paradigm. The current results clearly provide hintsfor using adaptive mesh approaches.^42,43 However, a full understanding of the efficiency of multilevel and Lagrangian approaches is essential before a study of the efficiency of parallel adaptive mesh methods for such a miscible flow problem can be undertaken.