A New Unified Stabilized Mixed Finite Element Method of the Stokes-Darcy Coupled Problem: Isotropic Discretization

Houédanou Koffi Wilfrid

Département de Mathématiques, Université d’Abomey-Calavi (UAC), Abomey-Calavi, Benin.

DOI: 10.4236/jamp.2021.97112                

Abstract

In this paper, we develop an a-priori error analysis of a new unified mixed finite element method for the coupling of fluid flow with porous media flow in R^N, N ∈ {2,3}, on isotropic meshes. Flows are governed by the Stokes and Darcy equations, respectively, and the corresponding transmission conditions are given by mass conservation, balance of normal forces, and the Beavers-Joseph-Saffman law. The approach utilizes a modification of the Darcy problem which allows us to apply a variant nonconforming Crouzeix-Raviart finite element to the whole coupled Stokes-Darcy problem. The well-posedness of the finite element scheme and its convergence analysis are derived. Finally, the numerical experiments are presented, which confirm the excellent stability and accuracy of our method.

Keywords

Coupled Stokes and Darcy Flows, Nonconforming Finite Element Method, Crouzeix-Raviart Element

1. Introduction

There are many serious problems currently facing the world in which the coupling between groundwater and surface water is important. These include questions such as predicting how pollution discharges into streams, lakes, and rivers its way into the water supply. This coupling is also important in technological applications involving filtration. We refer to the nice overview [1] and the references therein for its physical background, modeling, and standard numerical methods. One important issue in the modeling of the coupled Darcy-Stokes flow is the treatment of the interface condition, where the Stokes fluid meets the porous medium. In this paper, we only consider the so-called Beavers-Joseph-Saffman condition, which was experimentally derived by Beavers and Joseph in [2], modified by Saffman in [3], and later mathematically justified in [4] [5] [6] [7].

It is well known that the discretization of the velocity and the pressure, for both Stokes and Darcy problems and the coupled of them, has to be made in a compatible way in order to avoid instabilities. Since, usually, stable elements for the free fluid flow cannot been successfully applied to the porous medium flow, most of the finite element formulations developed for the Stokes-Darcy coupled problem are based on appropriate combinations of stable elements for the Stokes equations with stable elements for the Darcy equations. In [4] [6] [8] - [25], and in the references therein, we can find a large list of contributions devoted to numerically approximate the solution of this interaction problem, including conforming and nonconforming methods.

There are a lot of papers considering different finite element spaces in each flow region (see, for example, [21] [26] [27] and the references therein). In contrast to this, other articles use the same finite element spaces in both regions by, in general, introducing some penalizing terms (ref. for examples [10] [14] [22] and the references therein).

In [22], a conforming unified finite element has been proposed for the modified coupled Stokes-Darcy problem in a plane domain, which has simple and straightforward implementations. The authors apply the classical Mini-element to the whole coupled Stokes-Darcy coupled problem. An a-priori error analysis is performed with some numerical tests confirming the convergence rates.

In this article, we propose a modification of the Darcy problem which allows us to apply a variant nonconforming finite element to the whole coupled Stokes-Darcy problem. We use a variant nonconforming Crouzeix-Raviart finite element method that has so many advantages for the velocities and piecewise constant for the pressures in both the Stokes and Darcy regions, and apply a stabilization term penalizing the jumps over the element edges of the piecewise continuous velocities. We prove that the formulation satisfies the discrete inf-sup condition, obtaining as a result optimal accuracy with respect to solution regularity. Numerical experiments are also presented, which confirm the excellent stability and optimal performance of our method. The difference between our paper and the reference [22] is that our discretization is nonconforming in both the Stokes domain and Darcy domain (in $\Omega \subset {ℝ}^N$, $N=2$ or 3). As a result, additional terms are included in the priori error analysis that measures the non-conformity of the method. One essential difficulty in choosing the unified discretization is that, the Stokes side velocity is in $H^1$ while the Darcy side velocity is only in $H\left(\text{div}\right)$. Thus, we introduce a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element space (larger than the space $H_h$ used in [14] ). The choice of $H_h$ [see (34)] is more natural than the one introduced in [14] since the space $H_h$ approximates only $H\left(\text{div}\mathrm{,}{\Omega }_d\right)$ and not ${\left[H^1\left({\Omega }_d\right)\right]}^N$, while our a priori error analysis is only valid in this larger space.

The rest of the paper is organized as follows. In Section 2 we present the modified coupled Stokes-Darcy problem in $\Omega \subset {ℝ}^N$, $N=2$ or 3, notations and the weak formulation. Section 3 is devoted to the finite element discretization and the error estimation.

In Section 4, we present the results of numerical experiments to verify the predicted rates of convergence. Finally, we offer our conclusion and the further works in Section 5.

2. Preliminaries and Notation

2.1. Model Problem

We consider the model of a flow in a bounded domain $\Omega \subset {ℝ}^N$ ( $N=2$ or 3), consisting of a porous medium domain ${\Omega }_d$, where the flow is a Darcy flow, and an open region ${\Omega }_s=\Omega \backslash {\bar{\Omega }}_d$, where the flow is governed by the Stokes equations. The two regions are separated by an interface ${\Gamma }_I=\partial {\Omega }_d\cap \partial {\Omega }_s$. Let ${\Gamma }_l=\partial {\Omega }_l\backslash {\Gamma }_I$, $l=s,d$. Each interface and boundary is assumed to be polygonal ( $N=2$ ) or polyhedral ( $N=3$ ). We denote by $n_s$ (resp. $n_d$ ) the unit outward normal vector along $\partial {\Omega }_s$ (resp. $\partial {\Omega }_d$ ). Note that on the interface ${\Gamma }_I$, we have $n_s=-n_d$. Figure 1 and Figure 2 give a schematic representation of the geometry.

For any function v defined in $\Omega$, since its restriction to ${\Omega }_s$ or to ${\Omega }_d$ could play a different mathematical roles (for instance their traces on ${\Gamma }_I$ ), we will set $v_s=v_{|{\Omega }_s}$ and $v_d=v_{|{\Omega }_d}$.

Figure 1. A sketch of the geometry of the problem (case: $\partial {\Omega }_d\ne {\Gamma }_I$ ).

Figure 2. A sketch of the geometry of the problem (case $\partial {\Omega }_d={\Gamma }_I$ ).

In $\Omega$, we denote by $u$ the fluid velocity and by p the pressure. The motion of the fluid in ${\Omega }_s$ is described by the Stokes equations

$\{\begin{array}{ll}-2\mu \text{div}D\left(u\right)+\nabla p=f\hfill & \text{in}{\Omega }_s,\hfill \\ \text{div}u=g\hfill & \text{in}{\Omega }_s,\hfill \\ u=0\hfill & \text{on}{\Gamma }_s,\hfill \end{array}$ (1)

while in the porous medium ${\Omega }_d$, by Darcy’s law

$\{\begin{array}{ll}\mu K^{-1}u+\nabla p=f\hfill & \text{in}{\Omega }_d,\hfill \\ \text{div}u=g\hfill & \text{in}{\Omega }_d,\hfill \\ u\cdot n_d=0\hfill & \text{on}{\Gamma }_d,\hfill \end{array}$ (2)

Here, $\mu >0$ is the fluid viscosity, $D$ the deformation rate tensor defined by

$D{\left(\psi \right)}_{ij}:=\frac{1}{2}\left(\frac{\partial {\psi }_i}{\partial x_j}+\frac{\partial {\psi }_j}{\partial x_i}\right),1\le i,j\le N,$

and $K$ a symmetric and uniformly positive definite tensor representing the rock permeability and satisfying, for some constants $0<K_{*}\le K^{*}<+\infty$,

$K_{\mathrm{<mo>\; *\; </mo>}}{\xi }^{\text{T}}\xi \le {\xi }^{\text{T}}K\left(x\right)\xi \le K^{\mathrm{<mo>\; *\; </mo>}}{\xi }^{\text{T}}\xi \mathrm{,}\forall x\in {\Omega }_d\mathrm{,}\xi \in {ℝ}^N\mathrm{.}$

$f\in {\left[L^2\left(\Omega \right)\right]}^N$ is a term related to body forces and $g\in L^2\left(\Omega \right)$ a source or sink term satisfying the compatibility condition

${\displaystyle {\int }_{\Omega }}g\left(x\right)\text{d}x=0.$

Finally we consider the following interface conditions on ${\Gamma }_I$ :

$u_s\cdot n_s+u_d\cdot n_d=\mathrm{0,}$ (3)

$p_s-2\mu n_s\cdot D\left(u_s\right)\cdot n_s=p_d\mathrm{,}$ (4)

$\frac{\sqrt{{\kappa }_j}}{{\alpha }_1}2n_s\cdot D\left(u_s\right)\cdot {\tau }_j=-u_s\cdot {\tau }_j\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,}N-1.$ (5)

Here, Equation (3) represents mass conservation, Equation (4) the balance of normal forces, and Equation (5) the Beavers-Joseph-Saffman conditions. Moreover, ${\left\{{\tau }_j\right\}}_{j=1,\cdots ,N-1}$ denotes an orthonormal system of tangent vectors on ${\Gamma }_I$, ${\kappa }_j={\tau }_j\cdot K\cdot {\tau }_j$, and ${\alpha }_1$ is a parameter determined by experimental evidence.

Equations (1) to (5) consist of the model of the coupled Stokes and Darcy flows problem that we will study below.

2.2. New Weak Formulation

We begin this subsection by introducing some useful notations. If W is a bounded domain of ${ℝ}^N$ and m is a non negative integer, the Sobolev space $H^m\left(W\right)=W^{m\mathrm{,2}}\left(W\right)$ is defined in the usual way with the usual norm ${\Vert \cdot \Vert }_{m\mathrm{,}W}$ and semi-norm ${\left|\cdot \right|}_{m\mathrm{,}W}$. In particular, $H^0\left(W\right)=L^2\left(W\right)$ and we write ${\Vert \cdot \Vert }_W$ for ${\Vert \cdot \Vert }_{\mathrm{0,}W}$. Similarly we denote by ${\left(\cdot \mathrm{,}\cdot \right)}_W$ the $L^2\left(W\right)$ ${\left[L^2\left(W\right)\right]}^N$ or ${\left[L^2\left(W\right)\right]}^{N\times N}$ inner product. For shortness if W is equal to $\Omega$, we will drop the index $\Omega$, while for any $m\ge 0$, ${\Vert \cdot \Vert }_{m\mathrm{,}l}={\Vert \cdot \Vert }_{m\mathrm{,}{\Omega }_l}$, ${\left|\cdot \right|}_{m\mathrm{,}l}={\left|\cdot \right|}_{m\mathrm{,}{\Omega }_l}$ and ${\left(\mathrm{.,.}\right)}_l={\left(\mathrm{.,.}\right)}_{{\Omega }_l}$, for $l=s,d$. The space $H_0^m\left(\Omega \right)$ denotes the closure of $C_0^{\infty }\left(\Omega \right)$ in $H^m\left(\Omega \right)$. Let ${\left[H^m\left(\Omega \right)\right]}^N$ be the space of vector valued functions $v=\left(v_1\mathrm{,}\cdots \mathrm{,}{v}_N\right)$ with components $v_i$ in $H^m\left(\Omega \right)$. The norm and the semi-norm on ${\left[H^m\left(\Omega \right)\right]}^N$ are given by

${\Vert v\Vert }_{m,\Omega }:={\left({\displaystyle \underset{i=0}{\overset{N}{\sum }}}{\Vert v_i\Vert }_{m,\Omega }^2\right)}^{1/2}\text{and}{\left|v\right|}_{m,\Omega }:={\left({\displaystyle \underset{i=0}{\overset{N}{\sum }}}{\left|v_i\right|}_{m,\Omega }^2\right)}^{1/2}.$ (6)

For a connected open subset of the boundary $\Gamma \subset \partial {\Omega }_s\cup \partial {\Omega }_d$, we write ${\langle \mathrm{.,.}\rangle }_{\Gamma }$ for the $L^2\left(\Gamma \right)$ inner product (or duality pairing), that is, for scalar valued functions $\lambda$, $\eta$ one defines:

${\langle \lambda \mathrm{,}\eta \rangle }_{\Gamma }\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_{\Gamma }}\lambda \left(s\right)\eta \left(s\right)\text{d}s$ (7)

