Higher-order compositional modeling with Fickian diffusion in unstructured and anisotropic media

Advances in Water Resources 33 (2010) 951–968

Contents lists available at ScienceDirect

Advances in Water Resources j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / a d v wa t r e s

Higher-order compositional modeling with Fickian diffusion in unstructured and anisotropic media Joachim Moortgat a, Abbas Firoozabadi a,b,⁎ a b

Reservoir Engineering Research Institute, 385 Sherman Ave, Suite. 5, Palo Alto, CA, 94306, USA Yale University, New Haven, CT, USA

a r t i c l e

i n f o

Article history: Received 13 October 2009 Received in revised form 18 March 2010 Accepted 26 April 2010 Available online 12 May 2010 Keywords: Mixed hybrid finite element Discontinuous Galerkin Compositional modeling Porous media Heterogeneous media Slope limiter Fickian diffusion Unstructured grids

a b s t r a c t We present advances in compositional modeling of two-phase multi-component flow through highly complex porous media. Higher-order methods are used to approximate both mass transport and the velocity and pressure fields. We employ the Mixed Hybrid Finite Element (MHFE) method to simultaneously solve, to the same order, the pressure equation and Darcy's law for the velocity. The species balance equation is approximated by the discontinuous Galerkin (DG) approach, combined with a slope limiter. In this work we present an improved DG scheme where phase splitting is analyzed at all element vertices in the two-phase regions, rather than only as element averages. This approximation is higher-order than the commonly employed finite volume method and earlier DG approximations. The method reduces numerical dispersion, allowing for an accurate capture of shock fronts and lower dependence on mesh quality and orientation. Further new features are the extension to unstructured grids and support for arbitrary permeability tensors (allowing for both scalar heterogeneity, and shear anisotropy). The most important advancement in this work is the self-consistent modeling of two-phase multi-component Fickian diffusion. We present several numerical examples to illustrate the powerful features of our combined MHFE–dg method with respect to lower-order calculations, ranging from simple two component fluids to more challenging real problems regarding CO2 injection into a vertical domain saturated with a multicomponent petroleum fluid. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction Injection of CO2 in the subsurface for CO2 sequestration and for improved oil recovery has many attractive features. Capture and sequestration of large quantities of CO2 from coal-, oil-, and gas-fired power plants can reduce sustained global warming. CO2 injection in hydrocarbon reservoirs can improve oil recovery, providing economic benefits. One of the issues in subsurface modeling of CO2 relates to the density increase when it dissolves in water and petroleum fluids [3,59]. CO2 can also decrease the viscosity and allow higher liquid mobility. The presence of CO2 can cause vaporization of species from the petroleum fluids. These aspects add to the complexity of the numerical simulation of flows in permeable media in the subsurface. One important requirement for subsurface modeling of CO2 is the knowledge of the distribution of CO2 and the effect of heterogeneity (such as fractures) and anisotropy. Due to serious numerical dispersion and gridding, current models based on first-order finite difference and finite volume methods have low numerical accuracy. The use of higher-order methods has the potential to reduce the ⁎ Corresponding author. Yale University, New Haven, CT, USA. E-mail addresses: [email protected] (J. Moortgat), abbas.fi[email protected] (A. Firoozabadi). URL: http://www.rerinst.org (J. Moortgat). 0309-1708/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2010.04.012

numerical dispersion and grid orientation effect but has not been successful until recently. Part of the problem in higher-order methods is the need for closer attention to physical concepts: many physical non-linearities and instabilities can be dampened by the excessive numerical diffusion in low-order approach but need to be properly understood and distinguished from numerical artifacts in a higherorder treatment. In earlier papers, [34–36] Hoteit and Firoozabadi review the literature on high-order methods for numerical modeling of two-phase multicomponent flows and demonstrate powerful features of the combined discontinuous Galerkin and mixed finite element methods. Issues of the consistent mathematical formulation of upwinding in two-phase multicomponent flows are discussed in [46]. In the current paper, we report new key advances in numerical modeling of two-phase multicomponent flow in permeable media. We expand the previous work to unstructured grids and anisotropic media where numerical dispersion and grid orientation effects become even more pronounced. To reduce those two artificial effects, we propose to use a true discontinuous Galerkin method, where nodal values of species composition are captured both in single- and two-phase flows. In previous attempts [35,36], the formulation is based on the evaluation of cell average values in two-phase regions. Use of nodal phase compositions introduces new complexities in the use of slope limiters to prevent spurious oscillations. Our proposal increases the order of accuracy of the two-phase regions,

952

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

and therefore of the entire method when both single- and two-phase regions are present. Fickian diffusion has received relatively little attention in the literature on compositional modeling of oil reservoirs, partly due to the lack of experimental data for the full multi-component matrix of cross-species diffusion coefficients. Instead, correlations have been used that assume an effective diffusivity in which the diffusive flux of component i only depends on its own compositional gradient. This may be an oversimplification, even for three-components systems, as was already shown experimentally in [24]. In this work, we present new results when Fickian diffusion is accounted for. We calculate the diffusion coefficients from irreversible thermodynamics, as a function of temperature, pressure and phase composition, based on the work in [44]. Various discontinuous Galerkin methods have been developed for elliptic diffusion equations and convection–diffusion equations, such as those in [4–7,9,10,13,18, 40,42,54], with a clear review given in [58]. Historically, these methods have been developed in two seemingly disparate classes of ‘primal’ methods that use ‘internal penalties’ or bilinear forms at element edgesor ‘local DG’ (LDG) methods that enforce mass continuity across edges through the construction of certain ‘numerical fluxes’. In [2], however, a uniform framework was presented demonstrating that these two approaches are fully equivalent, and each bilinear form can be translated into a numerical flux, and vice versa. More recent work proves optimal convergence of the original LDG method on rectangular grids in [14], and has produced a new promising class of hybridized global DG methods with superconvergent properties [15–17,49,50]. The adoption of any one of these methods for the purpose of reservoir simulation, unfortunately, is not straightforward. The mathematical analysis of stability and convergence properties in all of the above references is performed under various simplified assumptions of, for instance, lower dimensionality, steady state, convection- or diffusion-dominated, periodic boundary conditions or constant diffusion coefficients, and do not necessarily hold up for multi-dimensional, multi-phase, multi-component convection–diffusion transport. The diffusive flux is particularly challenging because the matrix of diffusion coefficients is not positive definite, and its determinant may vary by an order of magnitude due to compositional effects. In [34] it was shown that the LDG method is stable for two-phase convective flow, and we find that this method also gives the most stable results when diffusion is included [13]. Diffusion is believed to be important in heterogeneous media, such as highly fractured reservoirs. Just like gravity, it can drive gas from the fractures into the rock matrix resulting in enhanced recovery of matrix oil [38]. In [47] we simulate laboratory experiments where CO2 is injected in a fractured core, and show that the recovery and composition in the production stream can only be accurately modeled by Fickian diffusion. In this paper, we show for the first time the importance of Fickian diffusion even in homogeneous media, where it has a significant effect in suppressing the formation of density-driven fingers. This is confirmed in our work on modeling laboratory experiments in unfractured cores (to be presented in a future publication), where again the experimental results can only be reproduced by accounting for Fickian diffusion. We do not consider capillarity, because in oil–gas systems the capillary pressure is often negligible due to the low interface tension at high reservoir pressures [30]. Capillarity is reduced even further in CO2 injection scenarios, because of the high solubility of CO2 in the oil and associated phase behavior effects. While the former reason applies to relatively homogeneous media, the latter effect may be true even in fractured media. This paper is structured along the following lines. We first present the problem formulation for multi-component two-phase flow in anisotropic permeable media in Section 2. We then cover the numerical solution of the mathematical problem using the mixed hybrid finite element (MHFE) method and discontinuous Galerkin (DG) method

(Section 3) and review details of the stability testing, phase-split calculations and the implementation of the slope limiter. Various numerical examples are presented in Section 4 to demonstrate improved accuracy and the efficiency and robustness of the proposed algorithm. The paper concludes with a short summary and conclusions in Section 5. 2. Mathematical model We summarize the system of equations and boundary and initial conditions needed to model the convective and diffusive flow of a multi-component two-phase fluid through porous media. In reservoir simulations, mass and momentum conservation are customarily expressed by the species balance equation and Darcy's law. ϕ

∂czi + ∇⋅Ui = fi ; ∂t

i = 1; …; nc ;

ð1Þ nc

ϑα = −λα K ð∇p−ρα gÞ; where ρα = cα ∑ xi;α Mi : i=1

ð2Þ

Here ϕ is the porosity, c is the overall molar density and K the permeability tensor. For each phase α = (o, g), i.e. oil or gas, ϑα is the phase flux, cα the molar density, λα(Sα) the mobility, Sα the saturation, ρα the mass density, g the gravitational acceleration and xi,α the mole fraction of component i. Furthermore, for each of the nc components i we have the total molar flux Ui, the overall mole fraction zi, the molar weight Mi and the sink and source term fi. In previous work, we have used a scalar permeability in modeling homogeneous and isotropic domains. Here we will extend this to include anisotropy and heterogeneity. Note that the relative permeability, krα, is incorporated in the mobility λα = krα / μα, with μα the phase viscosity. The total molar flux in Eq. (1) consists of the convective phase fluxes ϑα and the diffusive phase fluxes Ji,α: Ui =

∑

α = o;g


 cα xi;α ϑα + Sα Ji;α ; i = 1; …; nc

nc −1

Ji;α = −cα ∑ Dij;α ∇xj;α ; j=1

i = 1; …; nc −1;

nc

ð3aÞ

ð3bÞ

ð3cÞ

∑ Ji;α = 0:

j=1

The Fickian diffusion coefficients in the square (nc − 1) component matrices Dα = Dij,α are found using the unified model presented in 1 [44], and can be written as Dα = B− α Γα in terms of a matrix of thermodynamic factors Γα and a matrix Bα, both of the same square size and given by: Γij;α = xi;α

Bii;α =

∂ lnFi;α ∂xj;α

j

xj;α ;T;p

;

i; j = 1; …; nc −1

