A new unified stabilized mixed finite element method of the Stokes-Darcy coupled problem: Isotropic discretization

TL;DR

This paper introduces a new unified mixed finite element method for coupled Stokes-Darcy flows on isotropic meshes, providing stability, convergence analysis, and numerical validation for fluid-porous media interactions.

Contribution

It develops a novel finite element approach combining a modified Darcy problem with a nonconforming element for the coupled flow, along with theoretical analysis and numerical validation.

Findings

Method demonstrates excellent stability.

Convergence analysis confirms accuracy.

Numerical experiments validate theoretical results.

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 $\mathbb{R}^N$, $N\in\{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.

Tables

Table 1. Table 1. Structure of rigidity Matrix on 10 10 10 iterations.

Parameter $h$	${}_{}{}^o/_o^{}$ of Coefs. $\ne$ 0
1/16	79.668827
1/22	66.405634
1/28	41.751001
1/34	28.652204
1/40	20.873193
1/46	15.879941
1/52	12.485349
1/58	10.073294
1/66	08.298201
1/70	06.954061

Equations

$$\displaystyle\left\{\begin{array}[]{ccccccccc}-2\mu\mathop{\rm div}\nolimits\textbf{D}(\textbf{u})+\nabla p&=&\textbf{f}&\mbox{ in }&&\Omega_{s},&\\ \mathop{\rm div}\nolimits\textbf{u}&=&g&\mbox{ in }&&\Omega_{s},&\\ \textbf{u}&=&\textbf{0}&\mbox{ on }&&\Gamma_{s},&\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{cccccccccccc}\mu\textbf{K}^{-1}\mathbf{u}+\nabla p&=&\textbf{f}&\mbox{ in }&&\Omega_{d},&\\ \mathop{\rm div}\nolimits\textbf{u}&=&g&\mbox{ in }&&\Omega_{d},&\\ \textbf{u}\cdot\textbf{n}_{d}&=&0&\mbox{ on }&&\Gamma_{d}.&\end{array}\right.$$

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

$$K_{\ast }{\xi }^T\xi ⩽{\xi }^T\text{K}(x)\xi ⩽K^{\ast }{\xi }^T\xi ,\text{mbox}\forall x\in {\Omega }_d,\text{mbox}\xi \in {\mathbb{R}}^N.$$

$${\int }_{\Omega }g(x)dx=0.$$

$${\text{u}}_s\cdot {\text{n}}_s+{\mathbf{u}}_d\cdot {\text{n}}_d$$

$$p_s-2\mu {\text{n}}_s\cdot \text{D}({\text{u}}_s)\cdot {\text{n}}_s$$

$$\frac{\sqrt{{\kappa }_j}}{{\alpha }_1}2{\text{n}}_s\cdot \text{D}({\text{u}}_s)\cdot {\tau }_j$$

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

$$⟨\lambda ,\eta {⟩}_{\Gamma }$$

$$\text{H}(\mathrm{div},\Omega )$$

$$:=$$

$$\parallel \text{v}{\parallel }_{\text{H}}$$

$$Q=L_0^2(\Omega ):=\left\{q\in L^2(\Omega ):{\int }_{\Omega }q(x)dx=0\right\}.$$

$$({\text{f}}_s,{\text{v}}_s{)}_{{\Omega }_s}$$

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

$$({\text{f}}_s,\text{v}{)}_{{\Omega }_s}$$

$$-{\left(q,\text{mbox}div\text{mbox}{\text{u}}_s\right)}_{{\Omega }_s}$$

$${\left(\mu {\text{K}}^{-1}{\text{u}}_d,\text{v}\right)}_{{\Omega }_d}$$

$${\left(\text{mbox}div\text{mbox}{\text{u}}_d,q\right)}_{{\Omega }_d}$$

$$\left\{\begin{array}[]{ccc}\textbf{a}(\textbf{u},\textbf{v})+\textbf{b}(\textbf{v},p)&=&L(\textbf{v}),\hskip 28.45274pt\forall\textbf{v}\in\textbf{H},\\ \textbf{b}(\textbf{u},q)&=&G(q),\hskip 28.45274pt\forall q\in Q.\end{array}\right.$$

$$\text{a}(\text{u},\text{v})$$

$$\text{b}(\text{v},q)$$

$$L(\text{v}):=(\text{f},\text{v}{)}_{{\Omega }_s}+(\text{f},\text{v}{)}_{{\Omega }_d}\text{mbox}andG(q):=-(g,q{)}_{{\Omega }_s}-(g,q{)}_{{\Omega }_d}.$$

$${\int }_{{\Omega }_d}(\mathrm{div}{\text{u}}_d-g_d)\mathrm{div}\text{v}=0.$$

$${\left(\mu {\text{K}}^{-1}{\text{u}}_d,\text{v}\right)}_{{\Omega }_d}+{\left(\mathrm{div}{\text{u}}_d,\mathrm{div}\text{v}\right)}_{{\Omega }_d}-{\left(p_d,\mathrm{div}\text{v}\right)}_{{\Omega }_d}$$

$$+(p_d,{\text{v}}_d\cdot {\text{n}}_d{)}_{{\Gamma }_I}$$

$${\left(\text{mbox}div\text{mbox}{\text{u}}_d,q\right)}_{{\Omega }_d}$$

$$\left\{\begin{array}[]{ccc}\tilde{\textbf{a}}(\textbf{u},\textbf{v})+\textbf{b}(\textbf{v},p)&=&\tilde{L}(\textbf{v}),\hskip 28.45274pt\forall\textbf{v}\in\textbf{H},\\ \textbf{b}(\textbf{u},q)&=&G(q),\hskip 28.45274pt\forall q\in Q.\end{array}\right.$$

$$\tilde{\text{a}}(\text{u},\text{v})=2\mu (\text{D}(\text{u}),\text{D}(\text{v}){)}_s+\sum _{j=1}^{N-1}\frac{\mu {\alpha }_1}{\sqrt{{\kappa }_j}}{⟨{\text{u}}_s\cdot {\tau }_j,{\text{v}}_s\cdot {\tau }_j⟩}_{{\Gamma }_I}+\mu {\left({\text{K}}^{-1}\text{u},\text{v}\right)}_d+{\left(\mathrm{div}{\text{u}}_d,\mathrm{div}\text{v}\right)}_{{\Omega }_d}$$

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

$$\tilde{L}(\text{v}):=(\text{f},\text{v}{)}_{{\Omega }_s}+(\text{f},\text{v}{)}_{{\Omega }_d}+(\mathrm{div}\text{u},\mathrm{div}\text{v}{)}_{{\Omega }_D}\text{mbox}andG(q):=-(g,q{)}_{{\Omega }_s}-(g,q{)}_{{\Omega }_d}.$$

$$\parallel \text{u}{\parallel }_{\text{H}}+\parallel p{\parallel }_Q\le \tilde{C}\left(\parallel {\text{f}}_s{\parallel }_{{\Omega }_s}+\parallel {\text{f}}_d{\parallel }_{{\Omega }_d}+\parallel g_d{\parallel }_{{\Omega }_d}+\parallel g_s{\parallel }_{{\Omega }_s}\right).$$

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

Full text

**Keywords : ** Coupled Stokes and Darcy flows; Nonconforming finite element method; Crouzeix-Raviart element.

Mathematics Subject Classification [MSC]: 74S05, 74S10, 74S15, 74S20, 74S25, 7430.

1. 1 Introduction

2. 2 Preliminaries and notation

3. 2.1 Model problem

4. 2.2 New weak formulation

5. 3 A priori error analysis

6. 3.1 Finite element discretization

7. 3.2 A convergence analysis

8. 4 Numerical experiments

9. 5 Conclusion

10. 6 Acknowledgment

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 making its way into the water supply. This coupling is also important in technological applications involving filtration. We refer to the nice overview [9] 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 treatement 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 [4], modified by Saffman in [31], and later mathematically justified in [20, 21, 23, 28].

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 [19, 32, 27, 33, 25, 26, 30, 1, 12, 15, 20, 23, 22, 29, 6, 3, 13, 2, 34, 24, 18], 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, [14, 8, 13] 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 [2, 30, 27] and the references therein).

In [2], 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 conditions, obtaining as a result optimal accuracy with respect to solution regularity. Numerical experimants are also presented, which confirm the excellent stability and optimal performance of our method. The difference between our paper and the reference [2] is that our discretization is nonconforming in both the Stokes domain and Darcy domain (in $\Omega\subset\mathbb{R}^{N}$, $N=2\mbox{ or }3$). As a result, additional terms are included in the priori error analysis that measure 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(\mathop{\rm div}\nolimits)$. Thus, we introduce a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element space (larger than the space $\textbf{H}_{h}$ used in [30]). The choice of $\textbf{H}_{h}$ [see (38)] is more natural than the one introduced in [30] since the space $\textbf{H}_{h}$ approximates only $H(\mathop{\rm div}\nolimits,\Omega_{d})$ and not $[H^{1}(\Omega_{d})]^{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\mathbb{R}^{N}$, $N=2\mbox{ or }3$, notations and the weak formulation. Section 3 is devoted to the finite element discretization and the error estimation. Finally, in Section 4, we present the results of numerical experiments to verify the predicted rates of convergence.

2. Preliminaries and notation

2.1. Model problem

We consider the model of a flow in a bounded domain $\Omega\subset\mathbb{R}^{N}$ $(N=2\mbox{ or }3)$, consisting of a porous medium domain $\Omega_{d}$, where the flow is a Darcy flow, and an open region $\Omega_{s}=\Omega\smallsetminus\overline{\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}\smallsetminus\Gamma_{I}$, $l=s,d$. Each interface and boundary is assumed to be polygonal $(N=2)$ or polyhedral $(N=3)$. We denote by $\textbf{n}_{s}$ (resp. $\textbf{n}_{d}$) the unit outward normal vector along $\partial\Omega_{s}$ (resp. $\partial\Omega_{d}$). Note that on the interface $\Gamma_{I}$, we have $\textbf{n}_{s}=-\textbf{n}_{d}$. The Figures 2 and 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}}$.

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

[TABLE]

while in the porous medium $\Omega_{d}$, by Darcy’s law

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Finally we consider the following interface conditions on $\Gamma_{I}:$

[TABLE]

Here, Eq. (9) represents mass conservation, Eq. (10) the balance of normal forces, and Eq. (11) the Beavers-Joseph-Saffman conditions. Moreover, $\{\tau_{j}\}_{j=1,\ldots,N-1}$ denotes an orthonormal system of tangent vectors on $\Gamma_{I}$, $\kappa_{j}=\tau_{j}\cdot\textbf{K}\cdot\tau_{j}$, and $\alpha_{1}$ is a parameter determined by experimental evidence.

Eqs. (4) to (11) 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 $\mathbb{R}^{N}$ and $m$ is a non negative integer, the Sobolev space $H^{m}(W)=W^{m,2}(W)$ is defined in the usual way with the usual norm $\parallel\cdot\parallel_{m,W}$ and semi-norm $|\cdot|_{m,W}$. In particular, $H^{0}(W)=L^{2}(W)$ and we write $\parallel\cdot\parallel_{W}$ for $\parallel\cdot\parallel_{0,W}$. Similarly we denote by $(\cdot,\cdot)_{W}$ the $L^{2}(W)$ $[L^{2}(W)]^{N}$ or $[L^{2}(W)]^{N\times N}$ inner product. For shortness if $W$ is equal to $\Omega$, we will drop the index $\Omega$, while for any $m\geq 0$, $\parallel\cdot\parallel_{m,l}=\parallel\cdot\parallel_{m,\Omega_{l}}$, $|\cdot|_{m,l}=|\cdot|_{m,\Omega_{l}}$ and $(.,.)_{l}=(\cdot,\cdot)_{\Omega_{l}}$, for $l=s,d$. The space $H_{0}^{m}(\Omega)$ denotes the closure of $C_{0}^{\infty}(\Omega)$ in $H^{m}(\Omega)$. Let $[H^{m}(\Omega)]^{N}$ be the space of vector valued functions $\textbf{v}=(v_{1},\ldots,v_{N})$ with components $v_{i}$ in $H^{m}(\Omega)$. The norm and the seminorm on $[H^{m}(\Omega)]^{N}$ are given by

[TABLE]

For a connected open subset of the boundary $\Gamma\subset\partial\Omega_{s}\cup\partial\Omega_{d}$, we write $\langle.,.\rangle_{\Gamma}$ for the $L^{2}(\Gamma)$ inner product (or duality pairing), that is, for scalar valued functions $\lambda$, $\eta$ one defines:

[TABLE]

We also define the special vector-valued functions space

[TABLE]

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

[TABLE]

equipped with the norm

[TABLE]

and

[TABLE]

Multiplying the first equation of (4) by a test fonction $\bf v\in\bf H$ and the second one by $q\in Q$, integrating by parts over $\Omega_{s}$ the terms involving $\mathop{\rm div}\nolimits\textbf{D}(\textbf{u})$ and $\nabla p$, yield the variational form of Stokes equations:

[TABLE]

Using interface conditions $(\ref{cd2})$ and $(\ref{cd3})$ in (2.2), we obtain:

[TABLE]

We apply a similar treatment to the Darcy equations by testing the first equation of (8) with a smooth fonction $\textbf{v}\in\textbf{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:

[TABLE]

Now, incorporating the first boundary interface condition (9) and taking into account that the vector valued functions in H have (weakly) continuous normal components on $\Gamma_{I}$ (see [16, Theorem 2.5]), the mixed variational formulation of the coupled problem (4)-(11) can be stated as follows [27]: Find $(\textbf{u},p)\in\textbf{H}\times Q$ that satisfies

[TABLE]

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

[TABLE]

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

[TABLE]

It is easy to prove that a et b are continuous, b satisfies the continuous inf-sup condtion and a is coercive on the null space of b. It is also clear that $F$ and $G$ are continuous and bounded. Then, using the classical theory of mixed methods (see, e.g., [16, Theorem and Corollary 4.1 in Chapter I]) it follows the well-posedness of the continuous formulation (23) and so the following theorem holds [27]:

Theorem 2.1**.**

*If $\textbf{f}\in[L^{2}(\Omega)]^{N}$ and $g\in L_{0}^{2}(\Omega)$, there exists a unique solution $(\textbf{u},p)\in\textbf{H}\times Q$ to the problem (23).

Remark 2.1**.**

Note that if $g$ is of mean zero, (23) directly implies that (4), (8) and (9) hold ( the differential equations being understood in the distributional sense), while the interface conditions (10) and (11) are imposed in a weak sense. Also, we observe that the mixed variational formulation of the coupled problem (4)-(11) is equivalent to weak formulation (2.4) (and also (2.5) of [33]), with the particularity that, in our case, for any $\textbf{v}\in\textbf{H}$, we have that $\left<\textbf{v}_{s}-\textbf{v}_{d},\textbf{n}_{s}p_{s}\right>_{\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 [2]. Indeed, we observe that taking the second equation of Darcy’ problem (8) we can write, for any $\textbf{v}\in\textbf{H},$

[TABLE]

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

[TABLE]

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

In the same way that before, incorporating the boundary conditions (9) and remambering that, since $\textbf{v}\in\textbf{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 $(\textbf{u},p)\in\textbf{H}\times Q$ satisfying

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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 $(\textbf{u},p)\in\textbf{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{\textbf{a}}$ and the boundedness constants for $\tilde{\textbf{a}}$ and b, such that:

[TABLE]

We end this section with some notation. In $2D$, the $\mathop{\rm curl}\nolimits$ of a scalar function $w$ is given as usual by $\mathop{\rm curl}\nolimits w:=(\frac{\partial w}{\partial x_{2}},-\frac{\partial w}{\partial x_{1}})^{\top}$ while in $3D$, the $\mathop{\rm curl}\nolimits$ of a vector function w is given as usual by $\mathop{\rm curl}\nolimits\textbf{w}:=\nabla\times\textbf{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\leqslant cy$ and $c_{1}x\leqslant y\leqslant 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

[TABLE]

We assume that the family of triangulations is regular, in the sense that there exists $\sigma_{0}>0$ such that $\sigma_{h}\leqslant\sigma_{0}$, 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 Fig. 5, 5, 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}(T)$ (resp. ${\mathcal{N}}(T))$ the set of its edges $(N=2)$ or faces $(N=3)$ (resp. vertices) and set $\mathcal{E}_{h}=\displaystyle\bigcup_{T\in\mathcal{T}_{h}}\mathcal{E}(T)$, ${\mathcal{N}}_{h}=\displaystyle\bigcup_{T\in\mathcal{T}_{h}}{\mathcal{N}}(T)$. For $\mathcal{A}\subset\overline{\Omega}$ we define

[TABLE]

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

[TABLE]

where $\Omega_{s}^{+}=\Omega_{s}\cup\Gamma_{s}.$ Note that $\mathcal{E}_{h}(\Gamma_{I})$ is included in $\mathcal{E}_{h}(\partial\Omega_{d})$.

With every edges $E\in\mathcal{E}_{h}$, we associate a unit vector $\textbf{n}_{E}$ such that $\textbf{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 $\varphi$, we denote by $[\varphi]_{E}$ its jump across $E$ in the direction of $\textbf{n}_{E}$:

[TABLE]

For $i\in\{0,\cdots,N\}$, we set:

[TABLE]

The triplet $\{T,\mathbb{P}^{1}(T),\Sigma\}$ with $\Sigma=\{\sigma_{i}\}_{0\hskip 2.84544pt\leqslant i\hskip 2.84544pt\leqslant N}$ is finite element [10, Page 83]. The local basis functions are defined by:

[TABLE]

where for each $i\in\{0,\cdots,N\}$, $\lambda_{i}(T)$ is barycentric coordonates of $T\in\mathcal{T}_{h}$.

In classical reference element $\overline{T}$, the basis fonctions are given by:

[TABLE]

Based on the above notation, we introduce a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element space (larger than the space $\textbf{H}_{h}$ used in [30])

[TABLE]

and piecewise constant function space

[TABLE]

where $\mathbb{P}^{m}(T)$ 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 $\parallel\cdot\parallel$ while the norm on $\textbf{H}_{h}$ will be specified later on. The choice of $\textbf{H}_{h}$ is more natural than the one introduced in [30] since the space $\textbf{H}_{h}$ approximates only $H(\mathop{\rm div}\nolimits,\Omega_{d})$ and not $[H^{1}(\Omega_{d})]^{N}$, while our a priori error analysis is only valid in this larger space.

Let us introduce the discrete divergence operator $\mathop{\rm div}\nolimits_{h}\in\mathcal{L}(\textbf{H}_{h};Q_{h})\cap\mathcal{L}(\textbf{H};Q)$ by

[TABLE]

Then, we can introduce two bilinear forms

[TABLE]

and

[TABLE]

Then the finite element discretization of (27) is to find $(\textbf{u}_{h},p_{h})\in\textbf{H}_{h}\times Q_{h}$ such that

[TABLE]

This is the natural discretization of the modified weak formulation (27) except that the penalizing term $\textbf{J}(\textbf{u}_{h},\textbf{v}_{h})$ is added. This bilinear form $\textbf{J}(.,.)$ is defined by following the decomposition (44) of $\mathcal{E}_{h}$:

[TABLE]

where

[TABLE]

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 $\textbf{J}(\textbf{u},\textbf{v}).$

Remark 3.1**.**

The Eq. (43) have the matrix representation

[TABLE]

where U (resp. P) denote the coefficients of $\textbf{u}_{h}$ (resp. $p_{h}$) expanded with respect to a basis for $\textbf{H}_{h}$ (rep. $Q_{h}$).

We are now able to define the norm on $\textbf{H}_{h}$ (see [30]):

[TABLE]

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:

[TABLE]

From H$\ddot{o}$lder’s inequality, we derive the boundedness of $A_{h}(\cdot,\cdot)$ and $\textbf{b}_{h}(\cdot,\cdot)$:

Lemma 3.1**.**

(Continuity of forms) There holds:

[TABLE]

Theorem 3.1**.**

(Coercivity of $A_{h}$) There is an $\alpha>0$ such that:

[TABLE]

Proof.

Let $\textbf{v}_{h}\in\textbf{H}_{h}$. We have

[TABLE]

We introduce the local space

[TABLE]

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

[TABLE]

with the semi-norm

[TABLE]

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

[TABLE]

[TABLE]

$\bullet$ Estimate $\displaystyle\sum_{E\in\mathcal{E}_{h}(\Gamma_{s})}\int_{E}\left|\gamma_{\tau}(\textbf{v}_{h})\right|_{\mathbb{R}^{l}}$ ($l=1$ or $l=3$). We have by Cauchy-Schwarz inequality:

[TABLE]

Also, we have:

[TABLE]

Then,

[TABLE]

Hence we deduce

[TABLE]

$\bullet$ Now we estime the term $\displaystyle\sum_{E\in\mathcal{E}_{h}(\Gamma_{I})}\int_{E}\left|\gamma_{\tau}(\textbf{v}_{h})\right|_{\mathbb{R}^{l}}$. By Cauchy-Schwarz, we obtain:

[TABLE]

[TABLE]

Thus we deduce the estimation:

[TABLE]

Then,

[TABLE]

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

[TABLE]

Thus

[TABLE]

Hence,

[TABLE]

We have,

[TABLE]

The estimates (58), (59), (60) and (61), lead to (49). The proof is complete. ∎

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

[TABLE]

We define also the Crouzeix-Raviart interpolation operator $\textbf{r}_{h}:\textbf{W}\rightarrow\textbf{H}_{h}$ by:

[TABLE]

Lemma 3.2**.**

The operator $\textbf{r}_{h}$ is bounded: there is a constant $C_{5}>0$ depending on $\sigma$, $\mu$ and $N$ such that

[TABLE]

Proof.

The proof is similar to [30]. ∎

Then, we have the following result

Theorem 3.2**.**

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

[TABLE]

Proof.

We use Fortin argument i.e. for each $q_{h}\in Q_{h}$, we find $\textbf{v}_{h}\in\textbf{H}_{h}$ such that:

[TABLE]

Let $q_{h}\in Q_{h}\subset Q$. Then from [16, Corollary 2.4, Page 24], there exist vectoriel function $\textbf{v}\in[H_{0}^{1}(\Omega)]^{N}$ satisfying

[TABLE]

$[H_{0}^{1}(\Omega)]^{N}\subset W$, hence $\textbf{v}\in W$. We take $\textbf{v}_{h}=r_{h}\textbf{v}\in\textbf{H}_{h}$ and we have:

[TABLE]

Thus, we obtain

[TABLE]

Using the system $(\ref{IS1})$, we have:

[TABLE]

Also,

[TABLE]

From (70) et (71), we deduce:

[TABLE]

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

From Theorem 3.1 and Theorem 3.2 we have the following result:

Theorem 3.3**.**

There exists a unique solution $(\textbf{u}_{h},p_{h})\in\textbf{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 $\textbf{H}_{h}\nsubseteq\textbf{H}$, so the approximation error contains some extra consistency error terms. In fact, the abstract error estimates give the following result:

Lemma 3.3**.**

Let $(\textbf{u},p)\in\textbf{H}\times Q$ be the solution of problem (27) and $(\textbf{u}_{h},p_{h})\in\textbf{H}_{h}\times Q_{h}$ be the solution of the discrete problem (43). Then we have

[TABLE]

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

[TABLE]

Note that $\textbf{b}_{h}(\textbf{u},q_{h})=\textbf{b}(\textbf{u},q_{h})$, thus $E_{2h}=0$.

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

Assumption 3.1**.**

---

1. (1)

$\textbf{u}\in\textbf{H}$, $\textbf{u}_{s}\in[H^{2}(\Omega_{s})]^{N}$, $\textbf{u}_{d}\in[H^{2}(\Omega_{d})]^{N}$; 2. (2)

$p\in Q$, $p_{s}\in H^{1}(\Omega_{s})$, $p_{d}\in H^{1}(\Omega_{d})$.

We begin with the estimates for the terms: $\displaystyle\inf_{\textbf{v}_{h}\in\textbf{H}_{h}}\parallel\textbf{u}-\textbf{v}_{h}\parallel_{h}\mbox{ and }\displaystyle\inf_{q_{h}\in Q_{h}}\parallel p-q_{h}\parallel$

Lemma 3.4**.**

(Ref. [30]) There hold:

[TABLE]

Finally, let us consider the term $A_{h}(\textbf{u}_{h},\textbf{v}_{h})+\textbf{b}_{h}(\textbf{v}_{h},p_{h}-\tilde{L}(\textbf{v}_{h}).$ The smoothness assumption of u implies $\textbf{J}(\textbf{u},\textbf{v}_{h})=0$, thus $A_{h}(\textbf{u},\textbf{v}_{h})=\tilde{\textbf{a}}(\textbf{u},\textbf{v}_{h}),\forall\textbf{v}_{h}\in\textbf{H}_{h}$. Clearly,

[TABLE]

[TABLE]

Thus, we have

[TABLE]

where

[TABLE]

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}(T)$, denote by $P_{0}(E)$ the constant space of the restrictions to $E$ and $\pi_{E}$ the projection operator from $L^{(}E)$ on to $P_{0}(E)$ such that

[TABLE]

The operator $\pi_{E}$ has the property [7]:

[TABLE]

For any $\textbf{v}\in[L^{2}(E)]^{N}$, we let $\Pi_{E}\textbf{v}$ be the function in $[P_{0}(E)]^{N}$ such that

[TABLE]

Using inequality (80), we obtain

[TABLE]

Then we have the following lemma:

Lemma 3.5**.**

(Estimation the four face integrals) There holds:

[TABLE]

Proof.

1. (1)

Estimate (82): We begin with an estimate for the first term $R_{1}(\textbf{v}_{h})$. For any face $E\in\mathcal{E}_{h}(\Omega_{s}^{+})$, there exists at least one element $T\in\mathcal{T}_{h}^{s}$ such that $E\in\mathcal{E}(T)$. Then, from condition (79), H$\ddot{o}$der’s inequality and inequality (81), it follows that

[TABLE] 2. (2)

Estimate (83):

We have $\textbf{n}_{E}\cdot\textbf{D}(\textbf{u})_{|E}\in[L^{2}(E)]^{N}$, hence $\Pi_{E}(\textbf{n}_{E}\cdot\textbf{D}(\textbf{u}))\in P_{0}(E)^{N}$.

[TABLE]

Thus,

[TABLE]

Furthermore, summing on $E\in\mathcal{E}_{h}(\Omega_{s}^{+})$ faces, we obtain the estimate:

[TABLE] 3. (3)

For the terms $R_{3}(\textbf{v}_{h})$ and $R_{4}(\textbf{v}_{h})$, we use the same techniques as in the proof of the bounds for $R_{i}(\textbf{v}_{h})$, $i\in\{1,2\}$, and we obtain:

[TABLE]

The proof is complete. ∎

From Lemma 3.3, Lemma 3.4 and Lemma 3.5, now we derive the following convergence theorem:

Theorem 3.4**.**

Let the solution $(\textbf{u},p)$ of problem (27) satifies the smoothness assumption (Assumption 3.1). Let $(\textbf{u}_{h},p_{h})$ be the solution of the discrete problem (43). Then there exists a positive constant $C$ depending on $N$, $\mu$, $K_{*}$, $K^{*}$, $\alpha_{1}$ and $\sigma$ such that:

[TABLE]

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++ [11, 17] in isotropic coupled mesh of Fig. 13. The solutions have been represented by Mathematica software. For simplicity we choose each domain $\Omega_{l}$, $l\in\{s,d\}$ 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. $\overline{\Omega}=[0,1[\cup\{1\}\cup]1,2]$ like the show the Figure 10.

We consider the application $\phi:(x,y)\in\mathbb{R}^{2}\longmapsto\phi(x,y)=x^{2}(x-1)^{3}y^{2}(y-1)^{2}\in\mathbb{R}$ on the square $\Omega=]0,1[^{2}\cup]1,2[^{2}$. In $\Omega$, we define $\textbf{u}=(u_{1},u_{2})=\mathop{\rm curl}\nolimits\phi\left(-\frac{\partial\phi}{\partial y},\frac{\partial\phi}{\partial x}\right)$ and we obtain:

[TABLE]

We choose quadratic pressure $p\in L^{2}(\Omega)$ by

[TABLE]

Thus,

[TABLE]

The exact solution $(\textbf{u},p)$ satifies the following condition:

[TABLE]

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

[TABLE]

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

[TABLE]

Thus, $\textbf{f}_{s}(x,y)=(f_{1}(x,y),f_{2}(x,y))$ in $\Omega_{s}$ leads to

[TABLE]

and in $\Omega_{d}$, $\textbf{f}_{d}(x,y)=(k_{1}(x,y),k_{2}(x,y))$ is given by:

[TABLE]

5. Conclusion

In this contribution, we investigated a new mixed finite element method to solve the Stokes-Darcy fluid flow model without introducing any Lagrange multiplier. We proposed a modification of the Darcy problem which allows us to apply a slight variant nonconforming Crouzeix-Raviart element to the whole coupled Stokes-Darcy problem. The proposed method is probably one 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.

6. Acknowledgment

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.