We also define the special vector-valued functions space

$H\left(\text{div}\mathrm{,}\Omega \right)\mathrm{<mo>\; :\; </mo>}=\left\{v\in {\left[L^2\left(\Omega \right)\right]}^N\mathrm{<mo>\; :\; </mo>}\text{div}v\in L^2\left(\Omega \right)\right\}$ (8)

To give the variational formulation of our coupled problem we define the following two spaces for the velocity and the pressure:

$H\mathrm{<mo>\; :\; </mo>}=\left\{v\in H\left(\text{div}\mathrm{,}\Omega \right)\mathrm{<mo>\; :\; </mo>}{v}_s\in {\left[H^1\left({\Omega }_s\right)\right]}^N\mathrm{,}v=0\text{on}{\Gamma }_s\text{and}v\cdot n_d=0\text{on}{\Gamma }_d\right\}$

equipped with the norm

${\Vert v\Vert }_H\mathrm{<mo>\; :\; </mo>}={\left({\left|v\right|}_{\mathrm{1,}s}^2+{\Vert v\Vert }_d^2+{\Vert \text{div}v\Vert }_d^2\right)}^{1/2}\mathrm{,}$ (9)

and

$Q=L_0^2\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}=\left\{q\in L^2\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}{\displaystyle {\int }_{\Omega }}q\left(x\right)\text{d}x=0\right\}\mathrm{.}$ (10)

Multiplying the first equation of (1) by a test function $v\in H$ and the second one by $q\in Q$, integrating by parts over ${\Omega }_s$ the terms involving $\text{div}D\left(u\right)$ and $\nabla p$, yield the variational form of Stokes equations:

$\begin{array}{c}{\left(f_s,v_s\right)}_{{\Omega }_s}=2\mu {\left(D\left(u_s\right),D\left(v_s\right)\right)}_{{\Omega }_s}-{\left(p_s,\text{div}{v}_s\right)}_{{\Omega }_s}\\ +{\left(\left\{p_s-2\mu n_s\cdot D\left(u_s\right)\cdot n_s\right\},v_s\cdot n_s\right)}_{{\Gamma }_I}\\ +{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}}{\left(-2\mu n_s\cdot D\left(u_s\right)\cdot {\tau }_j,v_s\cdot {\tau }_j\right)}_{{\Gamma }_I}\end{array}$ (11)

$-{\left(q_s,\text{div}{u}_s\right)}_{{\Omega }_s}=-{\left(g_s,q_s\right)}_{{\Omega }_s}$ (12)

Using interface conditions (4) and (5) in (11), we obtain:

$\begin{array}{l}{\left(f_s,v\right)}_{{\Omega }_s}=2\mu {\left(D\left(u_s\right),D\left(v\right)\right)}_{{\Omega }_s}-{\left(p_s,\text{div}v\right)}_{{\Omega }_s}\\ +{\left(p_d,v_s\cdot n_s\right)}_{{\Gamma }_I}+{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}}\frac{\mu {\alpha }_1}{\sqrt{k_j}}{\left(u_s\cdot {\tau }_j,v_s\cdot {\tau }_j\right)}_{{\Gamma }_I}\forall v\in H\end{array}$ (13)

$-{\left(q,\text{div}{u}_s\right)}_{{\Omega }_s}=-{\left(g_s,q\right)}_{{\Omega }_s}\forall q\in Q$ (14)

We apply a similar treatment to the Darcy equations by testing the first equation of (2) with a smooth function $v\in H$ and the second on by $q\in Q$, integrating by parts over ${\Omega }_d$ the terms involving $\nabla p_d$, yield the variational form of Darcy equations:

${\left(\mu K^{-1}{u}_d,v\right)}_{{\Omega }_d}={\left(p_d,\text{div}v\right)}_{{\Omega }_d}+{\left(f_d,v\right)}_{{\Omega }_d}-{\left(p_d,v_d\cdot n_d\right)}_{{\Gamma }_I}\forall v\in H$ (15)

${\left(\text{div}{u}_d,q\right)}_{{\Omega }_d}={\left(g_d,q\right)}_{{\Omega }_d}\forall q\in Q$ (16)

Now, incorporating the first boundary interface condition (3) and taking into account that the vector valued functions in $H$ have (weakly) continuous normal components on ${\Gamma }_I$ (see [28], Theorem 2.5), the mixed variational formulation of the coupled problem (1)-(5) can be stated as follows [10]: Find $\left(u\mathrm{,}p\right)\in H\times Q$ that satisfies

$\{\begin{array}{l}a\left(u\mathrm{,}v\right)+b\left(v\mathrm{,}p\right)=L\left(v\right)\mathrm{,}\forall v\in H\mathrm{,}\\ b\left(u\mathrm{,}q\right)=G\left(q\right)\mathrm{,}\forall q\in Q\mathrm{.}\end{array}$ (17)

where the bilinear forms $a\left(\cdot \mathrm{,}\cdot \right)$ and $b\left(\cdot \mathrm{,}\cdot \right)$ are defined on $H\times H$ and $H\times Q$, respectively, as:

$a\left(u,v\right):=2\mu {\left(D\left(u\right),D\left(v\right)\right)}_s+{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}}\frac{\mu {\alpha }_1}{\sqrt{{\kappa }_j}}{\langle u_s\cdot {\tau }_j,v_s\cdot {\tau }_j\rangle }_{{\Gamma }_I}+\mu {\left(K^{-1}u,v\right)}_d$

$b\left(v,q\right):=-{\left(q,\text{div}v\right)}_{{\Omega }_s}-{\left(q,\text{div}v\right)}_{{\Omega }_d}$

By last, the linear forms L and G are defined as:

$L\left(v\right):={\left(f,v\right)}_{{\Omega }_s}+{\left(f,v\right)}_{{\Omega }_d}\text{and}G\left(q\right):=-{\left(g,q\right)}_{{\Omega }_s}-{\left(g,q\right)}_{{\Omega }_d}.$

Theorem 2.1. If $f\in {\left[L^2\left(\Omega \right)\right]}^N$ and $g\in L_0^2\left(\Omega \right)$ , there exists a unique solution $\left(u\mathrm{,}p\right)\in H\times Q$ to the problem (17).

Proof. We will establish the existence and uniqueness of a weak solution by using the classical theory of mixed methods so-called Brezzi conditions (see, e.g., [28], Theorem and Corollary 4.1 in Chapter I).

· It is easy to prove that $a$ and $b$ are continuous. It is also clear that F and G are continuous and bounded. In summary, the triangular inequality and Cauchy-Schwarz inequality lead to:

$\left|a\left(u\mathrm{,}v\right)\right|\lesssim {\Vert u\Vert }_H{\Vert v\Vert }_H\forall u\mathrm{,}v\in H\mathrm{,}$

$\left|b\left(v\mathrm{,}q\right)\right|\lesssim {\Vert v\Vert }_H\Vert q\Vert \forall v\in H\mathrm{,}\forall q\in Q\mathrm{,}$

$\left|L\left(v\right)\right|\lesssim {\Vert v\Vert }_H\forall v\in H\mathrm{,}$

$\left|G\left(q\right)\right|\lesssim \Vert q\Vert \forall q\in Q\mathrm{.}$

· Now we define the null space of $b$ i.e. $Z=\left\{v\in H,b\left(v\mathrm{,}q\right)=0\forall q\in L_0^2\left(\Omega \right)\right\}$. From the classical Korn’s inequality ( [8], Section 5), it follows that there is a positive constant C such that:

$a\left(v\mathrm{,}v\right)\ge C\left({\Vert \nabla v\Vert }_{{\Omega }_s}^2+{\Vert v\Vert }_{{\Omega }_d}^2\right)\mathrm{,}\forall v\in Z\mathrm{.}$ (18)

Let $q\in L^2\left(\Omega \right)$. Because $L^2\left(\Omega \right)=L_0^2\left(\Omega \right)\oplus ℝ$, then there exists unique $\left(q_1\mathrm{,}{q}_0\right)\in L_0^2\left(\Omega \right)\times ℝ$ such that $q=q_0+q_1$. Thus, we obtain:

$\begin{array}{c}b\left(v,q\right)=b\left(v,q_1\right)+b\left(v,q_0\right)=0+q_{0}b\left(v,1\right)\\ =-q_0{\displaystyle {\int }_{\Omega }}\nabla \cdot v=-q_0{\displaystyle {\int }_{{\Gamma }_s\cup {\Gamma }_d}}v\cdot n=0.\end{array}$

Hence, the coercivity of the bilinear form $a$,

$a\left(v,v\right)\ge C{\Vert v\Vert }_H^2\forall v\in Z$

follows.

· The second Brezzi condition that we need to verify is the inf-sup condition. Let $q\in Q=L_0^2\left(\Omega \right)$, then there exists $v\in {\left[H_0^1\left(\Omega \right)\right]}^N\backslash \left\{0\right\}\subset H$ such that,

$\nabla \cdot v=-q\text{in}\Omega ,$ (19)

with the estimation

${\Vert v\Vert }_H\lesssim {\Vert v\Vert }_{\mathrm{1,}\Omega }\lesssim \Vert q\Vert$ (20)

Thus, from (19) we have

$b\left(v\mathrm{,}q\right)=-{\displaystyle {\int }_{\Omega }}q\nabla \cdot v={\Vert q\Vert }^2\mathrm{.}$ (21)

Using the estimate (20), we obtain

$\frac{b\left(v\mathrm{,}q\right)}{{\Vert v\Vert }_H}\gtrsim \frac{{\Vert q\Vert }^2}{{\Vert v\Vert }_H}\gtrsim \Vert q\Vert \mathrm{.}$ (22)

Hence, the inf-sup condition,

$\underset{v\in H\mathrm{,}v\ne 0}{\sup }\frac{b\left(v\mathrm{,}q\right)}{{\Vert v\Vert }_H}\gtrsim \Vert q\Vert \forall q\in Q\mathrm{.}$ (23)

follows.

Remark 2.1. Note that if g is of mean zero, (17) directly implies that (1), (2) and (3) hold (the differential equations being understood in the distributional sense), while the interface conditions (4) and (5) are imposed in a weak sense. Also, we observe that the mixed variational formulation of the coupled problem (1)-(5) is equivalent to weak formulation (2.4) (and also (2.5) of [11] ), with the particularity that, in our case, for any $v\in H$ , we have that ${\langle v_s-v_d\mathrm{,}{n}_s{p}_s\rangle }_{{\Gamma }_I}=0$.

Now we introduce a modification to the Darcy equation, with the purpose in mind of the development of a unified discretization for the coupled problem, that is, the Stokes and Darcy parts be discretized using the same finite element spaces. The modification that we apply to the Darcy equation follows the idea (same argument) given in [22]. Indeed, we observe that taking the second equation of Darcy’s problem (2) we can write, for any $v\in H$,

${\displaystyle {\int }_{{\Omega }_d}}\left(\text{div}{u}_d-g_d\right)\text{div}v=0.$ (24)

Then, by adding this equation to the first equation of the variational form in (15), we get:

$\begin{array}{l}{\left(\mu K^{-1}{u}_d,v\right)}_{{\Omega }_d}+{\left(\text{div}{u}_d,\text{div}v\right)}_{{\Omega }_d}-{\left(p_d,\text{div}v\right)}_{{\Omega }_d}+{\left(p_d,v_d\cdot n_d\right)}_{{\Gamma }_I}\\ ={\left(f_d,v\right)}_{{\Omega }_d}+{\left(\text{div}v,g_d\right)}_{{\Omega }_d}\forall v\in H\end{array}$ (25)

${\left(\text{div}{u}_d,q\right)}_{{\Omega }_d}={\left(g_d,q\right)}_{{\Omega }_d}\forall q\in Q$ (26)

From now on, we work with this modified variational form of Darcy equations.

In the same way that before, incorporating the boundary conditions (3) and remambering that, since $v\in H$, it was (weakly) continuous normal components on ${\Gamma }_I$, the variational form of the modified Stokes-Darcy problem can be written as follows: Find $\left(u\mathrm{,}p\right)\in H\times Q$ satisfying

