Hybrid Mixed Methods Applied to Miscible Displacements with Adverse Mobility Ratio

Abstract

We propose stable and locally conservative hybrid mixed finite element methods to approximate the Darcy system and convection-diffusion problem, presented in a mixed form, to solve miscible displacements considering convective flows with adverse mobility ratio. The stability of the proposed formulations is achieved due to the choice of non-conforming Raviart-Thomas spaces combined to upwind scheme for the convection-dominated regimes, where the continuity conditions, between the elements, are weakly enforced by the introduction of Lagrange multipliers. Thus, the primal variables of both systems can be condensed in the element level leading a positive-definite global problem involving only the degrees of freedom associated with the multipliers. This approach, compared to the classical conforming Raviart-Thomas, present a reduction of the computational cost because, in both problems, the Lagrange multiplier is associated with a scalar field. In this context, a staggered algorithm is employed to decouple the Darcy problem from the convection-diffusion mixed system. However, both formulations are solved at the same time step, and the time discretization adopted for the convection-diffusion problem is the implicit backward Euler method. Numerical results show optimal convergence rates for all variables and the capacity to capture the formation and the propagation of the viscous fingering, as can be seen in the comparisons of the simulations of the Hele-Shaw cell with experimental results of the literature.

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. This work was also partially supported by CNPq and UFJF.

Keywords

1 Introduction

The miscible displacement of a higher viscosity fluid in a porous medium raises considerable attention to a variety of applications like hydrology, blood motion in vessels, industrial processes involving filtration such as sugar refining, carbon sink, oil recovery and groundwater exploration in porous media [6, 14, 30]. During the transport of an injected fluid in a reservoir, it is common the formation of fingers and channels in the flow, and this can happen even in a homogeneous media if the injected fluid is less viscous than the resident fluid. Injection of a fluid with different viscosity of the reservoir-resident fluid usually produces complex concentration patterns along with displacement [5, 6, 18, 28].

In the numerical point of view, since approximations are adopted to simulate the fluids displacement with different viscosities using meshes of elements with the cell diameter smaller than the length of a formed finger, this can lead to inaccurate representations of the interfaces of these fronts, due to the tendency that mixtures, with a higher mobility ratio M, have to generate highly branched fractal structures [5, 14, 29]. The mobility ratio relates the viscosity of the resident fluid $\nu _ {res}$ to the viscosity of the injected fluid $\nu _ {inj}$

$$\begin{aligned} M = \nu _{res}/\nu _{inj}, \end{aligned}$$

(1)

for $M>1$, or adverse mobility ratio, the injected fluid is less viscous than resident fluid, and in this case, the experiments predict that the displacement front is physically unstable where small perturbations can be forming multiple viscous fingerings [1, 18, 26, 27]. This unstable behavior is also true mathematically, as can be seen in [27] and in references therein. Therefore, this problem requires the employ of robust numerical methods that accurately solve this phenomenon.

The Partial Differential Equations (PDE) that govern the phenomenon of the displacement of the fluid mixtures consist of a system formed by Darcy’s problem and Transport problem. Some successful numerical approximations employing finite element methods to solve the Darcy problem can be found in [2, 8, 11, 16, 19, 25], for the Transport equation can be seen in [9, 10, 12, 17, 20, 23] and for the Darcy-Transport coupled system we mention the works [13, 15, 17, 21, 31].

The objective of this work is to propose an equivalent finite element method for Darcy and Transport problems. In this sense, we rewrite the Transport equation in a mixed form in terms of diffusive flux and concentration. Thus, we have two mixed systems. However, the use of finite element methods for mixed problems requires compatibility between the approximation spaces [3, 7]. Stable formulations have proven successful as can be seen in [8, 25]. The Raviart-Thomas spaces [25], referenced here by $\mathcal {RT}_k$, was developed for the mixed problems, like Darcy problem, is be able to simulate problems in heterogeneous medium with stability, mass conservation and optimal convergence rates. However, Arnold and Brezzi in [2] proposed an hybridization, employing Lagrange multipliers associated to the scalar field, that gives rise to a positive-definite system obtaining the same accuracy of the conforming Raviart-Thomas spaces but with fewer unknowns. An extension of this approach for the mixed convection-diffusion problems are developed and analyzed by Egger and Schöberl [12], where a formulation is proposed by the combination of upwind techniques used in discontinuous Galerkin methods for hyperbolic problems with conservative discretizations of mixed methods for elliptic problems.

