Convergence of a finite element method for degenerate two-phase flow in porous media

Let $\Omega\subset{\rm I\!R}^{d}$, $d=2$ or $3$, be a bounded connected Lipschitz domain with boundary $\partial\Omega$ and unit exterior normal ${\boldsymbol{n}}$, and let $T$ be a final time. With the last relation in (1), $s_{w}$ is the only unknown saturation; so we set $s=s_{w}$, and rewrite (1) almost everywhere in $\Omega\times]0,T[$ as

	$\displaystyle\partial_{t}(\varphi s)-\nabla\cdot(\eta_{w}\nabla p_{w})=f_{w}(s_{\mathrm{in}})\bar{q}-f_{w}(s)\underline{q}$		(2)
	$\displaystyle-\partial_{t}(\varphi s)-\nabla\cdot(\eta_{o}\nabla p_{o})=f_{o}(s_{\mathrm{in}})\bar{q}-f_{o}(s)\underline{q},$		(3)

complemented by a natural boundary condition almost everywhere on $\partial\Omega\times]0,T[$

$$\eta_{w}\nabla p_{w}\cdot{\bf n}=0,\quad\eta_{o}\nabla p_{o}\cdot{\bf n}=0,$$

(4)

and an initial condition almost everywhere in $\Omega$

$$s_{w}(\cdot,0)=s_{w}^{0},\quad 0\leq s_{w}^{0}\leq 1.$$

(5)

The fractional flows are related to the mobilities by

$$f_{w}=\frac{\eta_{w}}{\eta_{w}+\eta_{o}},\quad f_{o}=1-f_{w}.$$

(6)

Recall that the phase saturations sum up to 1 and the phase pressures are related by

$$p_{c}(s_{w})=p_{o}-p_{w}.$$

(7)

The first part of this work, up to Section 5.7, is done under the following basic assumptions:

Assumptions

• •

  The porosity $\varphi$ is piecewise constant in space, independent of time, positive, bounded and uniformly bounded away from zero.

• •

  The mobility of the wetting phase $\eta_{w}\geq 0$ is continuous and increasing. The mobility of the non-wetting phase $\eta_{o}\geq 0$ is continuous and decreasing. This implies that the function $f_{w}$ is increasing and the function $f_{o}$ is decreasing.

• •

  There is a positive constant $\eta_{\ast}$ such that

  $$\eta_{w}(s)+\eta_{o}(s)\geq\eta_{\ast},\quad\forall s\in[0,1].$$

  (8)

• •

  The capillary pressure $p_{c}$ is a continuous, strictly decreasing function in $W^{1,1}(0,1)$.

• •

  The flow rates at the injection and production wells, $\bar{q},\underline{q}\in L^{2}(\Omega\times]0,T[)$ satisfy

  $$\bar{q}\geq 0,\quad\underline{q}\geq 0,\quad\int_{\Omega}\bar{q}=\int_{\Omega}\underline{q}.$$

  (9)

• •

  The prescribed input saturation $s_{\mathrm{in}}$ satisfies almost everywhere in $\Omega\times]0,T[$

  $$0\leq s_{\mathrm{in}}\leq 1\color[rgb]{0,0,0}.$$

  (10)

Since $p_{c}$, $\eta_{\alpha}$, $f_{\alpha}$, $\alpha=w,o$ are bounded above and below, it is convenient to extend them continuously by constants to ${\rm I\!R}$.

Although the numerical scheme studied below does not discretize the global pressure, following [16], its convergence proof uses a number of auxiliary functions related to the global pressure. First, we introduce the primitive $g_{c}$ of $p_{c}$,

$$\forall x\in[0,1],\quad g_{c}(x)=\int_{x}^{1}p_{c}(s)ds.$$

(11)

Since $p_{c}$ is a continuous function on $[0,1]$, the function $g_{c}$ belongs to $\mathcal{C}^{1}([0,1])$. Next, we introduce the auxiliary pressures $p_{wg}$, $p_{wo}$, and $g$,

$$\forall x\in[0,1],\quad p_{wg}(x)=\int_{0}^{x}f_{o}(s)p_{c}^{\prime}(s)ds,\quad p_{og}(x)=\int_{0}^{x}f_{w}(s)p_{c}^{\prime}(s)ds,$$

(12)

$$\forall x\in[0,1],\quad g(x)=-\int_{0}^{x}\frac{\eta_{w}(s)\eta_{o}(s)}{\eta_{w}(s)+\eta_{o}(s)}p^{\prime}_{c}(s)ds.$$

(13)

Owing to (6),

$$\forall x\in[0,1],\quad p_{wg}(x)+p_{og}(x)=\int_{0}^{x}p_{c}^{\prime}(s)ds=p_{c}(x)-p_{c}(0).$$

(14)

Moreover, the derivative of $g$ satisfies formally the identities

$$\forall x\in[0,1],\quad\eta_{\alpha\color[rgb]{0,0,0}}(x)p^{\prime}_{\alpha g}(x)+g^{\prime}(x)=0,\quad\alpha=w,o.$$

(15)

By multiplying (2) and (3) with a smooth function $v$, say $v\in{\mathcal{C}}^{1}(\Omega\times[0,T])$ that vanishes at $t=T$, applying Green’s formula in time and space, and using the boundary and initial conditions (4) and (5), we formally derive a weak variational formulation

	$\displaystyle-\int_{0}^{T}\int_{\Omega}\varphi\,s\,\partial_{t}v+\int_{0}^{T}\int_{\Omega}\eta_{w}\nabla\,p_{w}\cdot\nabla\,v=$	$\displaystyle\int_{\Omega}\varphi\,s^{0}v(0)+\int_{0}^{T}\int_{\Omega}\big{(}f_{w}(s_{\mathrm{in}})\bar{q}-f_{w}(s)\underline{q}\big{)}v,$
	$\displaystyle\int_{0}^{T}\int_{\Omega}\varphi\,s\,\partial_{t}v+\int_{0}^{T}\int_{\Omega}\eta_{o}\nabla\,p_{o}\cdot\nabla\,v=$	$\displaystyle-\int_{\Omega}\varphi\,s^{0}v(0)+\int_{0}^{T}\int_{\Omega}\big{(}f_{o}(s_{\mathrm{in}})\bar{q}-f_{o}(s)\underline{q}\big{)}v.$

But in general, the pressures are not sufficiently smooth to make this formulation meaningful and following [7], by using (15), it is rewritten in terms of the artificial pressures,

$$\begin{split}-\int_{0}^{T}\int_{\Omega}\varphi\,s\,\partial_{t}v+\int_{0}^{T}\int_{\Omega}\big{(}\eta_{w}&\nabla(p_{w}+p_{wg}(s))+\nabla\,g(s)\big{)}\cdot\nabla\,v=\int_{\Omega}\varphi\,s^{0}v(0)\\ &+\int_{0}^{T}\int_{\Omega}\big{(}f_{w}(s_{\mathrm{in}})\bar{q}-f_{w}(s)\underline{q}\big{)}v,\\ \int_{0}^{T}\int_{\Omega}\varphi\,s\,\partial_{t}v+\int_{0}^{T}\int_{\Omega}\big{(}\eta_{o}&\nabla(p_{o}-p_{og}(s))-\nabla\,g(s)\big{)}\cdot\nabla\,v=-\int_{\Omega}\varphi\,s^{0}v(0)\\ &+\int_{0}^{T}\int_{\Omega}\big{(}f_{o}(s_{\mathrm{in}})\bar{q}-f_{o}(s)\underline{q}\big{)}v.\end{split}$$

(16)

The two formulas coincide when the pressures are slightly more regular. With the above assumptions, problem (16) has been analyzed in reference [1], where it is shown that it has a solution $s$ in $L^{\infty}(\Omega\times]0,T[)$ with $g(s)$ in $L^{2}(0,T;H^{1}(\Omega))$, $p_{\alpha}$, $\alpha=w,o$, in $L^{2}(\Omega\times]0,T[)$ with both $p_{w}+p_{wg}(s)$ and $p_{o}-p_{og}(s)$ in $L^{2}(0,T;H^{1}(\Omega))$.

The mesh ${\mathcal{T}}_{h}$ is a regular family of simplices $K$, with a constraint on the angle that will be used to enforce the maximum principle: each angle is not larger than $\pi/2$, see [6]. This is easily constructed in 2D. In 3D, since we only investigate convergence we can embed the domain in a triangulated box. Moreover, since the porosity $\varphi$ is a piecewise constant, to simplify we also assume that the mesh is such that $\varphi$ is a constant per element. The parameter $h$ denotes the mesh size i.e., the maximum diameter of the simplices. On this mesh, we consider the standard finite element space of order one

$$X_{h}=\{v_{h}\in{\mathcal{C}}^{0}(\bar{\Omega})\,;\,\forall K\in{\mathcal{T}}_{h},v_{h}|_{K}\in{\mathbb{P}}_{1}\}.$$

(17)

Thus the dimension of $X_{h}$ is the number of nodes, say $M$, of ${\mathcal{T}}_{h}$. Let $\phi_{i}$ be the Lagrange basis function, that is piecewise linear, and takes the value $1$ at node $i$ and the value 0 elsewhere. As usual, the Lagrange interpolation operator $I_{h}\in{\mathcal{L}}({\mathcal{C}}^{0}(\bar{\Omega});X_{h})$ is defined by

$$\forall v\in{\mathcal{C}}^{0}(\bar{\Omega}),\quad I_{h}(v)=\sum_{i=1}^{M}v_{i}\phi_{i},$$

(18)

where $v_{i}$ is the value of $v$ at the node of index $i$. It is easy to see that under the mesh condition, we have

$$\forall K,\quad\int_{K}\nabla\phi_{i}\cdot\nabla\phi_{j}\leq 0,\quad\forall i\neq j.$$

(19)

For a given node $i$, we denote by $\Delta_{i}$ the union of elements sharing the node $i$ and by $\mathcal{N}(i)$ the set of indices of all the nodes in $\Delta_{i}$. In the spirit of [19], we define

$$c_{ij}=\int_{\Delta_{i}\cap\Delta_{j}}|\nabla\phi_{i}\cdot\nabla\phi_{j}|,\quad\forall i,j.$$

(20)

Recall that the trapezoidal rule on a triangle or a tetrahedron $K$ is

$$\int_{K}f\approx\frac{1}{d+1}|K|\sum_{\ell=1}^{d+1}f_{i_{\ell}},$$

where $f_{i_{\ell}}$ is the value of the function $f$ at the $\ell^{th}$ node (vertex), with global number $i_{\ell}$, of $K$. For any region ${\mathcal{O}}$, the notation $|\mathcal{O}|$ means the measure (volume) of $\mathcal{O}$.

We define

$$m_{i}=\frac{1}{d+1}\sum_{K\in\Delta_{i}}|K|=\frac{1}{d+1}|\Delta_{i}|,$$

and taking into account the porosity $\varphi$, we define more generally

$$\tilde{m}_{i}(\varphi)=\frac{1}{d+1}\sum_{K\in\Delta_{i}}\varphi|_{K}|K|,$$

so that $m_{i}=\tilde{m}_{i}(1)$. It is well-known that the trapezoidal rule defines a norm on $X_{h}$, $\|\cdot\|_{h}$, uniformly equivalent to the $L^{2}$ norm. Let $U_{h}\in X_{h}$ and write

$$U_{h}=\sum_{i=1}^{M}U^{i}\phi_{i}.$$

The discrete $L^{2}$ norm associated with the trapezoidal rule is

$$\|U_{h}\|_{h}=\left(\sum_{i=1}^{M}m_{i}|U^{i}|^{2}\right)^{\frac{1}{2}}.$$

There exist positive constants $\underline{C}$ and $\overline{C}$, independent of $h$ and $M$, such that

$$\forall U_{h}\in X_{h},\quad\underline{C}\,\|U_{h}\|^{2}_{L^{2}({\Omega})}\leq\|U_{h}\|^{2}_{h}\leq\overline{C}\,\|U_{h}\|^{2}_{L^{2}({\Omega})}.$$

(21)

This is also true for other piecewise polynomial functions, but with possibly different constants. The scalar product associated with this norm is denoted by $(\cdot,\cdot)_{h}$,

$$\forall U_{h},V_{h}\in X_{h},\quad(U_{h},V_{h})_{h}=\sum_{i=1}^{M}m_{i}U^{i}V^{i}.$$

(22)

By analogy, we introduce the notation

$$\forall U_{h},V_{h}\in X_{h},\quad(U_{h},V_{h})^{\varphi}_{h}=\sum_{i=1}^{M}\tilde{m}_{i}(\varphi)U^{i}V^{i}.$$

(23)

The assumptions on the porosity $\varphi$ imply that (23) defines a weighted scalar product associated with the weighted norm $\|\cdot\|_{h}^{\varphi}$,

$$\forall U_{h}\in X_{h},\quad\|U_{h}\|_{h}^{\varphi}=\big{(}(U_{h},U_{h})^{\varphi}_{h}\big{)}^{\frac{1}{2}},$$

that satisfies the analogue of (21), with the same constants $\underline{C}$ and $\overline{C}$,

$$\forall U_{h}\in X_{h},\quad\underline{C}\,(\min_{\Omega}\varphi)\,\|U_{h}\|^{2}_{L^{2}({\Omega})}\leq\big{(}\|U_{h}\|_{h}^{\varphi}\big{)}^{2}\leq\overline{C}\,(\max_{\Omega}\varphi)\,\|U_{h}\|^{2}_{L^{2}({\Omega})}.$$

(24)

While discretizing the time derivative is fairly straightforward, discretizing the space derivatives is more delicate because we need a scheme that is consistent and satisfies the maximum principle for the saturation. For the moment, we freeze the time variable and focus on consistency in space. First, we recall a standard property of functions of $X_{h}$ on meshes satisfying (19).

Note that $c_{ij}$ vanishes when $j\notin{\mathcal{N}}(i)$. Therefore, when there is no ambiguity it is convenient to write the above double sums on $i$ and $j$ with $i$ and $j$ running from $1$ to $M$.

As an immediate consequence of Proposition 1, we have, by taking $V_{h}=U_{h}$,

$$\forall U_{h}\in X_{h},\quad\|\nabla\,U_{h}\|_{L^{2}({\Omega})}=\frac{1}{\sqrt{2}}\Big{(}\sum_{i,j=1}^{M}c_{ij}|U^{j}-U^{i}|^{2}\Big{)}^{\frac{1}{2}}.$$

(26)

Note that $d_{ij}=d_{ji}$ owing to (33). The first consequence of Proposition 2 is that the right-hand side of (29) is a consistent approximation of $(w,\nabla\,u\cdot\nabla\,v)$.

Now, if $w$ is in $W^{1,\infty}(\Omega)$, then again, standard finite element approximation shows that there exists a constant $C$, independent of $h$, $K\subset\Delta_{i}\cap\Delta_{j}$, and $w$, such that

$$\big{\|}w_{K}-w\big{\|}_{L^{\infty}(K)}\leq C\,h\,|w|_{W^{1,\infty}(K)}\leq C\,h\,|w|_{W^{1,\infty}(\Omega)}.$$

(35)

As a consequence, we will show that in the error formula (34), the average $w_{K}$ can be replaced by any value of $w$ in $K$. Since all $K$ in $\Delta_{i}\cap\Delta_{j}$ share the edge, say $e_{ij}$, whose end points are the nodes with indices $i$ and $j$, then we can pick the value of $w$ at any point, say $\tilde{W}^{i,j}$, of $e_{ij}$. At this stage, we choose this value freely, but we prescribe that it be symmetrical with respect to $i$ and $j$, i.e.,

$$\tilde{W}^{i,j}=\tilde{W}^{j,i}.$$

(36)

Then we have the following approximation result.

The above considerations show that

$$-\sum_{i,j=1}^{M}U^{i}c_{ij}\tilde{W}^{i,j}\big{(}V^{j}-V^{i}\big{)}\ \mbox{is a consistent approximation of}\ \int_{\Omega}w\nabla\,u\cdot\nabla\,v,$$

for any symmetric choice of $\tilde{W}^{i,j}$ in $e_{ij}$, the common edge of $\Delta_{i}\cap\Delta_{j}$. This will lead to the upwinded space discretization in the next subsection, see also [23]. Furthermore, for all real numbers $V^{i}$ and $\tilde{W}^{i,j}$ satisfying (36), $1\leq i,j\leq M$, the symmetry of $c_{ij}$ and anti-symmetry of $V^{j}-V^{i}$ imply

$$\sum_{i,j=1}^{M}c_{ij}\tilde{W}^{i,j}(V^{j}-V^{i})=0.$$

(40)

Let $\tau=\frac{T}{N}$ be the time step, $t_{n}=n\tau$, the discrete times, $0\leq n\leq N$. Regarding time, we shall use the standard $L^{2}$ projection $\rho_{\tau}$ defined on $]t_{n-1},t_{n}]$, for any function $f$ in $L^{1}(0,T)$, by