$\{\begin{array}{l}\tilde{a}\left(u\mathrm{,}v\right)+b\left(v\mathrm{,}p\right)=\tilde{L}\left(v\right)\mathrm{,}\forall v\in H\mathrm{,}\\ b\left(u\mathrm{,}q\right)=G\left(q\right)\mathrm{,}\forall q\in Q\mathrm{.}\end{array}$ (27)

where the bilinear forms $\tilde{a}\left(\cdot \mathrm{,}\cdot \right)$ and $b\left(\cdot \mathrm{,}\cdot \right)$ are defined on $H\times H$, $H\times Q$, respectively, as:

$\begin{array}{c}\tilde{a}\left(u\mathrm{,}v\right)=2\mu {\left(D\left(u\right),D\left(v\right)\right)}_s+{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}}\frac{\mu {\alpha }_1}{\sqrt{{\kappa }_j}}{\langle u_s\cdot {\tau }_j,v_s\cdot {\tau }_j\rangle }_{{\Gamma }_I}\\ +\mu {\left(K^{-1}u,v\right)}_d+{\left(\text{div}{u}_d,\text{div}v\right)}_{{\Omega }_d}\end{array}$

and

$b\left(v,q\right):=-{\left(q,\text{div}v\right)}_{{\Omega }_s}-{\left(q,\text{div}v\right)}_{{\Omega }_d}.$

By last, the linear forms $\tilde{L}$ and G are defined as:

$\tilde{L}\left(v\right):={\left(f,v\right)}_{{\Omega }_s}+{\left(f,v\right)}_{{\Omega }_d}+{\left(g,\text{div}v\right)}_d\text{and}G\left(q\right):=-{\left(g,q\right)}_{{\Omega }_s}-{\left(g,q\right)}_{{\Omega }_d}.$

Then, applying the classical theory of mixed methods it follows the well-posedness of the continuous formulation (27).

Theorem 2.2. There exists a unique $\left(u\mathrm{,}p\right)\in H\times Q$ solution to modified formulation (27). In addition, there exists a positive constant $\tilde{C}$ , depending on the continuous inf-sup condition constant for $b$ , the coercivity constant for $\tilde{a}$ and the boundedness constants for $\tilde{a}$ and $b$ , such that:

${\Vert u\Vert }_H+{\Vert p\Vert }_Q\le \tilde{C}\left({\Vert f_s\Vert }_{{\Omega }_s}+{\Vert f_d\Vert }_{{\Omega }_d}+{\Vert g_d\Vert }_{{\Omega }_d}+{\Vert g_s\Vert }_{{\Omega }_s}\right)\mathrm{.}$ (28)

We end this section with some notation. In 2D, the curl of a scalar function w is given as usual by $\text{curl}w\mathrm{<mo>\; :\; </mo>}={\left(\frac{\partial w}{\partial x_2}\mathrm{,}-\frac{\partial w}{\partial x_1}\right)}^{{\rm T}}$ while in 3D, the curl of a vector function $w$ is given as usual by $\text{curl}w\mathrm{<mo>\; :\; </mo>}=\nabla \times w$. Finally, let ${\mathbb{P}}^k$ be the space of polynomials of total degree not larger than k. In order to avoid excessive use of constants, the abbreviations $x\lesssim y$ and $x\sim y$ stand for $x\le cy$ and $c_{1}x\le y\le c_{2}x$, respectively, with positive constants independent of $x,y$ or ${\mathcal{T}}_h$.

3. A Priori Error Analysis

3.1. Finite Element Discretization

In this subsection, we will use a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element approximation for the velocity and piecewise constant approximation for the pressure.

Let ${\left\{{\mathcal{T}}_h\right\}}_{h>0}$ be a family of triangulations of $\Omega$ with nondegenerate elements (i.e. triangles for $N=2$ and tetrahedrons for $N=3$ ). For any $T\in {\mathcal{T}}_h$, we denote by $h_T$ the diameter of T and ${\rho }_T$ the diameter of the largest ball inscribed into T and set

$h=\underset{T\in {\mathcal{T}}_h}{\max }{h}_T\text{and}{\sigma }_h=\underset{T\in {\mathcal{T}}_h}{\max }\frac{h_T}{2r_T}$ (29)

We assume that the family of triangulations is regular, in the sense that there exists $\sigma >0$ such that ${\sigma }_h\le \sigma$, for all $h>0$. We also assume that the triangulation is conform with respect to the partition of $\Omega$ into ${\Omega }_s$ and ${\Omega }_d$, namely each $T\in {\mathcal{T}}_h$ is either in ${\Omega }_s$ or in ${\Omega }_d$ (see Figures 3-5).

Let ${\mathcal{T}}_h^s$ and ${\mathcal{T}}_h^d$ be the corresponding induced triangulations of ${\Omega }_s$ and ${\Omega }_d$. For any $T\in {\mathcal{T}}_h$, we denote by $\mathcal{E}\left(T\right)$ (resp. $\mathcal{N}\left(T\right)$ ) the set of its edges ( $N=2$ ) or faces ( $N=3$ ) (resp. vertices) and set ${\mathcal{E}}_h={\displaystyle \underset{T\in {\mathcal{T}}_h}{\cup }\mathcal{E}\left(T\right)}$, ${\mathcal{N}}_h={\displaystyle \underset{T\in {\mathcal{T}}_h}{\cup }\mathcal{N}\left(T\right)}$. For $\mathcal{A}\subset \bar{\Omega }$ we define

${\mathcal{E}}_h\left(\mathcal{A}\right)=\left\{E\in {\mathcal{E}}_h\mathrm{<mo>\; :\; </mo>}E\subset \mathcal{A}\right\}\mathrm{.}$

Notice that ${\mathcal{E}}_h$ can be split up in the form

${\mathcal{E}}_h={\mathcal{E}}_h\left({\Omega }_s^{+}\right)\cup {\mathcal{E}}_h\left({\Omega }_d\right)\cup {\mathcal{E}}_h\left(\partial {\Omega }_d\right)$ (30)

where ${\Omega }_s^{+}={\Omega }_s\cup {\Gamma }_s$. Note that ${\mathcal{E}}_h\left({\Gamma }_I\right)$ is included in ${\mathcal{E}}_h\left(\partial {\Omega }_d\right)$.

Figure 3. Isotropic element T in 2d.

Figure 4. Example of conforming mesh in 2d.

Figure 5. Example of nonconforming mesh in 2d.

With every edges $E\in {\mathcal{E}}_h$, we associate a unit vector $n_E$ such that $n_E$ is orthogonal to E and equals to the unit exterior normal vector to $\partial \Omega$ if $E\subset \partial \Omega$. For any $E\in {\mathcal{E}}_h$ and any piecewise continuous function $\phi$, we denote by ${\left[\phi \right]}_E$ its jump across E in the direction of $n_E$ :

${\left[\phi \right]}_E\left(x\right)\mathrm{<mo>\; :\; </mo>}=\{\begin{array}{ll}\underset{t\to 0+}{\lim }\phi \left(x+t{n}_E\right)-\underset{t\to 0+}{\lim }\phi \left(x-t{n}_E\right)\hfill & \text{for}\text{an}\text{interior}\text{edge/face}E\mathrm{,}\hfill \\ -\underset{t\to 0+}{\lim }\phi \left(x-t{n}_E\right)\hfill & \text{for}\text{a}\text{boundary}\text{edge/face}E\hfill \end{array}$

For $i\in \left\{\mathrm{0,}\cdots \mathrm{,}N\right\}$, we set (see Figure 6 in 2D for illustration):

${\sigma }_i\left(p\right)\mathrm{<mo>\; :\; </mo>}=\frac{1}{\left|E_i\right|}{\displaystyle {\int }_{E_i}}p\mathrm{,}\forall p\in {\mathbb{P}}^1\left(T\right)\mathrm{,}\text{where}{E}_i\in \mathcal{E}\left(T\right)$ (31)

The triplet $\left\{T\mathrm{,}{\mathbb{P}}^1\left(T\right)\mathrm{,}\Sigma \right\}$ with $\Sigma ={\left\{{\sigma }_i\right\}}_{0\le i\le N}$ is finite element [ [29], Page 83]. The local basis functions are defined by:

${\psi }_i\left(T\right)=1-N{\lambda }_i\left(T\right)\mathrm{,}i\in \left\{\mathrm{0,}\cdots \mathrm{,}N\right\}\mathrm{,}$ (32)

where for each $i\in \left\{\mathrm{0,}\cdots \mathrm{,}N\right\}$, ${\lambda }_i\left(T\right)$ is barycentric coordonates of $T\in {\mathcal{T}}_h$.

In classical reference element $\bar{T}$, the basis fonctions are given in ${ℝ}^2$ by:

$\{\begin{array}{l}{\bar{\psi }}_0\left(\bar{x}\mathrm{,}\bar{y}\right)=1-2\bar{y}\mathrm{,}\\ {\bar{\psi }}_1\left(\bar{x}\mathrm{,}\bar{y}\right)=-1+2\bar{x}+2\bar{y}\mathrm{,}\\ {\bar{\psi }}_2\left(\bar{x}\mathrm{,}\bar{y}\right)=1-2\bar{x}\mathrm{.}\end{array}$ (33)

Figures 7-9 represent the surface of these functions in space.

The measure of an element or edge/face is denoted by $\left|T\right|\mathrm{<mo>\; :\; </mo>}={\text{meas}}_N\left(T\right)$ and $\left|E\right|\mathrm{<mo>\; :\; </mo>}={\text{meas}}_{N-1}\left(E\right)$, respectively.

Based on the above notation, we introduce a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element space (larger than the space $H_h$ used in [14] )

$\begin{array}{l}{H}_h\mathrm{<mo>\; :\; </mo>}=\{v_h\mathrm{<mo>\; :\; </mo>}{v}_h{}_{\mathrm{<mo>\; |\; </mo>}T}\in {\left[{\mathbb{P}}^1\left(T\right)\right]}^N\forall T\in {\mathcal{T}}_h\mathrm{,}{\left({\left[v_h\right]}_E\mathrm{,}1\right)}_E=0\forall E\in {\mathcal{E}}_h\left({\Omega }_s^{+}\right)\mathrm{,}\\ {\left({\left[v_h\cdot n_E\right]}_E\mathrm{,1}\right)}_E=0\forall E\in {\mathcal{E}}_h\left({\Omega }_d\right)\cup {\mathcal{E}}_h\left(\partial {\Omega }_d\right)\}\end{array}$ (34)

and piecewise constant function space

$Q_h\mathrm{<mo>\; :\; </mo>}=\left\{q_h\in L_0^2\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}{q}_h{}_{\mathrm{<mo>\; |\; </mo>}T}\in {\mathbb{P}}^0\left(T\right)\forall T\in {\mathcal{T}}_h\right\}\mathrm{,}$

where ${\mathbb{P}}^m\left(T\right)$ is the space of the restrictions to T of all polynomials of degree less than or equal to m. The space $Q_h$ is equipped with the norm $\Vert \cdot \Vert$ while

Figure 6. ${\mathbb{P}}^1$ -nonconforming finite element T in 2d.

Figure 7. ${\bar{\psi }}_0$.

Figure 8. ${\bar{\psi }}_1$.

Figure 9. ${\bar{\psi }}_2$.

the norm on $H_h$ will be specified later on. The choice of $H_h$ is more natural than the one introduced in [14] since the space $H_h$ approximates only $H\left(\text{div}\mathrm{,}{\Omega }_d\right)$ and not ${\left[H^1\left({\Omega }_d\right)\right]}^N$, while our a priori error analysis is only valid in this larger space.

Let us introduce the discrete divergence operator ${\text{div}}_h\in \mathcal{L}\left(H_h\mathrm{<mo>\; ;\; </mo>}{Q}_h\right)\cap \mathcal{L}\left(H\mathrm{<mo>\; ;\; </mo>}Q\right)$ by

${\left({\text{div}}_h{v}_h\right)}_{\mathrm{<mo>\; |\; </mo>}T}=\text{div}\left(v_h{}_{\mathrm{<mo>\; |\; </mo>}T}\right)\mathrm{,}\forall T\in {\mathcal{T}}_h\mathrm{.}$ (35)

Then, we can introduce two bilinear forms

$\begin{array}{l}{\tilde{a}}_h\left(u\mathrm{,}v\right)\mathrm{<mo>\; :\; </mo>}=2\mu {\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}{\left(D\left(u\right)\mathrm{,}D\left(v\right)\right)}_T+{\displaystyle \underset{j\mathrm{<mo>\; =\; </mo>1}}{\overset{N-1}{\sum }}}\frac{\mu {\alpha }_1}{\sqrt{{\kappa }_j}}{\langle u_s\cdot {\tau }_j\mathrm{,}{v}_s\cdot {\tau }_j\rangle }_{{\Gamma }_I}\\ +\mu {\left(K^{-1}u\mathrm{,}v\right)}_{{\Omega }_d}+{\left({\text{div}}_{h}u\mathrm{,}{\text{div}}_{h}v\right)}_{{\Omega }_d}\mathrm{,}\forall u\mathrm{,}v\in H\cup H_h\end{array}$

and

$b_h\left(v\mathrm{,}q\right)\mathrm{<mo>\; :\; </mo>}=-{\left(q\mathrm{,}{\text{div}}_{h}v\right)}_{\Omega }\mathrm{,}\forall v\in H\cup H_h\mathrm{,}\forall q\in Q_h\mathrm{.}$

Then the finite element discretization of (27) is to find $\left(u_h\mathrm{,}{p}_h\right)\in H_h\times Q_h$ such that

$\{\begin{array}{l}{\tilde{a}}_h\left(u_h\mathrm{,}{v}_h\right)+b_h\left(v_h\mathrm{,}{p}_h\right)+J\left(u_h\mathrm{,}{v}_h\right)=\tilde{L}\left(v_h\right)\mathrm{,}\forall v_h\in H_h\mathrm{,}\\ b_h\left(u_h\mathrm{,}{q}_h\right)=G\left(q_h\right)\mathrm{,}\forall q_h\in Q_h\mathrm{.}\end{array}$ (36)

This is the natural discretization of the modified weak formulation (27) except that the penalizing term $J\left(u_h\mathrm{,}{v}_h\right)$ is added. This bilinear form $J\left(\mathrm{.,.}\right)$ is defined by following the decomposition (37) of ${\mathcal{E}}_h$ :

$J\left(u\mathrm{,}v\right)=J_{{\Omega }_s^{+}}\left(u\mathrm{,}v\right)+J_{{\Omega }_d}\left(u\mathrm{,}v\right)+J_{\partial {\Omega }_d}\left(u\mathrm{,}v\right)$ (37)

where

$J_{{\Omega }_s^{+}}\left(u\mathrm{,}v\right)\mathrm{<mo>\; :\; </mo>}=\left(1+2\mu \right){\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Omega }_s^{+}\right)}{\sum }}{h}_E^{-1}{\displaystyle {\int }_E}{\left[u\right]}_E\cdot {\left[v\right]}_E\text{d}s\mathrm{,}$