In this context, we propose an equivalent locally conservative numerical method for the Darcy and the convection-diffusion problems employing the Raviart-Thomas spaces in a non-conforming way in both formulations. In this case, the continuity is weakly imposed via Lagrange multipliers associated with the scalar field defined on the edges/faces of the elements. Thus, all interest variables (velocity, pressure, diffusive flux and concentration) degrees-of-freedom are condensed in the element level leading to a global problem involving only the degrees-of-freedom of the multiplier. After solving the global problem, the approximate solution of the multiplier is plugged into the local problems to recover the discontinuous approximation of the primal variables. This approach significantly reduces the computational cost compared to the use of conforming Raviart-Thomas spaces.

The methodology to solve this system is based on a staggered algorithm to decouple the hydrodynamics from the hyperbolic system, resulting in a scheme that uses a locally conservative hybrid mixed finite element method to approximate both problems. The two problems are solved in the same time scale, applying an implicit backward finite difference scheme in time to approximate the Transport equation. Moreover, the spurious oscillations characteristics of the convection-dominated regimes are mitigated through of an upwind scheme associated with the multiplier and concentration values, evaluated on the edges of the elements [12]. Numerical simulations are presented and demonstrate optimal convergence rates for all variables and a good capture of the viscous fingering instabilities on the interface between the miscible fluids compared to experimental simulations in Hele-Shaw cells with rectilinear flows.

This paper is organized as follows: In Sect. 2, we present the Darcy-transport model problem. In Sect. 3, notations and definitions required to present the hybrid methods are described. The stable and equivalent mixed hybrid method for the Darcy and transport is presented in Sect. 4. Section 5 is devoted to convergence study and numerical simulations, comparing the approximate solution with experimental results in Hele-Shaw cell, considering adverse mobility ratio and high Péclet number. And finally, in Sect. 6, we present the concluding remarks of this work.

2 Model Problem

The model problem is described by a PDE system composed by two subsystems, the Darcy problem, considering a incompressible flow and neglecting gravitational forces, and the Transport equation. Thus, let $\varOmega \subset \mathbb {R}^d$ be the domain, with $d = 2,3$, and the boundary $\partial \varOmega = \partial \varOmega _N \cup \partial \varOmega _D$ in time interval (0, T], the problem can be written as follows

Given the concentration c and the functions $\overline{p}$ and f, find the pair $[\mathbf {u}, p]$, such that:

$$\begin{aligned} \begin{array}{ll} \mathbf {u} = - \mathbb {K}\nabla p \quad &{}\text {in} \quad \varOmega ,\\ {\text {div}}( \mathbf {u} ) = f \quad &{}\text {in} \quad \varOmega ,\\ \mathbf {u} \cdot \mathbf {n} = 0 \quad &{}\text {on} \quad \partial \varOmega _N,\\ p = \overline{p} \quad &{}\text {on} \quad \partial \varOmega _D. \end{array} \end{aligned}$$

(2)

Given the Darcy velocity $\mathbf {u}$, porosity $\phi$ and the functions g, $c_0$ and $\overline{c}$, find c, such that:

$$\begin{aligned} \begin{array}{ll} \phi \dfrac{\partial c}{\partial t}+\mathbf {u} \cdot \nabla c-{\text {div}}(\mathbb {D} \nabla c)=g_{} &{}\text {in} \quad \varOmega \times (0, T], \\ c(\mathbf {x}, 0) = c_{0}(\mathbf {x}) \quad &{}\text {in} \quad \varOmega , \\ \mathbb {D} \nabla c \cdot \mathbf {n} = 0 \quad &{}\text {on} \quad \partial \varOmega _N \times (0, T], \\ c = \overline{c} \quad &{}\text {on} \quad \partial \varOmega _D \times (0, T]. \end{array} \end{aligned}$$

(3)

where the variables of interest, p, $\mathbf {u}$, and c are respectively the hydrostatic pressure, the Darcy’s velocity and the fluid concentration. The coefficients of the equations are $\phi = \phi (\mathbf {x})$ the effective porosity of the medium, $\mathbb {K} = \mathbb {K}(c,\mathbf {x}) = \frac{\mathbb {G}}{\nu }$ the hydraulic conductivity tensor, a proportionality coefficient that takes into account the characteristics of the medium, including size, distribution, form, and arrangement of the particles, and the viscosity of the fluids. Thus, $\mathbb {G} = \mathbb {G}(\mathbf {x})$ and $\nu = \nu (c)$ denotes the medium permeability and the fluid viscosity, respectively. Finally, $\mathbb {D} = \mathbb {D}(\mathbf {u})$ is the dispersion tensor that can be defined as

