The Influence of Velocity Field Approximationsin Tracer Injection Processes

ABSTRACT

Although the concentration is the most important variable in tracer injection processes, an efficient and accurate velocity field approximation is crucial to obtain a good physical behaviour for the problem. In this paper we analyse a Stabilized Dual Hybrid Mixed (SDHM) method to solve the Darcy’s system in the velocity and pressure variables that involves the conservation of mass and Darcy’s law. This approach is locally conservative, free of compromise between the finite element approximation spaces and capable of dealing with heterogeneous media with discontinuous properties. The tracer concentration is solved via a combination of the Streamline Upwind Petrov-Galerkin (SUPG) method in space with an implicit finite difference scheme in time. We also employ a semi-analytical approach (Abbaszadeh-Dehghani analytical solution) to integrate the transport equation. A numerical comparative study using the SDHM formulation, the Galerkin method and a post-processing technique to calculate the velocity field in combination with those concentration approximation methodologies are presented. In all comparisons, the SDHM formulation appears as the most efficient, accurate and almost free of spurious oscillations.

Keywords:
Miscible displacements; Hybridized method; Oil reservoir simulations

RESUMO

Embora a concentração seja a variável mais importante nos processos de injeção de traçadores, uma eficiente e precisa aproximação do campo de velocidades é crucial para obter um bom comportamento físico para o problema. Neste artigo, analisamos o método misto dual híbrido estabilizado (SDHM) para resolver o sistema de Darcy nas variáveis de velocidade e de pressão a partir da equação de conservação de massa e da lei de Darcy. Esta abordagem é localmente conservativa, livre de comprometimento entre os espaços de aproximação de elementos finitos e capaz de lidar com meios heterogêneos com propriedades descontínuas. A concentração do traçador é resolvida através de uma combinação do método Streamline Upwind Petrov-Galeklin (SUPG) no espaço com um método de diferenças finitas implícita no tempo. Também empregamos uma abordagem semi-analítica (solução analítica de Abbaszadeh-Dehghani) para integrar a equação de transporte. Um estudo comparativo numérico utilizando a formulação de SDHM, o método de Galerkin e uma técnica de pós-processamento para calcular o campo de velocidade em combinação com essas metodologias de aproximação da concentração são apresentados. Em todas as comparações, a formulação SDHM aparece como a mais eficiente, precisa e quase sem oscilações espúrias.

Palavras-chave:
Deslocamentos miscíveis; Métodos hibridizados; Simulações de reservatórios de óleos

1 INTRODUCTION

Information extracted from tracer breakthrough profiles at production wells plays an important role in reservoir engineering both in the characterization of reservoir heterogeneities as in the project of recovery techniques. Those profiles can be detected either experimentally or via the solution of a mathematical model, which describes the transport of substances through a porous medium. Although the concentration is the variable of most interest, approximation of the velocity field is crucial, since it is responsible for the flow displacement.

Standard Mixed Finite Element (MFE) methods have been extensively employed during the last decades to solve Darcy’s system^{27), (7), (8), (15), (13}. The main idea of these methods is simultaneously approximate pressure and velocity by using different spaces for each variable. This leads to a compatibility condition between the approximation spaces (the LBB condition⁶), and thus restricts the choice of stable finite element spaces. However, these choices are usually unstable with standard dual mixed formulation, as was illustrated in¹¹. A very common example of stable mixed method is given by the Raviart-Thomas spaces²⁷. To overcome the compatibility conditions, some stabilized mixed finite element methods were proposed in^{16), (19), (22), (4), (10), (11), (12), (5}. In general, these stabilized formulations use continuous Lagrangian finite element spaces and can be successfully employed in simulating Darcy flows in homogeneous porous media. However, continuous interpolations are not appropriate to heterogeneous porous media with discontinuous properties. This is due to the fact that on the interfaces of the discontinuities the normal component of Darcy velocity must be continuous (mass conservation) but the tangential component is discontinuous. Formulations based on continuous Lagrangian interpolation for velocity fail to represent the tangential discontinuity, producing inaccurate approximations and spurious oscillations. Therefore, based on hybridization techniques the Stabilized Dual Hybrid Mixed (SDHM) method that combines advantages of the Discontinuous Galerkin (DG) methods²⁸ with reduced computational was developed for calculating accurate velocity fields to miscible displacements in homogeneous and heterogeneous porous media^{23), (24), (25}. The SDHM formulation consists of a set of local problems defined at the element level coupled to a global system for the Lagrange multipliers. Stabilization terms are added to generate a stable and adjoint consistent formulation allowing greater flexibility in the choice of the approximation spaces. The Lagrange multiplier is identified as the pressure trace on the element interface, which is a natural choice.

In²⁴ the authors have shown that for regular solutions the SDHM method leads to optimum rates of convergence for the velocity and pressure fields, even when same order interpolations are employed. In addition, the mixture transport in homogeneous and heterogeneous media were recently analysed²⁵ by studying the influence of the mobility ratio and the permeability field variation. In these cases, the SDHM method was applied in the approximation of the velocity field coupled to the SUPG method in the calculation of the concentration field. The good performance of this methodology was verified through the presented numerical simulations.

In this paper, following the works mentioned in the above paragraph, we study the efficiency and robustness of the SDHM comparing with the usual Galerkin method and a post-processing technique^{21), (20}. To do this, we analyse the influence of the velocity field approximations on the behaviour of the prediction of tracer concentration at the producer well. The concentration approximation is obtained via the SUPG method combined with an implicit finite difference scheme, as in^{24), (25}. In order to validate the numerical results, a semi-analytical methodology, which consists of combining the analytical solution for the concentration given by Abbaszadeh-Dehghani¹ with an approximation for the velocity field is also provided.