$J_{{\Omega }_d}\left(u\mathrm{,}v\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Omega }_d\right)}{\sum }}{h}_E^{-1}{\displaystyle {\int }_E}{\left[u\right]}_E\cdot {\left[v\right]}_E\text{d}s\mathrm{,}$

and

$J_{\partial {\Omega }_d}\left(u\mathrm{,}v\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{E\in {\mathcal{E}}_h\left(\partial {\Omega }_d\right)}{\sum }}{h}_E^{-1}{\displaystyle {\int }_E}{\left[u\cdot n_E\right]}_E{\left[v\cdot n_E\right]}_E\text{d}s\mathrm{.}$

Here, $h_E$ is the length ( $N=2$ ) or diameter ( $N=3$ ) of E. Note that each element of ${\mathcal{E}}_h$ only contributes with one jump term in $J\left(u\mathrm{,}v\right)$.

Remark 3.1 The Equation (36) have the matrix representation

$M_{A}U+M_B^{\text{T}}P+M_{J}U=F$

$M_{B}U=G$

where $U$ (resp. $P$ ) denote the coefficients of $u_h$ (resp. $p_h$ ) expanded with respect to a basis for $H_h$ (resp. $Q_h$ ).

We are now able to define the norm on $H_h$ (see [14] ):

${\Vert v\Vert }_h\mathrm{<mo>\; :\; </mo>}={\left({\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}{\left|v\right|}_{\mathrm{1,}T}^2+{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}}{\langle v_s\cdot {\tau }_j\mathrm{,}{v}_s\cdot {\tau }_j\rangle }_{{\Gamma }_I}+{\Vert v\Vert }_{{\Omega }_d}^2+{\Vert {\text{div}}_{h}v\Vert }_{{\Omega }_d}^2+J\left(v,v\right)\right)}^{1/2}\mathrm{.}$

In the sequel, we will denote by $\alpha ,\beta$ and $C_i$ various constants independent of h. For the sake of convenience, we will define the bilinear form:

$A_h\left(u\mathrm{,}v\right)={\tilde{a}}_h\left(u\mathrm{,}v\right)+J\left(u\mathrm{,}v\right)\mathrm{.}$

From Hölder’s inequality, we derive the boundedness of $A_h\left(\cdot \mathrm{,}\cdot \right)$ and $b_h\left(\cdot \mathrm{,}\cdot \right)$ :

Lemma 3.1. (Continuity of forms) There holds:

$\left|\tilde{L}\left(v_h\right)\right|\le C_1{\Vert v_h\Vert }_h\mathrm{,}\forall v_h\in H_h\cup H\mathrm{,}$ (38)

$\left|G\left(q_h\right)\right|\le C_2\Vert q_h\Vert \mathrm{,}\forall q_h\in Q_h\mathrm{,}$ (39)

$\left|A_h\left(u_h\mathrm{,}{v}_h\right)\right|\le C_3{\Vert u_h\Vert }_h{\Vert v_h\Vert }_h\mathrm{,}\forall u_h\mathrm{,}{v}_h\in H_h\cup H\mathrm{,}$ (40)

$\left|b_h\left(v_h\mathrm{,}{q}_h\right)\right|\le C_4{\Vert v_h\Vert }_h\Vert q_h\Vert \mathrm{,}\forall v_h\in H_h\cup H\mathrm{,}\forall q_h\in Q_h\mathrm{.}$ (41)

Lemma 3.2. (Coercivity of A_h) There is an $\alpha >0$ such that:

$A_h\left(v_h\mathrm{,}{v}_h\right)\ge \alpha {\Vert v_h\Vert }_h^2\forall v_h\in H_h\mathrm{.}$ (42)

Proof. Let $v_h\in H_h$. We have

$\begin{array}{c}{A}_h\left(v_h\mathrm{,}{v}_h\right)=2\mu {\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}{\Vert D\left(v_h\right)\Vert }_T^2+\mu {\left(K^{-1}{v}_h\mathrm{,}{v}_h\right)}_{{\Omega }_d}+{\displaystyle \underset{j\mathrm{<mo>\; =\; </mo>1}}{\overset{N-1}{\sum }}}\frac{\mu {\alpha }_1}{{\kappa }_j}{\Vert v_h\cdot {\tau }_j\Vert }_{{\Gamma }_I}^2\\ +{\Vert {\text{div}}_{h}v\Vert }_{{\Omega }_d}^2+J_{{\Omega }_s^{+}}\left(v_h\mathrm{,}{v}_h\right)+J_{{\Omega }_d}\left(v_h\mathrm{,}{v}_h\right)+J_{\partial {\Omega }_d}\left(v_h\mathrm{,}{v}_h\right)\end{array}$

We introduce the local space

$H\left(\text{curl},T\right):=\{\begin{array}{ll}\left\{v\in {\left[L^2\left(T\right)\right]}^2:\text{curl}v\in L^2\left(T\right)\right\}\hfill & \text{if}N=2\hfill \\ \left\{v\in {\left[L^2\left(T\right)\right]}^3:\text{curl}v\in {\left[L^2\left(T\right)\right]}^3\right\}\hfill & \text{if}N=3\hfill \end{array}$

and for $\psi \in {\left[H^1\left(T\right)\right]}^N$, we define

${\gamma }_{\tau }\psi :=\{\begin{array}{ll}\psi \cdot {\tau }_{|\partial T}\hfill & \text{if}N=2,\hfill \\ \psi \times n_{|\partial T}\hfill & \text{if}N=3,\left(\tau \cdot n=0\text{on}\partial T\right),\hfill \end{array}$

with the semi-norm:

$\varphi \left(v_h\right)={\left|{\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}{\displaystyle {\int }_T}\text{curl}{v}_h\right|}_{{ℝ}^l}\mathrm{,}\left(\text{where}l=1\text{or}l=3\right)\mathrm{.}$ (43)

Using Young’s inequality and Green formula, we have:

$\begin{array}{c}\varphi \left(v_h\right)={\left|{\displaystyle {\int }_{{\Omega }_s}}\text{curl}{v}_h\right|}_{{ℝ}^l}={\left|{\displaystyle {\int }_{\partial {\Omega }_s}}{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}\\ ={\left|{\displaystyle {\int }_{{\Gamma }_s}}{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}+{\left|{\displaystyle {\int }_{{\Gamma }_I}}{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}\\ \lesssim {\displaystyle {\int }_{{\Gamma }_s}}{\left|{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}+{\displaystyle {\int }_{{\Gamma }_I}}{\left|{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}\mathrm{.}\end{array}$

· Estimate ${\displaystyle {\sum }_{E\in {\mathcal{E}}_h\left({\Gamma }_s\right)}}{\displaystyle {\int }_E}{\left|{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}$ ( $l=1$ or $l=3$ ). We have by Cauchy-Schwarz inequality:

$\begin{array}{c}{\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Gamma }_s\right)}{\sum }}{\displaystyle {\int }_E}{\left|{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}\le {\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Gamma }_s\right)}{\sum }}\left\{{\left({\displaystyle {\int }_E}{\left|{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}^2\right)}^{1/2}{\left|h_E\right|}^{1/2}\right\}\\ \le {\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Gamma }_s\right)}{\sum }}\left\{h_E^{-1/2}{\left({\displaystyle {\int }_E}{\left|{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}^2\right)}^{1/2}{h}_E\right\}\\ \le {\left({\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Gamma }_s\right)}{\sum }}{h}_E^{-1}{\displaystyle {\int }_E}{\left|{\left[v_h\right]}_E^2\right|}_{{ℝ}^N}\right)}^{1/2}{\left({\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Gamma }_s\right)}{\sum }}{h}_E^2\right)}^{1/2}\end{array}$

Also, we have:

${\left({\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Gamma }_s\right)}{\sum }}{h}_E^2\right)}^{1/2}\lesssim \mathrm{1,}$ (44)

Then,

${\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Gamma }_s\right)}{\sum }}{\displaystyle {\int }_E}{\left|{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}\lesssim {\left({\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Gamma }_s\right)}{\sum }}{h}_E^{-1}{\displaystyle {\int }_E}{\left|{\left[v_h\right]}_E\right|}_{{ℝ}^N}^2\right)}^{1/2}\mathrm{.}$ (45)

Hence we deduce