nc xi;α xk;α + ∑ ; i = 1; …; nc −1; Dinc ;α Dik;α

ð4Þ

ð5Þ

k=1 i≠k

Bij;α = −xi;α

! 1 1 − ; i; j = 1; …; nc −1; i≠j: Dij;α Dinc ;α

ð6Þ

Here, T is the temperature, Fi,α is the fugacity of component i in phase α and Dij;α are the Stefan–Maxwell diffusion coefficients for each binary pair i–j in phase α (with Dij;α = Dji;α ). In [44], the authors propose a correlation based on experimental data to find the Stefan– Maxwell diffusion coefficients in the infinite dilution limit, denoted as

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

D∞ ij,α. The Stefan–Maxwell diffusion coefficients can be derived from D∞ ij,α by the Vignes relation, generalized to multi-component mixtures in [43] as
x
x x = 2 nc
j;α i;α k;α ∞ ∞ ∞ ∞ Dij;α = Dij;α Dji;α ∏ Dik;α Dkj;α :

ð7Þ

field v). In the general cases that we will study here, the problem is highly coupled and non-linear. The system of equations is closed by a pressure equation —derived from volume balance considerations [1,60]– and the Peng–Robinson equation of state [52]:

k=1

ϕcf

i≠k

The phase molar densities and compositions are derived under the assumption of local thermodynamic equilibrium, which implies that for each component the fugacities in the two-phases are equal:

 Fi;o T; p; x1;o ; …; xnc −1;o = Fi;g T; p; x1;g ; …; xnc −1;g

i = 1; …; nc : ð8Þ

We recast Darcy's law in terms of the total flux: ϑ = ∑α ϑα = −λt K ð∇p−ρg Þ; with : λt = ∑α λα ; ρ = ∑α fα ρα and fα =

ð9Þ

λα ; λt

in terms of the total mobility λt and the fractional flow functions fα . The phase fluxes can be recovered — independent of the pressure — from the total flux by:

nc ∂p P + ∑ v i ∇⋅Ui = fi ; ∂t i=1

cα =

The phase mobilities (and thus fα) depend on the saturation through the relative permeability λα ∝ krα(S), where the saturation is defined in terms of the gas molar fraction ϖ as: ϖc c Sg = ; and So = ð1−ϖÞ ; with So + Sg = 1: cg co

ð16Þ

where cf and v̄i are the total compressibility and total partial molar volume, respectively, R is the universal gas constant and Zα depends on temperature, pressure and phase composition [29]. To define a proper physical model, the above equations have to be supplemented with the appropriate boundary and initial conditions on the computational domain Ω. We will denote the non-overlapping Dirichlet and Neumann boundary conditions by ΓD and ΓN and the boundary of Ω by Γ = ΓD ∪ ΓN. Initial conditions are provided for composition and pressure (and position x = (x,y)): 0

zi ðx; 0Þ = zi ðxÞ in Ω; 0

ð10Þ

ð15Þ

p ; Zα RT

pðx; 0Þ = p ðxÞ;