$$\rho_{\tau}(f)^{n}:=\rho_{\tau}(f)|_{]t_{n-1},t_{n}]}:=\frac{1}{\tau}\int_{t_{n-1\color[rgb]{0,0,0}}}^{t_{n\color[rgb]{0,0,0}}}f.$$

(41)

Regarding space, we shall use a standard element-by-element $L^{2}$ projection $\rho_{h}$ as well as a nodal approximation operator $r_{h}$ defined at each node ${\boldsymbol{x}}_{i}$ for any function $g\in L^{1}(\Omega)$ by

$$r_{h}(g)({\boldsymbol{x}}_{i})=\frac{1}{|\Delta_{i}|}\int_{\Delta_{i}}g,\quad 1\leq i\leq M,$$

(42)

and extended to $\Omega$ by $r_{h}(g)\in X_{h}$. The operator $\rho_{h}$ is defined for any $f$ in $L^{1}(\Omega)$ by $\rho_{h}(f)|_{K}=\rho_{K}(f)$ where, in any element $K$,

$$\rho_{K}(f)=\frac{1}{|K|}\int_{K}f.$$

(43)

The initial saturation $s_{w}^{0}$ is approximated by the operator $r_{h}$,

$$S_{h}^{0}=r_{h}(s_{w}^{0}).$$