$$\begin{aligned} \mathbb {D}(\mathbf {u})= \alpha _m\mathbb {I}+ \Vert \mathbf {u}\Vert \left[ \alpha _l \mathbb {E}(\mathbf {u}) + \alpha _t(\mathbb {I} - \mathbb {E}(\mathbf {u}))\right] , \quad \mathbb {E}(\mathbf {u}) = \dfrac{\mathbf {u} \otimes \mathbf {u}}{\Vert \mathbf {u}\Vert ^2}, \end{aligned}$$

with $\Vert \mathbf {u}\Vert ^2=\sum _{i=1}^d u_i^2$, $\otimes$ denoting the tensorial product, $\mathbb {I}$ the identity tensor, $\alpha _m$ being a molecular diffusion coefficient, $\alpha _l$ the longitudinal dispersion and $\alpha _t$ the transverse dispersion. In miscible displacement of a fluid through another in a reservoir, the dispersion is physically more important than the molecular diffusion [15, 24]. Thus, we assume the following properties $0< \alpha _m \le \alpha _l, \quad \alpha _l\ge \alpha _t > 0 \quad \text{ and } \quad 0< \phi \le 1.$

With the gravitational effect neglected, besides the mobility ratio, another dimensionless parameter determines the behavior of the model is the Péclet number, $\mathrm {P_e}=\Vert \mathbf {u}\Vert L / \Vert \mathbb {D}\Vert$, where L is the channel length. For miscible displacements in an petroleum reservoir, the viscosity of the fluid mixture may depend on the concentration of the injected fluid through a nonlinear function, the quarter-power viscosity law [30]

$$\begin{aligned} \nu (c) = \nu _{res} [1 - c + M^{ \frac{1}{4} } c]^{-4}, \quad c \in [0,1] \end{aligned}$$

(4)

where M defined in Eq. (1) denotes the mobility ratio. From the Eq. (4) we can observe that, for $M \ne 1$ the subsystems (2) and (3) become tightly coupled.

In order to generate equivalence between transport problem (3) and Darcy problem (2), we rewrite the transport problem in a mixed form including the diffusive flux $\varvec{\sigma } = -\mathbb {D} \nabla c$, which gives rise to the following problem

$$\begin{aligned} \begin{array}{ll} \varvec{\sigma }+ \mathbb {D} \nabla c=0 &{}\text {on} \quad \varOmega \times (0, T], \\ \phi \dfrac{\partial c}{\partial t}+\mathbf {u} \cdot \nabla c + {\text {div}}(\varvec{\sigma })=g &{}\text {on} \quad \varOmega \times (0, T], \end{array} \end{aligned}$$

(5)

with boundary and initial conditions given in (3).

3 Notations and Definitions

To introduce the equivalent stable hybrid formulation for Darcy (2) and Transport (3) systems, we first recall some notations and definitions. Therefore, let $H^m(\varOmega )$ the usual Sobolev space equipped with norm $\Vert \cdot \Vert _{m,\varOmega } = \Vert \cdot \Vert _{m}$ and semi-norm $|\cdot |_{m,\varOmega } = |\cdot |_{m}$, with $m\ge 0$. For $m=0$, we present $L^2(\varOmega )=H^0(\varOmega )$ as the space of square integrable functions and $H^1_0(\varOmega )$ the subspace of functions in $H^1(\varOmega )$ with zero trace on $\partial \varOmega$. In additional, we also define the Hilbert space associated to the divergence operator

$$\begin{aligned} {H}(\mathrm {div},\varOmega ) = {H}(\mathrm {div}) = \{\mathbf {w}\in [L^2(\varOmega )]^d, \, \mathrm {div}\mathbf {w} \in L^2(\varOmega )\}, \end{aligned}$$

with norm

$$\begin{aligned} \Vert \mathbf {w}\Vert _{{H}(\mathrm {div})}^2 = \Vert \mathbf {w}\Vert _{0}^2 + \Vert \mathbf {w}\Vert _{0}^2. \end{aligned}$$