${\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Gamma }_s\right)}{\sum }}{\displaystyle {\int }_E}{\left|{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}\lesssim {\left(J_{{\Omega }_s^{+}}\left(v_h\mathrm{,}{v}_h\right)\right)}^{1/2}\mathrm{.}$ (46)

· Now we estime the term ${\displaystyle {\sum }_{E\in {\mathcal{E}}_h\left({\Gamma }_I\right)}}{\displaystyle {\int }_E}{\left|{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}$. By Cauchy-Schwarz, we obtain:

${\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Gamma }_I\right)}{\sum }}{\displaystyle {\int }_E}{\left|{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}\le {\left({\displaystyle {\int }_{{\Gamma }_I}}{\left|{\gamma }_{\tau }\left(v_h\right)\right|}_{{ℝ}^l}^2\right)}^{1/2}{\left|{\Gamma }_I\right|}^{1/2}\lesssim {\tilde{a}}_h{\left(v_h\mathrm{,}{v}_h\right)}^{1/2}\mathrm{.}$

Thus we deduce the estimation:

${\left(\varphi \left(v_h\right)\right)}^2\lesssim J_{{\Omega }_s^{+}}\left(v_h\mathrm{,}{v}_h\right)+{\tilde{a}}_h\left(v_h\mathrm{,}{v}_h\right)\mathrm{.}$ (47)

Then,

$\begin{array}{l}{J}_{{\Omega }_s^{+}}\left(v_h\mathrm{,}{v}_h\right)+a_h\left(v_h\mathrm{,}{v}_h\right)+{\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}{\Vert D\left(v_h\right)\Vert }_T^2\\ \gtrsim {\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}{\Vert D\left(v_h\right)\Vert }_T^2+{\left(\varphi \left(v_h\right)\right)}^2+J_{{\Omega }_s^{+}}\left(v_h\mathrm{,}{v}_h\right)\mathrm{.}\end{array}$

We apply Korn’s discrete inequality [30] and we get:

$J_{{\Omega }_s^{+}}\left(v_h\mathrm{,}{v}_h\right)+a_h\left(v_h\mathrm{,}{v}_h\right)+{\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}{\Vert D\left(v_h\right)\Vert }_T^2\gtrsim {\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}{\Vert \nabla \left(v_h\right)\Vert }_T^2\mathrm{.}$ (48)

Thus

$J\left(v_h\mathrm{,}{v}_h\right)+{\tilde{a}}_h\left(v_h\mathrm{,}{v}_h\right)+{\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}{\Vert D\left(v_h\right)\Vert }_T^2\gtrsim {\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}{\Vert \nabla \left(v_h\right)\Vert }_T^2+J\left(v_h\mathrm{,}{v}_h\right)\mathrm{,}$

Hence,

$A_h\left(v_h\mathrm{,}{v}_h\right)\gtrsim {\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}{\Vert \nabla \left(v_h\right)\Vert }_T^2+J\left(v_h\mathrm{,}{v}_h\right).$ (49)

We have,

$A_h\left(v_h\mathrm{,}{v}_h\right)\ge {\displaystyle \underset{j=1}{\overset{N-1}{\sum }}}{\Vert v_h\cdot {\tau }_j\Vert }_{{\Gamma }_I}^2\mathrm{,}$ (50)

$A_h\left(v_h\mathrm{,}{v}_h\right)\ge {\Vert v_h\Vert }_{{\Omega }_d}^2$ (51)

$A_h\left(v_h\mathrm{,}{v}_h\right)\ge {\Vert {\text{div}}_h{v}_h\Vert }_{{\Omega }_d}^2$ (52)

The estimates (49), (50), (51) and (52), lead to (42). The proof is complete. o

In order to verify the discrete inf-sup condition, we define the space:

$W\mathrm{<mo>\; :\; </mo>}=\left\{v\in H\mathrm{<mo>\; :\; </mo>}{v}_{\mathrm{<mo>\; |\; </mo>}{\Omega }_d}\in {\left[H^1\left({\Omega }_d\right)\right]}^N\right\}\mathrm{.}$ (53)

We define also the Crouzeix-Raviart interpolation operator $r_h\mathrm{<mo>\; :\; </mo>}W\to H_h$ by:

${\displaystyle {\int }_E}{\left(r_{h}v\right)}_s\text{d}s={\displaystyle {\int }_E}{v}_s\text{d}s\mathrm{,}\forall E\in {\mathcal{E}}_h\left({\bar{\Omega }}_s\right)\mathrm{,}\forall v\in W\mathrm{,}$ (54)

${\displaystyle {\int }_E}{\left(r_{h}v\right)}_d\text{d}s={\displaystyle {\int }_E}{v}_d\text{d}s\mathrm{,}\forall E\in {\mathcal{E}}_h\left({\bar{\Omega }}_d\right)\mathrm{,}\forall v\in W\mathrm{.}$ (55)

Lemma 3.3. The operator $r_h$ is bounded: there is a constant $C_5>0$ depending on $\sigma$ , $\mu$ and N such that

${\Vert r_{h}v\Vert }_h\le C_5{\left({\Vert v\Vert }_{\mathrm{1,}s}^2+{\Vert v\Vert }_{\mathrm{1,}d}^2\right)}^{1/2}\mathrm{,}\forall v\in W\mathrm{.}$ (56)

Proof. The proof is similar to ( [14], Lemma 4.5, Page 2695). o

Then, we have the following result

Lemma 3.4. (Inf-Sup condition) There exists a positive constant $\beta$ depending on $\sigma$ , $\mu$ and N such that

$\underset{q_h\in Q_h}{\inf }\underset{v_h\in H_h}{\sup }\frac{b_h\left(v_h,q_h\right)}{{\Vert v_h\Vert }_h\Vert q_h\Vert }\ge \beta .$ (57)

Proof. We use Fortin argument [31], i.e. for each $q_h\in Q_h$, we find $v_h\in H_h$ such that:

$b_h\left(v_h,q_h\right)={\Vert q_h\Vert }_{\Omega }^2\text{and}{\Vert v_h\Vert }_h\lesssim {\Vert q_h\Vert }_{\Omega }\mathrm{.}$

Let $q_h\in Q_h\subset Q$. Then from ( [28], Corollary 2.4, Page 24), there exist vectoriel function $v\in {\left[H_0^1\left(\Omega \right)\right]}^N$ satisfying

$\{\begin{array}{l}\text{div}v=-q_h\mathrm{,}\text{in}\Omega \\ {\Vert v\Vert }_{\mathrm{1,}\Omega }\lesssim {\Vert q_h\Vert }_{\Omega }\mathrm{.}\end{array}$ (58)

${\left[H_0^1\left(\Omega \right)\right]}^N\subset W$, hence $v\in W$. We take $v_h=r_{h}v\in H_h$ and we have:

$\begin{array}{c}{b}_h\left(v-r_{h}v\mathrm{,}{q}_h\right)=-{\displaystyle \underset{T\in {\mathcal{T}}_h}{\sum }}{\displaystyle {\int }_T}{q}_h\text{div}\left(v-r_{h}v\right)\mathrm{,}\\ =-{\displaystyle \underset{T\in {\mathcal{T}}_h}{\sum }}{\displaystyle {\int }_{\partial T}}{q}_h{n}_T\cdot \left(v-r_{h}v\right)\\ =-{\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Omega }_s^{+}\right)}{\sum }}{\displaystyle {\int }_E}{q}_h{n}_E\cdot \left(v-r_{h}v\right)-{\displaystyle \underset{E\in {\mathcal{E}}_h\left({\bar{\Omega }}_d\right)}{\sum }}{\displaystyle {\int }_E}{q}_h{n}_E\cdot \left(v-r_{h}v\right)\\ =0\left(\text{fromtheidentities}\left(54\right)\text{and}\left(55\right)\right).\end{array}$

Thus, we obtain

$b_h\left(v\mathrm{,}{q}_h\right)=b_h\left(r_{h}v\mathrm{,}{q}_h\right)\mathrm{.}$

Using the system (58), we have:

$b_h\left(r_{h}v\mathrm{,}{q}_h\right)=-{\displaystyle {\int }_{\Omega }}{q}_h\text{div}\left(v\right)={\Vert q_h\Vert }^2\gtrsim {\Vert v\Vert }_{\mathrm{1,}\Omega }\Vert q_h\Vert \mathrm{.}$ (59)

Also,

${\Vert v_h\Vert }_h={\Vert r_{h}v\Vert }_h\lesssim {\Vert v\Vert }_{\mathrm{1,}\Omega }\mathrm{.}$ (60)

From (59) et (60), we deduce:

$b_h\left(r_{h}v\mathrm{,}{q}_h\right)\gtrsim {\Vert v_h\Vert }_h\Vert q_h\Vert \mathrm{,}\forall q_h\in Q_h\mathrm{.}$ (61)

The Inf-Sup condition holds and the proof is complete. o

From Lemma 3.2 and Lemma 3.4 we have the following result:

Theorem 3.1. There exists a unique solution $\left(u_h\mathrm{,}{p}_h\right)\in H_h\times Q_h$ to the problem (27).

3.2. A convergence Analysis

We now present an a priori analysis of the approximation error: The use of nonconforming finite element leads to $H_h⊈H$, so the approximation error contains some extra consistency error terms. In fact, the abstract error estimates from [14] give the following result:

Lemma 3.5. Let $\left(u\mathrm{,}p\right)\in H\times Q$ be the solution of problem (27) and $\left(u_h\mathrm{,}{p}_h\right)\in H_h\times Q_h$ be the solution of the discrete problem (36). Then we have

${\Vert u-u_h\Vert }_h+\Vert p-p_h\Vert \lesssim \underset{v_h\in H_h}{\inf }{\Vert u-v_h\Vert }_h+\underset{q_h\in Q_h}{\inf }\Vert p-q_h\Vert +E_{1h}+E_{2h}\mathrm{.}$ (62)

where $E_{1h}$ and $E_{2h}$ are the consistency error terms define by:

$E_{1h}=\underset{v_h\in H_h}{\sup }\frac{\left|A_h\left(u\mathrm{,}{v}_h\right)+b_h\left(v_h\mathrm{,}p\right)-{\left(f\mathrm{,}{v}_h\right)}_{\Omega }-{\left(g_d\mathrm{,}{\text{div}}_h{v}_h\right)}_{{\Omega }_d}\right|}{{\Vert v_h\Vert }_h}\mathrm{,}$ (63)

$E_{2h}=\underset{q_h\in Q_h}{\sup }\frac{\left|b_h\left(u\mathrm{,}{q}_h\right)+{\left(g\mathrm{,}{q}_h\right)}_{\Omega }\right|}{\Vert q_h\Vert }\mathrm{.}$ (64)

Note that $b_h\left(u\mathrm{,}{q}_h\right)=b\left(u\mathrm{,}{q}_h\right)$, thus $E_{2h}=0$.

For estiming the approximation error, we assume that the solution $\left(u\mathrm{,}p\right)$ of problem (27) satisfies the smoothness assumptions:

Assumption 3.1.

1) $u\in H$, $u_s\in {\left[H^2\left({\Omega }_s\right)\right]}^N$, $u_d\in {\left[H^2\left({\Omega }_d\right)\right]}^N$ ;

2) $p\in Q$, $p_s\in H^1\left({\Omega }_s\right)$, $p_d\in H^1\left({\Omega }_d\right)$.

We begin with an estimate for the terms ${\inf }_{v_h\in H_h}{\Vert u-v_h\Vert }_h$ and ${\inf }_{q_h\in Q_h}\Vert p-q_h\Vert$.

Lemma 3.6. (Ref. [14] ) There hold:

$\underset{v_h\in H_h}{\inf }{\Vert u-v_h\Vert }_h\lesssim h\left({\left|u\right|}_{\mathrm{2,}s}+{\left|u\right|}_{\mathrm{2,}d}\right)\mathrm{,}$ (65)

$\underset{q_h\in Q_h}{\inf }\Vert p-q_h\Vert \lesssim h\left({\left|p\right|}_{\mathrm{1,}s}+{\left|p\right|}_{\mathrm{1,}d}\right)\mathrm{.}$ (66)

Finally, let us consider the term $A_h\left(u_h\mathrm{,}{v}_h\right)+b_h\left(v_h\mathrm{,}{p}_h\right)-\tilde{L}\left(v_h\right)$. The smoothness assumption of $u$ implies $J\left(u\mathrm{,}{v}_h\right)=0$, thus $A_h\left(u\mathrm{,}{v}_h\right)=\tilde{a}\left(u\mathrm{,}{v}_h\right)\mathrm{,}\forall v_h\in H_h$. Clearly,