The outline of this work is as follows. Section 2 considers the mathematical model of the problem. The stabilized dual hybrid mixed (SDHM) method for the Darcy’s system is presented in Section 3. Tracer concentration approximations based on the SUPG approach and on the semi-analytical methodology are discussed in Section 4, as well as its combination with the usual Galerkin method and a post-processing technique. Numerical experiments are reported in Section 5. Finally, some conclusions are given in Section 6.

2 THE MODEL PROBLEM

The equations governing tracer flows in miscible displacements²⁶ can be described by the first order system of partial differential equations in the Darcy velocity u=u(x,t) and the pressure p=p(x) coming from the mass conservation of the mixture and the Darcy’s law:

$divu=fin\Omega \times \left(0,T\right),$ (2.1)

$u=-\mathbb{K}\nabla pin\Omega \times \left(0,T\right),$ (2.2)

combined with a convection dominated diffusion-convection equation, which expresses the conservation of mass of the injected fluid (the tracer concentration equation), c=c(x,t), gives by

$\varphi \frac{\partial c}{\partial t}+u·\nabla c-div\left(\mathbb{D}\nabla c\right)=gin\Omega \left(0,T\right),$ (2.3)

with the following initial and boundary conditions

$c\left(x,0\right)=c_0\left(x\right)in\Omega ,$ (2.4)

$\mathbb{D}\nabla c·n=0on\partial \Omega \times \left(0,T\right),$ (2.5)

$u·n=0on\partial \Omega \times \left(0,T\right).$ (2.6)

In (1-6) K=K(x) is the permeability tensor and Φ is the porosity. Functions f and g are the source terms, n is the exterior normal to ∂Ω and

$\mathbb{D}=\mathbb{D}\left(u\right)=\left({\alpha }_{mol}+{\alpha }_t\left|u\right|\right)\mathrm{𝕀}+\frac{{\alpha }_l-{\alpha }_t}{\left|u\right|}u\otimes u$ (2.7)

is the dispersion-diffusion tensor²⁶ where α_mol , α_l and α_t are the molecular diffusion, longitudinal and transversal dispersion coefficients, respectively, with the domain Ω ⊂ R ⁿ , ∂Ω its boundary and T ∈ R,T>0.

3 THE HYBRID METHOD FOR DARCY FLOW PROBLEM

Let T _h = {K} be a regular finite element mesh on the domain $\Omega \left(\Omega =\underset{\mathcal{K}}{\cup }\mathcal{K}\right)$. The set of all edges of all elements K is E _h = {e: e is an edge of K for all K ∈ T _h } with ${\mathcal{E}}_h^0$ denoting the set of the interior edges. To introduce the hybrid formulation we first multiply equations (1)-(2) by their respective weighting functions and integrate by parts on each element K, getting the following local weak form

${\left({\mathbb{K}}^{-1}u,v\right)}_{\mathcal{K}}-{\left(p,divv\right)}_{\mathcal{K}}+{\int }_{\partial \mathcal{K}}p\left(v·n_{\mathcal{K}}\right)ds=0,\forall v\in {\mathcal{U}}_{\mathcal{K}},$ (3.1)

$-{\left(divu,q\right)}_{\mathcal{K}}+{\left(f,q\right)}_{\mathcal{K}}=0,\forall q\in {\mathcal{L}}_{\mathcal{K}}.$ (3.2)

U_k = {v ∈ [L ²(K)]², div v ∈ L ²(K), ∀K ∈ T _h } and L _K = {q ∈ L ²(K), ∀K ∈ T _h }, are the local function spaces on each element K. For each element K ∈ T _h , (w, v)_K = ∫_K wvd _K denotes the usual L ²(K), inner product, ∂K is set of all edges of element K and n _K is the external normal to ∂K.

For a given $p=\overline{p}$ on ∂K, we can solve the set of local problems:

For each K ∈ T _h , find {u, p} ∈ U _K × L _K , ∀{v, q} ∈ U _K × L _K such that

${\left({\mathbb{K}}^{-1}u,v\right)}_{\mathcal{K}}-{\left(p,divv\right)}_{\mathcal{K}}-{\left(divu,q\right)}_{\mathcal{K}}=-{\left(f,q\right)}_{\mathcal{K}}-{\int }_{\partial \mathcal{K}}\overline{p}\left(v·n_{\mathcal{K}}\right)ds.$

Following the ideas of Arnold et al.², an approximation for the pressure trace, $\overline{p}$, can be obtained by solving a global problem associated with the dual hybrid mixed formulation. Defining the function spaces M = {μ ∈ L ²(e), ∀e ∈ E _h }, U = Π_K U _K and L = Π_K L _K the dual hybrid formulation consists in:

Find u∈U, p∈Q and λ∈M, such that

$\sum _{\mathcal{K}\in {\mathcal{T}}_h}\left[{\left({\mathbb{K}}^{-1}u,v\right)}_{\mathcal{K}}-{\left(p,divv\right)}_{\mathcal{K}}+{\int }_{\partial \mathcal{K}}\lambda \left(v·n_{\mathcal{K}}\right)ds\right]=0,\forall v\in \mathcal{U},$ (3.3)

$\sum _{\mathcal{K}\in {\mathcal{T}}_h}\left[-{\left(divu,q\right)}_{\mathcal{K}}+{\left(f,q\right)}_{\mathcal{K}}\right]=0,\forall q\in \mathcal{L},$ (3.4)

$\sum _{\mathcal{K}\in {\mathcal{T}}_h}{\int }_{\partial \mathcal{K}}\mu \left(u·n_{\mathcal{K}}\right)ds=0,\forall \mu \in \mathcal{M}.$ (3.5)