Let $\mathcal {T}_h$ be a regular finite element partition of the domain $\varOmega$, defined by

$$\begin{aligned} \mathcal {T}_h = \{K\} := \text {an union of all elements } K \end{aligned}$$

and let

$$\begin{aligned} \mathcal {E}_h = \{e; e \text { is an edge/face of } K \text { for all } K\in \mathcal {T}_h\} \end{aligned}$$

denotes the set of all edges/faces e of all elements K,

$$\begin{aligned} \mathcal {E}_h^0 = \{e \in \mathcal {E}_h; e \text { is an interior edge/face}\} \end{aligned}$$

is the set of interior egdes/faces, and

$$\begin{aligned} \mathcal {E}^{\partial }_h = \{e \in \mathcal {E}_h; e \subset \partial \varOmega \} \end{aligned}$$

the set of edges/faces of $\mathcal {E}_h$ on the boundary of $\varOmega$. We assume that the domain $\varOmega$ is a polygonal and $\mathcal {T}_h$ is a regular partition of $\varOmega$. Thus, there exists $c > 0$ such that $h\le c h_e$, where $h_e$ is the diameter of the edge/face $e \in \partial K$ and h, the mesh parameter, is the element diameter. For each element K we associate a unit outward normal vector $\mathbf {n}_K$.

The $\mathcal {RT}_k$ spaces [25] are constructed by mapping polynomials defined on the reference element $\hat{K}=[-1,1]^2$ to each element K of the mesh $\mathcal {T}_h$. We denote by $\mathbf {F}_{K}:\hat{K}\rightarrow K$ the invertible, bilinear map of the two domains in $\mathbb {R}^d$, $d=2,3$. A scalar-valued function $\hat{\varphi }$ on $\hat{K}$ transforms to a function $\varphi =P^0_{K}\hat{\varphi }$ on K by the composition

$$\begin{aligned} \varphi (\mathbf {x})=(P^0_{K}\hat{\varphi })(\mathbf {x})=\hat{\varphi }(\hat{\mathbf {x}}), \end{aligned}$$

(6)

with $\mathbf {x}=\mathbf {F}_{K}(\hat{\mathbf {x}})$. A vector-valued function $\hat{\mathbf {\varphi }}$ on $\hat{K}$ transforms to a function $\mathbf {\varphi }=P^1_{K}\mathbf {\hat{\varphi }}$ on K via the Piola transform

$$\begin{aligned} \mathbf {\varphi }(\mathbf {x})=(P^1_{K}\hat{\mathbf {\varphi }})(\mathbf {x}) =\frac{1}{J_{K}(\hat{\mathbf {x}})}[\text {D}\mathbf {F}_{K}(\hat{\mathbf {x}})]\hat{\mathbf {\varphi }}(\hat{\mathbf {x}}), \end{aligned}$$

(7)

where $\text {D}\mathbf {F}_{K}(\hat{\mathbf {x}})$ is the Jacobian matrix of the mapping $\mathbf {F}_{K}$ and $J_K(\hat{\mathbf {x}})$ is its determinant.

The (discontinuous) $\mathcal {RT}_k$ space of index k, $\mathcal {U}_h^k\times \mathcal {P}_h^k$, is defined to be

$$\begin{aligned} \mathcal {U}_h^k&= \left\{ \mathbf{v} _h \in [L^2(\varOmega )]^2; \mathbf{v} _h|_K \in P_K^1(\mathbb {Q}_{k+1, k}(\hat{K}) \times \mathbb {Q}_{k, k+1}(\hat{K})), \forall K \in \mathcal {T}_h \right\} , \end{aligned}$$

(8)

$$\begin{aligned} \mathcal {P}_h^k&= \left\{ q_h \in L^2(\varOmega ); q_h|_{K} \in P_K^0(\mathbb {Q}_{k}(\hat{K})), \forall K \in \mathcal {T}_h \right\} , \end{aligned}$$

(9)

where $\mathbb {Q}_{i,j}$ denotes the set of polynomial functions of degree up to i in x and up to j in y. We also define the following sets of functions on the mesh skeleton:

$$\begin{aligned} \mathcal {M}_h^k&= \left\{ \mu \in L^2(\mathcal {E}_h);~ \mu |_e \in p_k(e), ~\forall e \in \mathcal {E}_h^0,~ \mu |_e = \overline{\mu }, \forall e \in \mathcal {E}_h^\partial \cap \partial \varOmega _D \right\} , \end{aligned}$$