$\begin{array}{l}-\tilde{L}\left(v_h\right)=-{\left(f\mathrm{,}{v}_h\right)}_{\Omega }-{\left(g\mathrm{,}{\text{div}}_h{v}_h\right)}_{{\Omega }_d}\\ =-{\left(f\mathrm{,}{v}_h\right)}_{{\Omega }_s}-{\left(f\mathrm{,}{v}_h\right)}_{{\Omega }_d}-{\left(g\mathrm{,}{\text{div}}_h{v}_h\right)}_{{\Omega }_d}\\ ={\left(2\mu \text{div}D\left(u\right)-\nabla p\mathrm{,}{v}_h\right)}_{{\Omega }_s}-{\left(\mu K^{-1}u+\nabla p\mathrm{,}{v}_h\right)}_{{\Omega }_d}-{\left(\text{div}u\mathrm{,}{\text{div}}_h{v}_h\right)}_{{\Omega }_d}\\ =2\mu {\left(\text{div}D\left(u\right)\mathrm{,}{v}_h\right)}_{{\Omega }_s}-{\left(\nabla p\mathrm{,}{v}_h\right)}_{{\Omega }_s}-\mu {\left(K^{-1}u\mathrm{,}{v}_h\right)}_{{\Omega }_d}\\ -{\left(\nabla p\mathrm{,}{v}_h\right)}_{{\Omega }_d}-{\left(\text{div}u\mathrm{,}{\text{div}}_h{v}_h\right)}_{{\Omega }_d}\end{array}$

$\begin{array}{l}={\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}\left\{2\mu {\left(\text{div}D\left(u\right)\mathrm{,}{v}_h\right)}_T-{\left(\nabla p\mathrm{,}{v}_h\right)}_T\right\}\\ +{\displaystyle \underset{T\in {\mathcal{T}}_h^d}{\sum }}\left\{\mu {\left(K^{-1}u\mathrm{,}{v}_h\right)}_T-{\left(\nabla p\mathrm{,}{v}_h\right)}_T\right\}-{\left(\text{div}u\mathrm{,}{\text{div}}_h{v}_h\right)}_{{\Omega }_d}\\ ={\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}\{-2\mu {\left(D\left(u\right)\mathrm{,}D\left(v_h\right)\right)}_T+2\mu {\left(n_T\cdot D\left(u\right)\mathrm{,}{v}_h\right)}_{\partial T}+{\left(p\mathrm{,}\text{div}{v}_h\right)}_T\\ -{\left(v_h\cdot n_T\mathrm{,}p\right)}_{\partial T}\}+{\displaystyle \underset{T\in {\mathcal{T}}_h^d}{\sum }}\{-\mu {\left(K^{-1}u\mathrm{,}{v}_h\right)}_T+{\left(p\mathrm{,}\text{div}{v}_h\right)}_T\\ -{\left(v_h\cdot n_T\mathrm{,}p\right)}_{\partial T}\}-{\left(\text{div}u\mathrm{,}{\text{div}}_h{v}_h\right)}_{{\Omega }_d}\end{array}$

$\begin{array}{l}=-\left\{{\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}2\mu {\left(D\left(u\right)\mathrm{,}D\left(v_h\right)\right)}_T\right\}+\mu {\left(K^{-1}u\mathrm{,}{v}_h\right)}_{{\Omega }_d}-{\left(p\mathrm{,}{\text{div}}_h{v}_h\right)}_{\Omega }\\ -{\left(\text{div}u\mathrm{,}{\text{div}}_h{v}_h\right)}_{{\Omega }_d}+{\displaystyle \underset{T\in {\mathcal{T}}_h^s}{\sum }}\left\{2\mu {\left(n_T\cdot D\left(u\right)\mathrm{,}{v}_h\right)}_{\partial T}-{\left(v_h\cdot n_T\mathrm{,}p\right)}_{\partial T}\right\}\\ -{\displaystyle \underset{T\in {\mathcal{T}}_h^d}{\sum }}{\left(v_h\cdot n_T\mathrm{,}p\right)}_{\partial T}\end{array}$

$\begin{array}{l}=-{\tilde{a}}_h\left(u\mathrm{,}{v}_h\right)+{\displaystyle \underset{j\mathrm{<mo>\; =\; </mo>1}}{\overset{N-1}{\sum }}}\frac{\mu {\alpha }_1}{\sqrt{{\kappa }_j}}{\left(u_s\cdot {\tau }_j\mathrm{,}{v}_{h\mathrm{,}s}\cdot {\tau }_j\right)}_{{\Gamma }_I}-b_h\left(u\mathrm{,}{v}_h\right)\\ +2\mu {\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Omega }_s^{+}\right)}{\sum }}{\left(n_E\cdot D\left(u\right)\mathrm{,}{\left[v_h\right]}_E\right)}_E-{\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Omega }_d\right)\cup {\mathcal{E}}_h\left(\partial {\Omega }_d\right)}{\sum }}{\left({\left[v_h\cdot n_E\right]}_E\mathrm{,}{p}_d\right)}_E\\ -{\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Omega }_s^{+}\right)}{\sum }}{\left({\left[v_h\cdot n_E\right]}_E\mathrm{,}{p}_s\right)}_E\mathrm{.}\end{array}$

Thus, we have

$\begin{array}{l}\tilde{a}\left(u\mathrm{,}{v}_h\right)+J\left(u\mathrm{,}{v}_h\right)+b_h\left(u\mathrm{,}{v}_h\right)-\tilde{L}\left(v_h\right)\\ =R_1\left(v_h\right)+R_2\left(v_h\right)+R_3\left(v_h\right)+R_4\left(v_h\right)\mathrm{,}\end{array}$ (67)

where

$R_1\left(v_h\right)={\displaystyle \underset{j\mathrm{<mo>\; =\; </mo>1}}{\overset{N-1}{\sum }}}\frac{\mu {\alpha }_1}{\sqrt{{\kappa }_j}}{\left(u_s\cdot {\tau }_j\mathrm{,}{v}_{h\mathrm{,}s}\cdot {\tau }_j\right)}_{{\Gamma }_I},$

$R_2\left(v_h\right)=2\mu {\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Omega }_s^{+}\right)}{\sum }}{\left(n_E\cdot D\left(u\right)\mathrm{,}{\left[v_h\right]}_E\right)}_E,$

$R_3\left(v_h\right)={\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Omega }_d\right)\cup {\mathcal{E}}_h\left(\partial {\Omega }_d\right)}{\sum }}{\left({\left[v_h\cdot n_E\right]}_E\mathrm{,}{p}_d\right)}_E,$

$R_4\left(v_h\right)={\displaystyle \underset{E\in {\mathcal{E}}_h\left({\Omega }_s^{+}\right)}{\sum }}{\left({\left[v_h\cdot n_E\right]}_E\mathrm{,}{p}_s\right)}_E\mathrm{.}$

In order to evaluate the four face integrals, let us introduce two projections operators in the following.

For any $T\in {\mathcal{T}}_h$ and $E\in \mathcal{E}\left(T\right)$, denote by ${\mathbb{P}}^0\left(E\right)$ the constant space of the restrictions to E and ${\pi }_E$ the projection operator from $L^2\left(E\right)$ on to ${\mathbb{P}}^0\left(E\right)$ such that

${\displaystyle {\int }_E}{\pi }_{E}v={\displaystyle {\int }_E}v\text{d}s\mathrm{.}$ (68)

The operator ${\pi }_E$ has the property [32]:

${\Vert v-{\pi }_E\Vert }_{\mathrm{0,}E}\lesssim h_E^{1/2}{\left|v\right|}_{\mathrm{1,}T}\forall v\in H^1\left(T\right)\mathrm{.}$ (69)

For any $v\in {\left[L^2\left(E\right)\right]}^N$, we let ${\Pi }_{E}v$ be the function in ${\left[{\mathbb{P}}^0\left(E\right)\right]}^N$ such that

${\left({\Pi }_{E}v\right)}_i={\pi }_E{v}_i\mathrm{,1}\le i\le N\mathrm{.}$

Using inequality (69), we obtain

${\Vert v-{\Pi }_{E}v\Vert }_{\mathrm{0,}E}\lesssim h_E^{1/2}{\left|v\right|}_{\mathrm{1,}T}\forall v\in {\left[H^1\left(T\right)\right]}^N\mathrm{.}$ (70)

Then we have the following lemma:

Lemma 3.7. (Estimation the four face integrals) There holds:

$\left|R_1\left(v_h\right)\right|\le \underset{1\le j\le N-1}{\max }\left(\frac{\mu {\alpha }_1}{\sqrt{{\kappa }_j}}\right)h{\Vert u_s\Vert }_{\mathrm{1,}{\Omega }_s}{\Vert v_h\Vert }_h$ (71)

$\left|R_2\left(v_h\right)\right|\lesssim {\left|u\right|}_{\mathrm{2,}s}{\Vert v_h\Vert }_h$ (72)

$\left|R_3\left(v_h\right)\right|\lesssim h\left({\left|p\right|}_{\mathrm{1,}s}+{\left|p\right|}_{\mathrm{1,}d}\right){\Vert v_h\Vert }_h$ (73)

$\left|R_4\left(v_h\right)\right|\lesssim h\left({\left|p\right|}_{\mathrm{1,}d}\right){\Vert v_h\Vert }_h\mathrm{.}$ (74)

Proof.

1) Estimate (71): We begin with an estimate for the first term $R_1\left(v_h\right)$. For any face $E\in {\mathcal{E}}_h\left({\Omega }_s^{+}\right)$, there exists at least one element $T\in {\mathcal{T}}_h^s$ such that $E\in \mathcal{E}\left(T\right)$. Then, from condition (68), Höder’s inequality and inequality (70), it follows that

$\begin{array}{c}\left|R_1\left(v_h\right)\right|\le {\displaystyle \underset{j=1}{\overset{N-1}{\sum }}}\frac{\mu {\alpha }_1}{\sqrt{{\kappa }_j}}{\left({\displaystyle {\int }_{{\Gamma }_I}}{\left|u_s\cdot {\tau }_j\right|}^2\right)}^{1/2}{\left({\displaystyle {\int }_{{\Gamma }_I}}{\left|v_{h\mathrm{,}s}\cdot {\tau }_j\right|}^2\right)}^{1/2}\\ \le {\displaystyle \underset{j=1}{\overset{N-1}{\sum }}}\frac{\mu {\alpha }_1}{\sqrt{{\kappa }_j}}{\Vert u_s\Vert }_{\mathrm{1,}{\Omega }_s}{\left({\displaystyle {\int }_{{\Gamma }_I}}{\left|v_{h\mathrm{,}s}\cdot {\tau }_j\right|}^2\right)}^{\mathrm{1\; <mo>\; /\; </mo>2}}\\ \le \underset{1\le j\le N-1}{\max }\left(\frac{\mu {\alpha }_1}{\sqrt{{\kappa }_j}}\right)h{\Vert u_s\Vert }_{\mathrm{1,}{\Omega }_s}{\Vert v_h\Vert }_h\mathrm{.}\end{array}$

2) Estimate (72):

We have ${\left\{n_E\cdot D\left(u\right)\right\}}_{\mathrm{<mo>\; |\; </mo>}E}\in {\left[L^2\left(E\right)\right]}^N$, hence ${\Pi }_E\left(n_E\cdot D\left(u\right)\right)\in {\left[{\mathbb{P}}^0\left(E\right)\right]}^N$.

${\displaystyle {\int }_E}{\Pi }_E\left(n_E\cdot D\left(u\right)\right)\cdot {\left[v_h\right]}_E={\Pi }_E\left(n_E\cdot D\left(u\right)\right){\displaystyle {\int }_E}{\left[v_h\right]}_E=0.$ (75)

Thus,

$\begin{array}{c}{\displaystyle {\int }_E}{n}_E\cdot D\left(u\right)\cdot {\left[v_h\right]}_E={\displaystyle {\int }_E}\left(n_E\cdot D\left(u\right)-{\Pi }_E\left(n_E\cdot D\left(u\right)\right)\right)\cdot {\left[v_h\right]}_E\\ ={\displaystyle {\int }_E}\left(I-{\Pi }_E\right)\left(n_E\cdot D\left(u\right)\right)\cdot {\left[v_h\right]}_E\\ \lesssim {\Vert h_E^{1/2}\left(I-{\Pi }_E\right)\left(n_E\right)\cdot D\left(u\right)\Vert }_E{\Vert h_E^{-1/2}{\left[v_h\right]}_E\Vert }_E\\ \lesssim h_E{\left|D\left(u\right)\right|}_{\mathrm{1,}T}{h}_E^{-1/2}{\Vert {\left[v_h\right]}_E\Vert }_E\mathrm{.}\end{array}$