The Lagrange multiplier λ is identified with trace of the pressure on the all edges of the elements K, $\lambda =\overline{p}$, unlike the classical primal hybrid formulation of Raviart-Thomas²⁷ and the stabilized formulation of Ewing et al.¹⁷, where the multiplier is identified with the flux. The third equation, (3.5), weakly imposes the continuity of the normal component of the velocity field (flux continuity) and the flux boundary condition u·n=0 on ∂Ω.

3.1 The SDHM Formulation

To generate a stable and adjoint consistent formulation, allowing greater flexibility in the choice of the finite element approximation spaces for velocity and pressure fields and the Lagrange multiplier, we add to the system (3.3)-(3.5) inner stabilization terms associated with least square residual forms, coming from the mass balance, the Darcy’s law and the curl of Darcy’s law, as in³. Furthermore, we add a stabilization term for the multiplier according to¹⁰ obtaining the Stabilized Dual Hybrid Mixed (SDHM) formulation:

Find u∈U, p∈Q and λ∈M ∀v∈U, q∈Q and μ∈M, such that

$\sum _{\mathcal{K}\in {\mathcal{T}}_h}\left[{\left({\mathbb{K}}^{-1}u,v\right)}_{\mathcal{K}}-{\left(p,divv\right)}_{\mathcal{K}}+{\int }_{\partial \mathcal{K}}\lambda \left(v·n_{\mathcal{K}}\right)ds+\frac{1}{2}{\left({||{\mathbb{K}}^{-1}||}_{\infty }\left(divu-f\right),divv\right)}_{\mathcal{K}}-\frac{1}{2}{\left({\mathbb{K}}^{-1}u+\nabla p,v\right)}_{\mathcal{K}}+\frac{1}{2}{\left({||\mathbb{K}||}_{\infty }\left(\nabla \times {\mathbb{K}}^{-1}u\right),\nabla \times {\mathbb{K}}^{-1}v\right)}_{\mathcal{K}}\right]-0,$ (3.6)

$\sum _{\mathcal{K}\in {\mathcal{T}}_h}\left[-{\left(divu,q\right)}_{\mathcal{K}}+{\left(f,q\right)}_{\mathcal{K}}-\frac{1}{2}{\left({\mathbb{K}}^{-1}u+\nabla p,\mathbb{K}\nabla q\right)}_{\mathcal{K}}+{\int }_{\partial \mathcal{K}}{||\mathbb{K}||}_{\infty }\beta \left(p-\lambda \right)qds\right]=0,$ (3.7)

$\sum _{\mathcal{K}\in {\mathcal{T}}_h}\left[{\int }_{\partial \mathcal{K}}\mu \left(u·n_{\mathcal{K}}\right)ds+{\int }_{\partial \mathcal{K}}{||\mathbb{K}||}_{\infty }\beta \left(\lambda -p\right)\mu ds\right]=0.$ (3.8)

where U = Π_K (H ¹(K) × H ¹(K)), L = ΠKH ¹(K), $\mathcal{M}=\left\{\mu \in L^2\left(e\right),\forall e\in {\mathcal{E}}_h^0\right\}$, ∇× is the curl operator and $\beta =\frac{{\beta }_0}{h}$ is the stabilization parameter associated with the Lagrange multiplier λ, and β₀ ∈ R is independent of h.

3.2 Finite Element Approximations

Let ${\mathcal{U}}_h^m=\left\{v\in {\mathcal{U}}_{\delta }:{\overline{)v}}_{\mathcal{K}}\in {\mathcal{R}}^m\times {\mathcal{R}}^m\forall \mathcal{K}\in {\mathcal{T}}_h\right\},{\mathcal{L}}_h^l=\left\{q\in {\mathcal{L}}_{\delta }:{\overline{)q}}_{\mathcal{K}}\in {\mathcal{R}}^l\forall \mathcal{K}\in {\mathcal{T}}_h\right\}and{\mathcal{M}}_h^s=\left\{\mu \in \mathcal{M}:{\overline{)\mu }}_e\in {\mathcal{P}}^s\forall e\in {\mathcal{E}}_h^0\right\}$ be the discontinuous Lagrangian finite element spaces where R ^r is the polynomial set with degree less than or equal to r if K is a triangle, or less than or equal to r in each cartesian variable if K is a quadrilateral (r=l or m), and P ^s is the polynomial set of degree less than or equal to s on each edge e. We can now present a finite element approximation for the stabilized dual hybrid mixed formulation introduced in the last section. Considering that {u _h , p _h }, belonging to the broken function spaces, are defined independently on each element K ∈T _h , we observe that system (3.6)-(3.8) can be split into a set of local problems defined on each element K coupled to the global problem defined on E _h , as follow:

Local problems: Find $u_h\in {\mathcal{U}}_h^m,p_h\in {\mathcal{L}}_h^l$, for each $\mathcal{K}\in {\mathcal{T}}_h,\forall v_h\in {\mathcal{U}}_h^{m}and\forall q_h{\mathcal{L}}_h^l$ such that

$\begin{array}{c}\left({\mathbb{K}}^{-1}{u}_h,v_h\right)\mathcal{K}-\left(p_h,div{v}_h\right)\mathcal{K}+{\int }_{\partial \mathcal{K}}{\lambda }_h\left(v_h·n_{\mathcal{K}}\right)ds-\frac{1}{2}{\left({\mathbb{K}}^{-1}{u}_h+\nabla p_h,v_h\right)}_{\mathcal{K}}\\ +\frac{1}{2}{\left({||{\mathbb{K}}^{-1}||}_{\infty }\left(div{u}_h-f\right),div{v}_h\right)}_{\mathcal{K}}+\frac{1}{2}{\left({||\mathbb{K}||}_{\infty }\left(\nabla \times {\mathbb{K}}^{-1}{u}_h\right),\nabla \times {\mathbb{K}}^{-1}{v}_h\right)}_{\mathcal{K}}=0,\end{array}$ (3.9)