(10)

$$\begin{aligned} \bar{\mathcal {M}}_h^k&= \left\{ \mu \in L^2(\mathcal {E}_h);~ \mu |_e \in p_k(e), ~\forall e \in \mathcal {E}_h^0,~ \mu |_e = 0, \forall e \in \mathcal {E}_h^\partial \cap \partial \varOmega _D \right\} , \end{aligned}$$

(11)

where $p_k(e)$ denotes the set of polynomial functions of degree up to k on e, and $\overline{\mu }$ is the Dirichlet boundary condition function associated to $\overline{p}$ in the Darcy problem and $\overline{c}$ in the transport equation.

4 Equivalent Hybrid Mixed Method

From the definitions presented in the previous section, we can write the following hybrid mixed formulation for the Darcy and Transport systems. For this, we define the product spaces $\mathbf {V}_h^k = \mathcal {U}_h^k \times \mathcal {P}_{h}^{k} \times \mathcal {M}_h^k$ and $\bar{\mathbf {V}}_h^k = \mathcal {U}_h^k \times \mathcal {P}_{h}^{k} \times \bar{\mathcal {M}}_h^k$ and the variable sets

$$ \mathbf {X}^D_h = [\mathbf {u}_{h}, p_{h}, \lambda _h^p] \in \mathbf {V}_h^k \qquad \text {and} \qquad \mathbf {X}^T_h = [\varvec{\sigma }, c_{h}, \lambda _h^c] \in \mathbf {V}_h^k $$

concerning the variables related to the Darcy problem (2) and the variables related to transport problem respectively, we can show the following semi-discrete formulation to the Darcy problem

Given $c_h$, find $\mathbf {X}^D_h \in \mathbf {V}_h^k$, such that

$$\begin{aligned} A(\mathbf {X}^D_h, \mathbf {Y}_h) = F (\mathbf {Y}_h), \qquad \forall \mathbf {Y}_h \in \bar{\mathbf {V}}_h^k. \end{aligned}$$

(12)

On the other hand, the transport problem can be formulated by

Given $\mathbf {u}_h$, find $\mathbf {X}^T_h \in \mathbf {V}_h^k$, such that $\mathbf {Y}_h \in \bar{\mathbf {V}}_h^k$

$$\begin{aligned} \phi \dfrac{\partial c_h}{\partial t} + B(\mathbf {X}^T_h, \mathbf {Y}_h) = G(\mathbf {Y}_h), \qquad \forall \mathbf {Y}_h \in \bar{\mathbf {V}}_h^k. \end{aligned}$$

(13)

In this context, we define the following generalized bilinear form for Darcy and transport problems

$$\begin{aligned} A\left( \mathbf {X}_h, \mathbf {Y}_h\right) = \sum _{K \in \mathcal {T}_{h}}\bigg [&\int _{K} \mathbb {C} \mathbf {w}_{h} \cdot \mathbf {v}_{h} d \mathbf {x}-\int _{K} s_{h} \nabla \cdot \mathbf {v}_{h} d \mathbf {x} + \int _{\partial K} \lambda _{h}\left( \mathbf {w}_{h} \cdot \mathbf {n}_K \right) ds \nonumber \\&+\int _{\partial K} \mu _{h}\left( \mathbf {w}_{h} \cdot \mathbf {n}_K \right) ds -\int _{K} q_{h} \nabla \cdot \mathbf {w}_{h} d \mathbf {x} \bigg ], \end{aligned}$$

(14)

where $\mathbf {X}_h = [\mathbf {w}_{h}, s_h, \lambda _h]$ and $\mathbf {Y}_h = [\mathbf {v}_{h}, q_h, \mu _h]$. Taking $\mathbf {X}_h = \mathbf {X}_h^D$ and $\mathbb {C} = \mathbb {K}^{-1}$ this form can be adapted to Darcy’s problem and to $\mathbf {X}_h = \mathbf {X}_h^T$ and $\mathbb {C} = \mathbb {D}^{-1}$ can be adapted to the transport problem. The source term of (12) and (13) is obtained directly from the multiplication of the respective functions f and g by $q_h$ with a negative sign.

The convective term of the transport equation is stabilized by an upwind strategy proposed by [12], which can be defined as