ϑα = fα ðϑ−Gα Þ; ( λo ðρo −ρg ÞKg; α = g; Gα = λg ðρg −ρo ÞKg; α = o:

953

ð17aÞ

in Ω;

ð17bÞ

whereas the boundary conditions of the velocity field and pressure Eq. (9) are given in terms of the imposed injection rate and fixed pressure, qN and pD, respectively, as: N

N

cðx; t Þzi ðx; t Þϑðx; t Þ⋅n = qi ðx; t Þ on Γ : D

D

pðx; t Þ = p ðx; t Þ; on Γ ;

ð18aÞ ð18bÞ

ð11Þ 3. Numerical model

The compositions zi and xi,α are also constrained, by: nc

nc

nc

i=1

i=1

i=1

∑ xi;o = ∑ xi;g = ∑ zi = 1;

ð12aÞ

ϖxi;g + ð1−ϖÞxi;o = zi :

ð12bÞ

As the molar volume occupied by the two-phase fluid is not equal to the sum of the molar volumes of the individual fluids in single phase, the total molar density also is not simply the sum of the phase densities. Instead, from Eq. (11), we can express the total molar density in terms of the phase molar densities as c=

co cg : ϖco + ð1−ϖÞcg

ð13Þ

Definitions in Eq. (9)–(13) suggest an alternative form of the species balance equation Eq. (1) for i = i, …, nc − 1: " ∑ ϕ α

nc −1 ∇xj;α ∂Sα cα xi;α f ðS Þ + ∇⋅Sα cα xi;α α α ðϑ−Gα Þ− ∑ Dij;α Sα xi;α ∂t j=1

!# = fi ;

ð14Þ which more explicitly shows some of the non-linearities of the problem. As an example, we generally assume a power-law dependence for the relative permeability as krα = |krα | (Sα)p. Then, only for a linear dependence p = 1 and in the absence of gravity and diffusion, does the system Eq. (14) (plus Darcy's law) have the same order of complexity as the well studied single phase flow problems where ∂u/∂t + ∇ · uv = f and v ∝ ∇u (for some state variable u and velocity

We will solve the main equations presented in the previous section, Eqs. (1), (9), and (15) using an implicit pressure, explicit composition scheme. More specifically, Darcy's law is solved by the mixed hybrid finite element method in the lowest-order Raviart– Thomas approximation space [53] and the transport equation by the linear or bilinear discontinuous Galerkin method. Many details of this approach can be found in previous work in [34–36,46] and will not be reiterated here, unless for self-consistency. Rather, we focus on our generalization to unstructured grids, an improved bilinear method for rectangular grids, including Fickian diffusion and considering arbitrary permeability tensors. Some of the algebraic manipulations are deferred to Appendices A–C. We will mainly consider unstructured triangulations of the domain Ω, where mesh cells will be denoted by K with boundary ∂K made up of 3 edges E. NK is the total number of edges and NE the number of edges not belonging to ΓD. 3.1. Mixed hybrid finite element method Darcy's law Eq. (9) together with the pressure equation Eq. (15) are solved implicitly for the pressure using the MHFE method. The DG method may be more accurate for convection-dominated problems, but the MHFE method provides higher-order convergence for the flux variable than DG, with the same degree of polynomials. In the conforming finite element model, only the pressure is calculated directly, and the velocity field is found by differentiation. This leads to an element-wise constant velocity with no continuity imposed on the normal component through element edges, which cannot properly model the physics [48]. Both the mixed finite element method (MFE) and mixed hybridized (MHFE) models solve for the pressure and

954

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

velocity simultaneously and to the same order. However, the former leads to a saddle point problem, whereas the latter was developed to lead to a positive definite linear system by appending as extra degrees of freedom the averaged pressures at the element edges [8,11,12]. This positive definiteness of the linear system unfortunately breaks down for the more complicated (compressible and non-conservative) pressure equation (Eq. (15)). The MHFE method, however, still provides a more elegant and efficient mixed method compared to the fully equivalent non-hybridized method, as demonstrated in [12]. The essential part of the MHFE method is to approximate the convective and gravitational fluxes by a superposition of the linear basis vector functions of the lowest-order Raviart–Thomas space [53], with the normal components of the flux across element edges as its Kg coefficients, denoted by qK,E (and qα,K,E, qK,E): Kg

ϑðt; xÞ = ∑ qK;E ðt ÞwK;E ðxÞ; K·g = ∑ qK;E wK;E : E2∂K

ð19Þ

E2∂K

The advantage of expanding of K · g, rather than g, will be apparent shortly. To apply our model to unstructured triangular grids, we employ the basis vectors on the unit triangle:  w1 =

 x ; y−1

  x w2 = ; y

 w3 =

 x−1 : y

Since both the MHFE and DG methods are strictly local on each element K, we will only need basis (vector) functions on the unit triangle. Quantities on the corresponding triangle in the mesh are obtained by a simple affine mapping (see Appendix A). Our choice of node and edge numbering in each of these triangles is illustrated in Fig. 1. The approximation to Darcy's law then proceeds in the following steps: −1

1. multiply Eq. (10) by (λtK) wK,E (the total velocity concept is used to guarantee that λt is positive definite, whereas λα can vanish), 2. transform the equation into the weak form by integrating over K, 3. integrate by parts to turn the integral over the pressure gradient into a surface integral over p and a volume integral over the divergence of w and 4. use the divergence theorem on the latter integral. 5. Finally, expand ϑ and K · g as in Eq. (19) to find for each edge E: ∑

E′2∂K

qK;E′ −1 ∫ w ·K ·wK;E′ λK;t K K;E Kg

= pK −tpK;E + ∑ qK;E′ ρK ∫K wK;E ·K ′2

E ∂K

−1

·wK;E′ ;

ð20Þ

with the definitions Eq. (21a) below and the total mobility, permeability and density assumed element-wise constant. By defining a gravitational flux as K · g, only one integral has to be evaluated to describe both the mesh- and the physical geometrical properties, in terms of the basis vector fields and the permeability tensor, respectively (see Appendix A).

In Eq. (20) the element and edge averaged pressures are defined in terms of the surface of the element A and the length of the edge L as: 1 pK ≡ ∫K p; A

1 tpK;E ≡ ∫E p; L

qK;E′ = ∫E ϑ·nK;E′ ;

Kg

qK;E′ = ∫E ðK·gÞ·nK;E :

ð21bÞ

Multiplying Eq. (20) with the inverse of the integral over the permeability, results in an explicit expression for the normal components of the total flux through each edge E (where the coefficients can be found in Appendix A): qK;E = αK;E pK − ∑ βK;E;E′ tpK;E′ −γK;E ; E2∂K: E′2∂K

ð22Þ

By imposing the boundary condition Eq. (18a) and flux continuity through edges E = K ∩ K′, such that qK,E = − qK',E we can eliminate the flux qK, E by summing Eq. (22) over each two neighboring elements to obtain a linear system for the entire domain Ω: T

ð23Þ

R P−MTP = V;

where the elements of R, M and V are simply the coefficients in Eq. (22) summed over neighboring elements [34]. The pressure equation Eq. (15) is expressed in terms of the flux qK,E through the definition of the phase velocities. Since we will treat the pressure implicitly, we adopt the backward Euler scheme on the time derivative and write for the new time-step n + 1: DP

n + 1

˜ n + 1 = G: − RTP

ð24Þ

The matrices D, R̃ and G are given in Appendix A in terms of the coefficients in Eq. (22) (evaluated at the previous time-step) with additional dependence on the anisotropic permeability in G. Finally, we can use Eq. (23) to eliminate P n + 1 from Eq. (24) and obtain our final equation for the main variables in the MHFE model, the pressure traces at the element edges, TPn + 1:
 T −1 n M−R D R˜ TP

+ 1

T

=R D

−1

G−V;

ð25Þ

where we use the fact that D is diagonal and invertible [33]. Subsequently, the solutions TPn + 1 are used in Eq. (24) to find Pn + 1. Details of the upwinding, necessary in evaluating the divergence term in Eq. (15), are discussed in [46,55]. 3.2. Discontinuous Galerkin method To solve the transport equation, we will approximate the solution to Eq. (1) on each element by a superposition of linear (on triangles) or bilinear (on rectangles) functions. Eq. (1) is a second order PDE in terms of the compositions. The local discontinuous Galerkin method is a mixed method, in which the problem is computationally simplified by rewriting Eq. (1) instead as a system of first order PDE's ϕ


 ∂czi + ∑ ∇· cα xi;α ϑα + Sα Ji;α = 0; ∂t α

ξi;α + ∇xi;α = 0; nc −1

Ji;α −cα ∑ Dij;α ξi;α = 0; j=1

Fig. 1. Node and edge numbering on unit (left) and mesh (right) triangular elements.

ð21aÞ

ð26aÞ ð26bÞ ð26cÞ

in terms of the auxiliary variables ξi,α and Ji,α [13]. The advantage is that the above equation can be solved individually by the traditional DG techniques as developed for the convection equation. The auxiliary

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

variables should not be seen as increasing the problem size, and are only introduced as a matter of convenient. They can be eliminated when one prefers to solve a global system of equations, such as in [15]. Since the solution to Eq. (1) for each separate component is decoupled and, furthermore, the summation over the two phases α in Eq. (1) is trivial, we will introduce the more succinct notation z = czi, x = cα xi,α and Jdiff = Sα Ji,α and write the equations for phase α. We provide here some details specific to triangular elements, used in our extension to unstructured grids. The DG method consists of approximating the unknown variables in a given partial differential equation by a superposition in terms of the basis functions of the chosen approximation space. As the method is discontinuous and fully local per element, any order of basis functions can be used, unlike for continuous methods, and we can map each element to the unit triangle (see Fig. 1 and Appendix A). We will choose the linear function space with bases: φ1 = y;

φ2 = x; and φ3 = 1−x−y:

ð27Þ

In terms of these basis functions, the concentrations can be approximated by: 3

3

i=1

i=1

z = ∑ zK;i ðt Þφi ðx; yÞ; x = ∑ xK;i ðt Þφi ðx; yÞ;

ð28Þ

where zi,K and xi,K are the nodal values of z and x in element K (note that we evaluate Eq. (1) on the unit triangle, so the φi are the same for every element K). The next steps are similar to those in the MHFE method discussed in Section 3: 1. multiply Eq. (1) by appropriate weight functions, which are chosen to be the same φi, 2. integrate by parts over each element K to put the equation in its weak form, 3. use the divergence theorem and expand z, x and ϑ in terms of the bases ϕi and wi. For each element we obtain at each node i (summation over the index j is implied): ϕK

dzK;j ∫ φφ− dt K i j

xK;j ∑ qK;E ∫K φj wE ⋅∇φi +

ð29aÞ ð29bÞ

E2∂K

x˜E;j ∑

E2∂K

qK;E ∫ φφ + jE j E i j

∑ Jdiff;K;E ∫E φi = ∫K φi fK ;

ð29cÞ ð29dÞ

E2∂K

where ϕK and fK are assumed element-wise constant, x̃E,j is the upstream value of x at node j with respect to the phase velocity ϑα through edge E (discussed in more detail in Section 3.2.2), and Jdiff, K,E= Jdiff · nE is the averaged diffusive flux through edge E. We assume here that Jdiff,K,E is continuous across edges, and compute ξi,α in Eq. (26b) from the difference (gradient) in average phase compositions in the two elements neighboring edge E. Notice that Eq. (29a) is an ordinary differential equation since the spatial dependence is captured in the shape functions φ. The integrals in Eq. (29a) can be evaluated algebraically and the mass matrix in Eq. (29a) is invertible so that we can solve Eq. (29a) explicitly using the forward Euler derivative. Details are given in Appendix B for triangular meshes, and generalized further to bilinear DG on rectangular meshes in Appendix C. Higher-order approximations of the time derivative, such as proposed in [28,37,41] for non-compositional problems may also be

955

used. We also note that fully implicit methods are often employed in solving advection–diffusion problems to avoid the more stringent CFL stability constraint on the time-step in treating the parabolic diffusion term. This approach, however, has been developed mainly for single phase flow, and incompressible two-phase flow without mass transfer between the phases [25–27]. These methods cannot be generalized to compositional two-phase flow because the stability analyses (and accurate phase splitting calculations) are highly non-linear and require iterative solution methods that cannot be implemented in implicit schemes. Besides, the larger time-steps allowed by unconditionally stable implicit methods does not guarantee accuracy of the results for time-steps larger than the CFL condition, because the pressure equation is strongly coupled to mass transport in two-phase flow. Finally, a fully implicit scheme would require the simultaneous solution of a much larger system of coupled equations which may ultimately be more computationally expensive, given a desired accuracy, than an efficient explicit scheme updated at much smaller time-steps, especially since the latter can be massively parallelized. 3.2.1. Stability and phase-split calculations In this section we discuss how we obtain the phase molar density cα and molar fractions xi,α needed in the transport equation. First, a stability analysis is performed to determine whether the element is in single or in two-phase, based on the element pressure, temperature and composition zi from the previous time-step [29]. When the flow is in single phase, the xi,α are simply equal to the zi, cα = c, Sα and ϖ are either zero or one, and the thermodynamic equilibrium equation Eq. (8) is irrelevant. In this regime, the composition at the nodes zi can be updated self-consistently from its values at the previous time-step, using Eq. (29a). Single phase regions are treated the same as in [35]. When the element is in two-phase, phase-split calculations (referred to as flash) have to be performed to determine the xi,α. In [35,36], this was only done once per element to speed up the calculation, using the cell average pressure, temperature and composition. However, this effectively reduces the simulation of the two-phase regions to a lowerorder finite volume treatment, leading to increased numerical dispersion. In the examples presented in [35] of binary and tertiary mixtures, the two-phase region was limited to a relatively sharp shock front and the speed up in calculations may have justified the cruder treatment of the two-phase region. In this work, we will present more complicated problems such as CO2 injection into a multi-component petroleum fluid. In reservoirs, the two-phase region may cover large parts of the domain and the flash calculations and stability have to be done more carefully. More generally, the adoption of higher-order methods is precisely motivated by a need to accurately capture the shocks or phase boundaries, and a more accurate treatment of the two-phase regions may increase the order or accuracy of the entire method (a similar argument holds for too much slope-limiting at discontinuities). We propose that stability analysis and flash be performed at each of the nodes, using the nodal pressure, temperature and composition. First, we determine whether the node is in single or two-phase. When in twophase, we calculate the gas molar fraction ϖ and the compositions in each phase xi,α using a number of techniques that can be found in [29] and are briefly outlined at the end of this section. Subsequently, the phase molar densities cα at the nodes are calculated using the Peng– Robinson equation of state. For the saturation, compressibility, partial molar volume and phase mobilities, we only need element averaged values in the MHFE part of our algorithm, which requires accurate element-average values of cα and xi,α. A simple extrapolation or average of the nodal xi,α, however, is non-trivial, especially when some nodes are in single phase and others in two-phase. Therefore, to maximize accuracy and stability of the code, we perform an additional stability and flash calculation based on the cell averaged pressure, temperature and composition to find the element average cα and xi,α. These average values are only used in the MHFE update of fluxes and pressures, while the nodal values are used for the DG update of the overall compositions.

956

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

For an unstructured triangular grid, this means a maximum of four stability and flash calculations per element (less when using the procedure at the end of this section). In [34–36] a linear DG approximation was used in terms of only three degrees of freedom for the overall composition: the composition at the center of the element, and its gradients in the horizontal and vertical directions. The overall compositions on the edges were found by extrapolation, while the phase compositions were assumed element-wise-constant. To treat two-phase regions with the same accuracy as single phase regions, phase compositions are required at all nodes (or edge centers), even when using linear DG for the overall composition, with a fifth computation to find averaged properties that are used in the MHFE update. The method presented here was adopted in [46] for rectangular grids, where five flash calculations were performed for the average phase compositions and the phase compositions at the edge centers. However, the same three degrees of freedom as in [34–36] were used for the overall compositions. In two-phase fluids, a bilinear DG approximation with the same four degrees of freedom for both the overall and phase compositions may be more suitable, given that the computation of the phase compositions is far more CPU expensive that the algebraic DG update of the overall compositions. The increase in CPU time due to the full DG calculation of the phase behavior is moderate on large grids when the CPU time is dominated by the linear solver in the implicit pressure derivation. For the flash calculations, we have developed a very fast flash algorithm that proceeds in the following steps: • Stability analysis provides a first estimate of the ratios Ki = xi,g/xi,o. • The Ki from stability are used in the successive substitution iterative (SSI) procedure, which does not need accurate initial guesses, to provide more accurate values for Ki. • The Ki from SSI are used as initial values for the Newton method to calculate accurate Ki and xi,α in only very few iterations (in fact, we use a modified method that uses the logarithm of xi,o for more stable behavior at small xi,o). • If the Newton method fails (in rare cases where SSI does not provide accurate enough initial guesses for Ki), the slower but robust SSI is used until the end. Further innovations are used to avoid having to perform either stability or flash calculations or both in most elements: • When Ki values are available from the previous time-step, the element is assumed to be in two-phase (without performing stability). A flash calculation is then attempted directly and only if that fails do we return to stability analysis. • If no Ki values are available from the previous time-step because the element was in single phase, and all the surrounding elements were also in single phase, then we directly update the composition from the zi at the previous time-step. • Finally, within each element, the first Ki values to be calculated can be used as more accurate initial input for the flash calculations at the other nodes (although, depending on the grid refinement and the gradients in composition, the nodal values at the previous time-step can be more accurate initial guesses). The flash module is both stabilized and optimized for speed by automatically adapting the maximum allowed number of iterations nit,max for each flash method, to ensure that the resulting Ki and compositions are always within the expected range while keeping nit,max as small as possible. 3.2.2. Improvements in the DG scheme After digressing in the previous section on the phase behavior calculations, we continue to summarize additional progress we have made in the DG treatment of the transport equation. An advantageous side-effect of using nodal state variables on the triangular grids, rather than the edge based degrees of freedom that

were used on the structured rectangular meshes, is that the upwinding is improved: instead of transferring only the edge center upstream density and composition, upwinding is done at two nodes per edge, which also communicates information on the upstream gradients in density and composition. As a result, the sensitivity to mesh orientation significantly reduces, as will be demonstrated in the examples in Section 4.1. It is well known that for the DG method to be stable, a slope limiter has to be used to avoid spurious oscillations in the solution when it takes on extremal values at the nodes. In other words, the composition zK has to be reconstructed such that its nodal values zK,j lie between the extremal values of the (average) zK in all the elements sharing the node j. Instead of the simple standard minmod function to limit the slopes of zi in each element, we have implemented the more sophisticated and well tested slope limiter presented in [32]. This concludes the discussion of the improvements in our compositional modeling. A summary of the steps in the full algorithm is given in Appendix D. 4. Results In this section we will discuss several examples to demonstrate unique features and advances in our compositional modeling algorithm. Unless explicitly mentioned otherwise, all simulations are performed using the updated full discontinuous Galerkin treatment of the transport equation presented in Section 3.2.1. Different linear solvers may be used to solve Eq. (25). Direct methods are very general and robust and can be more efficient in cases where finding and computing a good preconditioner for iterative methods can be computationally expensive. However, iterative methods have to be used when the system size becomes very large. For the relatively small examples presented here, we have tested both the unsymmetric-pattern multifrontal method for sparse LU factorization [20–23] and the (parallelized) Pardiso interface [51,56,57]. The latter is easy to implement (and is included with the Intel compiler as part of the Math Kernel Library), and ranks highly in a comparison of ten widely used solvers in [31]. For a structurally symmetric system, such as in Eq. (25), it first computes a symmetric fillin reducing permutation followed by parallel numerical factorization and partial pivoting in the supernodes. An approximation of the solution is found by forward and backward substitution and iterative refinements. For the numerical examples discussed in this section, both solvers resulted in near identical CPU times (where Pardiso was used in serial mode). Indications of CPU times for the simulations are summarized in Table 1. The code was compiled with the Intel ifort compiler and the computations were performed on a 2.8 GHz Intel Core 2 Duo iMac with 2 GB SDRAM. Table 1 General indication of CPU times. CPU times FVM

DG-FVM

DG

Bilinear DG

Example 1 Mesh 1 Mesh 2 Diffusion 20 × 20 40 × 40 80 × 80

b1 min b1 min b1 min – – –

b 1 min 1 min b 1 min – – –

1 min 2 min 1 min – – –

– – – b1 min 2 min 27 min

Example Example Example Example Example Example Example

– – – – – – –

– – – – – – –

1.5 h 2 min 1.5 h 11 min ∼5.5 h ∼2 h ∼2.5 h

– – – – – ∼5 h ∼5.5 h

2a 2b 2c 3 4a 4b 5

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968 Table 2 Simulation parameters for binary methane, propane mixture.

Injected gas [mole fraction] Initial fluid [mole fraction] Acentric factor Critical temperature [K] Critical pressure [bar] Molar weight [g/mol] Critical volume [m3/kg] Volume shift parameter Binary interaction coefficient Initial pressure (single phase and two-phase) Initial temperature (single phase and two-phase) Porosity and permeability

C1

C3

1.0 0.0 0.011 189.7 45.8 16 0.0061 0.0

0.0 1.0 0.15 369.8 42.5 44 0.0045 0.0 0.037 69 bar 311 K 10 md

50 bar 397 K 0.2

For simplicity, we assume a relative permeability of the form krα = (Sα)2 in all the examples. 4.1. Example 1: comparison of the different order methods We will start by demonstrating the significant improvements in both reducing numerical diffusion and mesh dependence, achieved by using the full discontinuous Galerkin approximation (denoted by DG) to the species balance equation and (physical) Fickian diffusion. We compare our current method to both the finite volume approach (FVM), and the ‘mixed’ discontinuous Galerkin method in which flash calculations are performed only for element averaged quantities (DG–FVM). Note that the latter two methods only refer to the mass transport update, and that we use the MHFE method for the pressures and fluxes in all three methods. The results should, therefore, still be more accurate than those from a pure FVM. 4.1.1. Problem set-up To best isolate the characteristic features of the problem, we consider a simple two component system where methane (C1) gas is injected into the lower right corner of a horizontal 50 × 50 m2 domain that is initially saturated with liquid propane (C3). The production well is in the opposite (upper left) corner, and the initial pressure and temperature are chosen such that a two-phase region develops behind the shock front. The initial compositions and other relevant parameters are given in Table 2. We use the same code to simulate all three models, but in the FVM case we only calculate and update element averaged values for c, zi, cα and xi,α. In the mixed DG–FVM, we calculate nodal values for czi = cαxi,α in the single phase regions. However, we only perform one flash and stability calculation per element in the two-phase region to obtain cα and xi,α, which essentially reduces the method to a FVM. Finally, for the full DG method, the phase behavior is analyzed at all the nodes. To study mesh effects, we compare two different grids with the same number of elements (1572). The first mesh (mesh-1) is created using an optimized algorithm that produces high quality elements (acute triangles) and a high degree of structure. The second mesh (mesh-2) is subsequently obtained by randomly perturbing each vertex by up to roughly Δx/2, with Δx the average triangle edge lengths, resulting in many low quality elements and reduced regularity. In all the runs in this section, the injection rate is 42.5 m2/day at surface conditions,1 or 0.43 PV/yr. 4.1.2. Results Figs. 2–4 compare the three different order methods at 60% pore volume injection (PVI), based on the methane composition (in molar fraction contours). The effect of Fickian diffusion, compared to results neglecting diffusion (Fig. 2), is shown for the two different meshes in Figs. 3–4. 1 We denote our injection rates in m2/day as the simulations are performed in 2D. Equivalently, we could assume a 1 m thickness of the domain and a m3/day rate.

957

The physical system has an almost constant (plateau) composition in the two-phase region with sharp fronts at the single to 2-phase interfaces. This is illustrated schematically in Fig. 5, which shows a surface plot for the methane composition (also at at 60% PVI), computed using the DG method on a high quality 104 element triangular mesh (shown), including diffusion.2 The fronts are captured accurately on the coarse grids by the DG method and to a lesser extent by the DG–FVM in the runs both with and without Fickian diffusion and on both meshes. The FVM on the other hand suffers from excessive numerical diffusion and the two-phase plateau is all but smeared out in all three cases. The other improvements are related to mesh-sensitivity. It is well known that the upwinding performed in updating the transport equation introduces a slight dependence on the mesh geometry: values for the state variables are taken from upwind neighbors which is sensitive to the orientation of the element edges. When the mesh is highly structured, the slight enhancement in the resulting velocities along axes of symmetry can accumulate and lead to unphysical orientation of the velocity field along those axes. For a rectangular grid this leads to overestimating the flow along horizontal and vertical lines.3 On a high quality triangular grid the orientation is along the domain boundaries and diagonal lines and on low quality grids the solution can further be sensitive to mesh irregularities. In the higher-order DG method where upwinding communicates information on both the state variable and its gradient (as discussed in Section 3.2.2) this mesh dependence is greatly reduced. In fact, with Fickian diffusion the scheme is further stabilized and the mesh dependence is completely eliminated, as can be seen in the right panels of Figs. 3–4. By comparison, both the FVM and the DG–FVM methods show a very strong dependence of the solution on the orientation of the mesh diagonals, as compared to the ‘physical’ propagation direction which should be along the diagonal from the lower right to the upper left corners of the domain. In comparing the left two panels in Figs. 2–3, one can observe that the results are improved by including the Fickian diffusion. However, even when the symmetry is reduced by perturbing the element vertices (Fig. 4) the orientation effect is significant for the lower-order methods. In comparing the DG–FVM and DG methods, it is apparent that the front between the single phase liquid and two-phase regions is essentially identical for both methods. Conversely, at the miscibility front between the single gas phase and two-phase regions, the phase behavior is determined by the highly non-linear two-phase region, which in the DG–FVM approach is treated to lower order than in our full DG model. It is clear from the figures that convergence to the physical solution, i.e. without the strong grid orientation effects, is only attained in the full DG treatment. To further confirm this, we show in Fig. 6 the results of a bilinear DG simulations of the same example, including diffusion, on 20× 20, 40× 40 and 80× 80 rectangular grids. We find close agreement between the solution on the latter two grids and those in the right panel in Fig. 3. All the following examples will use the linear or bilinear DG methods only.

4.2. Example 2: heterogeneous media In this example we will present results for two-phase flow in heterogeneous media where the permeability tensor is assumed diagonal with Kyy = Kxx and we only vary the scalar magnitude of the permeability in different regions of the domain.

2 The plot is rotated with respect to the Figs. 3–4, to better illustrate the features in 3D. Injection is in the far corner and production at the front. 3 This is often missed or misinterpreted in the literature since the results on a structured grid still look very smooth. On unstructured grids, the occurrence of mesh orientation is more obvious.

958

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

Fig. 2. Example 1a: methane composition (contours of molar fraction) calculated with FVM (left), mixed DG-FVM (middle) and DG (right) without diffusion at 60% PVI. Also shown are the physical symmetry axis (solid line) pointing from the injection to the production well, and grid orientation (dashed line) on mesh-1.

Fig. 3. Example 1b: as in Fig. 2, but including Fickian diffusion on mesh-1.

4.2.1. Example 2a Consider the flow through a layered domain with 4 layers in which the permeability is 10 md, alternating with 5 layers where ∥ K ∥ = 103

md and each layer is 30 m wide, as in Fig. 7. The physical dimensions of the entire domain are 0.5 km × 0.27 km and the mesh has 3026 elements. Again we suppress additional asymmetry caused by the

Fig. 4. Example 1c: as in Fig. 3 but on the lower quality mesh-2.

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

959

full width of the left boundary and produced at the right boundary where the pressure is kept constant. The injection rate is 842.5 m2/day at surface conditions, injecting one pore volume in 3.5 years. In the right panel in Fig. 7 the gas saturation is shown at 30% PVI. Clearly, the flow in the highly permeable layers is much faster than in those with lower permeability. However, we expect this to change when we include capillarity, similar to the results presented in [37] for incompressible and immiscible two-phase flow. There, the contrast in the capillary pressure between the high- and low permeability layers facilitates cross-flow between the layers, which slows down the flow in the highly permeable layers to essentially the same rate as in the low permeability regions. To a lesser extent, diffusion can play the same role, as will be demonstrated in the next example.

Fig. 5. Example 1: Schematic illustration of the physical solution for the methane composition (calculated on a 100 × 100 triangular grid, including diffusion).

gravitational field by performing all runs in this section on horizontal domains. The fluids are the same as in the previous example (with parameters as in Table 2), but now methane is injected through the

4.2.2. Example 2b While the previous example did not include Fickian diffusion, here we do and we revisit one of the examples in [35] of a two-layer horizontal domain consisting of one 50 m × 10 m layer with permeability 0.1 md and an equal-size layer with 10 md. We choose the same conditions (C1 and C3, 311 K, 40 bar, 40 × 20 elements) such that the domain is partly in two-phase. Fig. 8 shows methane (molar) concentration at 40% PVI with and without diffusion. Without diffusion, breakthrough occurs already at this early time, since almost all the methane propagates through the layer with the highest permeability. When diffusion is accounted for, C1 behind the lower front (in the more permeable layer) can diffuse into the less permeable layer (and similarly for C3). As a result of the cross-flow, the front is sharper and breakthrough is delayed in heterogeneous

Fig. 6. Example 1b: as in Fig. 3 but using bilinear DG on rectangular meshes of 400 (A), 1600 (B) and 6400 (C) elements.

Fig. 7. Example 2a. (A) mesh and permeability for layers of 10 md (white) and 103 md (black). (B) gas saturation at 30% PVI. Injection is through the left face, and production at the right face, where the pressure is held constant.

960

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

Fig. 8. Example 2b: two layers of 0.1 md (top) and 10.0 md (bottom), with the boundary indicated by the horizontal dashed line at y = 10 m. Methane composition at 40% PVI is shown without Fickian diffusion (A) and with diffusion (B).

media. This effect is not always recognized in the literature, but can be exceedingly important in highly fractured media.

4.2.3. Example 2c In the last example of scalar heterogeneity, we use the same domain and mesh as Example 2a, and assign to each element a random permeability in the range of 10–103 md (the distribution is generated with the Fortran random number generator, and is not representative of spatially correlated Gaussian permeabilities found in geological media). The permeability is shown in Fig. 9 (left). Methane is injected in the lower left corner and produced at the upper right end, and the fluid is again partly in two-phase. We show in the right panel the gas saturation Sg (in pseudo-color) and velocity streamlines at 95% PVI. The example demonstrates that, even in complex heterogeneous media, we can robustly capture the two-phase shock front with very little numerical dispersion.

4.3. Example 3: anisotropic media While the heterogeneities discussed in the previous section affect flow paths considerably in fractured media, a second type of permeability is incorporated to study geologically sheared rock matrices. In this example we consider such an anisotropic domain in which the permeabilities are different in the horizontal, vertical and (or) diagonal directions. In particular, we choose two permeability tensors such that the first has a higher permeability aligned with the diagonal from the injection to the production well, while the second is orthogonal to that:  Kleft = 10 md ×

5 3

 3 ; 2

 Kright = 10 md ×

 5 −3 : −3 2

ð30Þ

The size of the domain, triangulation and location of injection and production wells are the same as in Example 2c as are the fluid

Fig. 9. Example 2c: a random permeability distribution throughout the domain in the range [1, 100] × K0 = 10 md (A). Gas saturation (contour plot), and streamlines of the velocity field at 95% PVI (B).

Fig. 10. Example 3a-b: methane composition (molar fraction) at 40% PVI in anisotropic media. (A) Example 3a: K = Kleft (Eq. (30)). (B) Example 3b: K = Kright (Eq. (30)).

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

961

Table 3 Relevant date for Examples 4 and 5 (Sections 4.4 and 4.5).

Injected gas (mole fraction) Initial fluid (mole fraction) Acentric factor Critical temperature [K] Critical pressure [bar] Molar weight [g/mol] Critical volume [m3/kg] Volume shift parameter Binary interaction coefficients CO2 N2 C1 C2 − C3 C4 − C5 C6 − C10 C11 − C24 C25+ p, T, porosity and permeability:

CO2

N2

C1

C2 − C3

C4 − C5

C6 − C10

C11 − C24

C25+

1.0 0.0086 0.24 304.1 73.8 44 0.0021 0.060

0.0 0.0028 0.04 126.2 33.9 28 0.0032 −0.289

0.0 0.4451 0.01 190.6 46.0 16 0.0062 −0.016

0.0 0.1207 0.12 327.8 46.5 35 0.0047 −0.095

0.0 0.0505 0.21 435.6 36.1 630 0.437 −0.060

0.0 0.1328 0.42 574.4 25.0 110 0.0043 0.047

0.0 0.1660 0.66 709.0 15.0 212 0.0044 0.149

0.0 0.0735 1.73 891.5 7.6 463 0.0042 0.495

0 0 0.15 0.15 0.15 0.15 0.15 0.08 250 bar

0 0.1 0.1 0.1 0.1 0.1 0.1 403.15 K

0 0.035 0.039 0.047 0.064 0.105 0.2

0 0 0 0 0 10 md

0 0 0 0

0 0 0

0 0

0

Fig. 11. Example 4a–b: CO2 composition at 100% PV injection from the bottom left (production at the top-right). (A) Example 4a: anisotropic permeability Kyy = 5 md, Kxx = 250 × Kyy = 1.25 d, Kxy = 0. (B) Example 4b: isotropic scalar permeability K = 10 md.

properties, except here the flow is in single phase and diffusion is included. Fig. 10 shows the injection of C1 (molar fraction contours) into a liquid C3-saturated permeable medium at 40% PVI for the two orientations of the shear. Even for this very moderate anisotropy, the effect on recovery is significant. In the left figure where the flow is aligned with the orientation of the shear, methane breakthrough occurs at 40% PVI, whereas with the flow orthogonal to the shear, breakthrough is delayed until ∼ 90% PVI, after nearly all the propane is produced. The simple examples presented so far are selected to isolate new features in our model. Next we proceed to more physically relevant problems were several of issues discussed so far interplay in a more complex fashion. In particular, we will analyze the injection of CO2 into a vertical (gravitational field) domain that is initially saturated with a petroleum fluid, described by 8 pure- and pseudo-components. The first example will consider injection from the bottom and production at the top in an anisotropic domain, whereas in the second example we will

study injection from the top and production in the bottom in a homogeneous and isotropic domain.

4.4. Example 4a–b: CO2 injection in vertical anisotropic and isotropic domain from the bottom The last examples are related to our practical work on actual reservoirs, where CO2 (and/or nitrogen) injection is being considered to enhance oil recovery. CO2 has several unique features that are highly advantageous, particularly in oil recovery from fractured reservoirs. When CO2 dissolves in oil, it can reduce the viscosity of the mixture, leading to enhanced flow rates. At the same time, swelling4 of the mixed oil occurs, causing more oil to flow out of the reservoir. Another unusual property of CO2 is that when it dissolves in oil, the density may increase 4 This means that when CO2 dissolves in oil, the resulting mixture occupies a larger volume that the oil alone.

962

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

Fig. 12. Example 4b: as in Fig. 11(B), but simulated using bilinear DG on a rectangular 200 × 40 element mesh.

[3,59]. The effect of local density increases on gravitational drainage is considered in the next example (Section 4.5). At the same time as enabling enhanced oil production, the capture and sequestration of CO2 is clearly desirable in the face of the threat of global warming. A proper understanding of this process and its merits in complex reservoirs is a subject of high interest. We will consider a domain that is initially saturated with a petroleum fluid comprised of three standard components (CO2, N2 and C1) and 5 pseudo-components. The grouping into pseudo-components is performed to keep our compositional modeling manageable in terms of computing time, and is based on amounts and volatility range of the original components. The properties of the pure components are obtained from the literature, while those of the pseudo-components are obtained by molar averages. The last component, representing all hydrocarbons beyond C25, is tuned to match the bubble point temperature of the oil. All fluid parameters used in our model are listed in Table 3, and we assume a connate water saturation of 20% and a residual oil saturation of 30%. The (vertical) domain in these examples has physical dimensions 0.5 × 0.1 km2 with 4610 computational elements and is strongly sheared in the horizontal direction, as is often found in geological formations, with permeability:  K = 5 md ×

 250 0 ; 0 1

ð31Þ

CO2 is injected at a low rate of 0.008 PV/yr in the bottom of the domain at (22, 0)m, slightly displaced from the vertical boundary such that the flow is not dominated by flow along the impermeable boundary. Similarly, production is at (457, 100)m. The initial pressure, and temperature (250 bar and 403.15 K) are such that most of the domain will be in two-phase after injection. For comparison, we repeat the computation for an isotropic domain with scalar absolute permeability of 10 md. The results are illustrated in Fig. 11 at 100% PVI (t = 125 yr). Due to the density stratification in the gravitational field, the flow is highly buoyant and mostly vertical before spreading out when it hits the impermeable top cap of the domain. The ‘drooping’ of oil–CO2 mixture, referred to as fingering, can be seen to occur at the gas–liquid interface in the isotropic medium (Fig. 10). To verify that these features are not artifacts of the irregularities in the unstructured mesh, and to study the numerical diffusion at the phase boundary, we compare in Fig. 12 to a simulation performed on a finer 200 × 40 rectangular grid and using bilinear DG. We find that the results are in good agreement with only slightly more numerical diffusion at the fronts on the coarser unstructured mesh. The gravitational fingers are an instability that are triggered by numerical errors (representing small perturbations of a meta-stable state in reality) and are not expected to start at the exact same location from one simulation to the next. The fingering will be discussed in more detail in the next example. When the rock matrix has an anisotropic permeability and is sheared in the horizontal direction, the CO2 permeates more of the

domain by spreading horizontally (Fig. 10). If the CO2 were denser than the oil, which applies to certain reservoirs, this would lead to a significant enhancement in oil recovery. In this particular example, the light gas still rises quickly to the top of the domain, where the high horizontal permeability allows fast spreading and very early breakthrough. The resulting low recovery, compared to the simulations for an isotropic domain on unstructured and structured meshes is shown in Fig. 13. 4.5. Example 5: CO2 injection in vertical isotropic domain from the top In our final example, we consider the most interesting problem of CO2 injection from the top and production at the bottom of the domain, where all other parameters are the same as in the previous example (Section 4.4). Due to the density stratification, the lighter CO2 first spreads horizontally before gravitationally draining towards the production well in the bottom left. More interestingly, CO2 has the unique property that when it dissolves in oil, the resulting mixture is denser than the original oil and will consequently slowly sink towards the bottom. At the same time, the opposite effect occurs: the lighter components in the oil evaporate at the gas–oil interface and can locally decrease the gas density (CO2 is relatively heavy). This results in ‘plumes’ of lighter gas that buoyantly rise to the top of the domain. These non-linear fingering processes have received little attention in the literature. In this example, we choose the low injection rate of 0.008 PV/yr. The overall CO2 composition is shown at 25% (top) and 100% (middle)

Fig. 13. Example 4a–b: oil recovery (liquid volume fraction) through bottom injection into a sheared domain with 250 × larger permeability in the horizontal direction, as compared to an isotropic domain.

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

PVI in Fig. 14. The middle figure clearly demonstrates that we can resolve the fingering process, and we can see both downward drooping of denser oil, and upward rising of lighter gas plumes. However, in these simulations we neglected Fickian diffusion. When we include diffusion, the components in the heavier mixture can diffuse to their surroundings and vice versa, and the instability is dampened. The bottom panel in Fig. 14 shows this result at 100% PVI, where the fingering is completely suppressed. Fingering also does not occur when we increase the injection rate by a factor three, which is enough for the convection to dominate the flow. Again, we compare to bilinear DG results on a finer 200 × 40 element rectangular mesh in Fig. 15. For this case, we do not see the downward fingers develop, which may be due to the higher accuracy of the simulation, in which the phase boundary may not be perturbed enough to trigger the onset of the instability (this may still happen at a later time). In a simulation that neglects Fickian diffusion, we do see the occurrence of gas plumes in the overall CO2 composition (top panel), and more clearly from the C2,g − C3,g composition in the gas phase (middle panel). When we do account for Fickian diffusion, the result is the same as on the unstructured grid and the instability is dampened. The effect of diffusion on oil recovery, however, is negligible in this case (Fig. 16). Diffusion is considered to be mainly important in heterogeneous media, such as sheared formations or fractured rock matrices [47], but as we have illustrated in these examples, diffusion may play an important role even in homogeneous media. We have seen that gravitational fingering may not occur in real systems when the effect of diffusion is accounted for as in our formulation. In all the literature on fingering that we are aware of, this diffusion effect is ignored.

963

5. Summary and conclusions We have presented several refinements and additions to higher-order numerical algorithms for compositional modeling. To summarize some of the features of our current model: we efficiently and robustly simulate compressible two-phase flow in both homogeneous and heterogeneous porous media, on structured or unstructured grids with full compositional modeling. Fickian diffusion is properly included by using a self-consistent model that calculates the diffusion coefficients from irreversible thermodynamics, and by including the cross-diffusion terms. Both momentum transport, in the form of Darcy's law, and mass flow are solved using higher-order methods to improve the reliability of our results. In particular, we analyze the pressure and velocity fields using the hybridized mixed finite element method, which has been shown to be superior to the finite difference and finite volume methods in heterogeneous media and has desirable features as compared to the classical MFE method. Mass transport is more accurately modeled by the full discontinuous Galerkin method, both in the single phase regions and in two-phase regions. In the latter, several stability and flash calculations are performed per element to improve the shock capturing and reduce the numerical dispersion, as demonstrated in Section 4. Furthermore, we have extended our model to allow for tensorial permeability to study anisotropic and heterogeneous media. Finally, we have implemented unstructured grids in our model. This allows careful analysis of sensitivity to mesh quality and orientation and, more importantly, is essential in the treatment of arbitrary discrete fractures, which we will present in a forthcoming paper. In Section 4 we have first presented simple examples to illustrate key features in isolated settings before proceeding to the complex problems

Fig. 14. Example 5: isotropic domain with permeability: Kxx = Kyy = 10 md, Kxy = 0. Injection is from the top right and production at the bottom left. (A) 25% PVI and (B) 100% PVI without diffusion. (C) with diffusion at 100% PVI.

964

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

Fig. 15. Example 5: structured 200 × 40 element mesh at 100% PVI. Top: Overall CO2 composition, and C2 − C3 composition in the gas phase (middle) without diffusion. Bottom: CO2 composition with diffusion at 100% PVI.

of CO2 injection into an anisotropic domain saturated with a petroleum fluid. We can resolve the onset of the fingering instability, due to the phase behavior properties of CO2 mixtures, when a layer of carbondioxide is sitting on top of the oil for a long enough time-scale. However, this only occurs when diffusion is neglected. Fickian diffusion allows the components in the mixture to diffuse from regions of higher densities to regions of lower density and vice versa, equilibrating the system and dampening the instability.

Appendix A. Mixed hybrid finite element approximation To perform the integrals in Eq. (20) on triangular elements, it is useful to map and evaluate the expressions onto the unit triangle. Using the numbering in Fig. 1, we map node 2 to the origin by the translation (− x2, − y2). Furthermore, we define in terms of the nodecoordinates in the triangle element: A = x1 −x2 ; B = x3 −x2 C = y1 −y2 ;

D = y3 −y2

such that the affine mapping from the unit triangle to the mesh triangle, and its inverse can be expressed as:  T=

 B ; D

A C

T

−1

=

1 ∥T∥



D −C

 −B ; A

ð32Þ

denoting the determinant of T as ∥ T ∥ = AD − BC. The inverse of the permeability tensor in the mesh element K̃ and the mapping onto the unit triangle K are: −1

K K˜ =

−1

1 2 Kxx Kyy −Kxy

 −Kxy ; Kxx

D Kyy + ðB + CÞDKxy + BCKxx ∥K∥ ∥T∥2


Fig. 16. Example 5: oil recovery with and without Fickian diffusion.

Kyy −K xy

ð33Þ

2

KK;xx =

−1 KK;xy



=−

;

 B2 + AD Kxy + BðDKxx + AKxx Þ ∥K∥ ∥T∥2

ð34Þ

;

ð35Þ

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968


−1

KK;xy = −

−1

KK;yy =

 C 2 + AD Kxy + C ðDKxx + AKxx Þ ∥K∥ ∥T∥

2

;

ð36Þ

2

BCKyy + ðB + C ÞAKxy + A Kxx ∥K∥ ∥T∥2

:

ð37Þ

Taking the product with the xK,j, summing over the fluxes and multiplying with the inverse of the mass matrix, we get (omitting the K subscripts): 0 1 2x1 ðq1 + q2 −3q3 Þ + x2 ðq1 + q2 −3q3 Þ + x3 ðq1 + q2 −q3 Þ term 2 = @ x1 ðq1 −3q2 + q3 Þ + 2x2 ðq1 −3q2 + q3 Þ + x3 ðq1 −q2 + q3 Þ A x1 ðq2 −3q1 + q3 Þ + x2 ðq2 −3q1 + q3 Þ + 2x3 ð−q1 + q2 −3q3 Þ

ð42Þ

Next, the integral can be evaluated on the unit triangle, using the basis vector fields wE, as defined in Section 3 by: BK;E;E′ =

1 1−x −1 ∫0 ∫0 wE ⋅KK ⋅wE′ dxdy:

ð38Þ

The result is then simply mapped back onto the actual mesh element by = jjT jjBK;E;E′ . The second integral involving the heterogeneous B K;E;E′ ˜ permeability, in Eq. (21b), can be evaluated explicitly over each edge E. In particular, we assume that K is element-wise constant, and average the scalar magnitude of the permeability over the two neighboring elements at E to ensure continuity in the pressure equation. The coefficients in Eq. (22) are given, as in [35], by: −1

βK;E;E′ = λK;t BK;E;E′ ; αK;E = ∑ βK;E;E′ ;

ð39aÞ

h i γK;E = ∫E Kg⋅nK;E ∑ λK;α ρK;α :

ð39bÞ

E′ ∈∂K

The integrals in the third term (Eq. (29c)) are evaluated as: 0 2 1 @ 1 ∫E1 φi φj = 24 0

1 2 0

1 0 0A 0

ð43aÞ

0 2 1 @ ∫E2 φi φj = 0 24 1

0 0 0

1 1 0A 2

ð43bÞ

0 0 1 @ ∫E3 φi φj = 0 24 0

0 2 1

1 0 1A 2

ð43cÞ

Summing over the edges, taking the product with the upwind nodal composition and multiplying with M− 1, we find: 0

α

1 ð5 x˜11 + x˜12 Þq1 + ð5 x˜21 + x˜23 Þq2 −3ð x˜32 + x˜33 Þq3 @ term 3 = ð x˜11 + 5 x˜12 Þq1 −3ð x˜21 + x˜23 Þq2 + ð5 x˜32 + x˜33 Þq3 A −3ð x˜11 + x˜12 Þq1 + ð x˜21 + 5 x˜23 Þq2 + ð5 x˜32 + x˜33 Þq3

ð44Þ

and those in Eq. (24), in terms of ΦK = ϕKcK,f A / Δt by: nc

where x̃E,j is the upwind value of x at node j of edge E. Finally, term 4 is simply the sum of the averaged diffusive fluxes through the three edges of the element: ∑ E ∈ ∂K Jdiff, K, E. The updating scheme for the saturations at the new time-step n + 1 can now be written schematically as:

R˜ K;E′ = ∑ vK;i ∑ ∫E βK;E;E′ ∑ cα xi;α fα ; α

E2∂K

i=1

DK = ΦK + ∑ R˜K;E′ ; n

E′ 2∂K nc

GK = ΦK pK − ∑ vK;i ∑ ∑α ∫E cα xi;α λα ρα Kg⋅n:

n + 1

zK;j

E2∂K

i=1

Appendix B. Linear discontinuous Galerkin approximation on triangular grids We evaluate explicitly the integrals in the DG treatment of the transport equation Eq. (29a). First, the mass matrix in Eq. (29a) is: 0

M = ∫ K φi φj =

2 1 @ 1 24 1

1 2 1

1 1 A: 2

ð40Þ

0 1 2 1 −3 1 @ ∫K φi wE1 ·∇φj = 1 2 −3 A 24 1 1 −2 0

2 −3 1 @ 1 −2 24 1 −3

0 −2 1 @ ∫K φi wE3 ·∇φj = −3 24 −3

1 2 1

n

= zK; j +

 n 2Δt ∑ ðterm 2 + term 3 + term 4Þ + fK ϕK A α

ð45Þ

where fK is the element-wise constant sink/source and diffusion term and the factor 2/A is the inverse mapping of the unit- to the mesh triangle, with A the mesh element surface area. The forward Euler method is used for the derivative, and we indicate that Eq. (42), Eq. (44) and the diffusion fluxes should be summed over the two phases. Appendix C. Bilinear discontinuous Galerkin approximation on rectangular grids

1

Next we evaluate the integrals in Eq. (29b):

∫K φi wE2 ·∇φj =

965

ð41aÞ

1

1 1A 2 1 1 1A 2

ð41bÞ

ð41cÞ

In this Appendix, we briefly outline the generalization of the integrals for the linear DG approximation in the previous Appendix to bilinear DG on rectangular grids. We choose the scalar shape functions φj and lowest-order Raviart–Thomas basis vector functions on an element with surface K = lx × ly, as: φ1 ðxÞ =

   ly 1 lx −x −y ; lx ly 2 2

ð46aÞ

φ2 ðxÞ =

   ly 1 lx +x −y ; lx ly 2 2

ð46bÞ

φ3 ðxÞ =

   ly 1 lx +x +y ; lx ly 2 2

ð46cÞ

φ4 ðxÞ =

   ly 1 lx −x +y ; lx ly 2 2

ð46dÞ

966

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

0

and

1 −2 −1 1 2 B 1 −1 −2 2 1 C C; ∫T φj wT ·∇φk = B 6 @ −2 −4 4 2 A −4 −2 2 4

ð50cÞ

0 1 l 1 @x + x A wK;R ðxÞ = 2 ; lx ly 0

ð47aÞ

0 1 l 1 @ x− x A wK;L ðxÞ = 2 ; lx ly 0

ð47bÞ

1 4 2 −2 −4 1 B 2 4 −4 −2 C C: ∫B φj wB ·∇φk = B 6 @ 1 2 −2 −1 A 2 1 −1 −2

ð47cÞ

The terms are summed over the two phases and four edges and we approximate the time derivative by the forward Euler method, as in Eq. (45).

0 wK;T ðxÞ =

1

0

1 @ l A; lx ly y + y 2 0

wK;B ðxÞ =

0

Appendix D. Algorithm

1

0

1 @ l A; lx ly y− y 2

ð47dÞ

in terms of the edges E = R(ight), L(eft), T(op), B(ottom), which are numbered as E = (1, 2, 3, 4). The relevant integral matrices in Eq. (29a) are given next. The inverse of the mass matrix is 0 4 −2
−1 B −2 4 B ∫ K φj φk = 4@ 1 −2 −2 1

1 1 −2 −2 1C C: 4 −2 A −2 4

ð48Þ

Note that the inner product of Eq. (48) with the identity vector I = (1,1,1,1) is 4I, and that ∫Kϕk = 1/4 ∀k = 1, …, 4 such that the right hand side of Eq. (29a) again reduces to fKI when the sink/source term f (x, t) is taken element-wise-constant fK. The bilinear form of the contour integrals in Eq. (29a) are given by 0

0 1B 0 B ∫R φj φk = @ 6 0 0

0 2 1 0

0 1 2 0

1 0 0C C; 0A 0

ð49aÞ

0 2 1B 0 B ∫L φj φk = @ 6 0 1

0 0 0 0

0 0 0 0

1 1 0C C; 0A 2

ð49bÞ

0 0 1B 0 B ∫T φj φk = @ 6 0 0

0 0 0 0

0 0 2 1

1 0 0C C; 1A 2

ð49cÞ

0 2 1B 1 B ∫B φj φk = @ 6 0 0

1 2 0 0

0 0 0 0

1 0 0C C: 0A 0

ð49dÞ

Finally, the volume integrals are given by: 0

−2 2 1 B −4 4 ∫R φj wR ·∇φk = B 6 @ 1 −1 −1 1 0

4 1B2 ∫L φj wL ·∇φk = B 6@1 2

ð50dÞ

1 1 −1 2 −2 C C; −2 2A 2 −2

1 −4 −2 2 −2 −1 1 C C; −1 −2 2 A −2 −4 4

ð50aÞ

ð50bÞ

For clarity we summarize the structure of our numerical algorithm. First the problem is initialized by the following steps: • Discretize the domain, either as structured rectangular grids or an unstructured or structured triangulation. In the latter case, optimize the numbering of nodes and edges for the linear solver. • Post-process the geometry: for each edge and node find all the neighboring elements, calculate the element volumes and outward pointing normals to the element boundaries. • Read the properties of the multi-component fluid from an input file: binary interaction parameters, critical temperature, pressure and volume, molar weight, acentric factor, volume shift, viscosity, molar volume at the boiling point, critical density, normal boiling point temperature and finally, initial composition. These data will be used in the calculation of both phase behavior and the diffusion coefficients. • Read the initial and boundary conditions for the overall composition, temperature (generally assumed isothermal), pressure, injection rate or production rate, porosity, permeability, number and location of injectors and producers, residual oil and gas saturations and either tabulated or formula parameters for the relative permeability's dependence on saturation. • Perform a stability analysis to determine which regions in the domain are in single and which in two-phase and do a flash calculation to find the initial molar phase densities cα, compositions xi,α and ϖ, as well as the derived quantities such as saturation Sα, etc. • Compute the total partial molar volumes and total compressibility for the pressure equation. • Update the injection wells. • Calculate the mass matrix as in Eq. (38) (note that this has to be done only once). After the initialization we start the adaptive time evolution by repeating the following steps until either a maximum time or injected pore volume is reached. Generally, we end the simulation when 100% of the pore volume is injected. 1. From the temperature, pressure and phase compositions, calculate the Fickian diffusion coefficients Dij at user-specified timeintervals (Dij is only weakly time dependent). 2. Construct and solve the linear system Eq. (25) for the pressure traces tpK,E on each element edge and from Eq. (23) find the element averaged pressures pK. 3. Compute qK,E, the normal components of the total convection flux through edges E, from Eq. (22). 4. From the total flux compute the phase fluxes ϑα (or rather qα,K,E) using Eq. (10) with the saturation and viscosity from the previous time-step. Determine the upwind direction for each mesh element. 5. Evaluate the diffusion flux as in Eq. (3b)–Eq. (3c) using the diffusion coefficients calculated in the first step, and the molar density and composition gradients from the previous time-step. 6. With the computed flux, check the Courant–Friedrichs–Lewy (CFL) condition on the time-step [19] to keep the system stable

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

7.

8.

9. 10. 11.

12. 13. 14.

and obey the maximum principle in the diffusion equation [39]. If the CFL condition is violated, decrease the time-step and return to step 2 with the original data from the previous time-step. Update explicitly (using the forward Euler scheme) the species balance equation for the overall composition zi as in Eq. (29a) or Eq. (45), with the new phase convective and diffusive (as a source) fluxes and the cα xi,α from the previous time-step. Upwinding is done for cα xi,α in the convection integral. Perform a minimal reconstruction of the nodal overall composition such that zi does not take on extremal values with respect to the surrounding elements, to avoid spurious oscillations and instability of the scheme. Adapt the time-step based on the CFL condition, mass balance and maximum temporal gradients in pressure. Perform the stability analysis and flash calculations to find new phase molar densities and compositions. Update all the derived quantities, such as mole fraction, saturation, viscosity, compressibility and partial molar volume. The viscosities are calculated from the composition using the method in [45]. Update the source term. Check mass balance. At specified intervals, expressed in terms of the fraction of the pore volume injected, write output data.

References [1] Acs G, Doleschall S, Farkas E. General purpose compositional model. SPE J 1985;25 (4):543–53. [2] Arnold DN, Brezzi F, Cockburn B, Marini LD. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J Numer Anal 2002;39:1749–79. [3] Ashcroft S, Ben Isa M. Effect of dissolved gases on the densities of hydrocarbons. J Chem Eng Data 1997;42(6):1244–8. [4] Babuska I, Zlámal M. Nonconforming elements in the finite element method with penalty. SIAM J Numer Anal 1973;10:863–75. [5] Bassi F, Rebay S. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations. J Comput Phys 1997;131:267–79. [6] Bassi F, Rebay S, Mariotti G, Pedinotti S, Savini M. A high-order accurate discontinuous finite element method for inviscid and viscous turbomachinery flows. Fluid Dynamics and Thermodynamics. Antwerpen, Belgium: Technologisch Instituut; 1997. p. 99–108. [7] Baumann CE, Oden JT. A discontinuous hp finite element method for convection– diffusion problems. Comput Methods Appl Mech Engrg 1999;175:311–41. [8] Brezzi F, Fortin M. Mixed and hybrid finite element methods. New York: SpringerVerlag; 1991. [9] Brezzi F, Manzini G, Marini D, Pietra P, Russo A. Discontinuous finite elements for diffusion problems. Atti Convegno in onore di F. Brioschi (Milan, 1997). Milan, Italy: Istituto Lombardo, Accademia di Scienze e Lettere; 1999. p. 197–217. [10] Brezzi F, Manzini G, Marini D, Pietra P, Russo A. Discontinuous Galerkin approximations for elliptic problems. Numer Meth Partial Differ Equ 2000;16:365–78. [11] Chavent G, Cohen G, Jaffré J, Eymard R, Guérillot DR, Weill L. Discontinuous and mixed finite elements for two-phase incompressible flow. SPE Reservoir Eng 1990;5(4):567–75. [12] Chavent G, Roberts JE. A unified physical presentation of mixed, mixed-hybrid finiteelements and standard finite-difference approximations for the determination of velocities in water-flow problems. Adv Water Resour 1991;14(6):329–48. [13] Cockburn B, Dawson C. Some extensions of the local discontinuous Galerkin method for convection–diffusion equations in multidimensions. In: Whiteman JR, editor. Proceedings of the Conference on the Mathematics of Finite Elements and Applications: MAFELAP X. New York: Elsevier; 2000. p. 225–38. [14] Cockburn B, Dong B, Guzman J. Optimal convergence of the original DG method for the transport-reaction equation on special meshes. SIAM J Numer Anal 2008;46(3): 1250–65. [15] Cockburn B, Dong B, Guzman J, Restelli M, Sacco R. A hybridizable discontinuous Galerkin method for steady-state convection–diffusion–reaction problems. SIAM J Sci Comput 2009;31(5):3827–46. [16] Cockburn B, Gopalakrishnan J, Lazarov R. Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems. SIAM J Numer Anal 2009;47(2):1319–65. [17] Cockburn B, Guzman J, Soon S-C, Stolarski HK. An analysis of the embedded discontinuous Galerkin method for second-order elliptic problems. SIAM J Numer Anal 2009;47(4):2686–707. [18] Cockburn B, Shu C-W. The local discontinuous Galerkin method for timedependent convection–diffusion systems. SIAM J Numer Anal 1998;35:2440–63. [19] Courant R, Friedrichs K, Lewy H. Partial differential equations of mathematical physics. Math Ann 1928;100:32–74.

967

[20] Davis TA. Algorithm 832: Umfpack v4.3—an unsymmetric-pattern multifrontal method. ACM Trans Math Softw 2004;30(2):196–9. [21] Davis TA. A column pre-ordering strategy for the unsymmetric-pattern multifrontal method. ACM Trans Math Softw 2004;30(2):165–95. [22] Davis TA, Duff IS. An unsymmetric-pattern multifrontal method for sparse LU factorization. SIAM J Matrix Anal Appl 1997;18(1):140–58. [23] Davis TA, Duff IS. A combined unifrontal/multifrontal method for unsymmetric sparse matrices. ACM Trans Math Softw 1999;25(1):1–20. [24] Duncan JB, Toor HL. An experimental study of three component gas diffusion. AIChE J 1962;8(1):38–41. [25] Epshteyn Y, Rivière B. On the solution of incompressible two-phase flow by a p-version discontinuous Galerkin method. Commun Numer Meth Eng 2006;22(7):741–51. [26] Epshteyn Y, Rivière B. Fully implicit discontinuous finite element methods for two-phase flow. Appl Numer Math 2007;57(4):383–401. [27] Epshteyn Y, Rivière B. Analysis of hp discontinuous Galerkin methods for incompressible two-phase flow. J Comput Appl Math 2009;225(2):487–509. [28] Farthing MW, Kees CE, Miller CT. Mixed finite element methods and higher-order temporal approximations. Adv Water Resour 2002;25(1):85–101. [29] Firoozabadi A. Thermodynamics of hydrocarbon reservoirs. USA: McGraw-Hill Professional; 1999. [30] Firoozabadi A, Katz D, Soroosh H, Sajjadian VA. Surface tension of reservoir crudeoil/gas systems recognizing the asphalt in the heavy fraction. SPE Reservoir Eng J 1988;3(1):265–72. [31] Gould N, Hu Y, Scott J. A numerical evaluation of sparse direct solvers for the solution of large sparse, symmetric linear systems of equations. Technical report RAL-TR-2005-005, Council for the Central Laboratory of the Research Councils; 2005. p. 1–32. ftp://ftp.numerical.rl.ac.uk/pub/reports/ghsRAL200505.pdf. [32] Hoteit H, Ackerer P, Mose R, Erhel J, Philippe B. New two-dimensional slope limiters for discontinuous Galerkin methods on arbitrary meshes. Int J Numer Meth Eng 2004;61(14):2566–93. [33] Hoteit H, Erhel J, Mose R, Philippe B, Ackerer P. Numerical reliability for mixed methods applied to flow problems in porous media. SO Comput Geosc 2002;6(2):161–94. [34] Hoteit H, Firoozabadi A. Multicomponent fluid flow by discontinuous Galerkin and mixed methods in unfractured and fractured media. Water Resour Res 2005;41(11): W11412. [35] Hoteit H, Firoozabadi A. Compositional modeling by the combined discontinuous Galerkin and mixed methods. SPE J 2006;11(1):19–34. [36] Hoteit H, Firoozabadi A. Compositional modeling of discrete-fractured media without transfer functions by the discontinuous Galerkin and mixed methods. SPE J 2006;11(3):341–52. [37] Hoteit H, Firoozabadi A. Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures. Adv Water Resour 2008;31(1): 56–73. [38] Hoteit H, Firoozabadi A. Numerical modeling of diffusion in fractured media for gas-injection and -recycling schemes. SPE J 2009;14(2):323–37. [39] Hoteit H, Mose R, Philippe B, Ackerer P, Erhel J. The maximum principle violations of the mixed-hybrid finite-element method applied to diffusion equations. Int J Numer Methods Eng 2002;55(12):1373–90. [40] Douglas J, Dupont T. Lecture Notes in Phys. Ch. interior penalty procedures for elliptic and parabolic Galerkin methodsBerlin: Springer-Verlag; 1976. [41] Kees CE, Miller CT. Higher order time integration methods for two-phase flow. Adv Water Resour 2002;25(2):159–77. [42] Kirby RM, Karniadakis GE. Selecting the numerical flux in discontinuous Galerkin methods for diffusion problems. J ScientiÞc Comput 2005;22(1):385–411. [43] Kooijman H, Taylor R. Estimation of diffusion coefficients in multicomponent liquid systems. Ind Eng Chem Res 1991;30:1217–22. [44] Leahy-Dios A, Firoozabadi A. Unified model for nonideal multicomponent molecular diflusion coefficients. AlChE J 2007;53(11):2932–9. [45] Lohrenz J, Bray BG, Clark CR. Calculating viscosities of reservoir fluids from their compositions. J of Pet Technol 1964;16(10):1171–6. [46] Mikysta J, Firoozabadi A. Implementation of higher-order methods for robust and efficient compositional simulation. J Comput Phys 2010;229:2898–913. [47] Moortgat J, Firoozabadi A, Moravvej Farshi M. A new approach to compositional modeling of CO2 injection in fractured media compared to experimental data. SPE 2009;124918-MS:1–11. [48] Mose R, Siegel P, Ackerer P, Chavent G. Application of the mixed hybrid finiteelement approximation in a groundwater-flow model – luxury or necessity. Water Resour Res 1994;30(11):3001–12. [49] Nguyen NC, Peraire J, Cockburn B. An implicit high-order hybridizable discontinuous Galerkin method for linear convection–diffusion equations. J Comput Phys MAY 20 2009;228(9):3232–54. [50] Nguyen NC, Peraire J, Cockburn B. An implicit high-order hybridizable discontinuous Galerkin method for nonlinear convection–diffusion equations. J Comput Phys DEC 10 2009;228(23):8841–55. [51] Pardiso. PARDISO Solver Project. http://www.pardiso-project.org2009. [52] Peng D-Y, Robinson DB. A new two-constant equation of state. Ind Eng Chem Fundam 1976;15(1):59–64. [53] Raviart PA, Thomas JM. Primal hybrid finite-element methods for second-order elliptic equations. Math Comput 1977;31(138):391–413. [54] Rivière B, Wheeler M, Girault V. Improved energy estimates for interior penalty, constrained and discontinuous Galerkin methods for elliptic problems I. SIAM J Numer Anal 1999;3:337–60. [55] Sammon PH. An analysis of upstream differencing. SPE Reservoir Eng 1988;3(3): 1053–6. [56] Schenk O, Bollhoefer M, Roemer R. On large-scale diagonalization techniques for the Anderson model of localization. SIAM Rev 2008;50:91–112.

968

J. Moortgat, A. Firoozabadi / Advances in Water Resources 33 (2010) 951–968

[57] Schenk O, Waechter A, Hagemann M. Matching-based preprocessing algorithms to the solution of saddle-point problems in large-scale nonconvex interior-point optimization. J Comput Optim Appl 2007;36(2–3):321–41. [58] Shu CW. Different formulations of the discontinuous Galerkin method for the viscous terms. In: Shi Z-C, Mu M, Xue W, Zou J, editors. Advances in scientic computing. Mascou: Science Press; 2001. p. 144–55.

[59] Simon R, Rosman A, Zana E. Phase-behavior properties of CO2-reservoir oil systems. SPE J 1978;18(1):20–6. [60] Watts JW. A compositional formulation of the pressure and saturation equations. SPE Reservoir Eng 1986;1(3):243–52.