Furthermore, summing on $E\in {\mathcal{E}}_h\left({\Omega }_s^{+}\right)$ faces, we obtain the estimate:

$\left|R_2\left(v_h\right)\right|\lesssim h{\left|u\right|}_{\mathrm{2,}s}{\Vert v_h\Vert }_h\mathrm{.}$ (76)

3) For the terms $R_3\left(v_h\right)$ and $R_4\left(v_h\right)$, we use the same techniques as in the proof of the bounds for $R_i\left(v_h\right)$, $i\in \left\{\mathrm{1,2}\right\}$, and we obtain:

$\left|R_3\left(v_h\right)\right|\lesssim h\left({\left|p\right|}_{\mathrm{1,}s}+{\left|p\right|}_{\mathrm{1,}d}\right){\Vert v_h\Vert }_h,$

$\left|R_4\left(v_h\right)\right|\lesssim h\left({\left|p\right|}_{\mathrm{1,}d}\right){\Vert v_h\Vert }_h\mathrm{.}$

The proof is complete. o

From Lemma 3.5, Lemma 3.6 and Lemma 3.7, now we derive the following convergence theorem:

Theorem 3.2 Let the solution $\left(u\mathrm{,}p\right)$ of problem (27) satifies the smoothness assumption (Assumption 3.1). Let $\left(u_h\mathrm{,}{p}_h\right)$ be the solution of the discrete problem (36). Then there exists a positive constant C depending on $N,\mu ,K_{\mathrm{<mo>\; *\; </mo>}},K^{\mathrm{<mo>\; *\; </mo>}},{\alpha }_1$ and $\sigma$ such that:

${\Vert u-u_h\Vert }_h+\Vert p-p_h\Vert \le Ch\left({\left|u\right|}_{\mathrm{2,}s}+{\left|u\right|}_{\mathrm{2,}d}+{\left|p\right|}_{\mathrm{1,}s}+{\left|p\right|}_{\mathrm{1,}d}\right)\mathrm{.}$ (77)

4. Numerical Experiments

In this section we present one test case to verify the predicted rates of convergence. The numerical simulations have been performed on the finite element code FreeFem++ [33] [34] in isotropic coupled mesh of Figure 10.

The solutions have been represented by Mathematica software. For simplicity we choose each domain ${\Omega }_l$, $l\in \left\{s\mathrm{,}d\right\}$ as the unit square, ${\alpha }_1=\mu =1$, and the permeability tensor $K$ is taken to be the identity. The interface ${\Gamma }_I$, is the line $x=1$, i.e. $\bar{\Omega }={\bar{\Omega }}_s\cup {\bar{\Omega }}_d$ like the show Figure 11.

We consider the application $\varphi \mathrm{<mo>\; :\; </mo>}\left(x\mathrm{,}y\right)\in {ℝ}^2\mapsto \varphi \left(x\mathrm{,}y\right)=x^2{\left(x-1\right)}^3{y}^2{\left(y-1\right)}^2\in ℝ$ on the open domain $\Omega ={\Omega }_s\cup {\Omega }_d$. In $\Omega$, we define $u=\left(u_1,u_2\right)=\text{curl}\varphi =\left(-\frac{\partial \varphi }{\partial y},\frac{\partial \varphi }{\partial x}\right)$ and we obtain:

$u_1\left(x\mathrm{,}y\right)\mathrm{<mo>\; :\; </mo>}=-2{\left(-1+x\right)}^3{x}^2\left(-1+y\right)y\left(-1+2y\right)$ (78)

$u_2\left(x\mathrm{,}y\right)\mathrm{<mo>\; :\; </mo>}={\left(-1+x\right)}^{2}x\left(-2+5x\right){\left(-1+y\right)}^2{y}^2$ (79)

We choose quadratic pressure $p\in L^2\left(\Omega \right)$ by

$p\left(x\mathrm{,}y\right)=x^2-2xy+\frac{y^2}{2}-1.$ (80)

Thus,

Figure 10. Isotropic mesh ${\mathcal{T}}_{1/400}$ on coupled domain $\Omega \subset {ℝ}^2$.

Figure 11. The domain Ω in 2d.

${\displaystyle {\int }_{\Omega }}p\left(x,y\right)\text{d}x\text{d}y=0\text{and}\nabla p=\left(2x-2y,-2x+y\right).$ (81)

The exact solution $\left(u\mathrm{,}p\right)$ satisfies the following condition:

$\text{div}u=0=g\text{in}\Omega ,$ (82)

$u=0\text{on}\partial \Omega ,$ (83)

and the Beavers-Joseph-Saffman interface conditions on ${\Gamma }_I$ [ ${\Gamma }_I:x=1$ ]:

$u_s\cdot n_s+u_d\cdot n_d=0\text{on}{\Gamma }_I,$ (84)

$p_s-2\mu n_s\cdot D\left(u_s\right)\cdot n_s=p_d\text{on}{\Gamma }_I,$ (85)

$\frac{\sqrt{k_j}}{{\alpha }_1}2n_s\cdot D\left(u_s\right)\cdot {\tau }_j=-u_s\cdot {\tau }_j\text{on}{\Gamma }_I,j=1,\cdots ,N-1.$ (86)

Furthermore, we obtain the right-hand terms $f$ define by

$\{\begin{array}{l}{f}_s=-2\mu \text{div}D\left(u\right)+\nabla p\text{in}{\Omega }_s\mathrm{,}\\ f_d=\mu K^{-1}u+\nabla p\text{in}{\Omega }_d.\end{array}$ (87)

Thus, $f_s\left(x\mathrm{,}y\right)=\left(f_1\left(x\mathrm{,}y\right)\mathrm{,}{f}_2\left(x\mathrm{,}y\right)\right)$ in ${\Omega }_s$ leads to

$\begin{array}{c}{f}_1\left(x,y\right)=4\left(-1+x\right)\left(-1+2y\right)(-6x^3+3x^4+\left(-1+y\right)y-8x\left(-1+y\right)y\\ +x^2\left(3+10\left(-1+y\right)y\right))+2x-2y,\end{array}$

$\begin{array}{c}{f}_2\left(x,y\right)=-2(9{\left(-1+y\right)}^2{y}^2-12x^3\left(1+6\left(-1+y\right)y\right)+5x^4\left(1+6\left(-1+y\right)y\right)\\ -2x\left(1+6\left(-1+y\right)y\left(1+3\left(-1+y\right)y\right)\right)\\ +x^2\left(9+6\left(-1+y\right)y\left(9+5\left(-1+y\right)y\right)\right))-2x+y,\end{array}$

and in ${\Omega }_d$, $f_d\left(x\mathrm{,}y\right)=\left(k_1\left(x\mathrm{,}y\right)\mathrm{,}{k}_2\left(x\mathrm{,}y\right)\right)$ is given by:

$k_1\left(x\mathrm{,}y\right)={\left(-1+x\right)}^{2}x\left(-2+5x\right){\left(-1+y\right)}^2{y}^2+2x-2y\mathrm{,}$

$k_2\left(x\mathrm{,}y\right)={\left(-1+x\right)}^{2}x\left(-2+5x\right){\left(-1+y\right)}^2{y}^2-2x+y\mathrm{.}$

Figure 12 and Figure 13 give an illustration of isotropic and anisotropic meshes.

We study now the convergence of the speed and the pressure in each subdomain on a triangulation ${\mathcal{T}}_{\mathrm{1\; <mo>\; /\; </mo>200}}$ for the non-conforming finite elements of Crouzeix-Raviart. The results obtained are consistent, as shown in Figures 14-17. The order of convergence in norm $L^2$ approximately equal to 2 for the speed in each domain ${\Omega }_l$, $l\in \left\{s\mathrm{,}d\right\}$ and of order 1 for the pressure. Then, we represent on the same triangulation, the curves of isovalours of each component of the speed in ${\Omega }_l$, $l\in \left\{s\mathrm{,}d\right\}$ (see Figures 18-21).

Finally, we study the structure of the stiffness matrix in $\Omega$ on a uniform triangulation. We find that, when the discretization step h becomes more and more infinitely small, the matrix becomes more and more sparse (cf. Table 1), which partly shows the effectiveness of the method digital used.

The exact solutions $u=\left(u_1\mathrm{,}{u}_2\right)$ and p are represented by Figures 22-24, while the second members $f_1$, $f_2$, $k_1$ and $k_2$ are represented respectively by Figures 25-28.

Figure 12. Example of isotropic mesh ${\mathcal{T}}_{\mathrm{1\; <mo>\; /\; </mo>200}}$ in 2d.

Figure 13. Example of anisotropic mesh in 2d.

Figure 14. Error for the velocity ${\Vert u-u_h\Vert }_h$ in ${\Omega }_s$ (log/log plot).

Figure 15. Error for the pressure $\Vert p-p_h\Vert$ in ${\Omega }_s$ (log/log plot).

Figure 16. Error for the velocity ${\Vert u-u_h\Vert }_h$ in ${\Omega }_d$ (log/log plot).

Figure 17. Error for the pressure $\Vert p-p_h\Vert$ in ${\Omega }_d$ (log/log plot).

Figure 18. The isovalue of the first velocity component $u_1$ in ${\Omega }_s$.

Figure 19. The isovalue of the second velocity component $u_2$ in ${\Omega }_s$.

Figure 20. The isovalue of the first velocity component $u_1$ in ${\Omega }_d$.

Figure 21. The isovalue of the second velocity component $u_2$ in ${\Omega }_d$.

Table 1. Structure of rigidity Matrix on 10 iterations.

Figure 22. Component $u_1$ in $\Omega$.

5. Summary

· In this contribution, we have investigated a new mixed finite element method to solve the Stokes-Darcy fluid flow model without introducing any Lagrange multiplier. We have proposed a modification of the Darcy problem which allows us to apply a slight variant nonconforming Crouzeix-Raviart element to

Figure 23. Component $u_2$ in $\Omega$.

Figure 24. Pressure p in $\Omega$.

Figure 25. Right-hand term $f_1$ in ${\Omega }_s$.

Figure 26. Right-hand term $f_2$ in ${\Omega }_s$.

Figure 27. Right-hand term $k_1$ in ${\Omega }_d$.

Figure 28. Right-hand term $k_2$ in ${\Omega }_d$.

the whole coupled Stokes-Darcy problem. The proposed method is probably the cheapest method for Discontinuous Galerkin (DG) approximation of the coupled system, has optimal accuracy with respect to solution regularity, and has simple and straightforward implementations. Numerical experiments have been also presented, which confirm the excellent stability and accuracy of our method.

· Further works: Further it is well known that an internal layer appears at the interface ${\Gamma }_I$ as the permeability tensor degenerates, in that case anisotropic meshes have to be used in this layer (see for instance [8] [35] ). Hence we intend to extend our results to such anisotropic meshes.

6. Nomenclatures

· $\Omega \subset {ℝ}^N\mathrm{,}N\in \left\{\mathrm{2,3}\right\}$ bounded domain

· ${\Omega }_d$ : the porous medium domain

· ${\Omega }_s=\Omega \backslash {\bar{\Omega }}_d$

· ${\Gamma }_I=\partial {\Omega }_s\cap \partial {\Omega }_d$

· ${\Gamma }_l=\partial {\Omega }_l\backslash {\Gamma }_I,l=s,d$

· $n_s$ (resp. $n_d$ ) the unit outward normal vector along $\partial {\Omega }_s$ (resp. $\partial {\Omega }_d$ )

· $u$ : the fluid velocity

· p: the fluid pressure

· In 2D, the curl of a scalar function w is given as usual by

$\text{curl}w\mathrm{<mo>\; :\; </mo>}={\left(\frac{\partial w}{\partial x_2}\mathrm{,}-\frac{\partial w}{\partial x_1}\right)}^{{\rm T}}$

· In 3D, the curl of a vector function $w=\left(w_1\mathrm{,}{w}_2\mathrm{,}{w}_3\right)$ is given as usual by $\text{curl}w\mathrm{<mo>\; :\; </mo>}=\nabla \times w$ namely,

$\text{curl}w\mathrm{<mo>\; :\; </mo>}=\left(\frac{\partial w_3}{\partial x_2}-\frac{\partial w_2}{\partial x_3}\mathrm{,}\frac{\partial w_1}{\partial x_3}-\frac{\partial w_3}{\partial x_1}\mathrm{,}\frac{\partial w_2}{\partial x_1}-\frac{\partial w_1}{\partial x_2}\right)$

· ${\mathbb{P}}^k$ : the space of polynomials of total degree not larger than k

· ${\mathcal{T}}_h$ : triangulation of $\Omega$

· ${\mathcal{T}}_h^l$ : the corresponding induced triangulation of ${\Omega }_l$, $l\in \left\{s\mathrm{,}d\right\}$

· For any $T\in {\mathcal{T}}_h$, $h_T$ is the diameter of T and ${\rho }_T=2r_T$ is the diameter of the largest ball inscribed into T

· $h\mathrm{<mo>\; :\; </mo>}=\underset{T\in {\mathcal{T}}_h}{\max }{h}_T$ and ${\sigma }_h\mathrm{<mo>\; :\; </mo>}=\underset{T\in {\mathcal{T}}_h}{\max }\frac{h_T}{{\rho }_T}$

· ${\mathcal{E}}_h$ : the set of all the edges or faces of the triangulation

· $\mathcal{E}\left(T\right)$ : the set of all the edges ( $N=2$ ) or faces ( $N=3$ ) of a element T

· ${\mathcal{E}}_h\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{T\in {\mathcal{T}}_h}{\cup }\mathcal{E}(\; T\; )}$