$$\begin{aligned} \displaystyle A_{conv}(\mathbf {X}^T_h, \mathbf {Y}_h) = \int _{K} c_{h} \mathbf {u}_h \cdot \nabla q_{h} \mathrm {d} \mathbf {x} +\int _{\partial K} \mathbf {u}_h \cdot \mathbf {n}_K \left\{ \lambda _{h}^c / c_{h}\right\} ( \mu _{h}-q_{h}) ds \end{aligned}$$

(15)

with

$$\begin{aligned} \{\lambda / c\}:=\left\{ \begin{array}{ll}{\lambda ,} &{} {e \subset \partial K^{\mathrm {in}}} \\ {c,} &{} {e \subset \partial K^{\mathrm {out}}}\end{array}\right. , \end{aligned}$$

where the element outflow boundary is defined as $\partial K^{\text{ out }}:=\{e \in \partial K : \mathbf {u} \cdot \mathbf {n}_K>0\}$ and the element inflow boundary as $\partial K^{\mathrm {in}}= \partial K \backslash \partial K^{\mathrm {out}}$.

Therefore, the form $A(\cdot ,\cdot )$ of the formulation (12) is given by the bilinear form $A\left( \cdot , \cdot \right)$ defined in (14) taking $\mathbf {X}_h = \mathbf {X}_h^D$ and $\mathbb {C} = \mathbb {K}^{-1}$. On the other hand, the compact bilinear form $B(\cdot ,\cdot )$ of the formulation (13) is defined as

$$\begin{aligned} B(\mathbf {X}^T_h, \mathbf {Y}_h) = A\left( \mathbf {X}_h^T, \mathbf {Y}_h\right) + A_{conv} \left( \mathbf {X}_h^T, \mathbf {Y}_h\right) , \end{aligned}$$

(16)

with $\mathbb {C} = \mathbb {D}^{-1}$.

It is important to emphasize that according to the numerical analysis of the presented hybrid formulations, using Raviart-Thomas spaces, a priori error estimates in the $L^2(\varOmega )$-norm ensures optimal convergence rates for the velocity, pressure and concentration. For more details see [2, 12].

4.1 Time Discretization and Resolution Algorithm

Setting the time step $\varDelta t > 0$, such that $N=T/\varDelta t$ and $t_n = n \varDelta t$ with $n=0,1,2, ..., N$ and the time interval $I_h = \{0 = t_0, t_1, ..., t_N = T\}$ which defines a partition of $I = (0,T]$, we can discretize the time derivative term of the semi-discrete formulation (13) employing implicit backward finite difference scheme as follows

Given $\mathbf {u}_h^{n+1}$ and $c_h^{0}$, find $\mathbf {X}^{T,n+1}_h \in \mathbf {V}_h^k\times I_h$, such that for all $\mathbf {Y}_h \in \bar{\mathbf {V}}_h^k$

$$\begin{aligned} \phi \dfrac{c_h^{n+1} - c_h^{n}}{\varDelta t} + B(\mathbf {X}^{T,n+1}_h, \mathbf {Y}_h) = G(\mathbf {Y}_h), \quad \text {com} \quad n=0,1,2, ..., N \end{aligned}$$

(17)

The resolution methodology is focused on decoupling Darcy and transport problems. Thus, given the initial condition $c_h^{0}$, the viscosity is evaluated using the Eq. (4) and used to solve the Darcy problem (12), once the velocity $\mathbf {u}_h$ is computed the transport problem is solved using the formulation (17). This resolution algorithm is repeated until it reaches the final time T.

To reduce the computational cost of the problem resolution at each time step, the static condensation technique is employed. Thus, element-level problems are condensed in favor of the Lagrange multiplier, generating a global system with degrees of freedom associated with the multiplier only, which in this case is a scalar. Also, the new system of equations is positive-definite, allowing for using simpler and more robust solvers. After a global system resolution, the variables $\mathbf {u}_h$, $p_h$, ${\varvec{\sigma }}_h$ and $c_h$ are retrieved in the element level [15, 16]. In this work, the deal.II library [4] is used to solve these problems.

5 Numerical Results

In this section, the proposed hybrid method (12)–(17) is tested through convergence studies and numerical simulations that are compared with experimental results presented by Malhotra et al. [18] in adverse mobility rate scenarios.

5.1 Convergence Study

Here, we study the convergence rates of the proposed hybrid mixed formulations in a domain $\varOmega =[0,1]^{2}$ and in the time interval $I = [0,0.2]$, adopting the sources $f = 0$ e $g = 2 \pi ^{2} \sin (\pi (x+y-2 t))$ and the parameters $\mathbb {G} = \mathbb {I}$, $\nu (c)=0.2 c - 0.5$, $\alpha _{mol} = 0.01$, $\alpha _l = \alpha _t = 0$ and $\phi =1.0$, is possible to derive the following analytical solution [22]

$$\begin{aligned} \mathbf {u} = [1,1]^T ,\quad p = \frac{2}{10\pi } \cos (\pi (x+y-2 t))+\dfrac{1}{2}(x+y) ,\quad c = \sin (\pi (x+y-2 t)). \end{aligned}$$

(18)

The initial and boundary conditions are determined by the evaluation of the exact solution on time $t=0$ and on the boundary of the domain $\partial \varOmega _D$, respectively.

For the h-convergence study, we adopt meshes with 4, 16, 64, 256, 1024, 4096 quadrilateral elements with the same polynomial order $k=0,1,2$ for the Lagrange multiplier and the primal variables (velocity, pressure, diffusive flux and concentration) and the time step $\varDelta t = h^{k+1}$ to reduce the effects of the error associated to the time discretization. Hence, the spatial error governs the overall error. The results can be seen in Table 1, where it is possible to observe optimal convergence rates for all variables, i.e., order $k+1$.

Table 1. h-Convergence study from the approximations $c_h$, $\mathbf {u}_h$ and $p_h$.

Fig. 1.

Hele-Shaw cell experiment (left) developed by Malhotra et al. [18] compared with the approximate solution (right) at time $t=22,48,63,68,80,90$ s, adopting $1500 \times 100$ elements, $M = 50$, $\mathrm {P_e} = 10 234$ and time step $\varDelta t = 0.1$ s.

5.2 Numerical Simulations in Hele-Shaw Cells

The following results are performed in a Hele-Shaw cell where the flow channel is 84 cm long and 5 cm wide, as described in the work of Malhotra et al. [18]. In this experimental work, Malhotra and collaborators injecting water with dye into glycerol solutions to quantify the growth of the mixing zone in miscible viscous fingering.

In this example, dyed water was injected into $79\%$ glycerol solution at a rate of 4.69 ml/min with a mobility ratio $M=50$ and a Péclet number $\mathrm {P_e} = 10 234$ gives rise to the profiles presented in the left side of the Fig. 1. In this context, adopting the same experimental data in a mesh of $1500 \times 100$ elements with polynomial order $k=1$, time step $\varDelta t = 0.1$ s, $\mathbb {G} = \mathbb {I}$, $\alpha _{mol} = 1.53\times 10^{-7}$ m$^2$/s, $\alpha _{l} = \alpha _{mol}$ and $\alpha _{t} = 1.53\times 10^{-8}$ m$^2$/s and $\phi =1.0$, we develop numerical simulations to compare with experimental results at time $t=22,48,63,68,80,90$ s, as can be seen in the right side of the Fig. 1.

As can be seen in Fig. 1, the numerical results show a good performance of the proposed method in capturing both the formation of viscous fingers and the propagation of the water front. Moreover, the upwind scheme employed is capable of stabilizing the convection-dominated regime caused by the reduction of the resident fluid viscosity adopted in the experiment that generates a Péclet number of the $10^4$ order.

6 Conclusions

In this work, we proposed an equivalent stable hybrid mixed method adopting non-conforming Raviart-Thomas spaces for the Darcy and the Transport systems to solve miscible displacements with adverse mobility ratio. For the convection-diffusion equation, an upwind scheme was employed to stabilizing the convection-dominated regimes, and the implicit backward Euler approach was used to the time discretization. The continuity was weakly imposed by the Lagrange multiplier associated with the pressure field for the Darcy problem and the concentration for the Transport mixed problem. This approach gives rise to a positive-definite global matrix with reduced computational cost compared to classical conforming Raviart-Thomas formulations. To solve this coupled problem, we employed a staggered approach to decoupling the systems using the same time scale for both problems. The numerical studies confirmed optimal convergence rates for velocity, pressure and concentration. In addition, simulations adopting an adverse mobility ratio, comparing the approximate solution with experimental results in Hele-Shaw cell, demonstrated that the proposed hybrid formulations are capable of capturing the formation and propagation of the viscous fingering even in cases of high Péclet number.