$-{\left(div{u}_h,q_h\right)}_{\mathcal{K}}+{\left(f,q_h\right)}_{\mathcal{K}}-\frac{1}{2}{\left({\mathbb{K}}^{-1}{u}_h+\nabla p_h,\mathbb{K}\nabla q_h\right)}_{\mathcal{K}}+{\int }_{\partial \mathcal{K}}{||\mathbb{K}||}_{\infty }\beta \left(p_h-{\lambda }_h\right)q_{h}ds=0,$ (3.10)

Global Problem: Find ${\lambda }_h\in {\mathcal{M}}_h^s,\forall {\mu }_h\in {\mathcal{M}}_h^s$ such that

$\sum _{\mathcal{K}\in {\mathcal{T}}_h}\left[{\int }_{\partial \mathcal{K}}{\mu }_h\left(u_h·n_{\mathcal{K}}\right)ds+{\int }_{\partial \mathcal{K}}{||\mathbb{K}||}_{\infty }\beta \left({\lambda }_h-p_h\right){\mu }_{h}ds\right]=0.$ (3.11)

We observe that the SDHM method is consistent and provides optimal rates of convergence²⁴, ensuring flexibility in the choice of the approximation spaces and interpolation functions, unlike the classical dual formulation which the stability is restricted to appropriate choices of the finite elements spaces, such as Raviart-Thomas²⁷ and Brezzi-Douglas-Marini (BDM)⁷ that requiring compromises between the approximation spaces. Moreover, the SDHM method is locally conservative for equal order approximations of all fields (l=m=s) and stable for any value of the edge stabilization parameter β, including β=0^{24), (23), (25}. The choose of the multiplier as the trace of pressure is crucial to assure that the local problems (3.9)-(3.10) are solvable for the variables ${\overline{)\left\{u_h,p_h\right\}}}_{\mathcal{K}}\in {\mathcal{U}}_h^m\times {\mathcal{L}}_h^l$ as a function of the multiplier λ_h .

3.2.1 Solver Strategy

Here we define on each element the operators

$\begin{array}{c}{a}_{\mathcal{K}}\left(\left[u_h,p_h\right];\left[v_h,q_h\right]\right)={\left({\mathbb{K}}^{-1}{u}_h,v_h\right)}_{\mathcal{K}}-{\left(p_h,div{v}_h\right)}_{\mathcal{K}}+\frac{1}{2}{||{\mathbb{K}}^{-1}||}_{\infty }{\left(div{u}_h,div{v}_h\right)}_{\mathcal{K}}\\ -{\left(div{u}_h,q_h\right)}_{\mathcal{K}}-\frac{1}{2}{\left({\mathbb{K}}^{-1}{u}_h+\nabla p_h,v_h+\mathbb{K}\nabla q_h\right)}_{\mathcal{K}}\\ +\frac{1}{2}\left({||\mathbb{K}||}_{\infty }\nabla \times \left({\mathbb{K}}^{-1}{u}_h\right),\nabla \left(\times {\mathbb{K}}^{-1}{v}_h\right)\right)\mathcal{K},\\ b_{\mathcal{K}}\left({\lambda }_h;\left[v_h,q_h\right]\right)={\int }_{\partial \mathcal{K}}{\lambda }_h\left(v_h·n\right)ds+{\int }_{\partial \mathcal{K}}\beta {||\mathbb{K}||}_{\infty }{\lambda }_h{q}_h{d}_s,\\ c_{\mathcal{K}}\left({\lambda }_h,{\mu }_h\right)={\int }_{\partial \mathcal{K}}{||\mathbb{K}||}_{\infty }\beta {\lambda }_h{\mu }_{h}ds,\end{array}$

and the functional

$f_{\mathcal{K}}\left(\left[v_h,q_h\right]\right)={\left(f,q_h\right)}_{\mathcal{K}}+\frac{1}{2}{||{\mathbb{K}}^{-1}||}_{\infty }{\left(f,div{v}_h\right)}_{\mathcal{K}}.$

Then, the SDHM method is now reformulated as

Find $u_h\in {\mathcal{U}}_h^m,p_h\in {\mathcal{L}}_h^l$, for each $\mathcal{K}\in {\mathcal{T}}_{h}and{\lambda }_h\in {\mathcal{M}}_h^s$ such that

$a_{\mathcal{K}}\left(\left[u_h,p_h\right];\left[v_h,q_h\right]\right)+b_{\mathcal{K}}\left({\lambda }_h;\left[v_h,q_h\right]\right)=f_{\mathcal{K}}\left(\left[v_h,q_h\right]\right),$ (3.12)

$\sum _{\mathcal{K}\in {\mathcal{T}}_h}{b}_{\mathcal{K}}^T\left(\left[u_h,p_h\right],{\mu }_h\right)+\sum _{\mathcal{K}\in {\mathcal{T}}_h}{c}_{\mathcal{K}}\left({\lambda }_h,{\mu }_h\right)=0,$ (3.13)

$\forall \left[v_h,q_h\right]\in {\mathcal{U}}_h^m\times {\mathcal{L}}_h^{l}and\forall {\mu }_h\in {\mathcal{M}}_h^s$. Considering A_K , B _K and C _K matrices generated respectively by the local operators a _K (⋅, ⋅), b _K (⋅, ⋅) and c _K (⋅, ⋅) and F _K the vector given by f _K (⋅), we can rewrite (3..12)-(3.13) in the following matrix form,

$A_{\mathcal{K}}U+B_{\mathcal{K}}\Lambda =F_{\mathcal{K}},\forall \mathcal{K}\in {\mathcal{T}}_h$ (3.14)