(44)

The input saturation $s_{\mathrm{in}}$ is approximated in space and time by

$$s_{\mathrm{in,h,\tau}}=\rho_{\tau}(r_{h}(s_{\mathrm{in}})).$$

(45)

Clearly, (10) implies in space and time

$$0\leq s_{\mathrm{in,h,\color[rgb]{1,0,0}\tau\color[rgb]{0,0,0}}}\leq 1.$$

In order to preserve (9), the functions $\bar{q}$ and $\underline{q}$ are approximated by the functions $\bar{q}_{h,\tau}$ and $\underline{q}_{h,\tau}$ defined with $r_{h}$ and corrected as follows:

$$\bar{q}_{h,\tau}=\rho_{\tau}\left(r_{h}(\bar{q})-\frac{1}{|\Omega|}\int_{\Omega}(r_{h}(\bar{q})-\bar{q})\right),\quad\underline{q}_{h,\tau}=\rho_{\tau}\left(r_{h}(\underline{q})-\frac{1}{|\Omega|}\int_{\Omega}(r_{h}(\underline{q})-\underline{q})\right).$$

(46)

Since $\bar{q}_{h,\tau}$ and $\underline{q}_{h,\tau}$ are piecewise linears in space, they are exactly integrated by the trapezoidal rule and we easily derive from (9) and (46) that we have for all $n$,