· $\mathcal{N}\left(T\right)$ : the set of all the vertices of a element T

· ${\mathcal{N}}_h:={\displaystyle \underset{T\in {\mathcal{T}}_h}{\cup }\mathcal{N}(\; T\; )}$

· For $\mathcal{A}\subset \bar{\Omega }$, ${\mathcal{E}}_h\left(\mathcal{A}\right)\mathrm{<mo>\; :\; </mo>}=\left\{E\in {\mathcal{E}}_h\mathrm{<mo>\; :\; </mo>}E\subset \mathcal{A}\right\}$

· For $E\in {\mathcal{E}}_h$, we associate a unit vector $n_E$ such that $n_E$ is orthogonal to E and equals to the unit exterior normal vector to $\partial \Omega$

· For $E\in {\mathcal{E}}_h$, ${\left[\varphi \right]}_E$ is the jump across E in the direction of $n_E$

· In order to avoid excessive use of constants, the abbreviations $x\lesssim y$ and $x\sim y$ stand for $x\le cy$ and $c_{1}x\le y\le c_{2}x$, respectively, with positive constants independent of $x,y$ or ${\mathcal{T}}_h$.

Acknowledgements

The author thanks Professor Emmanuel Creusé (University of Lille 1, France) for having sent us useful documents and for fruitful discussions concerning the numerical tests. We are thankful to the editor and the anonymous reviewers for many valuable suggestions to improve this paper. This article has been posted since April 05, 2019 according to the link https://arxiv.org/abs/1908.01892.

Data Availability

There are no data underlying the findings in this paper to be shared.

Funding

This work has received no funding.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Discacciati, M. and Quarteroni, A. (2009) Navier-Stokes/Darcy Coupling: Modeling, Analysis, and Numerical Approximation. Revista Matemática Complutense, 22, 315-426. https://doi.org/10.5209/rev_REMA.2009.v22.n2.16263
[2]	Beavers, G. and Joseph, D. (1967) Boundary Conditions at a Naturally Permeable Wall. Journal of Fluid Mechanics, 30, 197-207. https://doi.org/10.1017/S0022112067001375
[3]	Saffman, P. (1971) On the Boundary Condition at the Interface of a Porous Medium. Studies in Applied Mathematics, 1, 93-101. https://doi.org/10.1002/sapm197150293
[4]	Jäger, W. and Mikelić, A. (1996) On the Boundary Conditions of the Contact Interface between a Porous Medium and a Free Fluid. Annali della Scuola Normale Superiore. Classe di Scienze, 23, 403-465.
[5]	Jäger, W. and Mikelić, A. (2000) On the Interface Boundary Condition of Beavers, Joseph and Saffman. SIAM Journal on Applied Mathematics, 60, 1111-1127. https://doi.org/10.1137/S003613999833678X
[6]	Jäger, W., Mikelić, A. and Neuss, N. (2001) Asymptotic Analysis of the Laminar Visous Flow over a Porous Bed. SIAM Journal on Scientific Computing, 22, 2006-2028. https://doi.org/10.1137/S1064827599360339
[7]	Payne, L. and Straughan, B. (1998) Analysis of the Boundary Condition at the Interface between a Viscous Fluid and a Porous Medium and Related Modeling Questions. Journal de Mathématiques Pures et Appliquées, 77, 317-354. https://doi.org/10.1016/S0021-7824(98)80102-5
[8]	Houédanou, K.W. and Ahounou, B. (2016) A Posteriori Error Estimation for the Stokes-Darcy Coupled Problem on Anisotropic Discretization. Mathematical Methods in the Applied Sciences, 40, 3741-3774. https://doi.org/10.1002/mma.4261
[9]	Vassilev, D. and Yotov, I. (2009) Coupling Stokes-Darcy Flow with Transport. SIAM Journal on Scientific Computing, 5, 3661-3684. https://doi.org/10.1137/080732146
[10]	Nicaise, S., Ahounou, B. and Houédanou, W. (2015) A Residual-Based Posteriori Error Estimates for a Nonconforming Finite Element Discretization of the Stokes-Darcy Coupled Problem: Isotropic Discretization. Afrika Matematika, 27, 701-729. https://doi.org/10.1007/s13370-015-0370-3
[11]	William, L.J., Friedhelm, S. and Ivan, Y. (2002) Coupling Fluid Flow with Porous Media Flow. SIAM Journal on Numerical Analysis, 40, 2195-2218. https://doi.org/10.1137/S0036142901392766
[12]	Mardal, K.-A., Tai, X. and Winther, R. (2002) A Robust Finite Element Method for Darcy-Stokes Flow. SIAM Journal on Numerical Analysis, 40, 1605-1631. https://doi.org/10.1137/S0036142901383910
[13]	Mu, M. and Xu, J. (2007) A Two-Grid Method of a Mixed Stokes-Darcy Model for Coupling Fluid Flow with Porous Media Flow. SIAM Journal on Numerical Analysis, 45, 1801-1813. https://doi.org/10.1137/050637820
[14]	Rui, H. and Zhang, R. (2009) A Unified Stabilized Mixed Finite Element Method for Coupling Stokes and Darcy Flows. Computer Methods in Applied Mechanics and Engineering, 198, 2692-2699. https://doi.org/10.1016/j.cma.2009.03.011
[15]	Arbogast, T. and Brunson, D. (2007) A Computational Method for Approximating A Darcy-Stokes System Governing a Vuggy Porous Medium. Computational Geosciences, 11, 207-218. https://doi.org/10.1007/s10596-007-9043-0
[16]	Galvis, J. and Sarkis, M. (2007) Nonconforming Mortar Discretization Analysis for the Coupling Stokes-Darcy Equations. Electronic Transactions on Numerical Analysis, 26, 350-384.
[17]	Gatica, G.-N., Oyarzùa, R. and Sayas, F.-J. (2011) Convergence of a Family of Galerkin Discretizations for the Stokes-Darcy Coupled Problem. Numerical Methods for Partial Differential Equations, 27, 721-748.
[18]	Rivière, B. and Yotov, I. (2005) Locally Conservative Coupling of Stokes and Darcy Flows. SIAM Journal on Numerical Analysis, 42, 1959-1977. https://doi.org/10.1137/S0036142903427640
[19]	Chen, W. and Wang, Y. (2014) A Posteriori Error Estimate for H(div) Conforming Mixed Finite Element for the Coupled Darcy-Stokes System. Journal of Computational and Applied Mathematics, 255, 502-516. https://doi.org/10.1016/j.cam.2013.05.021
[20]	Babuška, I. and Gatica, G. (2010) A Residual-Based a Posteriori Error Estimator for the Stokes-Darcy Coupled Problem. SIAM Journal on Numerical Analysis, 48, 498-523. https://doi.org/10.1137/080727646
[21]	Gatica, G., Oyarzùa, R. and Sayas, F.-J. (2011) A Residual-Based a Posteriori Error Estimator for a Fully-Mixed Formulation of the Stokes-Darcy Coupled Problem. Computer Methods in Applied Mechanics and Engineering, 200, 1877-1891. https://doi.org/10.1016/j.cma.2011.02.009
[22]	Armentano, M.G. and Stockdale, M.L. (2018) A Unified Mixed Finite Element Approximations of the Stokes-Darcy Coupled Problem. Computers and Mathematics with Applications, 77, 2568-2584. https://doi.org/10.1016/j.camwa.2018.12.032
[23]	Yu, J., Mahbub, M.A.A., Shi, F. and Zheng, H. (2018) Stabilized Finite Element Method for the Stationary Mixed Stokes-Darcy Problem. Advances in Difference Equations, 2018, Article No. 346. https://doi.org/10.1186/s13662-018-1809-2
[24]	Li, R., Li, J., Chen, Z. and Gao, Y. (2015) A Stabilized Finite Element Method Based on Two Local Gauss Integrations for a Coupled Stokes-Darcy Problem. Journal of Computational and Applied Mathematics, 292, 92-104. https://doi.org/10.1016/j.cam.2015.06.014
[25]	Houédanou, K.W., Adetola, J. and Ahounou, B. (2017) Residual-Based a Posteriori Error Estimates for a Conforming Finite Element Discretization of the Navier-Stokes/Darcy Coupled Problem. Journal of Pure and Applied Mathematics: Advances and Applications, 18, 37-73. https://doi.org/10.18642/jpamaa_7100121867
[26]	Gatica, G.N., Meddahi, S. and Oyarzùa, R. (2009) A Conforming Mixed Finite Element Method for the Coupling of Fluid Flow with Porous Media Flow. IMA Journal of Numerical Analysis, 29, 86-108. https://doi.org/10.1093/imanum/drm049
[27]	Cui, M. and Yan, N. (2011) A Posteriori Error Estimate for the Stokes-Darcy System. Mathematical Methods in the Applied Sciences, 34, 1050-1064. https://doi.org/10.1002/mma.1422
[28]	Girault, V. and Raviart, P.-A. (1986) Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms. Springer Series in Computational Mathematics, Volume 5, Springer, Berlin. https://doi.org/10.1007/978-3-642-61623-5
[29]	Ern, A. (2005) Aide-mémoire eléments finis. Dunod, Paris.
[30]	Brenner, S. (2003) Korn’s Inequalities for Piecewise H1 Vector Fields. Mathematics of Computation, 73, 1067-1087. https://doi.org/10.1090/S0025-5718-03-01579-5
[31]	Brezzi, F. and Fortin, M. (1991) Mixed and Hybrid Finite Element Methods. Springer-Verlag, New York. https://doi.org/10.1007/978-1-4612-3172-1
[32]	Crouzeix, M. and Raviart, P. (1973) Conforming and Nonconforming Finite Element Methods for Solving the Stationary Stokes Equations I. ESAIM: Mathematical Modelling and Numerical Analysis, 7, 33-75. https://doi.org/10.1051/m2an/197307R300331
[33]	Hecht, F. (1998) The Mesh Adapting Software: BAMG. INRIA Report. http://wwwrocq.inria.fr/gamma/cdrom/www/bamg/eng.htm
[34]	Frederic, H. and Olivier, P. (2012) Freefem++. http://www.freefem.org
[35]	Houédanou, K.W. (2015) Analyse d’erreur a-posteriori pour quelques méthodes d’éléments finis mixtes pour le problème de transmission Stokes-Darcy: Discrétisations isotrope et anisotrope. Thèse de Doctorat, Université d’Abomey-Calavi, Godomey, 210 p. http://hal.archives-ouvertes.fr/tel-01373344