$\sum _{\mathcal{K}\in {\mathcal{T}}_h}{B}_{\mathcal{K}}^{T}U+\sum _{\mathcal{K}\in {\mathcal{T}}_h}{C}_{\mathcal{K}}\Lambda =0,$ (3.15)

Using the static condensation strategy, and given that the local matrix A _K is positive definite, we have from (3.14) that U = {u _h , p _h } can be expressing in terms of Λ = {λ_h } as

$U=A_{\mathcal{K}}^{-1}\left(F_{\mathcal{K}}-B_{\mathcal{K}}\Lambda \right).$ (3.16)

Then, replacing (3.16) in (3.15), we obtain the global system only for Λ = {λ_h }

$\sum _{\mathcal{K}\in {\mathcal{T}}_h}\left(C_{\mathcal{K}}-B_{\mathcal{K}}^T{A}_{\mathcal{K}}^{-1}{B}_{\mathcal{K}}\right)\Lambda =-\sum _{\mathcal{K}\in {\mathcal{T}}_h}{B}_{\mathcal{K}}^T{A}_{\mathcal{K}}^{-1}{F}_{\mathcal{K}}.$ (3.17)

After solving (3.17), the vector U is obtained from (3.16). A great advantage of this methodology is the size reduction of the overall system, now involving only the degrees of freedom associated with the multipliers Λ, leading to a reduced computational cost, since the time needed to solve all local problems is negligible compared to the time to solve the global system.

4 TRACER INJECTION APPROXIMATIONS

4.1 Fully discrete SUPG approximation

A common approach to transient problems is based on fully discrete formulations obtained by combining finite difference approximations in time with finite element methods in space. According to the Rothe method (or horizontal method of lines)¹⁸ of first discretizing in time and then in space on each discrete time level, we choose the partition I _Δ = {0 = t ⁰ < t ¹ < ... < t ^N = T} of the interval I=[0,T] with Δt _n = t ⁿ - t ^{n - 1} and Δt = max _n Δt ⁿ . Then, using a backward finite difference scheme to approximate the time derivative in equation (2.3), we have the following sequentially implicit time-stepping algorithm: for n=0,1,…,N-1, given f, g and c ₀(x), find c ^{n + 1} satisfying

$\varphi \frac{c^{n+1}-c^n}{\Delta t}+u^n·\nabla c^{n+1}-div\left(\mathbb{D}\left(u^n\right)\nabla c^{n+1}\right)=g^{n+1}in\Omega ,$ (4.1)

$\mathbb{D}\left(u^n\right)\nabla c^{n+1}·n=0on\partial \Omega$ (4.2)

with

$u^n=-\mathbb{K}\nabla p^{n}in\Omega ,$ (4.3)

$div{u}^n=f^{n}on\Omega ,$ (4.4)

and u ⁿ ⋅ n = 0 on ∂Ω and ∫_Ω p ⁿ dx = 0, t ∈ (0, T). The concentration is achieved at time n+1 and the velocity and pressure are given at time n. A complete numerical analysis, demonstrating existence and uniqueness of solution for the above semi-discrete system can be found in^{21), (20}. Note that this sequentially implicit method can be written in predictor-corrector form and the original system becomes partially uncoupled and linearized.

For the tracer injection processes, the velocity and pressure approximations are calculated just once at the beginning of the process, u ⁿ ≡ u∀n=0,1,…,N-1 and, then the concentration is obtained for all instants. Hereafter, to simplify the notation, the superscript n in the velocity and pressure approximations are dropped out.

We combine the semi-discrete approximation (4.1)-(4.2) with a stabilized finite element method in space (the SUPG method), and introduce the following fully discrete approximation for the concentration equation^{24), (23), (25}:

For time levels n=0,1,2,…, given u _h and c ⁰ find $c_h^{n+1}\in X_h^r$, such that

$B\left(u_h;c_h^{n+1},{\eta }_h\right)=F\left(c_h^n;{\eta }_h\right),\forall {\eta }_h\in X_h^r,$ (4.5)

with $\left(c_h^0,{\eta }_h\right)=\left(c^0,{\eta }_h\right),\forall {\eta }_h\in X_h^r$, where $X_h^r$ is a continuous Lagrangian finite element space of degree r≥1; and

$\begin{array}{c}B\left(u_h;c_h^{n+1},{\eta }_h\right)=\left(c_h^{n+1},{\eta }_h\right)+\Delta t\left(u_h·\nabla c_h^{n+1},{\eta }_h\right)+\Delta t\left(D\left(u_h\right)\nabla c_h^{n+1},\nabla {\eta }_h\right)\\ +\sum _{\mathcal{K}}\left(c_h^{n+1}+\Delta t{u}_h·\nabla c_h^{n+1},{\delta }_K{u}_h·\nabla {\eta }_h\right)\mathcal{K}\\ +\sum _{\mathcal{K}}\left(-\Delta tdiv\left(D\left(uh\right)\nabla c_h^{n+1}\right),{\delta }_K{u}_h·\nabla {\eta }_h\right)\mathcal{K},\end{array}$ (4.6)

$F\left(c_h^n,{\eta }_h\right)=\left(c_h^n+\Delta t{g}^{n+1},{\eta }_h\right)+\sum _{\mathcal{K}}{\left(c_h^n+\Delta t{g}^{n+1},{\delta }_K{u}_h·\nabla {\eta }_h\right)}_{\mathcal{K}}.$ (4.7)

The stabilization parameter δ_K is defined on each K ∈ T _h by^{9), (21), (20} as