$$\big{(}\bar{q}_{h}^{n},1\big{)}_{h}=\big{(}\underline{q}_{h}^{n},1\big{)}_{h}.$$

(47)

The set of primary unknowns is the discrete wetting phase saturation and the discrete wetting phase pressure, $S_{h}^{n}$ and $P_{w,h}^{n}$, defined pointwise at time $t_{n}$ by:

$$S_{h}^{n}=\sum_{i=1}^{M}S^{n,i}\phi_{i},\quad P_{w,h}^{n}=\sum_{i=1}^{M}P_{w}^{n,i}\phi_{i},\quad 1\leq n\leq N.$$

Then the discrete non-wetting phase pressure $P_{o,h}^{n}$ defined by

$$P_{o,h}^{n}=\sum_{i=1}^{M}P_{o}^{n,i}\phi_{i},\quad 1\leq n\leq N,$$

is a secondary unknown. The upwind scheme we propose for discretizing (2)–(3) is inspired by the finite volume scheme introduced and analyzed by Eymard al in [16]. For each time step $n$, $1\leq n\leq N$, the lines of the discrete equations are
$$\frac{\tilde{m}_{i}(\varphi)}{\tau}(S^{n,i}-S^{n-1,i})-\sum_{j=1}^{M}c_{ij}\eta_{w}(S^{n,ij}_{w})(P_{w}^{n,j}-P_{w}^{n,i})=m_{i}\left(f_{w}(s_{\mathrm{in}}^{n,i})\bar{q}^{n,i}-f_{w}(S^{n,i})\underline{q}^{n,i}\right),$$ (48) $$-\frac{\tilde{m}_{i}(\varphi)}{\tau}(S^{n,i}-S^{n-1,i})-\sum_{j=1}^{M}c_{ij}\eta_{o}(S^{n,ij}_{o})(P_{o}^{n,j}-P_{o}^{n,i})=m_{i}\left(f_{o}(s_{\mathrm{in}}^{n,i})\bar{q}^{n,i}-f_{o}(S^{n,i})\underline{q}^{n,i}\right),$$ (49) $$P_{o}^{n,i}-P_{w}^{n,i}=p_{c}(S^{n,i}),\quad 1\leq i\leq M,$$ (50) $$\sum_{i=1}^{M}m_{i}P_{w}^{n,i}=0.$$ (51) Here $i$ runs from $1$ to $M-1$ in (48) and from $1$ to $M$ in (49); the upwind values $S^{n,ij}_{w},S^{n,ij}_{o}$ are defined by