${\delta }_{\mathcal{K}}=\left\{\begin{array}{ll}\frac{h_{\mathcal{K}}}{2{||u_h||}_{L^{\infty }}\left(\mathcal{K}\right)},& P{e}_{\mathcal{K}}\ge 1,\\ 0,& 0<P{e}_{\mathcal{K}}<1,\end{array}\right.$ (4.8)

where

$P{e}_{\mathcal{K}}=\frac{m_{\mathcal{K}}{||u_h||}_{L^{\infty }}\left(\mathcal{K}\right)h_{\mathcal{K}}}{D_{\mathcal{K}}^2/d_{\mathcal{K}}}$ (4.9)

is the mesh-dependent Péclet number, $d_{\mathcal{K}}={\alpha }_{mol}+{\alpha }_t\underset{x\in \mathcal{K}}{inf}\left|u_h\left(x\right)\right|,m_{\mathcal{K}}=\frac{2}{3}min\left(\frac{1}{2},c_{inv}\right),D_{\mathcal{K}}=\sqrt{2{\left({\alpha }_{mol}+{\alpha }_l{||u_h||}_{L^{\infty }\left(\mathcal{K}\right)}\right)}^2+2{\left(3{\alpha }_l-2{\alpha }_t\right)}^2{||\nabla u_h||}_{L^{\infty }\left(\mathcal{K}\right)}^2{h^2}_{\mathcal{K}}{c}_{inv}}$, where c _inv is the typical inverse constant of finite element spaces. The terms in the right-hand side of (4.6)-(4.7) multiplied by δ_K are responsible for the additional stability of this method⁹. In the numerical experiments this proposal will be identified as the SUPG approach.

4.2 A semi-analytical methodology

Following¹⁴, we present a semi-analytical methodology which combines an analytical solution for the concentration given by¹ with a finite element approximation for the velocity field.

According to¹, neglecting transversal flux dispersion (α_mol = α_t = 0), the tracer concentration c(t) at the producer well is expressed by

$c\left(t\right)=\frac{\partial V_{tr}}{2q_t{\left(\pi {\alpha }_l\right)}^{1/2}}\sum _{n=1}^{N_{sl}}\frac{1}{1/2}exp\left(-\frac{{\left(t_{btn}-t\right)}^2}{4{\alpha }_l{l}_{\mathit{n}}}\right),$ (4.10)

where N _sl is the number of streamlines arriving at the production well, ĉ is the injection tracer concentration, V _tr is the injected tracer volume, q _t is the total injection rate, α_l is the longitudinal dispersion coefficient, t _btn is the breakthrough time for the streamline n given by $t_{btn}={\int }_{x_{iw}}^{x_{pw}}\frac{ds}{v_n}and{I}_n={\int }_{x_{iw}}^{x_{pw}}\frac{ds}{{v_n}^2}$ where v _n is the Darcy velocity on the streamline n, with x _iw and x _pw the coordinates of the injection and production wells, respectively.

In¹ the Darcy’s velocity v _n is also obtained analytically for particular tracer injection problems in homogeneous porous media. However, from the definitions of t _btn and I _n it is clear that the semi-analytical methodology (5.10) can be naturally reformulated replacing v _n by an approximation of the velocity field. In Section 5 we will observe, through numerical experiments, the influence of the velocity approximations on some tracer transport simulations in homogeneous and heterogeneous porous media solving the Darcy system (2.1)-(2.2) via the SDHM method, the classical Galerkin method and a post-processing technique.

Classical Galerkin method. Substituting (2) into (1) the Galerkin approximation of the elliptic sub-system reads: Find $p_h\in {\mathrm{\mathbb{N}}}_h^r$, such that

$\left({\mathbb{K}}^{-1}\nabla p_h,\nabla \phi \right)=\left(f,{\phi }_h\right),\forall {\phi }_h\in {\mathrm{\mathbb{N}}}_h^r,$ (4.11)

with ${\mathrm{\mathbb{N}}}_h^r$a continuous Lagrangian finite element space. From the pressure approximation (4.11) and Darcy’s law (2.2), it is natural to calculate the velocity field giving by

$u_G=-\mathbb{K}\nabla p_h.$ (4.12)

This approach generates a discontinuous velocity field at the element interfaces which do not satisfy the boundary condition U _G ⋅ n and, in addition, has a sub-optimal convergence rate.

Post-processing. To approximate the velocity field with improved accuracy we employed a post-processing technique²⁰ based on a variational formulation of Darcy’s law combined with the residual of the balance equation:

For p _h given by (4.11), find $u_{pp}\in S{S}_h^q=\left\{v_h\in {\mathrm{\mathbb{N}}}_h^r\times {\mathrm{\mathbb{N}}}_h^r,v_h·n=0on\partial \Omega \right\}$such that

$\left({\mathbb{K}}^{-1}{u}_{pp}+\nabla p_h,v_h\right)+\tau \left(div{u}_{pp}-f,div{v}_h\right)=0,\forall v_h\in S{S}_h^q.$ (4.13)

For homogeneous media with regular solution this post-processing presents a gain of O(h ^1/2)²⁰ compare to the classical Galerkin approximation (4.12). However, it does not produce accurate solution when applied to heterogeneous porous media as we will observe in the numerical experiments.

5 NUMERICAL SIMULATIONS

To illustrate the performance of the methodology proposed in Sections 3 and 4 we present the results of tracer injection simulations in a quarter of a repeated five-spot pattern in two dimensions consisting of a square domain (unit thickness) with side L=1000 ft. The injector well is located at the lower-left corner (x=y=0) and the producer well at the upper-right corner (x=y=L). The coefficients have been chosen as α_mol = 0.0, α_l = 1.0 ft, α_t = 0.0 ft and the porosity Φ =0.1. The tracer slug is 0.25% of the porous volume. These data were taken from²¹.

In all simulations we use uniform meshes of 80 × 80 bilinear quadrilateral elements to approximate concentration as well as velocity and pressure. Equal order interpolation functions (m=l=r=s=1) are employed to all variables. Discontinuous Lagrangian interpolation to the SDHM method and continuous Lagrangian interpolation to the Galerkin method and post-processing technique are used. The post-processing stabilization parameter is fixed as τ =1.0 and the stabilization parameter associated with the Lagrange multiplier in SDHM is β₀ = 0. The semi-analytical methodology parameters are ĉ=1.0, V _tr = 1000 ft ², q _t = 200 ft ²/day and N _sl = 399.

In order to compare the differents approaches described in Section 3, we take into account three scenarios: Case 1, Case 2 and Case 3. In Case 1 the porous medium is homogeneous with the permeability constant, K=1000 mD and the other scenarios consider heterogeneous porous media with subregions having different permeability values (Figure 1).

Figure 1:
Heterogeneous scenarios: (a) Case 2 and (b) Case 3.

5.1 Homogeneous Porous Medium

Figure 2 exhibits comparisons between the results of the SUPG approach and the semi-analytical methodology for the concentration with velocity approximations obtained by the SDHM formulation, the post-processing technique (PP) and the Galerkin method. For both approximations we see that the SDHM method has the closest profile to the analytical solution obtained by Abbaszadeh-Dehghani¹ (peak at 0.01) as showed in Figure 2-(b). The conservative property of the SDHM method is illustrated in Figure 3. Therefore, we can conclude that the SDHM formulation is the more stable and provides accurate results in all approaches considered here for the transport equation.

Figure 2:
Case 1 - Time history of the tracer concentration to different velocity field approximations.

Figure 3:
Case 1 - Tracer volume of the tracer concentration to different velocity field approximations.

Next, we study the influence of the Péclet number (Pe) on the tracer injection processes. Figure 4 shows the time history of the tracer concentration in the producer well with the SUPG approach for Péclet numbers corresponding to 100, 50, 25 and 12.5. As the Péclet number decrease the peak concentration value tends to the Abbaszadeh-Dehghani exact solution. For high Péclet numbers (Pe ≥50) spurious oscillations associated with the use of the SUPG approach are identified. These oscillations become less pronounced when we combine the SUPG approach with the SDHM method, as can be best viewed in the zoom regions plotted in Figure 4 (right side).

Figure 4:
Case 2- Tracer concentration maps. (a) Galerkin method, (b) post-processing technique and (c) the SDHM method.

5.2 Heterogeneous Porous Media

In this subsection, we consider two problems defined in the heterogeneous domains plotted in Figure 1 (Case 2 and Case 3). Figures 5 and 6 show the tracer concentration maps at t=1000 days when we use the SUPG approach combined with the Galerkin method, the post-processing technique and the SDHM method to Δt=5 days. Similar transport behaviours are observed for all velocity approximations.

Figure 5:
Case 3 - Tracer concentration maps. (a) Galerkin method, (b) post-processing technique and (c) the SDHM method.

Figure 6:
Case 2 - Time history of the tracer concentration to different velocity field approximations.

The time history of the tracer concentration in the producer well for Case 2 and Case 3 using the SUPG approach are showed in Figures 7 and 8 for the concentration with velocity approximations obtained by the SDHM formulation, the post-processing technique (PP) and the Galerkin method. Figure 7 exhibits two concentration peaks due to the influence of the lower permeability region, which act like a barrier to the flow. The SDHM formulation produces the highest peak according to the physical expected behavior. Note that the tracer concentration reaches the producer well and, the peak concentration is achieved, in an earlier time compared with that observed in Case 1 (see Figure 2). This is due to the fact that there are preferential flow paths where the velocity flow is increased as a result of the presence of a subregion of lower permeability.

Figure 7:
Case 3 - Time history of the tracer concentration to different velocity field approximations.

Figure 8:
Time history of the tracer concentration. Comparisons among Abbaszadeh-Dehghani solution with the semi-analytical methodology combined with the (a) Post-Processing technique and the (b) SDHM method.

Next in Figure 8 we observe a similar behavior to the homogeneous scenario (Case 1) but with higher concentration values. This can be better understand in Case 1 and Case 3 results plot in Figure 9, where we compare the time history of the tracer concentration obtained by the semi-analytical solution combined with the post-procesing technique and the SDHM formulation. The subregions with differents permeabilities have the same domain width. Therefore, they act as a delay to the flow. Thus, the concentration reaches the producer well later than the Case 1.

Figure 9:
Case 3 - Tracer volume for the tracer concentration. Comparison among the velocity field approximations: (a) and (b) SDHM, Post-Processing (PP) and Galerkin.

Figure 10:
Case 3 - Tracer volume for the tracer concentration. Comparison among the velocity field approximations: (a) and (b) SDHM, Post-Processing (PP) and Galerkin.

6 CONCLUSION

In this paper the tracer injection process into homogeneous and heterogeneous media was solved using two different methodologies for the concentration approximation combined with three formulations for the velocity field. When the SDHM formulation is employed in the calculation of velocity, combined with the SUPG approach their results lead to more accurate approximations for the concentration. This combination is able to capture the expected discontinuity properties of the solution and, consequently, a proper physical solution. The numerical experiments performed, illustrated the flexibility and robustness of this formulation.

REFERENCES

• ¹
  M. Abbaszadeh-Dehghani & W.E. Brigham. Technical report. Stanford University, (1982).
• ²
  D.N. Arnold & F. Brezzi. Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates. RAIRO MMAN, 19(7) (1985), 7-32.
• ³
  D.N. Arnold, F. Brezzi, B. Cockburn & L.D. Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM Journal on Numerical Analysis, 39(5) (2002), 1749-1779.
• ⁴
  G. Barrenechea, L.P. Franca & F. Valentin. A Petrov-Galerkin enriched method: A mass conservative finite element method for the Darcy equation. Computer Methods in Applied Mechanics and Engineering, 196 (2007), 2449-2464.
• ⁵
  T.P. Barrios, J.M. Cascón & M. González. A posteriori error analysis of an augmented mixed finite element method for Darcy flow. Computer Methods in Applied Mechanics and Engineering, 283 (2015), 909-922.
• ⁶
  F. Brezzi. On the existence, uniqueness and approximation of saddle point problems arising from Lagrange multipliers. Analyse numérique/Numerical Analysis (RAIRO), 8(R-2) (1974), 129-151.
• ⁷
  F. Brezzi, J. Douglas & L.D. Marini. Two families of mixed finite elements for second order elliptic problems. Numerische Mathematik, 47 (1985), 217-235.
• ⁸
  F. Brezzi & M. Fortin. “Mixed and Hybrid Finite Element Methods”. Springer-Verlag (1991).
• ⁹
  A.N. Brooks & T.J.R. Hughes. Streamline Upwind Petrov-Galerkin Formulations for Convection-Dominated flows with Particular emphasis on the Incompressible Navier-Stokes Equations. Comput. Methods Appl. Mech. Engrg., 32 (1982), 199-259.
• ¹⁰
  M.R. Correa & A.F.D. Loula. Stabilized velocity post-processings for Darcy flow in heterogenous porous media. Communications in Numerical Methods in Engineering, 23 (2007), 461-489.
• ¹¹
  M.R. Correa & A.F.D. Loula. Unconditionally stable mixed finite element methods for Darcy flow. Computer Methods in Applied Mechanics and Engineering, 197 (2008), 1525-1540.
• ¹²
  I.H.A. da Igreja. “Métodos de elementos finitos híbridos estabilizados para escoamentos de Stokes, Darcy e Stokes-Darcy acoplados”. Ph.D. thesis, Laboratório Nacional de Computação Científica, Petrópolis, Brasil (2015). URL http://tede.lncc.br/handle/tede/225
  » http://tede.lncc.br/handle/tede/225
• ¹³
  B.L. Darlow, R.E. Ewing & M.F. Wheeler. Mixed finite element method for miscible displacement problems in porous media. SPE Journal, (1984), 391-398.
• ¹⁴
  A. Datta-Gupta & M.J. King. A semianalytic approach to tracer flow modeling in heterogeneous permeable media. Advances in Water Resources, 18(1) (1995), 9-24.
• ¹⁵
  J.J. Douglas, R.E. Ewing & M.F. Wheeler. The approximation of the pressure by a Mixed-Method in the simulation of miscible displacement. R.A.I.R.O. Analyse Numerique, 17 (1983), 17-33.
• ¹⁶
  J.J. Douglas & J. Wang. An absolutely stabilized finite element method for the Stokes problem. Math. Comput., 52(186) (1989), 495-508.
• ¹⁷
  R.E. Ewing, J. Wang & Y. Yang. A stabilized discontinuous finite element method for elliptic problems. Numerical Linear Algebra with Applications, 10 (2003), 83-104.
• ¹⁸
  I. Harari. Stability of semidiscrete formulations for parabolic problems at small time steps. Comput. Methods Appl. Mech. Engrg., 193 (2004), 1491-1516.
• ¹⁹
  A.F.D. Loula, F.A. Rochinha & M.A. Murad. Higher-order gradient post-processings for second-order elliptic problems. Computer Methods in Applied Mechanics and Engineering, 128 (1995), 361-381.
• ²⁰
  S.M.C. Malta & A.F.D. Loula. Numerical analysis of finite element method for miscible displacements in porous media. Numerical Methods in Partial Differential Equations, 14 (1998), 519-548.
• ²¹
  S.M.C. Malta, A.F.D. Loula & E.L.M. Garcia. Numerical analysis of a stabilized finite element method for tracer injection simulations. Comput. Methods Appl. Mech. Engrg., 187 (2000), 119-136.
• ²²
  A. Masud & T.J.R. Hughes. A stabilized finite element method for Darcy flow. Computer. Methods Appl. Mech. Engrg, 191 (2002), 4341-4370.
• ²³
  Y.R. Núñez. “Métodos de elementos finitos híbridos aplicados a escoamentos miscíveis em meios porosos heterogêneos”. Ph.D. thesis, Laboratório Nacional de Computação Científica, Petrópolis, Brasil (2014). URL http://tede.lncc.br/handle/tede/203.
• ²⁴
  Y.R. Núñez, C.O. Faria, A.F.D. Loula & S.M.C. Malta. A mixed-hybrid finite element method applied to tracer injection processes. International Journal of Modeling and Simulation for the Petroleum Industry, 6(1) (2012), 51-59.
• ²⁵
  Y.R. Núñez, C.O. Faria, A.F.D. Loula & S.M.C. Malta. A hybrid finite element method applied to miscible displacements in heterogeneous porous media. Rev. Int. de Métodos Numér. Cálc. Diseño Ing., 33(1-2) (2017), 45-51.
• ²⁶
  D.W. Peaceman. “Fundamental of Numerical Reservoir Simulation”. Elsevier, Amsterdam (1977).
• ²⁷
  P.A. Raviart & J.M. Thomas. A Mixed Finite Element Method for Second Order Elliptic Problems. In I. Galligani & E. Magenes (editors), “Lecture Notes in Math”, volume 606. New York, Springer-Verlag (1977).
• ²⁸
  B. Rivière. “Discontinuous Galerkin Methods For Solving Elliptic And Parabolic Equations: Theory and Implementation”. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA (2008). Tend.

Publication Dates

• Publication in this collection
  May-Aug 2018

History

• Received
  14 Dec 2017
• Accepted
  19 Mar 2018