$$S_{w}^{n,ij}=\left\{\begin{array}[]{c}S^{n,i}\quad\mbox{if}\quad P_{w}^{n,i}>P_{w}^{n,j}\\ S^{n,j}\quad\mbox{if}\quad P_{w}^{n,i}<P_{w}^{n,j}\\ \max(S^{n,i},S^{n,j})\quad\mbox{if}\quad P_{w}^{n,i}=P_{w}^{n,j}\end{array}\right.$$

(52)

$$S_{o}^{n,ij}=\left\{\begin{array}[]{c}S^{n,i}\quad\mbox{if}\quad P_{o}^{n,i}>P_{o}^{n,j}\\ S^{n,j}\quad\mbox{if}\quad P_{o}^{n,i}<P_{o}^{n,j}\\ \min(S^{n,i},S^{n,j})\quad\mbox{if}\quad P_{o}^{n,i}=P_{o}^{n,j}\end{array}\right.$$

(53)

We observe that

$$S_{w}^{n,ij}=S_{w}^{n,ji},\quad S_{o}^{n,ij}=S_{o}^{n,ji},$$

so that, if we interpret in (48) (respectively, (49)) $\eta_{w}(S^{n,ij}_{w})$ (respectively, $\eta_{o}(S^{n,ij}_{o})$) as $\tilde{W}^{i,j}$, then (36) and hence (40) hold.