#  A posteriori error estimates for Darcy-Forchheimer’s problem

Toni Sayah† ,  Georges Semaan† and  Faouzi Triki‡

###### Abstract.

This work deals with the _a posteriori_ error estimates for the Darcy-
Forchheimer problem. We introduce the corresponding variational formulation
and discretize it by using the finite-element method. A posteriori error
estimate with two types of computable error indicators is showed. The first
one is linked to the linearization and the second one to the discretization.
Finally, numerical computations are performed to show the effectiveness of the
error indicators.

August 28, 2024.
† Laboratoire de "Mathématiques et Applications", Unité de Recherche
"Mathématoqies et Modélisation", CAR, Faculté des sciences de Beyrouth,
Université Saint-Joseph de Beyrouth, Liban.
Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>
‡ Laboratoire Jean Kuntzmann, UMR CNRS 5224, Université Grenoble-Alpes, 700
Avenue Centrale, 38401 Saint-Martin-d’Hères, France. E-mails:
<EMAIL_ADDRESS>FT was supported by the grant
ANR-17-CE40-0029 of the French National Research Agency ANR (project
MultiOnde)

##  1\. Introduction

Let $\Omega$ be a bounded open domain of ${\rm IR}^{d}$ ($d=2,3$) with
Lipschitz-continuous boundary $\Gamma=\partial\Omega$. We consider the Darcy-
Forchheimer equation

$\displaystyle\frac{\mu}{\rho}K^{-1}{\bf u}+\frac{\beta}{\rho}|{\bf u}|{\bf
u}+\nabla p={\bf f}\;\;\;\;\;\mbox{in}\;\;\;\Omega,$ (1.1)

with the divergence constraint

$\displaystyle\operatorname{div}{\bf
u}=b\;\;\;\;\;\mbox{in}\;\;\;\Omega{\color[rgb]{0,0,1},}$ (1.2)

and the boundary condition

$\displaystyle{\bf u}\cdot{\bf n}=g\;\;\;\;\;\mbox{on}\;\;\partial\Omega,$
(1.3)

where ${\bf u}$ represents the velocity, $p$ represents the pressure, ${\bf
n}$ is the unit exterior normal vector to $\Gamma$, $|.|$ denotes the
Euclidean norm and $|{\bf u}|^{2}={\bf u}\cdot{\bf u}$. The parameters $\rho$,
$\mu$ and $\beta$ represent the density of the fluid, its viscosity and its
dynamic viscosity, respectively. $\beta$ is also referred as Forchheimer
number which is a scalar positive constant. $K$ is the permeability tensor,
assumed to be uniformly positive definite and bounded such that there exists
two positive real numbers $K_{m}$ and $K_{M}$ verifying

$0<K_{m}\,{\bf{x}}\cdot{\bf{x}}\leq(K^{-1}({\bf{x}}){\bf{x}})\cdot{\bf{x}}\leq
K_{M}\,{\bf{x}}\cdot{\bf{x}}.$ (1.4)

For the compatibility, we suppose that $b$ and $g$ verify the following
compatibility condition:

$\displaystyle\int_{\Omega}b({\bf{x}})d{\bf{x}}=\displaystyle\int_{\Gamma}\mbox{g}(s)\,ds.$

We denote by Problem $(P)$ the system of equations ((1.1), (1.2), (1.3)).

Darcy’s law describes the creeping flow of Newtonian fluids in porous media.
It is represented by Equation (1.1) without the non-linear term
$\displaystyle\frac{\beta}{\rho}|{\bf u}|{\bf u}$ and is valid by experiment
under the condition that the creeping velocity is low and the porosity and
permeability are small enough by Darcy in $1856$ [3]. Forchheimer [25]
conducted flow experiments in sandpacks and recognized that when the velocity
is higher and the porosity is nonuniform, Darcy’s equation can be replaced by
Equation (1.1). A theoretical derivation of Forchheimer’s law can be found in
[32].
Multiple works approximated the Darcy-Forchheimer equation by using finite
element methods. Girault and Wheeler [17] approximated the velocity by
piecewise constants and the pressure by Crouzeix-Raviart element. They also
proposed and studied an alternating directions iterative method to solve the
system of nonlinear equations obtained by finite element discretizaton. Lopez
and al. [28] carried out numerical tests of the methods studied in [17] in
order to corroborate the results presented there. In [31], the authors
proposed and studied a mixed element approximation: the Raviart-Thomas mixed
element, Brezzi-Douglas-Marini mixed element. Salas J. et al. [33] presented a
theoretical study of the mixed finite element space, as proposed in [28]. In
[34], we considered the Darcy-Forchheimer problem coupled with the convection-
diffusion-reaction problem. We established existence of solutions by using a
Galerkin method and we proved uniqueness. Then, we introduced and analyzed a
numerical scheme based on the finite element method and we derived an optimal
a priori error estimates for the proposed numerical scheme. In this work, we
introduce a numerical iterative scheme, show the corresponding convergence,
establish the corresponding a posteriori error estimates and show
corresponding numerical investigations.

The present work investigates _a posteriori_ error estimates for the finite
element discretization of Darcy-Forchheimer problem.These estimates can be
used to evaluate the solution errors of the discrete problem without requiring
any a priori information on the exact solution. Indeed, the _a posteriori_
analysis controls the overall discretization error of a problem by providing
error indicators that are easy to compute. Once these error indicators are
constructed, their efficiency can be proven by bounding each indicator by the
local error. _A posteriori_ analysis was first introduced by I.
Babu$\check{s}$ka [5], developed by R. Verfürth [37], and has been the object
of a large number of publications. A posteriori error estimations have been
studied for several types of partial differential equations such that the
Stokes or Navier-Stokes equation (see for instance [37, 10, 4, 1, 12, 13,
19]), the Maxwell and Lamé equations [30, 8]. Many works have been established
for the Darcy flow, see for instance [2, 11, 14, 29]. In [15], Chen and Wang
establish optimal a poteriori error estimates for the $H({\rm div},\Omega)$
conforming mixed finite element method applied to the coupled Darcy-Stokes
system in two dimensions. For the Darcy equations with pressure dependent
viscosity, we refer to [26] and the references therein. In [20], we establish
a posteriori error estimates for Darcy’s problem coupled with the heat
equation. In this paper, we study the posteriori error estimates corresponding
to Problem $(P)$ and show corresponding numerical results.

This paper is organised as follow:

  * •

Section 2 describes the problem and the weak formulation.

  * •

Section 3 describes the discretization and studies the iterative scheme.

  * •

Section 4 is devoted to the a posteriori error estimates.

  * •

Section 5 shows numerical investigations.

##  2\. Variational formulation

In this section, we begin by introducing several definitions and notations in
order to write the weak formulation corresponding to Problem $(P)$.
Let $\alpha=(\alpha_{1},\alpha_{2},\dots\alpha_{d})$ be a $d$-uplet of non
negative integers, set $|\alpha|=\displaystyle\sum_{i=1}^{d}\alpha_{i}$, and
define the partial derivative $\partial^{\alpha}$ by

$\partial^{\alpha}=\displaystyle\frac{\partial^{|\alpha|}}{\partial
x_{1}^{\alpha_{1}}\partial x_{2}^{\alpha_{2}}\dots\partial
x_{d}^{\alpha_{d}}}.$

Then, for any positive integer $m$ and any number $q\geq 1$, we recall the
classical Sobolev space

$W^{m,q}(\Omega)=\\{v\in L^{q}(\Omega);\,\forall\,|\alpha|\leq
m,\;\partial^{\alpha}v\in L^{q}(\Omega)\\},$ (2.1)

equipped with the seminorm

$|v|_{W^{m,q}(\Omega)}=\big{\\{}\sum_{|\alpha|=m}\int_{\Omega}|\partial^{\alpha}v|^{q}\,d{\bf{x}}\,\big{\\}}^{\frac{1}{q}}$
(2.2)

and the norm

$\|v\|_{W^{m,q}(\Omega)}=\big{\\{}\sum_{0\leq k\leq
m}|v|_{W^{k,q}(\Omega)}^{q}\big{\\}}^{\frac{1}{q}}.$ (2.3)

When $q=2$, this space is the Hilbert space $H^{m}(\Omega)$. In particular,
the scalar product of $L^{2}(\Omega)$ is denoted by $(.,.)$. Furthermore, we
recall the following standard space for Darcy’s equations:

$L^{2}_{0}(\Omega)=\\{v\in
L^{2}(\Omega);\displaystyle\int_{\Omega}v\,d{\bf{x}}=0\\}.$

The definitions of these spaces are extended straightforwardly to vectors,
with the same notation, but with the following modification for the norms in
the non-Hilbert case. Let ${\bf v}$ be a vector valued function; we set

$\|{\bf v}\|_{L^{q}(\Omega)}=\big{(}\int_{\Omega}|{\bf
v}|^{q}\,d{\bf{x}}\,\big{)}^{\frac{1}{q}},$ (2.4)

where $|.|$ denotes the Euclidean vector norm.

We introduce the spaces:

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus 2.0ptX=L^{3}(\Omega)^{d},\\\
M=W^{1,\frac{3}{2}}(\Omega)\cap L^{2}_{0}(\Omega),\\\ {H=\left\\{{\bf v}\in
X;\operatorname{div}{\bf v}\in L^{3d/(d+3)}(\Omega)\right\\},}\\\ {V=\\{{\bf
v}\in X;\forall q\in M,\displaystyle\int_{\Omega}\nabla q\cdot{\bf
v}\,d{\bf{x}}=0\\}.}\end{array}$

The space $H$ endowed with the graph norm is complete. Moreover
$D(\overline{\Omega})$ is dense in $H$, and ${\bf v}\cdot{\bf n}$ belongs to
$W^{1,3/2}(\Gamma)$ for all $v\in H$ [35]. We have the following inf-sup
condition (see [17] for the proof)

$\underset{q\in M}{\inf}\;\;\underset{{\bf v}\in
X}{\sup}\;\;\displaystyle\frac{\displaystyle\int_{\Omega}\nabla q\cdot{\bf
v}\,d{\bf{x}}}{||\nabla q||_{L^{3/2}(\Omega)}||{\bf v}||_{L^{3}(\Omega)}}=1.$
(2.5)

The velocity ${\bf u}$ and pressure $p$ of Problem $(P)$ are considered
respectively in $H$ and $M$, while $b$ and $g$ are assumed to be respectively
in $L^{3d/(d+3)}(\Omega)$ and $L^{3(d-1)/d}(\Gamma)$ (see [17] for details).

In order to write the variational formulation associated to Problem $(P)$, we
introduce the mapping ${\bf v}\longrightarrow\mathcal{A}({\bf v})$ defined by:

$\begin{array}[]{rccl}\mathcal{A}:&L^{3}(\Omega)^{d}&\mapsto&L^{\frac{3}{2}}(\Omega)^{d}\\\
&{\bf v}&\mapsto&\mathcal{A}({\bf v})=\displaystyle\frac{\mu}{\rho}K^{-1}{\bf
v}+\frac{\beta}{\rho}|{\bf v}|{\bf v}.\end{array}$

We refer to [17, 24] for the following properties of $\mathcal{A}$.

######  Property 2.1.

We have the following monotonicity properties:

  1. (1)

for all ${\bf v}\in L^{3}(\Omega)^{d}$,

$\displaystyle\frac{\mu}{\rho}\int_{\Omega}K^{-1}{\bf v}\cdot{\bf
v}\,d{\bf{x}}\geq\displaystyle\frac{\mu K_{m}}{\rho}\left\|{\bf
v}\right\|_{L^{2}(\Omega)}^{2}$ (2.6)

  2. (2)

for all ${\bf v},{\bf w}\in L^{3}(\Omega)^{d}$,

$\displaystyle\frac{\beta}{\rho}\int_{\Omega}(|{\bf v}|{\bf v}-|{\bf w}|{\bf
w})({\bf v}-{\bf w})\,d{\bf{x}}\geq c_{m}\left\|{\bf v}-{\bf
w}\right\|_{L^{3}(\Omega)}^{3}.$ (2.7)

where $c_{m}$ is a strictly positive constant.

Recall the following Green’s formula [35]

$\forall q\in M,\;\forall v\in H,\quad\displaystyle\int_{\Omega}\nabla
q\cdot{\bf
v}\,d{\bf{x}}=\displaystyle-\int_{\Omega}q\operatorname{div}vd{\bf{x}}+\int_{\Gamma}q{\bf
v}\cdot{\bf n}ds.$ (2.8)

Using the Green formula (2.8) one can show that Problem $(P)$ is equivalent to
the following variational formulation: Find $({\bf u},p)\in X\times M$ such
that

$\left\\{\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus 2.0pt\forall{\bf
v}\in X,\quad\displaystyle\int_{\Omega}\mathcal{A}({\bf u})\cdot{\bf
v}\,d{\bf{x}}+\int_{\Omega}\nabla p\cdot{\bf
v}\,d{\bf{x}}=\displaystyle\int_{\Omega}{\bf f}\cdot{\bf v}\,d{\bf{x}},\\\
\forall q\in M,\quad\displaystyle\int_{\Omega}\nabla q\cdot{\bf
u}\,d{\bf{x}}=-\int_{\Omega}b\,q\,d{\bf{x}}+\int_{\Gamma}g\,q\,ds.\end{array}\right.$
(2.9)

It can be demonstrated that, for all $b\in L^{3d/(d+3)}(\Omega)$ and $g\in
L^{3(d-1)/d}(\Gamma)$, there is only one ${\bf u}_{l}\in{L^{3}(\Omega)/V}$
such that

$\forall q\in M,\,\displaystyle\int_{\Omega}\nabla q\cdot{\bf
u}_{l}\,d{\bf{x}}=-\int_{\Omega}b\,q\,d{\bf{x}}+\int_{\Gamma}g\,q\,ds,$ (2.10)

and

$\displaystyle||{\bf u}_{l}||_{L^{3}(\Omega)^{d}/V}\leq
C\big{(}||b||_{L^{3d/(d+3)}(\Omega)}+||g||_{L^{3(d-1)/d}(\Gamma)}\big{)},$

where $C>0$ is a constant depending only on $\Omega$ and $d$.

Using The Banach-Necas-Babu$\check{s}$ka Theorem (see for instance [6],
Theorem A.1 in [36]) one can show that (2.9), and equivalently the problem
$(P)$ admits a unique solution $({\bf u},p)\in X\times M$ satisfying the
relations

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus 2.0pt||{\bf
u}||_{L^{3}(\Omega)}\leq C\Big{(}||{\bf u}_{l}||_{L^{3}(\Omega)}^{3}+||{\bf
u}_{l}||_{L^{2}(\Omega)}^{2}+||{\bf f}||_{L^{2}(\Omega)}^{2}\Big{)}^{1/3},\\\
||\nabla p||_{L^{3/2}(\Omega)}\leq C\Big{(}||{\bf
u}||_{L^{3/2}(\Omega)}+||{\bf u}||_{L^{3}(\Omega)}^{2}+||{\bf
f}||_{L^{2}(\Omega)}\Big{)}.\end{array}$ (2.11)

For further details of the weak formulation corresponding to Problem $(P)$
presented in this section, We refer to [17].

##  3\. Discretization and a posteriori estimation

In this section, we introduce an iterative numerical scheme to solve the weak
variational formulation (2.9) of the Problem $(P)$, and study its convergence.
We establish the corresponding a posteriori error estimate in a second step.

From now on, we assume that $\Omega$ is a polygon when $d=2$ or polyhedron
when $d=3$, so it can be completely meshed. For the space discretization, we
consider a regular (see Ciarlet [16]) family of triangulations
$(\mathcal{T}_{h})_{h}$ of $\Omega$ which is a set of closed non degenerate
triangles for $d=2$ or tetrahedra for $d=3$, called elements, satisfying,

  * •

for each $h$, $\overline{\Omega}$ is the union of all elements of
$\mathcal{T}_{h}$;

  * •

the intersection of two distinct elements of $\mathcal{T}_{h}$ is either
empty, a common vertex, or an entire common edge (or face when $d=3$);

  * •

the ratio of the diameter $h_{\kappa}$ of an element $\kappa$ in
$\mathcal{T}_{h}$ to the diameter $\rho_{\kappa}$ of its inscribed circle when
$d=2$ or ball when $d=3$ is bounded by a constant independent of $h$: there
exists a positive constant $\sigma$ independent of $h$ such that,

$\displaystyle\max_{\kappa\in\mathcal{T}_{h}}\frac{h_{\kappa}}{\rho_{\kappa}}\leq\sigma.$
(3.1)

As usual, $h$ denotes the maximal diameter of all elements of
$\mathcal{T}_{h}$. To define the finite element functions, let $r$ be a non
negative integer. For each $\kappa$ in $\mathcal{T}_{h}$, we denote by ${\rm
IP}_{r}(\kappa)$ the space of restrictions to $\kappa$ of polynomials in $d$
variables and total degree at most $r$, with a similar notation on the faces
or edges of $\kappa$. For every edge (when $d=2$) or face (when $d=3$) $e$ of
the mesh $\mathcal{T}_{h}$, we denote by $h_{e}$ the diameter of $e$.

We shall use the following inverse inequality: for any numbers $p,q\geq 2$,
for any dimension $d$, and for any non negative integer $r$, there exist
constants $C_{I}(p)$ and $C_{J}(q)$ such that for any polynomial function
$v_{h}$ of degree $r$ on an element $\kappa$ or an edge (when $d=2$) or face
(when $d=3$) $e$ of the mesh $\mathcal{T}_{h}$,

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\|v_{h}\|_{L^{p}(\kappa)}\leq
C_{I}(p)h_{\kappa}^{\frac{d}{p}-\frac{d}{2}}\|v_{h}\|_{L^{2}(\kappa)}\\\
\mbox{and}\\\ \|v_{h}\|_{L^{q}(e)}\leq
C_{J}(q)h_{e}^{\frac{d-1}{q}-\frac{d-1}{2}}\|v_{h}\|_{L^{2}(e)}.\end{array}$
(3.2)

$C_{I}$ and $C_{J}$ depend on the regularity parameter $\sigma$ of (3.1). In
all the rest of the paper and for the simplicity, we will omit the dependency
of $C_{I}$ and $C_{J}$ with $p$ and $q$.

Let $X_{h}\subset X$ and $M_{h}\subset M$ be the discrete spaces corresponding
to the velocity and the pressure given by:

$\begin{split}Z_{h}=&\\{q_{h}\in
C^{0}(\bar{\Omega});\,\forall\,\kappa\in\mathcal{T}_{h},\;q_{h}|_{\kappa}\in{\rm
IP}_{1}\\},\\\ X_{h}=&\\{{\bf v}_{h}\in
L^{2}(\bar{\Omega})^{d};\,\forall\,\kappa\in\mathcal{T}_{h},\;{\bf
v}_{h}|_{\kappa}\in{\rm IP}_{0}^{d}\\},\\\ M_{h}=&Z_{h}\cap
L^{2}_{0}(\Omega).\end{split}$ (3.3)

They satisfy the following inf-sup condition (see [33]):

$\forall\,q_{h}\in M_{h},\;\sup_{{\bf v}_{h}\in
X_{h}}\displaystyle\frac{\displaystyle\int_{\Omega}\nabla q_{h}\cdot{\bf
v}_{h}\,d{\bf{x}}\,}{\|{\bf v}_{h}\|_{X}}\geq\beta_{1}\|q_{h}\|_{M_{h}},$
(3.4)

where $\beta_{1}$ is a positive constant independent of $h$.

Problem (2.9) can be discretized as following:

$\left\\{\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus 2.0pt\forall{\bf
v}_{h}\in X_{h},\quad\displaystyle\int_{\Omega}\mathcal{A}({\bf
u}_{h})\cdot{\bf v}_{h}\,d{\bf{x}}+\int_{\Omega}\nabla p_{h}\cdot{\bf
v}_{h}\,d{\bf{x}}=\displaystyle\int_{\Omega}{\bf f}\cdot{\bf
v}_{h}\,d{\bf{x}},\\\ \forall q_{h}\in
M_{h},\quad\displaystyle\int_{\Omega}\nabla q_{h}\cdot{\bf
u}_{h}\,d{\bf{x}}=-\int_{\Omega}bq_{h}\,d{\bf{x}}+\int_{\Gamma}gq_{h}\,ds.\end{array}\right.$
(3.5)

It is shown in [33] that there exists a unique ${\bf u}_{h,l}\in X_{h}$ such
that

$\forall q_{h}\in M_{h},\quad\displaystyle\int_{\Omega}\nabla q_{h}\cdot{\bf
u}_{h,l}\,d{\bf{x}}=\displaystyle-\int_{\Omega}b\,q_{h}d{\bf{x}}+\int_{\Gamma}g\,q_{h}ds,$
(3.6)

and ${\bf u}_{h,l}$ satisfies the following bound,

$||{\bf u}_{h,l}||_{L^{3}(\Omega)^{d}}\leq
C_{l}\big{(}||b||_{L^{\frac{3d}{3+d}}(\Omega)}+||g||_{L^{\frac{3(d-1)}{d}}(\Gamma)}\big{)}.$
(3.7)

It is also shown in [33] that Problem (3.5) admits a unique solution $({\bf
u}_{h},p_{h})\in X_{h}\times M_{h}$ satisfying exactly similar bounds as
(2.11). We refer to [33] for the proof of the following a priori error
estimates:

######  Theorem 3.1.

The solutions $({\bf u},p)$ of (2.9) and $({\bf u}_{h},p_{h})$ of
$\eqref{V1h}$ verify the following a priori error:
If $({\bf u},p)\in W^{1,4}(\Omega)^{d}\times W^{2,3/2}(\Omega)$, then there
exists strictly positive constants $C_{u},C_{p}$ independent of $h$ such that

$\begin{array}[]{ll}||{\bf u}-{\bf u}_{h}||_{L^{2}(\Omega)}\leq C_{u}h,\\\
||\nabla(p-p_{h})||_{3/2,\Omega}\leq C_{p}h.\end{array}$ (3.8)

###  3.1. Iterative algorithm:

To compute the solution of the non-linear problem (3.5), we introduce the
following iterative algorithm: for a given initial guess ${\bf u}_{h}^{0}\in
X_{h}$ and having ${\bf u}_{h}^{i}$ at each iteration $i$, we compute $({\bf
u}_{h}^{i+1},p^{i+1}_{h})$ solution of

$\left\\{\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus 2.0pt\forall{\bf
v}_{h}\in X_{h},\quad\displaystyle\int_{\Omega}\alpha({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})\cdot{\bf v}_{h}\,d{\bf{x}}+\frac{\mu}{\rho}\int_{\Omega}K^{-1}{\bf
u}^{i+1}_{h}\cdot{\bf v}_{h}\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}|{\bf
u}^{i}_{h}|{\bf u}^{i+1}_{h}\cdot{\bf v}_{h}\,d{\bf{x}}\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt\hskip 142.26378pt+\displaystyle\int_{\Omega}\nabla
p^{i+1}_{h}\cdot{\bf v}_{h}\,d{\bf{x}}=\displaystyle\int_{\Omega}{\bf
f}\cdot{\bf v}_{h}\,d{\bf{x}},\\\ \forall q_{h}\in
M_{h},\quad\displaystyle\int_{\Omega}\nabla q_{h}\cdot{\bf
u}^{i+1}_{h}\,d{\bf{x}}=-\int_{\Omega}bq_{h}\,d{\bf{x}}+\int_{\Gamma}gq_{h}\,ds,\end{array}\right.$
(3.9)

We next address the convergence of Scheme (3.9). The analysis of convergence
has two principal steps. The first is to derive an uniform bound to the
numerical iterative solution with respect to the iteration index. Then the
obtained bound is used to study the convergence of a linear explicit scheme.

######  Theorem 3.2.

Problem $(\ref{V1hi})$ admits a unique solution $({\bf
u}_{h}^{i+1},p_{h}^{i+1})\in X_{h}\times M_{h}$. Furthermore, if the initial
value ${\bf u}_{h}^{0}$ satisfies the condition

${||{\bf u}_{h}^{0}||_{L^{2}(\Omega)}\leq\displaystyle L_{1}({\bf f},{\bf
u}_{h,l},\alpha),}$ (3.10)

where

$\displaystyle L_{1}({\bf f},{\bf
u}_{h,l},\alpha)=\displaystyle\left(\ell_{0}+\ell_{1}\alpha\right)^{\frac{1}{2}},$
(3.11)

with

$\displaystyle\ell_{0}={\left(\frac{2\rho}{\mu
K_{m}}\right)^{2}}\Big{(}\displaystyle(\frac{3\rho}{2\mu
K_{m}}+\frac{1}{2})||{\bf f}||^{2}_{L^{2}(\Omega)}+(\frac{1}{2}+\frac{3\mu
K_{M}^{2}}{2\rho K_{m}})||{\bf
u}_{h,l}||^{2}_{L^{2}(\Omega)}+\frac{4\beta}{3\rho}||{\bf
u}_{h,l}||^{3}_{L^{3}(\Omega)}\Big{)},$ (3.12)
$\displaystyle\ell_{1}={\left(\frac{2\rho}{\mu K_{m}}\right)^{2}}||{\bf
u}_{h,l}||^{2}_{L^{2}(\Omega)},$ (3.13)

and if $\alpha\geq\alpha^{\star}$ with
$\alpha^{\star}=4(\gamma_{1}+\sqrt{\gamma_{1}^{2}+\gamma_{2}})^{2}$ where

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus
2.0pt\gamma_{1}&=&\displaystyle\frac{8\beta K_{M}}{3\rho
K_{m}}C_{I}^{3}h^{-\frac{d}{2}}\sqrt{\ell_{1}},\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\gamma_{2}&=&\displaystyle\frac{3\beta^{2}}{2\rho\mu
K_{m}}{C_{I}^{4}}h^{-\frac{2d}{3}}||{\bf
u}_{h,l}||^{2}_{L^{3}(\Omega)}+\frac{8\beta}{3\mu
K_{m}}{C_{I}^{3}}h^{-\frac{d}{2}}\|{\bf f}\|_{L^{2}(\Omega)^{2}}\\\
&&\displaystyle+\frac{8\beta K_{M}}{3\rho
K_{m}}{C_{I}^{3}}h^{-\frac{d}{2}}\sqrt{\ell_{0}}+\frac{8\beta^{2}}{3\rho^{2}}{C_{I}^{6}}h^{-d}\max(\frac{2\rho\ell_{0}}{\mu
K_{m}},\ell_{1}),\end{array}$

then the solution of Problem $(\ref{V1hi})$ satisfies the estimates

${||{\bf u}_{h}^{i+1}||_{L^{2}(\Omega)}\leq\displaystyle L_{1}({\bf f},{\bf
u}_{h,l},\alpha),}$ (3.14)

and

${||{\bf u}_{h}^{i+1}||^{3}_{L^{3}(\Omega)}\leq\displaystyle\frac{3\mu
K_{m}}{2\beta}L_{1}^{2}({\bf f},{\bf u}_{h,l},\alpha).}$ (3.15)

Proof. To prove the existence and uniqueness of the solution of Problem
$(\ref{V1hi})$ which is a square linear system in finite dimension, it
suffices to show the uniqueness. For a given ${\bf u}_{h}^{i}$, let $({\bf
u}_{h1}^{i+1},p_{h1}^{i+1})$ and $({\bf u}_{h2}^{i+1},p_{h2}^{i+1})$ two
different solutions of Problem $(\ref{V1hi})$ and let ${\bf w}_{h}={\bf
u}_{h1}^{i+1}-{\bf u}_{h2}^{i+1}$ and $\xi_{h}=p_{h1}^{i+1}-p_{h2}^{i+1}$,
then $({\bf w}_{h},\xi_{h})$ is the solution of the following problem:

$\left\\{\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus 2.0pt\forall{\bf
v}_{h}\in X_{h},\quad\displaystyle\int_{\Omega}\alpha{\bf w}_{h}\cdot{\bf
v}_{h}\,d{\bf{x}}+\frac{\mu}{\rho}\int_{\Omega}K^{-1}{\bf w}_{h}\cdot{\bf
v}_{h}\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}|{\bf u}^{i}_{h}|{\bf
w}_{h}\cdot{\bf
v}_{h}\,d{\bf{x}}+\displaystyle\int_{\Omega}\nabla\xi_{h}\cdot{\bf
v}_{h}\,d{\bf{x}}=0,\\\ \forall q_{h}\in
M_{h},\quad\displaystyle\int_{\Omega}\nabla q_{h}\cdot{\bf
w}_{h}\,d{\bf{x}}=0.\end{array}\right.$

By taking $({\bf v}_{h},q_{h})=({\bf w}_{h},\xi_{h})$ and by remarking that
$\displaystyle\frac{\beta}{\rho}\int_{\Omega}|{\bf u}^{i}_{h}||{\bf
w}_{h}|^{2}\,d{\bf{x}}\geq 0$, we obtain by using the properties of $K^{-1}$
the following bound:

$(\alpha+\displaystyle\frac{K_{m}\mu}{\rho})||{\bf
w}_{h}||^{2}_{L^{2}(\Omega)}\leq 0.$

Thus, we deduce that ${\bf w}_{h}=0$. The inf-sup condition (3.4) deduces that
$\xi_{h}=0$ and then, we get the uniqueness of the solution of Problem
$(\ref{V1hi})$.
To prove the bound (3.14), we need first to bound the error $\|{\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i}\|_{L^{2}(\Omega)}$ with respect to the previous
value ${\bf u}_{h}^{i}$. The second equation of Problem $(\ref{V1hi})$ allows
us to deduce the relation

$\forall q_{h}\in M_{h},\quad\displaystyle\int_{\Omega}\nabla q_{h}\cdot({\bf
u}^{i+1}_{h}-{\bf u}_{h}^{i})\,d{\bf{x}}=0.$ (3.16)

Then, the first equation of (3.9) with ${\bf v}_{h}={\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}$ gives

$\displaystyle\alpha\|{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}\|^{2}_{L^{2}(\Omega)}+\frac{\mu}{\rho}\int_{\Omega}K^{-1}{\bf
u}_{h}^{i+1}\cdot({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}|{\bf u}_{h}^{i}|{\bf
u}_{h}^{i+1}\cdot({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})\,d{\bf{x}}=\displaystyle\int_{\Omega}{\bf f}\cdot({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})\,d{\bf{x}}.$

By inserting $\pm{\bf u}_{h}^{i}$ in the second and the third terms of the
last equation we get

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\alpha\|{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}\|^{2}_{L^{2}(\Omega)}+\frac{\mu}{\rho}\int_{\Omega}K^{-1}|{\bf
u}_{h}^{i+1}-{\bf
u}_{h}^{i}|^{2}\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}|{\bf
u}_{h}^{i}||{\bf u}_{h}^{i+1}-{\bf u}_{h}^{i}|^{2}\,d{\bf{x}}\\\ \hskip
56.9055pt=\displaystyle\int_{\Omega}{\bf f}\cdot({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})\,d{\bf{x}}-\frac{\mu}{\rho}\int_{\Omega}K^{-1}{\bf
u}_{h}^{i}\cdot({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})\,d{\bf{x}}-\frac{\beta}{\rho}\int_{\Omega}|{\bf u}_{h}^{i}|{\bf
u}_{h}^{i}\cdot({\bf u}_{h}^{i+1}-{\bf u}_{h}^{i})\,d{\bf{x}}.\end{array}$

By using the properties of $K$, the Cauchy-Schwartz inequality and Relation
(3.2), we get, by remarking that the third term of the last equation is non-
negative, the bound

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\alpha\|{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}\|^{2}_{L^{2}(\Omega)}+\frac{\mu K_{m}}{\rho}\|{\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i}\|^{2}_{L^{2}(\Omega)}\leq\|{\bf
f}\|_{L^{2}(\Omega)}\|{\bf u}_{h}^{i+1}-{\bf u}_{h}^{i}\|_{L^{2}(\Omega)}\\\
\hskip 113.81102pt+\displaystyle\frac{\mu K_{M}}{\rho}\|{\bf
u}_{h}^{i}\|_{L^{2}(\Omega)}\|{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}\|_{L^{2}(\Omega)}+\frac{\beta}{\rho}{C_{I}^{3}}h^{-\frac{d}{2}}\|{\bf
u}_{h}^{i}\|^{2}_{L^{2}(\Omega)}\|{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}\|_{L^{2}(\Omega)}.\end{array}$

We simplify by $\|{\bf u}_{h}^{i+1}-{\bf u}_{h}^{i}\|_{L^{2}(\Omega)}$ to
obtain

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle(\alpha+\frac{\mu K_{m}}{2\rho})\|{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}\|_{L^{2}(\Omega)}\leq\displaystyle\|{\bf
f}\|_{L^{2}(\Omega)}+\frac{\mu K_{M}}{\rho}\|{\bf
u}_{h}^{i}\|_{L^{2}(\Omega)}+\frac{\beta}{\rho}{C_{I}^{3}}h^{-\frac{d}{2}}\|{\bf
u}_{h}^{i}\|^{2}_{L^{2}(\Omega)},\end{array}$

and then we get the following bound

$\displaystyle\|{\bf u}_{h}^{i+1}-{\bf u}_{h}^{i}\|_{L^{2}(\Omega)}\leq
L_{2}({\bf f},\|{\bf u}_{h}^{i}\|_{L^{2}(\Omega)}),$ (3.17)

where

$L_{2}({\bf f},\eta)=\displaystyle\frac{1}{\alpha+\frac{\mu
K_{m}}{2\rho}}\Big{(}\|{\bf f}\|_{L^{2}(\Omega)}+\frac{\mu
K_{M}}{\rho}\eta+\frac{\beta}{\rho}{C_{I}^{3}}h^{-\frac{d}{2}}\eta^{2}\Big{)},\;\;\eta\in\mathbb{R}_{+}.$

Then, we are in position to show the relation (3.14). Property (3.6) allows us
to deduce that the term ${\bf u}^{i+1}_{h,0}={\bf u}_{h}^{i+1}-{\bf u}_{h,l}$
is in $X_{h}$, and verifies

$\forall q_{h}\in M_{h},\quad\displaystyle\int_{\Omega}\nabla q_{h}\cdot{\bf
u}^{i+1}_{h,0}\,d{\bf{x}}=0.$ (3.18)

We consider the first equation of (3.9) with ${\bf v}_{h}={\bf
u}^{i+1}_{h,0}={\bf u}_{h}^{i+1}-{\bf u}_{h,l}$ and we obtain:

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\alpha\int_{\Omega}({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})\cdot{\bf
u}_{h}^{i+1}\,d{\bf{x}}+\frac{\mu}{\rho}\int_{\Omega}K^{-1}|{\bf
u}_{h}^{i+1}|^{2}\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}|{\bf
u}_{h}^{i+1}|^{3}\,d{\bf{x}}=\displaystyle\int_{\Omega}{\bf f}\cdot({\bf
u}_{h}^{i+1}-{\bf u}_{hl})\,d{\bf{x}}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\hskip 28.45274pt\displaystyle+\alpha\int_{\Omega}({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})\cdot{\bf
u}_{h,l}\,d{\bf{x}}+\frac{\mu}{\rho}\int_{\Omega}K^{-1}{\bf
u}_{h}^{i+1}\cdot{\bf
u}_{h,l}\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}(|{\bf u}_{h}^{i+1}|-|{\bf
u}^{i}_{h}|)|{\bf u}_{h}^{i+1}|^{2}\,d{\bf{x}}\\\ \hskip
56.9055pt\displaystyle+\frac{\beta}{\rho}\int_{\Omega}(|{\bf u}^{i}_{h}|-|{\bf
u}^{i+1}_{h}|){\bf u}_{h}^{i+1}\cdot{\bf
u}_{h,l}\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}|{\bf u}^{i+1}_{h}|\,{\bf
u}_{h}^{i+1}\cdot{\bf u}_{h,l}\,d{\bf{x}}.\end{array}$

By using the properties of $K$, the Cauchy-Schwartz inequality and the
relation
$a^{2}b\leq\frac{1}{3}\big{(}\frac{1}{\delta^{3}}b^{3}+2\delta^{3/2}a^{3}\big{)}$
(for any positive real numbers $a$ and $b$), we get:

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{\alpha}{2}||{\bf
u}_{h}^{i+1}||^{2}_{L^{2}(\Omega)}-\frac{\alpha}{2}||{\bf
u}_{h}^{i}||^{2}_{L^{2}(\Omega)}+\frac{\alpha}{2}||{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}||^{2}_{L^{2}(\Omega)}+\frac{\mu}{\rho}K_{m}||{\bf
u}_{h}^{i+1}||^{2}_{L^{2}(\Omega)}+\frac{\beta}{\rho}||{\bf
u}_{h}^{i+1}||^{3}_{L^{3}(\Omega)}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\hskip 28.45274pt\leq\displaystyle||{\bf f}||_{L^{2}(\Omega)}||{\bf
u}_{h}^{i+1}||_{L^{2}(\Omega)}+||{\bf f}||_{L^{2}(\Omega)}||{\bf
u}_{h,l}||_{L^{2}(\Omega)}+\alpha||{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}||_{L^{2}(\Omega)}||{\bf u}_{h,l}||_{L^{2}(\Omega)}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\hskip
42.67912pt+\displaystyle\frac{\mu}{\rho}K_{M}||{\bf
u}_{h}^{i+1}||_{L^{2}(\Omega)}||{\bf
u}_{h,l}||_{L^{2}(\Omega)}+\frac{\beta}{\rho}||{\bf u}^{i}_{h}-{\bf
u}^{i+1}_{h}||_{L^{3}(\Omega)}||{\bf u}_{h}^{i+1}||_{L^{3}(\Omega)}||{\bf
u}_{h,l}||_{L^{3}(\Omega)}\\\ \hskip
56.9055pt\displaystyle+\frac{\beta}{\rho}||{\bf u}^{i}_{h}-{\bf
u}^{i+1}_{h}||_{L^{3}(\Omega)}||{\bf
u}_{h}^{i+1}||^{2}_{L^{3}(\Omega)}+\frac{\beta}{\rho}||{\bf
u}^{i+1}_{h}||^{2}_{L^{3}(\Omega)}||{\bf
u}_{h,l}||_{L^{3}(\Omega)}.\end{array}$

We deduce by using the relation (3.2), that for any positive numbers
$\varepsilon_{i},i=1\dots,4$ and $\delta_{j},j=1,2$, we have the following
bound:

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{\alpha}{2}||{\bf
u}_{h}^{i+1}||^{2}_{L^{2}(\Omega)}-\frac{\alpha}{2}||{\bf
u}_{h}^{i}||^{2}_{L^{2}(\Omega)}+\frac{\alpha}{2}||{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}||^{2}_{L^{2}(\Omega)}+\frac{\mu}{\rho}K_{m}||{\bf
u}_{h}^{i+1}||^{2}_{L^{2}(\Omega)}+\frac{\beta}{\rho}||{\bf
u}_{h}^{i+1}||^{3}_{L^{3}(\Omega)}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\hskip 28.45274pt\leq\displaystyle\frac{1}{2\varepsilon_{1}}||{\bf
f}||^{2}_{L^{2}(\Omega)}+\frac{1}{2}\varepsilon_{1}||{\bf
u}_{h}^{i+1}||^{2}_{L^{2}(\Omega)}+\frac{1}{2}||{\bf
f}||^{2}_{L^{2}(\Omega)}+\frac{1}{2}||{\bf u}_{h,l}||^{2}_{L^{2}(\Omega)}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\hskip
42.67912pt\displaystyle+\frac{\alpha}{2\varepsilon_{2}}||{\bf
u}_{h}^{i+1}-{\bf
u}_{h}^{i}||^{2}_{L^{2}(\Omega)}+\frac{\alpha}{2}\varepsilon_{2}||{\bf
u}_{h,l}||^{2}_{L^{2}(\Omega)}+\frac{\mu^{2}}{2\rho^{2}\varepsilon_{3}}K_{M}^{2}||{\bf
u}_{h,l}||^{2}_{L^{2}(\Omega)}+\frac{1}{2}\varepsilon_{3}||{\bf
u}_{h}^{i+1}||^{2}_{L^{2}(\Omega)}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\hskip
42.67912pt\displaystyle+\frac{\beta^{2}}{2\rho^{2}\varepsilon_{4}}{C_{I}^{4}}h^{-\frac{2d}{3}}||{\bf
u}_{h,l}||^{2}_{L^{3}(\Omega)}||{\bf u}^{i}_{h}-{\bf
u}^{i+1}_{h}||^{2}_{L^{2}(\Omega)}+\frac{1}{2}\varepsilon_{4}||{\bf
u}_{h}^{i+1}||^{2}_{L^{2}(\Omega)}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\hskip
42.67912pt\displaystyle+\frac{\beta}{3\rho}\big{(}(\frac{1}{\delta_{1}})^{3}{C_{I}^{3}h^{-d/2}}||{\bf
u}^{i+1}_{h}-{\bf u}_{h}^{i}||^{3}_{L^{2}(\Omega)}+2\delta_{1}^{3/2}||{\bf
u}_{h}^{i+1}||^{3}_{L^{3}(\Omega)}\big{)}\\\ \hskip
42.67912pt\displaystyle+\frac{\beta}{3\rho}\big{(}(\frac{1}{\delta_{2}})^{3}||{\bf
u}_{h,l}||^{3}_{L^{3}(\Omega)}+2\delta_{2}^{3/2}||{\bf
u}^{i+1}_{h}||^{3}_{L^{3}(\Omega)}\big{)}.\end{array}$

We choose
$\varepsilon_{1}=\varepsilon_{3}=\varepsilon_{4}=\displaystyle\frac{\mu
K_{m}}{3\rho}$, $\varepsilon_{2}=2$,
$\delta_{1}=\delta_{2}=\displaystyle\big{(}\frac{1}{2}\big{)}^{2/3}$ and we
denote

$C_{1}(\|{\bf
u}_{h}^{i}\|_{L^{2}(\Omega)})=\displaystyle\frac{\alpha}{4}-\frac{3\beta^{2}}{2\rho\mu
K_{m}}{C_{I}^{4}}h^{-\frac{2d}{3}}||{\bf
u}_{h,l}||^{2}_{L^{3}(\Omega)}-\frac{4\beta}{3\rho}{C_{I}^{3}h^{-d/2}L_{2}({\bf
f},\|{\bf u}_{h}^{i}\|_{L^{2}(\Omega)}}),$

which is not necessarily positive at this level.
By using the bound (3.17) , we get the following bound

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{\alpha}{2}||{\bf
u}_{h}^{i+1}||^{2}_{L^{2}(\Omega)}-\frac{\alpha}{2}||{\bf
u}_{h}^{i}||^{2}_{L^{2}(\Omega)}+C_{1}(\|{\bf
u}_{h}^{i}\|_{L^{2}(\Omega)})||{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}||^{2}_{L^{2}(\Omega)}+\frac{\mu K_{m}}{2\rho}||{\bf
u}_{h}^{i+1}||^{2}_{L^{2}(\Omega)}+\frac{\beta}{3\rho}||{\bf
u}_{h}^{i+1}||^{3}_{L^{3}(\Omega)}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\hskip 28.45274pt\leq\displaystyle(\frac{3\rho}{2\mu
K_{m}}+\frac{1}{2})||{\bf
f}||^{2}_{L^{2}(\Omega)}+(\frac{1}{2}+\alpha+{\frac{3\mu K_{M}^{2}}{2\rho
K_{m}}})||{\bf u}_{h,l}||^{2}_{L^{2}(\Omega)}+\frac{4\beta}{3\rho}||{\bf
u}_{h,l}||^{3}_{L^{3}(\Omega)}\\\ \hskip 28.45274pt\leq\displaystyle{\frac{\mu
K_{m}}{2\rho}L_{1}^{2}({\bf f},{\bf u}_{h,l},\alpha).}\end{array}$ (3.19)

We now prove Estimate (3.14) by induction on $i$ under some conditions on
$\alpha$. Starting with the relation (3.10), we suppose that we have

$||{\bf u}_{h}^{i}||_{L^{2}(\Omega)}\leq\displaystyle L_{1}({\bf f},{\bf
u}_{h,l},\alpha).$ (3.20)

We are in one of the following two situations :

  * •

We have $\displaystyle||{\bf u}_{h}^{i+1}||_{L^{2}(\Omega)}\leq||{\bf
u}_{h}^{i}||_{L^{2}(\Omega)}$. We obviously deduce the bound

$||{\bf u}_{h}^{i+1}||_{L^{2}(\Omega)}\leq\displaystyle L_{1}({\bf f},{\bf
u}_{h,l},\alpha),$

from the induction hypothesis.

  * •

We have $\displaystyle||{\bf u}_{h}^{i+1}||_{L^{2}(\Omega)}\geq||{\bf
u}_{h}^{i}||_{L^{2}(\Omega)}$. By using the induction condition (3.20) and the
fact that the function $L_{2}$ is increasing with respect to $\eta$, we chose

$\displaystyle\displaystyle\frac{\alpha}{4}\geq\displaystyle\frac{3\beta^{2}}{2\rho\mu
K_{m}}{C_{I}^{4}}h^{-\frac{2d}{3}}||{\bf
u}_{h,l}||^{2}_{L^{3}(\Omega)}+\frac{4\beta}{3\rho}{C_{I}^{3}}h^{-d/2}{L_{2}({\bf
f},L_{1}({\bf f},{\bf u}_{h,l},\alpha))}$ (3.21)

and we get

$\displaystyle\displaystyle\frac{\alpha}{4}\geq\displaystyle\frac{3\beta^{2}}{2\rho\mu
K_{m}}{C_{I}^{4}}h^{-\frac{2d}{3}}||{\bf
u}_{h,l}||^{2}_{L^{3}(\Omega)}+\frac{4\beta}{3\rho}{C_{I}^{3}}h^{-d/2}{L_{2}({\bf
f},||{\bf u}_{h}^{i}||_{L^{2}(\Omega)})},$ (3.22)

which gives $C_{1}(||{\bf u}_{h}^{i}||_{L^{2}(\Omega)})\geq 0$ and then

$||{\bf u}_{h}^{i+1}||_{L^{2}(\Omega)}\leq\displaystyle L_{1}({\bf f},{\bf
u}_{h,l},\alpha).$

whence we deduce the relation (3.14). The bound (3.15) is a simple consequence
of Equation (3.19) and relation (3.14).

We now focus on the inequality (3.21). It is easy to show that,
$\forall\eta\in\mathbb{R}_{+}$,

$\displaystyle L_{2}({\bf f},\eta)\leq\displaystyle\frac{2\rho}{\mu
K_{m}}\|{\bf
f}\|_{L^{2}(\Omega)^{2}}+2\eta+\frac{\beta}{\rho}C_{I}^{3}h^{-\frac{d}{2}}\frac{\eta^{2}}{\alpha+\frac{\mu
K_{m}}{2\rho}},$ (3.23)

and then to get by using the definition of $L_{1}$,

$L_{2}({\bf f},L_{1}({\bf f},{\bf
u}_{h,l},\alpha))\leq\displaystyle\frac{2\rho}{\mu K_{m}}\|{\bf
f}\|_{L^{2}(\Omega)^{2}}+\frac{2K_{M}}{K_{m}}(\sqrt{\ell_{0}}+\sqrt{\ell_{1}}\sqrt{\alpha})+\frac{2\beta}{\rho}C_{I}^{3}h^{-\frac{d}{2}}\max(\frac{2\rho\ell_{0}}{\mu
K_{m}},\ell_{1}).$

Relation (3.21) and the last inequality allow us to obtain

$\displaystyle\frac{1}{4}\alpha-\frac{3\beta^{2}}{2\rho\mu
K_{m}}{C_{I}^{4}}h^{-\frac{2d}{3}}||{\bf
u}_{h,l}||^{2}_{L^{3}(\Omega)}-\frac{4\beta}{3\rho}{C_{I}^{3}}h^{-\frac{d}{2}}{L_{2}({\bf
f},L_{1}({\bf f},{\bf u}_{h,l},\alpha))\geq\phi(\alpha)},$ (3.24)

with

$\phi(\alpha)=\frac{1}{4}\alpha-\gamma_{1}\sqrt{\alpha}-\gamma_{2},$ (3.25)

where

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus
2.0pt\gamma_{1}&=&\displaystyle\frac{8\beta K_{M}}{3\rho
K_{m}}C_{I}^{3}h^{-\frac{d}{2}}\sqrt{\ell_{1}},\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\gamma_{2}&=&\displaystyle\frac{3\beta^{2}}{2\rho\mu
K_{m}}{C_{I}^{4}}h^{-\frac{2d}{3}}||{\bf
u}_{h,l}||^{2}_{L^{3}(\Omega)}+\frac{8\beta}{3\mu
K_{m}}{C_{I}^{3}}h^{-\frac{d}{2}}\|{\bf f}\|_{L^{2}(\Omega)^{2}},\\\
&&\displaystyle+\frac{8\beta K_{M}}{3\rho
K_{m}}{C_{I}^{3}}h^{-\frac{d}{2}}\sqrt{\ell_{0}}+\frac{8\beta^{2}}{3\rho^{2}}{C_{I}^{6}}h^{-d}\max(\frac{2\rho\ell_{0}}{\mu
K_{m}},\ell_{1}).\end{array}$

Therefore $\phi(\alpha)$ is a polynomial of second degree with respect to
$\sqrt{\alpha}$ and admits the only positive root
$\sqrt{\alpha^{*}}=2(\gamma_{1}+\sqrt{\gamma_{1}^{2}+\gamma_{2}})$. Thus we
get that $\phi(\alpha)\geq 0$ for all $\alpha\geq\alpha^{*}.$ $\hfill\Box$

######  Remark 3.3.

Notice that Theorem 3.2 indicates that $\alpha^{\star}$ blows up as $h^{-d}$,
and as $\beta^{2}$ when respectively $h$ tends to zero, and $\beta$ approaches
$+\infty$.

The next result provide the convergence of the solution $({\bf
u}_{h}^{i},p_{h}^{i})$ of Problem $(\ref{V1hi})$ in $L^{2}(\Omega)^{d}\times
L^{2}(\Omega)$ to the unique solution $({\bf u}_{h},p_{h})$ of the Problem
$(\ref{V1h})$.

######  Theorem 3.4.

Assume that there exists $\beta_{0}>0$ such that, for every element
$\kappa\in{\mathcal{T}}_{h}$, we have

$h_{\kappa}\geq\beta_{0}h,$

(which means that the family of triangulations is uniformly regular). Under
the assumptions of Theorem 3.2, and if $\alpha$ satisfies the condition

$\displaystyle\alpha>\displaystyle\max(\alpha^{\star},\alpha^{\star\star})$
(3.26)

$\alpha^{\star\star}$ is the largest positive zero of the polynomial function:

$\tilde{\phi}(\alpha)=-\tilde{\gamma}_{0}-\tilde{\gamma}_{1}\alpha-\tilde{\gamma}_{2}\alpha^{2}+\frac{1}{8}\alpha^{3},$
(3.27)

where $\tilde{\gamma}_{i}(h,{\bf f},{\bf u}_{h,l}),i=2,3,4$ are strictly
positive constants given by

$\displaystyle\tilde{\gamma}_{0}$ $\displaystyle=$
$\displaystyle\ell_{0}^{2}\tilde{C}h^{-2d},$ $\displaystyle\tilde{\gamma}_{1}$
$\displaystyle=$ $\displaystyle 2\ell_{0}\ell_{1}\tilde{C}h^{-2d},$
$\displaystyle\tilde{\gamma}_{2}$ $\displaystyle=$
$\displaystyle\ell_{1}^{2}\tilde{C}h^{-2d},$ $\displaystyle\tilde{C}$
$\displaystyle=$ $\displaystyle\frac{9\beta^{4}C_{I}^{12}}{32\mu
K_{m}\rho^{3}}$

then the sequence of solutions $({\bf u}_{h}^{i},p_{h}^{i})$ of Problem
$(\ref{V1hi})$ converges in $L^{2}(\Omega)^{d}\times L^{2}(\Omega)$ to the
solution $({\bf u}_{h},p_{h})$ of Problem $(\ref{V1h})$.

Proof. We take the difference between the equations (3.9) and (3.5) with ${\bf
v}_{h}={\bf u}^{i+1}_{h}-{\bf u}_{h}$ and we obtain the equation

$\begin{array}[]{rr}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{\alpha}{2}||{\bf u}_{h}^{i+1}-{\bf
u}_{h}||^{2}_{L^{2}(\Omega)}-\frac{\alpha}{2}||{\bf u}_{h}^{i}-{\bf
u}_{h}||^{2}_{L^{2}(\Omega)}+\frac{\alpha}{2}||{\bf u}_{h}^{i+1}-{\bf
u}^{i}_{h}||^{2}_{L^{2}(\Omega)}+\frac{\mu}{\rho}\int_{\Omega}K^{-1}({\bf
u}_{h}^{i+1}-{\bf u}_{h})^{2}\,d{\bf{x}}\\\
+\displaystyle\frac{\beta}{\rho}(|{\bf u}_{h}^{i}|{\bf u}_{h}^{i+1}-|{\bf
u}_{h}|{\bf u}_{h},{\bf u}_{h}^{i+1}-{\bf u}_{h})=0.\end{array}$

The last term in the previous equation, denoted by $T$, can be decomposed as

$T=\displaystyle\frac{\beta}{\rho}((|{\bf u}_{h}^{i}|-|{\bf u}^{i+1}_{h}|){\bf
u}_{h}^{i+1},{\bf u}_{h}^{i+1}-{\bf u}_{h})+\frac{\beta}{\rho}(|{\bf
u}_{h}^{i+1}|{\bf u}_{h}^{i+1}-|{\bf u}_{h}|{\bf u}_{h},{\bf u}_{h}^{i+1}-{\bf
u}_{h})=T_{1}+T_{2}.$

We denote by $T_{1}$ and $T_{2}$, respectively the first and the second terms
in the right-hand side of the last equation. Using (2.7), we have $T_{2}\geq
0$. Then we derive by using (3.15), (3.2) and (1.4),

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{\alpha}{2}||{\bf u}_{h}^{i+1}-{\bf
u}_{h}||^{2}_{L^{2}(\Omega)}-\frac{\alpha}{2}||{\bf u}_{h}^{i}-{\bf
u}_{h}||^{2}_{L^{2}(\Omega)}+\frac{\alpha}{2}||{\bf u}_{h}^{i+1}-{\bf
u}^{i}_{h}||^{2}_{L^{2}(\Omega)}+\frac{K_{m}\mu}{\rho}\|{\bf u}_{h}^{i+1}-{\bf
u}_{h}\|^{2}_{L^{2}(\Omega)}+T_{2}\leq|T_{1}|\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\hskip 28.45274pt\leq\displaystyle\frac{\beta}{\rho}\int_{\Omega}|{\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i}|\,|{\bf u}_{h}^{i+1}|\,|{\bf u}_{h}^{i+1}-{\bf
u}_{h}|d{\textbf{x}}\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\hskip
28.45274pt\leq\displaystyle\frac{\beta}{\rho}||{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}||_{L^{3}(\Omega)}\;||{\bf u}_{h}^{i+1}||_{L^{3}(\Omega)}\,||{\bf
u}_{h}^{i+1}-{\bf u}_{h}||_{L^{3}(\Omega)}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\hskip
28.45274pt\leq\displaystyle\frac{\beta}{\rho}C_{I}^{2}h^{-\frac{d}{3}}||{\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i}||_{L^{2}(\Omega)}\;||{\bf
u}_{h}^{i+1}||_{L^{3}(\Omega)}\,||{\bf u}_{h}^{i+1}-{\bf
u}_{h}||_{L^{2}(\Omega)}\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\hskip
28.45274pt\leq\displaystyle\frac{\beta}{\rho}\big{(}\frac{3\mu
K_{m}}{2\beta}\big{)}^{1/3}C_{I}^{2}h^{-\frac{d}{3}}L^{2/3}_{1}({\bf f},{\bf
u}_{h,l},\alpha)||{\bf u}_{h}^{i+1}-{\bf u}_{h}^{i}||_{L^{2}(\Omega)}\,||{\bf
u}_{h}^{i+1}-{\bf u}_{h}||_{L^{2}(\Omega)}.\end{array}$

We denote by $C=\displaystyle\frac{\beta}{\rho}\big{(}\frac{3\mu
K_{m}}{2\beta}\big{)}^{1/3}C_{I}^{2}L^{2/3}_{1}({\bf f},{\bf u}_{h,l},\alpha)$
and we use the inequality
$ab\leq\displaystyle\frac{1}{2\varepsilon}a^{2}+\frac{\varepsilon}{2}b^{2}$
(with $\varepsilon=\displaystyle\frac{K_{m}\mu}{\rho}$) to obtain the
following bound

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{\alpha}{2}||{\bf u}_{h}^{i+1}-{\bf
u}_{h}||^{2}_{L^{2}(\Omega)}-\frac{\alpha}{2}||{\bf u}_{h}^{i}-{\bf
u}_{h}||^{2}_{L^{2}(\Omega)}+\frac{\alpha}{2}||{\bf u}_{h}^{i+1}-{\bf
u}^{i}_{h}||^{2}_{L^{2}(\Omega)}+\frac{K_{m}\mu}{2\rho}{||{\bf
u}_{h}^{i+1}-{\bf u}_{h}||^{2}_{L^{2}(\Omega)}}\\\ \hskip
56.9055pt\leq\displaystyle\frac{\rho C^{2}}{2K_{m}\mu}h^{-\frac{2d}{3}}||{\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i}||^{2}_{L^{2}(\Omega)}.\end{array}$

We choose

$\displaystyle\frac{\alpha}{2}>\displaystyle\frac{\rho
C^{2}}{2K_{m}\mu}h^{-\frac{2d}{3}},$ (3.28)

denote by $C_{1}=\displaystyle\frac{1}{2}(\alpha-\displaystyle\frac{\rho
C^{2}}{K_{m}\mu}h^{-\frac{2d}{3}})$ and obtain

$\displaystyle\frac{\alpha}{2}||{\bf u}_{h}^{i+1}-{\bf
u}_{h}||^{2}_{L^{2}(\Omega)}-\frac{\alpha}{2}||{\bf u}_{h}^{i}-{\bf
u}_{h}||^{2}_{L^{2}(\Omega)}+C_{1}||{\bf u}_{h}^{i+1}-{\bf
u}^{i}_{h}||^{2}_{L^{2}(\Omega)}+\frac{\mu K_{m}}{2\rho}{||{\bf
u}_{h}^{i+1}-{\bf u}_{h}||^{2}_{L^{2}(\Omega)}}\leq 0.$ (3.29)

We clearly have that (3.28) is satisfied if

$\tilde{\phi}(\alpha)>0.$ (3.30)

Simple calculation shows that $\tilde{\phi}^{\prime}(\alpha)$ is a polynomial
function of degree two that is positive for large $\alpha$ and posses a unique
positive zero. Since $\tilde{\phi}(0)<0$, we deduce that
$\tilde{\phi}(\alpha)$ vanishes at least once on $\mathbb{R}_{+}$. By denoting
$\alpha^{\star\star}$ the largest root of $\tilde{\phi}(\alpha)$, we deduce
that $\tilde{\phi}(\alpha)\geq 0$ for $\alpha\geq\alpha^{\star\star}$.

We deduce from (3.29) that, for all $i\geq 1$, we have (if $||{\bf
u}_{h}^{i}-{\bf u}_{h}||_{L^{2}(\Omega)}\neq 0$)

$||{\bf u}_{h}^{i+1}-{\bf u}_{h}||_{L^{2}(\Omega)}<||{\bf u}_{h}^{i}-{\bf
u}_{h}||_{L^{2}(\Omega)},$

and we deduce the convergence of the sequence $({\bf u}_{h}^{i+1}-{\bf
u}_{h})$ in $L^{2}(\Omega)^{d}$ and then the convergence of the sequence ${\bf
u}_{h}^{i}$ in $L^{2}(\Omega)^{d}$. By taking the limit of (3.29) and
remarking that $T_{2}\geq 0$, we get

$\underset{{\small i\rightarrow+\infty}}{\lim}\Big{(}{||{\bf u}_{h}^{i+1}-{\bf
u}_{h}||^{2}_{L^{2}(\Omega)}}\Big{)}\leq 0.$

We deduce then that ${||{\bf u}_{h}^{i+1}-{\bf u}_{h}||_{L^{2}(\Omega)}}$
converges to $0$ and ${\bf u}_{h}^{i+1}$ converges to ${\bf u}_{h}$ in
$L^{2}(\Omega)^{d}$.
For the convergence of the pressure, we take the difference between the
equations (3.9) and (3.5) and we obtain for all ${\bf v}_{h}\in X_{h}$ the
equation

$\begin{array}[]{rr}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\int_{\Omega}\nabla(p_{h}^{i+1}-p_{h}){\bf
v}_{h}\,d{\bf{x}}=\displaystyle-\alpha\int_{\Omega}({\bf u}_{h}^{i+1}-{\bf
u}^{i}_{h}){\bf v}_{h}\,d{\bf{x}}+\frac{\mu}{\rho}\int_{\Omega}K^{-1}({\bf
u}_{h}-{\bf u}_{h}^{i+1}){\bf v}_{h}\,d{\bf{x}}\\\ \hskip
56.9055pt\displaystyle+\frac{\beta}{\rho}((|{\bf u}_{h}|-|{\bf
u}_{h}^{i}|){\bf u}_{h},{\bf v}_{h})+\frac{\beta}{\rho}(|{\bf u}^{i}_{h}|({\bf
u}_{h}-{\bf u}_{h}^{i+1}),{\bf v}_{h}).\end{array}$

We get by using the inverse inequality (3.2) the following:

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{\displaystyle\Big{|}\int_{\Omega}\nabla(p_{h}^{i+1}-p_{h}){\bf
v}_{h}\,d{\bf{x}}\Big{|}}{||{\bf
v}_{h}||_{L^{3}(\Omega)}}&\leq&\displaystyle(\alpha||{\bf u}_{h}^{i}-{\bf
u}_{h}^{i+1}||_{L^{2}(\Omega)}+\frac{\mu K_{m}}{\rho}||{\bf u}_{h}-{\bf
u}_{h}^{i+1}||_{L^{2}(\Omega)})\frac{||{\bf v}_{h}||_{L^{2}(\Omega)}}{||{\bf
v}_{h}||_{L^{3}(\Omega)}}\\\
&&\displaystyle{+\frac{\beta}{\rho}C_{I}h^{-\frac{d}{6}}||{\bf u}_{h}-{\bf
u}_{h}^{i}||_{L^{2}(\Omega)}(||{\bf u}_{h}||_{L^{3}(\Omega)}+||{\bf
u}_{h}^{i}||_{L^{3}(\Omega)}).}\end{array}$

Owning the inf-sup condition (3.4), we deduce the relation

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle||\nabla(p_{h}^{i+1}-p_{h})||_{L^{3/2}(\Omega)}\leq\displaystyle\frac{1}{\beta_{u}}\Big{(}\alpha{|\Omega|^{1/6}}||{\bf
u}_{h}^{i}-{\bf u}_{h}^{i+1}||_{L^{2}(\Omega)}+{|\Omega|^{1/6}}\frac{\mu
K_{m}}{\rho}||{\bf u}_{h}-{\bf u}_{h}^{i+1}||_{L^{2}(\Omega)}\\\ \hskip
56.9055pt\displaystyle{+\frac{\beta}{\rho}C_{I}h^{-\frac{d}{6}}||{\bf
u}_{h}-{\bf u}_{h}^{i}||_{L^{2}(\Omega)}(||{\bf
u}_{h}||_{L^{3}(\Omega)}+||{\bf
u}_{h}^{i}||_{L^{3}(\Omega)})\Big{)}.}\end{array}$

Thus for a given mesh (given $h$), the strong convergence of ${\bf u}_{h}^{i}$
to ${\bf u}_{h}$ in $L^{2}(\Omega)^{d}$ implies the strong convergence of
$\nabla p_{h}^{i}$ to $\nabla p_{h}$ in $L^{\frac{3}{2}}(\Omega)$.
Furthermore, the fact that $p_{h}^{i}$ and $p_{h}$ are in the discrete space
of ${\rm IP}_{1}$ finite elements $M_{h}\subset L^{2}_{0}(\Omega)$ allows us
to deduce the strong convergence of $p_{h}^{i}$ to $p_{h}$ in $L^{2}(\Omega)$
.$\hfill\Box$

######  Remark 3.5.

One can show that

$\displaystyle\alpha^{\star\star}\geq\max\left(8\ell_{1}^{2}C_{h},4(\ell_{0}\ell_{1}C_{h})^{\frac{1}{2}},2(\ell_{0}^{2}C_{h})^{\frac{1}{3}}\right),$

where
$C_{h}=\left(\frac{\beta^{2}C_{I}^{4}}{2K_{m}\mu\rho}\right)^{3}h^{-2d}$.
Consequently $\alpha^{\star\star}$ tends to infinity when $h$ becomes small or
$\beta$ approaches $+\infty$.

######  Remark 3.6.

The condition $(\ref{cond1f})$ supposes that the initial value of the
algorithm is small related to the data $f$. We can for example consider
$u_{h}^{0}=0$.

######  Remark 3.7.

Theorems 3.4 and 3.2 require the conditions
$\alpha>\max(\alpha^{\star},\alpha^{\star\star})$ and (3.26) to get the
convergence of the numerical scheme (3.9). This conditions cannot be computed
easily in practice, especially for the numerical investigations. In fact, This
result of convergence states that for a given mesh (for a given $h$), the
iterative solution $(u_{h}^{i},p_{h}^{i})$ converges to $(u_{h},p_{h})$ when
$i\rightarrow+\infty$.

###  3.2. A posteriori error estimates

As usual, for a posteriori  error estimates, we introduce the following
notations. We denote by

  * •

$\Gamma_{h}^{i}$ the set of edges (when $d=2$) or faces (when $d=3$) of
$\kappa$ that are not contained in $\partial\Omega$.

  * •

$\Gamma_{h}^{b}$ the set of edges (when $d=2$) or faces (when $d=3$) of
$\kappa$ which are contained in $\partial\Omega$.

For every element $\kappa$ in $\mathcal{T}_{h}$, we denote by $w_{\kappa}$ the
union of elements $K$ of $\mathcal{T}_{h}$ such that $\kappa\cap K\neq\phi$.
Furthermore, for every edge (when $d=2$) or face (when $d=3$) $e$ of the mesh
$\mathcal{T}_{h}$, we denote by

  * •

$\omega_{e}$ the union of elements of $\mathcal{T}_{h}$ adjacent to $e$.

  * •

$[\cdot]_{e}$ the jump through $e\in\Gamma_{h}^{i}$.

From now on, to simplify, we set $d=3$. The extension to two dimensions is
straightforward and simpler. We suppose also that $b\in L^{3}(\Omega)$ and
$g\in L^{3}(\Gamma)$.

In this and the next sections, the a posteriori error estimates are
established when the solution is slightly smoother.

Let $R_{h}$ be a Clément-type interpolation operator [18]. We have the
following error estimate: for all $\kappa$ in $\mathcal{T}_{h}$, for all $e$
in $\partial\kappa$ and for all $q\in W^{1,3/2}(\Omega)$,

$\|q-R_{h}q\|_{L^{3/2}(\kappa)}\leq
C_{\kappa}h_{\kappa}|q|_{W^{1,3/2}(w_{\kappa})}$ (3.31)

and

$\|q-R_{h}q\|_{L^{3/2}(e)}\leq C_{e}h_{e}^{1/3}|q|_{W^{1,3/2}(w_{e})},$ (3.32)

where $C_{e}$ and $C_{\kappa}$ are positive constants independent of $h$.

######  Remark 3.8.

Relations (3.31) and (3.32) are a direct consequence of the following two
properties:

  * •

for all integers $\ell$, $0\leq\ell\leq 2$, and for all $p$, $0\leq
p\leq+\infty$, there exists a constant $C$, independent of $h_{\kappa}$, such
that for all $\kappa\in\mathcal{T}_{h}$ and all function $q\in
W^{\ell,p}(w_{\kappa})$, the following inequalities hold (see for instance
[7], Theorem 1):

$||q-R_{h}q||_{L^{p}(\kappa)}\leq
Ch_{\kappa}^{\ell}|q|_{W^{\ell,p}(w_{\kappa})}$

and when $\ell\geq 1$

$|q-R_{h}q|_{W^{1,p}(\kappa)}\leq
Ch_{\kappa}^{\ell-1}|q|_{W^{\ell,p}(w_{\kappa})}.$

  * •

Let $s\in]0,1[$ and $p>\displaystyle\frac{1}{s}$ with $p\in[1,+\infty[$ or
$s=1$ with $p\in[1,+\infty]$. Then there exists $c$, uniform with respect to
the mesh such that the following trace inequality (see [22], Lemma 7.2) holds
for all $q\in W^{s,p}(\kappa)$ and all $\kappa\in\mathcal{T}_{h}$:

$||q||_{L^{p}(e)}\leq
c\big{(}h_{\kappa}^{-\frac{1}{p}}||q||_{L^{p}(\kappa)}+h_{\kappa}^{s-\frac{1}{p}}|q|_{W^{s,p}(\kappa)}\big{)}.$

####  3.2.1. Upper error bound

In order to establish upper bounds, we introduce, on every edge ($d=2$) or
face ($d=3$) $e$ of the mesh, the function

$\phi_{h,1}^{e}=\left\\{\begin{array}[]{ccl}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{1}{2}\,[\textbf{u}_{h}^{i+1}\cdot\textbf{n}]_{e}\quad\mbox{if
}e\in\Gamma_{h}^{i},\\\ \textbf{u}_{h}^{i+1}\cdot\textbf{n}-g_{h}\quad\mbox{if
}e\in\Gamma_{h}^{b},\end{array}\right.$ (3.33)

where $g_{h}$ is an approximation of $g$ which is constant on each $e$.

A standard calculation shows that the solutions of the problems (2.9) and
(3.9) verify for all $({\bf v},q)\in X\times M$ and $({\bf v}_{h},q_{h})\in
X_{h}\times M_{h}$:

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{\mu}{\rho}\int_{\Omega}K^{-1}({\bf u}-{\bf
u}^{i+1}_{h})\cdot{\bf v}\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}(|{\bf
u}|{\bf u}-|{\bf u}^{i}_{h}|{\bf u}^{i+1}_{h})\cdot{\bf
v}\,d{\bf{x}}+\displaystyle\int_{\Omega}\nabla(p-p_{h}^{i+1})\cdot\textbf{v}\,d{\bf{x}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\hskip
28.45274pt=\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\big{[}\int_{\kappa}(-\nabla
p_{h}^{i+1}-{{\alpha({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})}}-{{\frac{\mu}{\rho}K^{-1}{\bf
u}^{i+1}_{h}}}-\frac{\beta}{\rho}|{\bf u}^{i}_{h}|{\bf u}^{i+1}_{h}+{\bf
f}_{h})\cdot({\bf v}-{\bf v}_{h})\,d{\bf{x}}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\hskip 142.26378pt+\displaystyle\int_{\kappa}({\bf f}-{\bf
f}_{h})\cdot{\bf v}\,d{\bf{x}}\,+\alpha\displaystyle\int_{\kappa}({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})\cdot{\bf v}\,d{\bf{x}}\big{]},\end{array}$
(3.34)

and by using the second equations of the systems (2.9) and (3.9), the fact
that ${\rm div}\,{\bf u}_{h}^{i+1}=0$ (as in each element $\kappa$ we have:
${\bf u}_{h}^{i+1}\in{\rm IP}_{0}^{d}(\kappa)$), and by applying the
integration by parts to the term
$\displaystyle\int_{\Omega}\nabla(q-q_{h})\cdot{\bf u}_{h}^{i+1}d{\bf{x}}$ and
using the definition of $\phi_{h,1}^{e}$ we get,

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\int_{\Omega}\nabla
q\cdot(\textbf{u}-\textbf{u}_{h}^{i+1})\,d{\bf{x}}&=&\displaystyle\int_{\Omega}\nabla
q\cdot\textbf{u}\,d{\bf{x}}-\int_{\Omega}\nabla
q_{h}\cdot\textbf{u}_{h}^{i+1}\,d{\bf{x}}-\int_{\Omega}\nabla(q-q_{h})\cdot\textbf{u}_{h}^{i+1}\,d{\bf{x}}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt&=&\displaystyle-\sum_{\kappa\in\mathcal{T}_{h}}\Big{[}\int_{\kappa}(q-q_{h})(b-b_{h})\,d{\bf{x}}\,+\int_{\kappa}b_{h}(q-q_{h})\,d{\bf{x}}\\\
&&\qquad\qquad\displaystyle+\displaystyle\sum_{e\in\partial\kappa}\int_{e}\phi_{h,1}^{e}(q-q_{h})\,ds\,\Big{]}-\sum_{e\in\Gamma_{h}^{b}}\int_{e}(g_{h}-g)(q-q_{h})\,ds,\end{array}$
(3.35)

where ${\bf f}_{h}$ (resp. $b_{h}$) is an approximation of ${\bf f}$ (resp.
$b$) which is constant on each element $\kappa$ of $\mathcal{T}_{h}$.
From the error equations (3.34) and (3.35), we deduce the following error
indicators for each $\kappa\in\mathcal{T}_{h}$,

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\eta_{\kappa,i}^{(L)}&=&||{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}||_{L^{2}(\kappa)},\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\eta_{\kappa,i}^{(D_{1})}&=&\displaystyle\|-\nabla
p_{h}^{i+1}-{{\alpha({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})-\frac{\mu}{\rho}K^{-1}{\bf u}^{i+1}_{h}}}-\frac{\beta}{\rho}|{\bf
u}^{i}_{h}|{\bf u}^{i+1}_{h}+{\bf f}_{h}\,\|_{L^{2}(\kappa)},\\\ \vskip 6.0pt
plus 2.0pt minus
2.0pt\displaystyle\eta_{\kappa,i}^{(D_{2})}&=&h_{\kappa}\|b_{h}\|_{L^{3}(\kappa)}+\displaystyle\sum_{e\in\partial\kappa}h_{e}^{\frac{1}{3}}\|\phi_{h,1}^{e}\|_{L^{3}(e)}.\end{array}$
(3.36)

The term $h_{\kappa}\parallel b_{h}\parallel_{L^{3}(\kappa)}$ which appears in
$\eta_{\kappa,i}^{(D_{2})}$ is an indicator since it represents the quantity
$h_{\kappa}(\parallel b_{h}-\operatorname{div}({\bf
u}_{h}^{i+1})\parallel_{L^{3}(\kappa)})$ as $\operatorname{div}({\bf
u}_{h}^{i+1})=0$ in each element $\kappa$.

In order to establish an a posteriori error estimate, we need to bound the
numerical solution ${\bf u}_{h}^{i+1}$ in $L^{6}(\Omega)$ which is the subject
of the next lemma.

######  Lemma 3.9.

Let $d=3$ and let the mesh satisfy (3.1). Under the assumptions of Theorems
3.1, 3.2 and 3.4, and if the exact velocity ${\bf u}\in W^{1.6}(\Omega)^{d}$,
there exists an integer $i_{0}$ depending on $h$ such that for all $i\geq
i_{0}$, the solution $({\bf u}_{h}^{i+1},p_{h}^{i+1})$ of $\eqref{V1hi}$
verifies the following bound:

$\|{\bf u}^{i+1}_{h}\|_{L^{6}(\Omega)}\leq\widehat{C}({\bf u},p),$ (3.37)

where $\widehat{C}$ is a constant depending on the exact solution $({\bf
u},p)$ of (2.9).

###### Proof.

We consider the case $d=3$. Let $({\bf u},p)$ be the solution of (2.9), $({\bf
u}_{h},p_{h})$ the solution of $\eqref{V1h}$ and $({\bf
u}_{h}^{i+1},p_{h}^{i+1})$ the solution of $\eqref{V1hi}$.
By using (3.2) (for $p=6$), the properties of the operator $R_{h}$ and the a
priori error estimate (3.8), the term $\|{\bf u}^{i+1}_{h}\|_{L^{6}(\Omega)}$
can be bounded as following:

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus 2.0pt\|{\bf
u}^{i+1}_{h}\|_{L^{6}(\Omega)}&\leq&\|{\bf u}^{i+1}_{h}-{\bf
u}_{h}\|_{L^{6}(\Omega)}+\|{\bf u}_{h}-R_{h}({\bf
u})\|_{L^{6}(\Omega)}+\|R_{h}({\bf u})-{\bf u}\|_{L^{6}(\Omega)}+\|{\bf
u}\|_{L^{6}(\Omega)}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt&\leq&C\big{(}h^{-1}(\|{\bf u}^{i+1}_{h}-{\bf
u}_{h}\|_{L^{2}(\Omega)}+\|{\bf u}_{h}-R_{h}({\bf
u})\|_{L^{2}(\Omega)})+h\|{\bf u}\|_{W^{1,6}(\Omega)^{6}}\big{)}+\|{\bf
u}\|_{L^{6}(\Omega)}\\\ &\leq&C\big{(}h^{-1}(\|{\bf u}^{i+1}_{h}-{\bf
u}_{h}\|_{L^{2}(\Omega)}+\|{\bf u}_{h}-{\bf u}\|_{L^{2}(\Omega)}+\|{\bf
u}-R_{h}({\bf u})\|_{L^{2}(\Omega)}\big{)}+h\|{\bf
u}\|_{W^{1,6}(\Omega)}\big{)}\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\vskip
6.0pt plus 2.0pt minus 2.0pt&&\displaystyle+\|{\bf u}\|_{L^{6}(\Omega)}\\\
&\leq&C(h^{-1}\|{\bf u}^{i+1}_{h}-{\bf u}_{h}\|_{L^{2}(\Omega)}+C_{u}({\bf
u},p)+|{\bf u}|_{W^{1,2}(\Omega)}+h\|{\bf u}\|_{W^{1,6}(\Omega)}+\|{\bf
u}\|_{L^{6}(\Omega)}).\end{array}$ (3.38)

As ${\bf u}^{i+1}_{h}$ converges to ${\bf u}_{h}$ in $L^{2}(\Omega)$, then
there exists an integer $i_{0}$ depending on $h$ such that for all $i\geq
i_{0}$ we have

$\|{\bf u}^{i+1}_{h}-{\bf u}_{h}\|_{L^{2}(\Omega)}\leq h.$ (3.39)

Then, Equation (3.38) gives by using (3.39) the following bound for all $i\geq
i_{0}$:

$\|{\bf u}^{i+1}_{h}\|_{L^{6}(\Omega)}\leq\widehat{C}({\bf u},p).$ (3.40)

∎

Our main goal is to get an upper bound of the error between the exact solution
$({\bf u},p)$ of (2.9) and the numerical solution $({\bf
u}_{h}^{i+1},p_{h}^{i+1})$ of $\eqref{V1hi}$. To get this desired result, we
start by the following lemma which can be proved by using the inf-sup
condition (2.5).

######  Lemma 3.10.

There exists a velocity ${\bf v}_{r}$ in $X$ that solves the following
equation : $\forall q\in M$

$\displaystyle\int_{\Omega}\nabla q\cdot{\bf v}_{r}\,d{\bf{x}}$
$\displaystyle=\sum_{\kappa\in\mathcal{T}_{h}}\big{[}\int_{\kappa}(q-R_{h}(q))(-b+b_{h})\,d{\bf{x}}-\int_{\kappa}b_{h}(q-R_{h}(q))\,d{\bf{x}}$
(3.41)
$\displaystyle-\sum_{e\in\partial\kappa}\int_{e}\phi_{h,1}^{e}(q-R_{h}(q))ds\big{]}-\sum_{e\in\Gamma_{h}^{b}}\int_{e}(g_{h}-g)(q-R_{h}(q))ds$

satisfying the following bound

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus 2.0pt\|{\bf
v}_{r}\|_{L^{3}(\Omega)}&\leq&\displaystyle\hat{C}_{2}\big{(}\sum_{\kappa\in\mathcal{T}_{h}}\big{[}(\eta_{\kappa,i}^{(D_{2})})+h_{\kappa}||b-b_{h}||_{L^{3}(\kappa)}\big{]}+\sum_{e\in\Gamma_{h}^{b}}h_{e}^{1/3}\|g_{h}-g\|_{L^{3}(e)}\big{)}.\end{array}$
(3.42)

###### Proof.

Equation (3.35) with $q_{h}=R_{h}(q)$ and the inf-sup condition (2.5) imply
that there exists a ${\bf v}_{r}\in X$ such that (3.41) is verified and
satisfying

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus
2.0pt\|\textbf{v}_{r}\|_{L^{3}(\Omega)}&\leq&\underset{q\in
M}{\sup}\frac{1}{|\nabla
q|_{L^{3/2}(\Omega)}}\Big{|}\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\big{[}\int_{\kappa}(q-R_{h}(q))(-b+b_{h})\,d{\bf{x}}\,-\int_{\kappa}b_{h}(q-R_{h}(q))\,d{\bf{x}},\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt&&\displaystyle-\displaystyle\sum_{e\in\partial\kappa}\int_{e}\phi_{h,1}^{e}(q-R_{h}(q))\,ds\,\big{]}-\sum_{e\in\Gamma_{h}^{b}}\int_{e}(g_{h}-g)(q-R_{h}(q))\,ds\Big{|}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt&&\leq\underset{q\in
M}{\sup}\frac{1}{|\nabla
q|_{L^{3/2}(\Omega)}}\Big{|}\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\|q-R_{h}(q)\|_{L^{3/2}(\kappa)}\,\|b-b_{h}\|_{L^{3}(\kappa)}+\|q-R_{h}(q)\|_{L^{3/2}(\kappa)}\,\|b_{h}\|_{L^{3}(\kappa)}\\\
&&\displaystyle+\displaystyle\sum_{e\in\partial\kappa}\|\phi_{h,1}^{e}\|_{L^{3}(e)}\,\|q-R_{h}(q)\|_{L^{3/2}(e)}\big{]}+\sum_{e\in\Gamma_{h}^{b}}\|g_{h}-g\|_{L^{3}(e)}\|q-R_{h}(q)\|_{L^{3/2}(e)}\Big{|}.\end{array}$
(3.43)

Thus, from the properties of the operator $R_{h}$, the regularity of
$\mathcal{T}_{h}$ and the following Holder’s inequality ($p=3/2$, $q=3$)

$\displaystyle\sum_{k=1}^{n}a_{k}b_{k}\leq\displaystyle\big{(}\sum_{k=1}^{n}a_{k}^{p}\big{)}^{1/p}\big{(}\sum_{k=1}^{n}b_{k}^{q}\big{)}^{1/q},$
(3.44)

we infer after cubing the last relation:

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus
2.0pt\|\textbf{v}_{r}\|^{3}_{L^{3}(\Omega)}&\leq&\displaystyle
C_{2}\big{(}\sum_{\kappa\in\mathcal{T}_{h}}\big{[}(\eta_{\kappa,i}^{(D_{2})})^{3}+h_{\kappa}^{3}||b-b_{h}||^{3}_{L^{3}(\kappa)}\big{]}+\sum_{e\in\Gamma_{h}^{b}}h_{e}\|g_{h}-g\|^{3}_{L^{3}(e)}\big{)}.\end{array}$

Finally, we obtain (3.42) by taking the cubic root of the previous inequality.
∎

######  Theorem 3.11.

Under the assumptions of Lemma 3.9, there exists an integer $i_{0}$ depending
on $h$ such that for all $i\geq i_{0}$, the solutions $({\bf u},p)$ of (2.9)
and $({\bf u}_{h}^{i+1},p_{h}^{i+1})$ of $\eqref{V1hi}$ verify the following
error inequalities:

$\begin{split}\|{\bf v}_{r}\|_{L^{3}(\Omega)}+\|{\bf u}-{\bf
u}_{h}^{i+1}\|_{L^{2}(\Omega)}+\left\|{\bf u}-{\bf
u}_{h}^{i+1}\right\|_{L^{3}(\Omega)^{d}}^{3/2}\leq
C\displaystyle\Big{[}\sum_{\kappa\in\mathcal{T}_{h}}\Big{(}\eta_{\kappa,i}^{(D_{1})}+\eta_{\kappa,i}^{(D_{2})}+\eta_{\kappa,i}^{(L)}\Big{)}\\\
\hskip
56.9055pt+\sum_{\kappa\in\mathcal{T}_{h}}\Big{(}\|\textbf{f}-\textbf{f}_{h}\|_{L^{2}(\kappa)}+h_{k}\|b-b_{h}\|_{L^{3}(\kappa)}+\sum_{e\in\Gamma_{h}^{b}}h_{e}^{\frac{1}{3}}||g_{h}-g||_{L^{3}(e)}\Big{)}\Big{]},\end{split}$
(3.45)

where $\mathbf{z_{0}}={\bf u}-{\bf u}_{h}^{i+1}-{\bf v}_{r}$ and $C$ is a
constant depending on ${\bf u}$, ${\bf f}_{h}$ (resp. $b_{h}$) is an
approximation of ${\bf f}$ (resp. $b$) which is constant on each element
$\kappa$ of $\mathcal{T}_{h}$, and $g_{h}$ is an approximation of $g$ which is
constant on each face $e$ of $\mathcal{T}_{h}\cap\Gamma$.

###### Proof.

The velocity error equation (3.34) can be written as

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{\mu}{\rho}\int_{\Omega}K^{-1}({\bf u}-{\bf
u}^{i+1}_{h})\cdot{\bf v}\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}(|{\bf
u}|{\bf u}-|{\bf u}^{i+1}_{h}|{\bf u}^{i+1}_{h})\cdot{\bf
v}\,d{\bf{x}}+\displaystyle\int_{\Omega}\nabla(p-p_{h}^{i+1})\cdot\textbf{v}\,d{\bf{x}}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle=-\frac{\beta}{\rho}\int_{\Omega}((|{\bf u}^{i+1}_{h}|-|{\bf
u}^{i}_{h}|){\bf u}^{i+1}_{h})\cdot{\bf
v}\,d{\bf{x}}+\displaystyle\alpha\int_{\Omega}({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})\cdot{\bf v}\,d{\bf{x}}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt+\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\big{[}\int_{\kappa}(-\nabla
p_{h}^{i+1}-{{\frac{\mu}{\rho}K^{-1}{\bf u}_{h}^{i+1}}}-\alpha({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})-\frac{\beta}{\rho}|{\bf u}_{h}^{i}|{\bf
u}_{h}^{i+1}+{\bf f}_{h})\cdot({\bf v}-{\bf
v}_{h})\,d{\bf{x}}+\int_{\kappa}({\bf f}-{\bf f}_{h})\cdot{\bf
v}\,d{\bf{x}}\big{]}.\end{array}$ (3.46)

Now, to simplify we set
$\textbf{z}_{0}=\textbf{u}-\textbf{u}_{h}^{i+1}-\textbf{v}_{r}$ and we test
(3.46) with $\textbf{v}=\textbf{z}_{0}$ and $\textbf{v}_{h}={\bf 0}$.
By construction, (3.41) and (3.35) imply that we have with $q_{h}=R_{h}q$

$\mbox{for all }q\in M,\qquad\qquad\int_{\Omega}\nabla
q\cdot\textbf{z}_{0}\,d{\bf{x}}=0.$ (3.47)

Hence (3.46) reduces to

$\displaystyle\frac{\mu}{\rho}$
$\displaystyle\int_{\Omega}K^{-1}\mathbf{z_{0}}\cdot\mathbf{z_{0}}\,d{\bf{x}}+\frac{\mu}{\rho}\int_{\Omega}K^{-1}{\bf
v}_{r}\cdot\mathbf{z_{0}}\,d{\bf{x}}$ (3.48)
$\displaystyle+\frac{\beta}{\rho}\int_{\Omega}(|{\bf u}|{\bf u}-|{\bf
u}_{h}^{i+1}|{\bf u}_{h}^{i+1}|)({\bf u}-{\bf
u}_{h}^{i+1})\,d{\bf{x}}-\frac{\beta}{\rho}\int_{\Omega}(|{\bf u}|{\bf
u}-|{\bf u}_{h}^{i+1}|{\bf u}_{h}^{i+1})\cdot{\bf v}_{r}\,d{\bf{x}}$
$\displaystyle=\alpha\int_{\Omega}({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})\cdot\mathbf{z_{0}}\,d{\bf{x}}-\frac{\beta}{\rho}\int_{\Omega}\big{(}(|{\bf
u}_{h}^{i+1}|-|{\bf u}_{h}^{i}|){\bf
u}_{h}^{i+1}\big{)}\cdot\mathbf{z_{0}}\,d{\bf{x}}$
$\displaystyle+\sum_{\kappa\in\mathcal{T}_{h}}\Big{[}\int_{\kappa}\big{(}-\nabla
p_{h}^{i+1}-\alpha({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})-\frac{\mu}{\rho}K^{-1}{\bf u}_{h}^{i+1}-\frac{\beta}{\rho}|{\bf
u}_{h}^{i}|{\bf u}_{h}^{i+1}+{\bf
f}_{h}\big{)}\cdot\mathbf{z_{0}}\,d{\bf{x}}+\int_{\kappa}\big{(}{\bf f}-{\bf
f}_{h}\big{)}\cdot\mathbf{z_{0}}\,d{\bf{x}}\Big{]}.$

We decompose the fourth term of the left hand side as following

$\displaystyle\frac{\beta}{\rho}\int_{\Omega}(|{\bf u}|{\bf u}-|{\bf
u}^{i+1}_{h}|{\bf u}^{i+1}_{h})\cdot{\bf
v}_{r}\,d{\bf{x}}=\displaystyle\frac{\beta}{\rho}\int_{\Omega}|{\bf u}|({\bf
u}-{\bf u}^{i+1}_{h})\cdot{\bf
v}_{r}\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}(|{\bf u}|-|{\bf
u}^{i+1}_{h}|){\bf u}^{i+1}_{h}\cdot{\bf v}_{r}\,d{\bf{x}}.$

Equation (3.48) gives by inserting $\pm{\bf u}$ in the second term of the
right hand side and by using Property 2.1, the following inequality

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{\mu}{\rho}\int_{\Omega}K^{-1}\textbf{z}_{0}\cdot\textbf{z}_{0}\,d{\bf{x}}+c_{m}\left\|{\bf
u}-{\bf u}_{h}^{i+1}\right\|_{L^{3}(\Omega)^{d}}^{3}\leq\\\ \vskip 6.0pt plus
2.0pt minus
2.0pt\displaystyle\frac{\mu}{\rho}\int_{\Omega}|K^{-1}\textbf{v}_{r}|\,|\textbf{z}_{0}|\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}|{\bf
u}|\,|{\bf u}-{\bf u}^{i+1}_{h}|\,|{\bf
v}_{r}|\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}|{\bf u}-{\bf
u}^{i+1}_{h}|\,|{\bf u}^{i+1}_{h}|\,|{\bf v}_{r}|\,d{\bf{x}}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\displaystyle+\alpha\int_{\Omega}|{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}|\,|\textbf{z}_{0}|\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}|{\bf
u}^{i+1}_{h}-{\bf u}^{i}_{h}|\,|{\bf u}^{i+1}_{h}-{\bf
u}|\,|\textbf{z}_{0}|\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\Omega}|{\bf
u}^{i+1}_{h}-{\bf u}^{i}_{h}|\,|{\bf u}|\,|\textbf{z}_{0}|\,d{\bf{x}}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt+\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\big{[}\int_{\kappa}|-\nabla
p_{h}^{i+1}-\frac{\mu}{\rho}K^{-1}{\bf u}^{i+1}_{h}-\alpha({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})-\frac{\beta}{\rho}|{\bf u}^{i}_{h}|{\bf
u}^{i+1}_{h}+{\bf
f}_{h}|\,\,|\textbf{z}_{0}|\,d{\bf{x}}\,+\displaystyle\int_{\kappa}|{\bf
f}-{\bf f}_{h}|\,|\textbf{z}_{0}|\,d{\bf{x}}\,\big{]}.\end{array}$ (3.49)

By using the relation ${\bf u}-{\bf u}_{h}^{i+1}=\textbf{z}_{0}+{\bf v}_{r}$
and taking into account that $W^{1,6}(\Omega)\subset L^{\infty}(\Omega)$, the
last bound allows us to obtain the following inequality

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus 2.0pt\displaystyle\frac{\mu
K_{m}}{\rho}\|\textbf{z}_{0}\|^{2}_{L^{2}(\Omega)}+c_{m}\left\|{\bf u}-{\bf
u}_{h}^{i+1}\right\|_{L^{3}(\Omega)^{d}}^{3}\leq\displaystyle\frac{\mu
K_{M}}{\rho}\|\textbf{v}_{r}\|_{L^{2}(\Omega)}\,\|\textbf{z}_{0}\|_{L^{2}(\Omega)}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\displaystyle+\frac{\beta}{\rho}\|{\bf
u}\|_{L^{6}(\Omega)}(\|\textbf{z}_{0}\|_{L^{2}(\Omega)}+\|{\bf
v}_{r}\|_{L^{2}(\Omega)})\|{\bf v}_{r}\|_{L^{3}(\Omega)}\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt+\displaystyle\frac{\beta}{\rho}(\|{\bf
v}_{r}\|_{L^{2}(\Omega)}+\|\textbf{z}_{0}\|_{L^{2}(\Omega)})\|{\bf
u}^{i+1}_{h}\|_{L^{6}(\Omega)}\,\|{\bf
v}_{r}\|_{L^{3}(\Omega)}\displaystyle+\displaystyle\alpha\|{\bf
u}_{h}^{i+1}-{\bf
u}_{h}^{i}\|_{L^{2}(\Omega)}\,\|\textbf{z}_{0}\|_{L^{2}(\Omega)}\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\displaystyle+\frac{\beta}{\rho}\|{\bf
u}^{i+1}_{h}-{\bf u}^{i}_{h}\|_{L^{\infty}(\Omega)}\,(\|{\bf
v}_{r}\|_{L^{2}(\Omega)}+\|\textbf{z}_{0}\|_{L^{2}(\Omega)})\,\|\textbf{z}_{0}\|_{L^{2}(\Omega)}+\frac{\beta}{\rho}\|{\bf
u}^{i+1}_{h}-{\bf u}^{i}_{h}\|_{L^{2}(\Omega)}\,\|{\bf
u}\|_{L^{\infty}(\Omega)}\,\|\textbf{z}_{0}\|_{L^{2}(\Omega)}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt+\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\|-\nabla
p_{h}^{i+1}-\frac{\mu}{\rho}K^{-1}{\bf u}^{i+1}_{h}-\alpha({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})-\frac{\beta}{\rho}|{\bf u}^{i}_{h}|{\bf
u}^{i+1}_{h}+{\bf
f}_{h}\|_{L^{2}(\kappa)}\,\|\textbf{z}_{0}\|_{L^{2}(\kappa)}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt+\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\|{\bf
f}-{\bf
f}_{h}\|_{L^{2}(\kappa)}\,\|\textbf{z}_{0}\|_{L^{2}(\kappa)}.\end{array}$
(3.50)

Following Lemma 3.9, there exists an integer $i_{0}$ depending on $h$ such
that for all $i\geq i_{0}$, ${\bf u}_{h}^{i+1}$ is bounded in $L^{6}(\Omega)$.
Furthermore, we use following inverse inequality

$\|{\bf u}_{h}^{i+1}-{\bf u}_{h}^{i}\|_{L^{\infty}(\Omega)}\leq h^{-d/2}\|{\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i}\|_{L^{2}(\Omega)}$

and the convergence of the sequence ${\bf u}_{h}^{i}$ to ${\bf u}_{h}$ in
$L^{2}(\Omega)^{d}$ to deduce that we can choose $i_{0}$ sufficiently large
such that for all $i\geq i_{0}$ we have $\|{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}\|_{L^{\infty}(\Omega)}\leq\displaystyle\frac{\mu K_{m}}{2\beta}$.
Thus, Equation 3.50 gives:

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus 2.0pt\displaystyle\frac{\mu
K_{m}}{2\rho}\|\textbf{z}_{0}\|^{2}_{L^{2}(\Omega)}+c_{m}\left\|{\bf u}-{\bf
u}_{h}^{i+1}\right\|_{L^{3}(\Omega)^{d}}^{3}\leq\displaystyle\frac{\mu
K_{M}}{\rho}\|\textbf{v}_{r}\|_{L^{2}(\Omega)}\,\|\textbf{z}_{0}\|_{L^{2}(\Omega)}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\displaystyle+\frac{\beta}{\rho}\|{\bf
u}\|_{L^{6}(\Omega)}(\|\textbf{z}_{0}\|_{L^{2}(\Omega)}+\|{\bf
v}_{r}\|_{L^{2}(\Omega)})\|{\bf v}_{r}\|_{L^{3}(\Omega)}\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt+\displaystyle\frac{\beta}{\rho}(\|{\bf
v}_{r}\|_{L^{2}(\Omega)}+\|\textbf{z}_{0}\|_{L^{2}(\Omega)})\|{\bf
u}^{i+1}_{h}\|_{L^{6}(\Omega)}\,\|{\bf
v}_{r}\|_{L^{3}(\Omega)}\displaystyle+\displaystyle\alpha\|{\bf
u}_{h}^{i+1}-{\bf
u}_{h}^{i}\|_{L^{2}(\Omega)}\,\|\textbf{z}_{0}\|_{L^{2}(\Omega)}\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\displaystyle+\frac{\mu K_{m}}{2\rho}\,\|{\bf
v}_{r}\|_{L^{2}(\Omega)}\,\|\textbf{z}_{0}\|_{L^{2}(\Omega)}+\frac{\beta}{\rho}\|{\bf
u}^{i+1}_{h}-{\bf u}^{i}_{h}\|_{L^{2}(\Omega)}\,\|{\bf
u}\|_{L^{\infty}(\Omega)}\,\|\textbf{z}_{0}\|_{L^{2}(\Omega)}+\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\|{\bf
f}-{\bf f}_{h}\|_{L^{2}(\kappa)}\,\|\textbf{z}_{0}\|_{L^{2}(\kappa)}\\\ \vskip
6.0pt plus 2.0pt minus
2.0pt+\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\|-\nabla
p_{h}^{i+1}-{\frac{\mu}{\rho}}K^{-1}{\bf u}^{i+1}_{h}-\alpha({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})-\frac{\beta}{\rho}|{\bf u}^{i}_{h}|{\bf
u}^{i+1}_{h}+{\bf
f}_{h}\|_{L^{2}(\kappa)}\,\|\textbf{z}_{0}\|_{L^{2}(\kappa)}.\end{array}$
(3.51)

We use Lemma 3.9 and the decomposition
$ab\leq\displaystyle\frac{1}{2\varepsilon}a^{2}+\frac{1}{2}\varepsilon b^{2}$
for all the terms containing $b=\|\textbf{z}_{0}\|_{L^{2}(\Omega)}$ in the
right hand side of Equation (3.51), with $\varepsilon$ sufficiently small such
that all the terms of $\|\textbf{z}_{0}\|_{L^{2}(\Omega)}$ in the right hand
side will be absorbed by the term $\displaystyle\frac{\mu
K_{m}}{2\rho}\|\textbf{z}_{0}\|^{2}_{L^{2}(\Omega)}$ of the left hand side of
(3.51). We then get, by using the inequality $\|{\bf
v}_{r}\|_{L^{2}(\Omega)}\leq|\Omega|^{1/6}\|{\bf v}_{r}\|_{L^{3}(\Omega)}$ and
by taking the square root of the inequality, the bound

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\|\textbf{z}_{0}\|_{L^{2}(\Omega)}+\left\|{\bf u}-{\bf
u}_{h}^{i+1}\right\|_{L^{3}(\Omega)^{d}}^{3/2}\leq\displaystyle\overline{C}\Big{(}\|{\bf
v}_{r}\|_{L^{3}(\Omega)}+\sum_{\kappa\in\mathcal{T}_{h}}\|{\bf f}-{\bf
f}_{h}\|_{L^{2}(\kappa)}+\sum_{\kappa\in\mathcal{T}_{h}}\|{\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i}\|_{L^{2}(\kappa)}\\\
\qquad\qquad\qquad\qquad+\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\|-\nabla
p_{h}^{i+1}-{\frac{\mu}{\rho}}K^{-1}{\bf u}^{i+1}_{h}-\alpha({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})-\frac{\beta}{\rho}|{\bf u}^{i}_{h}|{\bf
u}^{i+1}_{h}+{\bf f}_{h}\|_{L^{2}(\kappa)}\Big{)},\end{array}$ (3.52)

where $\overline{C}$ is a constant depending on $({\bf u},p)$. By using
Relations (3.42) and (3.52), and the following inequality:

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus 2.0pt\|{\bf u}-{\bf
u}_{h}^{i+1}\|_{L^{2}(\Omega)}&\leq&\|\textbf{z}_{0}\|_{L^{2}(\Omega)}+\|{\bf
v}_{r}\|_{L^{2}(\Omega)}\\\
&\leq&\|\textbf{z}_{0}\|_{L^{2}(\Omega)}+|\Omega|^{1/6}\|{\bf
v}_{r}\|_{L^{3}(\Omega)}\end{array}$

we get the desired result. ∎

######  Remark 3.12.

Theorem 3.11 gives an upper bound for the error ${\bf u}-{\bf
u}_{h}^{i+1}=\textbf{z}_{0}+{\bf v}_{r}$ in $L^{2}(\Omega)^{d}$ and an upper
bound of ${\bf v}_{r}$ in $X$ which is a part of ${\bf u}-{\bf u}_{h}^{i+1}$.
Furthermore, it gives an upper bound of $\|{\bf u}-{\bf
u}_{h}^{i+1}\|_{L^{3}(\Omega)}$ but unfortunately, with the indicators to the
power of $2/3$.

In the next theorem, we will bound the error between the gradient of the exact
and numerical pressures with respect to the indicators in
$L^{3/2}(\Omega)^{d}$.

######  Theorem 3.13.

Under the assumptions of Theorem 3.11 and we assume that ${\bf f}\in
L^{2}(\Omega)^{d}$, then there exists a positive real number $i_{1}$ depending
on $h$ such that for all $i\geq i_{1}$ we have the following bound between the
solutions $({\bf u},p)$ of (2.9) and $({\bf u}_{h}^{i+1},p_{h}^{i+1})$ of
$\eqref{V1hi}$:

$\begin{array}[]{rcl}||\nabla(p-p_{h}^{i+1})||_{L^{3/2}(\Omega)}&\leq&C\displaystyle\Big{[}\sum_{\kappa\in\mathcal{T}_{h}}\Big{(}\eta_{\kappa,i}^{(D_{1})}+\eta_{\kappa,i}^{(D_{2})}+\eta_{\kappa,i}^{(L)}\Big{)}\\\
&&+\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\Big{(}\|\textbf{f}-\textbf{f}_{h}\|_{L^{2}(\kappa)}+h_{k}\|b-b_{h}\|_{L^{3}(\kappa)}+\sum_{e\in\Gamma_{h}^{b}}h_{e}^{\frac{1}{3}}||g_{h}-g||_{L^{3}(e)}\Big{)}\Big{]},\end{array}$
(3.53)

where $C$ is a constant depending on $({\bf u},p)$.

###### Proof.

Let $({\bf u},p)$ and $({\bf u}_{h}^{i+1},p_{h}^{i+1})$ the solutions of (2.9)
and (3.9). We test Equation (3.34) with $\textbf{v}_{h}={\bf 0}$ to get,

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\int_{\Omega}\nabla(p-p_{h}^{i+1})\cdot\textbf{v}\,d{\bf{x}}=\displaystyle-\frac{\mu}{\rho}\int_{\Omega}K^{-1}({\bf
u}-{\bf u}^{i+1}_{h})\cdot{\bf
v}\,d{\bf{x}}-\frac{\beta}{\rho}\int_{\Omega}(|{\bf u}|{\bf u}-|{\bf
u}^{i}_{h}|{\bf u}^{i+1}_{h})\cdot{\bf v}\,d{\bf{x}}\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt\hskip
28.45274pt+\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\big{[}\int_{\kappa}(-\nabla
p_{h}^{i+1}-\alpha({\bf u}_{h}^{i}-{\bf
u}_{h}^{i+1})-\frac{\mu}{\rho}K^{-1}{\bf u}^{i+1}_{h}-\frac{\beta}{\rho}|{\bf
u}^{i}_{h}|{\bf u}^{i+1}_{h}+{\bf f}_{h})\cdot{\bf v}\,d{\bf{x}}\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\hskip 142.26378pt+\displaystyle\int_{\kappa}({\bf
f}-{\bf f}_{h})\cdot{\bf v}\,d{\bf{x}}\,+\alpha\displaystyle\int_{\kappa}({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})\cdot{\bf v}\,d{\bf{x}}\big{]}.\end{array}$
(3.54)

By using the Cauchy-Schwartz inequality, we get after dividing the previous
inequality by $||{\bf v}||_{L^{3}(\Omega)}$:

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\frac{\displaystyle\Big{|}\int_{\Omega}\nabla(p_{h}^{i+1}-p){\bf
v}\,d{\bf{x}}\Big{|}}{||{\bf v}||_{L^{3}(\Omega)}}\leq C\displaystyle(||{\bf
u}-{\bf u}_{h}^{i+1}||_{L^{2}(\Omega)}+||{\bf u}_{h}^{i}-{\bf
u}_{h}^{i+1}||_{L^{2}(\Omega)})\frac{||{\bf v}||_{L^{2}(\Omega)}}{||{\bf
v}||_{L^{3}(\Omega)}}\\\
\displaystyle+C_{1}\Big{(}\sum_{\kappa\in\mathcal{T}_{h}}||{\bf f}-{\bf
f}_{h}||^{2}_{L^{2}(\kappa)}\Big{)}^{1/2}\frac{||{\bf
v}||_{L^{2}(\Omega)}}{||{\bf
v}||_{L^{3}(\Omega)}}+\frac{\beta}{\rho}\Big{|}\int_{\Omega}(|{\bf u}|{\bf
u}-|{\bf u}^{i}_{h}|{\bf u}^{i+1}_{h})\cdot{\bf
v}\,d{\bf{x}}\Big{|}\frac{1}{||{\bf v}||_{L^{3}(\Omega)}}\\\ +\displaystyle
C_{2}\Big{(}\sum_{\kappa\in\mathcal{T}_{h}}||-\nabla p_{h}^{i+1}-\alpha({\bf
u}_{h}^{i}-{\bf u}_{h}^{i+1})-\frac{\mu}{\rho}K^{-1}{\bf
u}^{i+1}_{h}-\frac{\beta}{\rho}|{\bf u}^{i}_{h}|{\bf u}^{i+1}_{h}+{\bf
f}_{h}||^{2}_{L^{2}(\kappa)}\Big{)}^{1/2}\frac{||{\bf
v}||_{L^{2}(\Omega)}}{||{\bf v}||_{L^{3}(\Omega)}}.\end{array}$ (3.55)

By using the relation $||{\bf v}||_{L^{2}(\Omega)}\leq|\Omega|^{1/6}||{\bf
v}||_{L^{3}(\Omega)}$, all the term of the right hand side of the previous
bound can be treated as in the previous theorem except the third one which can
be bounded as following:

$\begin{array}[]{rcl}\Big{|}(|{\bf u}|{\bf u}-|{\bf u}^{i}_{h}|{\bf
u}^{i+1}_{h},{\bf v})\Big{|}&\leq&\Big{|}((|{\bf u}|-|{\bf u}_{h}^{i}|){\bf
u},{\bf v})\Big{|}+\Big{|}(|{\bf u}_{h}^{i}|({\bf u}-{\bf u}_{h}^{i+1}),{\bf
v})\Big{|}\\\ &\leq&\big{(}||{\bf u}-{\bf u}^{i}_{h}||_{L^{2}(\Omega)}||{\bf
u}||_{L^{6}(\Omega)}+||{\bf u}^{i}_{h}||_{L^{6}(\Omega}||{\bf u}-{\bf
u}_{h}^{i+1}||_{L^{2}(\Omega)}\big{)}||{\bf v}||_{L^{3}(\Omega)}\end{array}$
(3.56)

We consider Relation (3.55). We use the following triangle inequality

$||{\bf u}-{\bf u}^{i}_{h}||_{L^{2}(\Omega)}\leq||{\bf u}-{\bf
u}^{i+1}_{h}||_{L^{2}(\Omega)}+||{\bf u}^{i+1}_{h}-{\bf
u}^{i}_{h}||_{L^{2}(\Omega)},$

the fact that the term $||{\bf u}^{i}_{h}||_{L^{6}(\Omega}$ is bounded, the
inf-sup condition (2.5) and Theorem 3.11, to get the desired error bound on
the pressure given by Equation (3.53). ∎

######  Remark 3.14.

The bounds (3.45) and (3.53) constitute our a posteriori error estimates where
we bound the error between the exact solution $({\bf u},p)$ of (2.9) and the
numerical solution $({\bf u}_{h}^{i+1},p_{h}^{i+1})$ of (3.9) with respect to
the indicators $\eta_{\kappa,i}^{(L)}$, $\eta_{\kappa,i}^{(D_{1})}$ and
$\eta_{\kappa,i}^{(D_{2})}$. But to get the bounds of the indicators which are
the subject of the next subsection (Section 3.2.2), we need to add the
following theorem where we add a supplementary bound giving an error bound of
the exact and numerical solutions.

######  Theorem 3.15.

Under the assumptions of Lemma 3.9, there exists an integer $i_{0}$ depending
on $h$ such that for all $i\geq i_{0}$, the solutions $({\bf u},p)$ of (2.9)
and $({\bf u}_{h}^{i+1},p_{h}^{i+1})$ of $\eqref{V1hi}$ verify the following
error inequalities:

$\begin{array}[]{ll}\|\displaystyle\frac{\beta}{\rho}(|{\bf u}|{\bf u}-|{\bf
u}_{h}^{i}|{\bf u}_{h}^{i+1})+\nabla(p-p_{h}^{i+1})\|_{L^{2}(\Omega)}\leq
C\displaystyle\Big{[}\sum_{\kappa\in\mathcal{T}_{h}}\Big{(}\eta_{\kappa,i}^{(D_{1})}+\eta_{\kappa,i}^{(D_{2})}+\eta_{\kappa,i}^{(L)}\Big{)}\\\
\hskip
85.35826pt\displaystyle+\sum_{\kappa\in\mathcal{T}_{h}}\Big{(}\|\textbf{f}-\textbf{f}_{h}\|_{L^{2}(\kappa)}+h_{k}\|b-b_{h}\|_{L^{3}(\kappa)}+\sum_{e\in\Gamma_{h}^{b}}h_{e}^{\frac{1}{3}}||g_{h}-g||_{L^{3}(e)}\Big{)}\Big{]},\end{array}$
(3.57)

where $C$ is a constant depending on the exact solution $({\bf u},p)$ of
(2.9).

###### Proof.

Let ${\bf u}\in W^{1,6}(\Omega)^{d}$. Then, Equation (1.1) allows us to get
that the pressure is such that $\nabla p\in L^{2}(\Omega)^{d}$. Thus, the
velocity error equation (3.34) is valid for all ${\bf v}\in L^{2}(\Omega)^{d}$
and can be written as

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\int_{\Omega}\big{(}\frac{\beta}{\rho}(|{\bf u}|{\bf
u}-|{\bf u}^{i}_{h}|{\bf u}^{i+1}_{h})+\nabla(p-p_{h}^{i+1})\big{)}\cdot{\bf
v}\,d{\bf{x}}=\displaystyle-\frac{\mu}{\rho}\int_{\Omega}K^{-1}({\bf u}-{\bf
u}^{i+1}_{h})\cdot{\bf v}\,d{\bf{x}}+\alpha\int_{\Omega}({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})\cdot{\bf v}\,d{\bf{x}}\\\ \vskip 6.0pt plus
2.0pt minus
2.0pt+\displaystyle\sum_{\kappa\in\mathcal{T}_{h}}\big{[}\int_{\kappa}(-\nabla
p_{h}^{i+1}-{{\frac{\mu}{\rho}K^{-1}{\bf u}_{h}^{i+1}}}-\alpha({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})-\frac{\beta}{\rho}|{\bf u}_{h}^{i}|{\bf
u}_{h}^{i+1}+{\bf f}_{h})\cdot({\bf v}-{\bf
v}_{h})\,d{\bf{x}}+\int_{\kappa}({\bf f}-{\bf f}_{h})\cdot{\bf
v}\,d{\bf{x}}\big{]}.\end{array}$ (3.58)

By taking ${\bf v}_{h}={\bf 0}$ and ${\bf
v}=\displaystyle\big{(}\frac{\beta}{\rho}(|{\bf u}|{\bf u}-|{\bf
u}^{i}_{h}|{\bf u}^{i+1}_{h})+\nabla(p-p_{h}^{i+1})\big{)}$, applying the
Cauchy-Schwartz and simplifying by $||{\bf v}||_{L^{2}(\Omega)}$, we get the
result after using Theorem 3.11. ∎

######  Remark 3.16.

Finally, the bounds (3.45), (3.53) and (3.57) constitute our a posteriori
error estimates.

####  3.2.2. Bounds of the indicators

In order to establish the efficiency of the a posteriori error estimates, we
recall the following properties (see R. Verfürth,[38], Chapter 1). For an
element $\kappa$ of $\mathcal{T}_{h}$, we consider the bubble function
$\psi_{\kappa}$ (resp. $\psi_{e}$ for the face $e$) which is equal to the
product of the $d+1$ barycentric coordinates associated with the vertices of
$\kappa$ (resp. of the $d$ barycentric coordinates associated with the
vertices of $e$). We also consider a lifting operator ${\mathcal{L}}_{e}$
defined on polynomials on $e$ vanishing on $\partial e$ into polynomials on
the at most two elements $\kappa$ containing $e$ and vanishing on
$\partial\kappa\setminus e$, which is constructed by affine transformation
from a fixed operator on the reference element.

######  Property 3.17.

Denoting by $Pr(\kappa)$ the space of polynomials of degree smaller than $r$
on $\kappa$. The following properties hold:

$\forall v\in
P_{r}(\kappa),\qquad\begin{cases}c||v||_{0,\kappa}\leq||v\psi^{1/2}_{\kappa}||_{0,\kappa}\leq
c^{\prime}||v||_{0,\kappa},&\\\ |v|_{1,\kappa}\leq
ch_{\kappa}^{-1}||v||_{0,\kappa}.&\end{cases}$ (3.59)

######  Property 3.18.

Denoting by $Pr(e)$ the space of polynomials of degree smaller than $r$ on
$e$, we have

$\forall\;v\in P_{r}(e),\qquad c\|v\|_{0,e}\leq\|v\psi_{e}^{1/2}\|_{0,e}\leq
c^{\prime}\|v\|_{0,e},$

and, for all polynomials $v$ in $Pr(e)$ vanishing on $\partial e$, if $\kappa$
is an element which contains $e$,

$\|{\mathcal{L}}_{e}v\|_{0,\kappa}+h_{e}\mid{\mathcal{L}}_{e}v\mid_{1,\kappa}\leq
ch^{1/2}_{e}\|v\|_{0,e}.$

We have the following bounds of the indicators:

######  Theorem 3.19.

Let $d=3$, $({\bf u},p)$ and $({\bf u}_{h}^{i+1},p_{h}^{i+1})$ the solutions
of (2.9) and of $\eqref{V1hi}$. We have the following bounds of the
indicators: for each element $\kappa\in\mathcal{T}_{h}$,

$\eta^{(L)}_{\kappa,i}\leq||{\bf u}-{\bf u}_{h}^{i}||_{L^{2}(\kappa)}+||{\bf
u}-{\bf u}_{h}^{i+1}||_{L^{2}(\kappa)},$ (3.60)

and

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus
2.0pt\eta^{(D_{2})}_{\kappa,i}&\leq&C\displaystyle\Big{(}\|{\bf
v}_{r}\|_{L^{3}(w_{\kappa})}+h_{k}\|b-b_{h}\|_{L^{3}(w_{\kappa})}+\sum_{e\in\partial\kappa}h_{e}^{\frac{1}{3}}||g_{h}-g||_{L^{3}(e)}\Big{)},\end{array}$
(3.61)

where $C$ is a constant independent of the mesh step but depends on the exact
solution $({\bf u},p)$.

###### Proof.

The bound (3.60) is a simple consequence of the definition of
$\eta_{\kappa,i}^{(L)}$ and the triangle inequality.
In order to prove (3.61), we consider first Equation (3.35) with $q_{h}=0$ and

$q=q_{\kappa}=\left\\{\begin{array}[]{lcl}b_{h}\psi_{\kappa}&\mbox{on
}\kappa,\\\ 0&\mbox{on }\Omega\backslash\kappa,\end{array}\right.$

where $\psi_{\kappa}$ is the bubble functions on a given element
$\kappa\in\mathcal{T}_{h}$. We obtain by using Relation (3.47) the following
equation:

$\displaystyle\int_{\kappa}b_{h}^{2}\psi_{\kappa}\,d{\bf{x}}=-\displaystyle\int_{\kappa}\nabla(b_{h}\psi_{\kappa})\cdot{\bf
v}_{r}\,d{\bf{x}}-\int_{\kappa}(b_{h}\psi_{\kappa})(b-b_{h})\,d{\bf{x}}.\,$
(3.62)

Then we use Property 3.17, the Cauchy-Schwartz inequality and the relation
$||{\bf v}||_{L^{2}(\kappa)}\leq|\kappa|^{1/6}||{\bf v}||_{L^{3}(\kappa)}$,
and by multiplying by $h_{\kappa}$ to get:

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus
2.0pth_{\kappa}||b_{h}||_{L^{2}(\kappa)}&\leq&ch_{\kappa}(h_{\kappa}^{-1}||{\bf
v}_{r}||_{L^{2}(\kappa)}+||b-b_{h}||_{L^{2}(\kappa)})\\\
&\leq&c_{1}(h_{\kappa}^{\frac{1}{2}}||{\bf
v}_{r}||_{L^{3}(\kappa)}+h^{\frac{3}{2}}_{\kappa}||b-b_{h}||_{L^{3}(\kappa)}).\end{array}$
(3.63)

Then we get by using the inverse inequality (3.2) with $p=3$,

$h_{\kappa}||b_{h}||_{L^{3}(\kappa)}\leq c_{2}(||{\bf
v}_{r}||_{L^{3}(\kappa)}+h_{\kappa}||b-b_{h}||_{L^{3}(\kappa)}),$ (3.64)

which is the part of the indicator $\eta_{\kappa,i}^{(D_{2})}$ corresponding
to $b_{h}$.
Again, we consider Equation (3.35) with $q_{h}=0$ and

$q=q_{e}=\left\\{\begin{array}[]{lcl}\mathcal{L}_{e,\kappa}\big{(}\phi_{h,1}^{e}\psi_{e}\big{)}&\mbox{on
}\\{\kappa,\kappa^{\prime}\\},\\\ 0&\mbox{ on
}\Omega\backslash(\kappa\cup\kappa^{\prime}),\\\ \end{array}\right.$

where $\psi_{e}$ is the bubble function of $e$ and $\kappa^{\prime}$ denotes
the other element of $\mathcal{T}_{h}$ that share $e$ with $\kappa$. We get
the following equation:

$\displaystyle\int_{e}\phi_{h,1}^{e}q\,ds=\displaystyle-\int_{\kappa\cup\kappa^{\prime}}q(b-b_{h})\,d{\bf{x}}\,-\int_{\kappa\cup\kappa^{\prime}}b_{h}q\,d{\bf{x}}-\int_{e}(g_{h}-g)q\,ds-\int_{\kappa\cup\kappa^{\prime}}\nabla
q\cdot{\bf v}_{r}\,d{\bf{x}}.$

Properties 3.17 and 3.18 allow us to get the following bound:

$||\phi_{h,1}^{e}||_{L^{2}(e)}\leq
C\big{(}h_{e}^{\frac{1}{2}}||b-b_{h}||_{L^{2}(\kappa\cup\kappa^{\prime})}+h_{e}^{\frac{1}{2}}||b_{h}||_{L^{2}(\kappa\cup\kappa^{\prime})}+||g-g_{h}||_{L^{2}(e)}+h_{e}^{-\frac{1}{2}}||{\bf
v}_{r}||_{L^{2}(\kappa\cup\kappa^{\prime})}\big{)}.$

By using again the inverse inequality 3.2 and the relation $\|{\bf
v}_{h}\|_{L^{2}(\kappa)}\leq|\kappa|^{1/6}\|{\bf v}_{h}\|_{L^{3}(\kappa)}$, we
obtain the bound:

$h_{e}^{\frac{1}{3}}||\phi_{h,1}^{e}||_{L^{3}(e)}\leq
C_{1}\big{(}h_{e}||b-b_{h}||_{L^{3}(\kappa\cup\kappa^{\prime})}+h_{e}||b_{h}||_{L^{3}(\kappa\cup\kappa^{\prime})}+h_{e}^{\frac{1}{3}}||g-g_{h}||_{L^{3}(e)}+||{\bf
v}_{r}||_{L^{3}(\kappa\cup\kappa^{\prime})}\big{)}.$ (3.65)

Hence, we bound the part of $\eta_{\kappa,i}^{(D_{2})}$ corresponding to
$\phi_{h,1}^{e}$. Relations (3.64) and (3.65) give Relation (3.61). ∎

######  Theorem 3.20.

Let $d=3$ and let the mesh satisfy (3.1). Under the assumptions of Lemma 3.9,
we have the following bound: for each element $\kappa\in\mathcal{T}_{h}$,

$\begin{array}[]{rcl}\vskip 6.0pt plus 2.0pt minus
2.0pt\eta^{(D_{1})}_{\kappa,i}&\leq&C\displaystyle\Big{(}\eta_{\kappa,i}^{(L)}+\|{\bf
u}-{\bf u}_{h}^{i+1}\|_{L^{2}(w_{\kappa})}+\|\frac{\beta}{\rho}\big{(}|{\bf
u}|{\bf u}-|{\bf u}_{h}^{i+1}|{\bf
u}_{h}^{i+1}\big{)}+\nabla(p-p_{h}^{i+1})\|_{L^{2}(\kappa)}\\\ &&\hskip
28.45274pt\displaystyle+\|K^{-1}-K_{h}^{-1}\|_{L^{3}(w_{\kappa})}+\|\textbf{f}-\textbf{f}_{h}\|_{L^{2}(w_{\kappa})}+h_{k}\|b-b_{h}\|_{L^{3}(w_{\kappa})}+\sum_{e\in\partial\kappa}h_{e}^{\frac{1}{3}}||g_{h}-g||_{L^{3}(e)}\Big{)},\end{array}$
(3.66)

where $C$ is a constant independent of the mesh step but depends on the exact
solution $({\bf u},p)$ and $K_{h}^{-1}$ is an approximation of $K^{-1}$ which
is a constant tensor in each triangle.

###### Proof.

Let us now prove Relation (3.66). We consider Equation (3.34) with ${\bf
v}_{h}=0$ and

${\bf v}={\bf
v}_{\kappa}=\left\\{\begin{array}[]{lcl}\displaystyle\big{(}-\nabla
p_{h}^{i+1}-\alpha({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})-\frac{\mu}{\rho}K_{h}^{-1}{\bf
u}^{i+1}_{h}-\frac{\beta}{\rho}|{\bf u}^{i}_{h}|{\bf u}^{i+1}_{h}+{\bf
f}_{h}\big{)}\psi_{\kappa}&&\mbox{ on }\kappa,\\\ 0&&\mbox{ on
}\Omega\backslash\kappa,\\\ \end{array}\right.$

where $K_{h}^{-1}$ is an approximation of $K^{-1}$ which is a constant tensor
in each triangle. We obtain the following equation:

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus
2.0pt\displaystyle\int_{\kappa}|(-\nabla p_{h}^{i+1}-\alpha({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})-\frac{\mu}{\rho}K_{h}^{-1}{\bf
u}^{i+1}_{h}-\frac{\beta}{\rho}|{\bf u}^{i}_{h}|{\bf u}^{i+1}_{h}+{\bf
f}_{h})\psi_{\kappa}^{1/2}|^{2}d{\bf{x}}=\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\hskip
56.9055pt\displaystyle\frac{\mu}{\rho}\int_{\kappa}(K^{-1}-K_{h}^{-1}){\bf
u}^{i+1}_{h}\cdot{\bf v}\,d{\bf{x}}+\frac{\mu}{\rho}\int_{\kappa}K^{-1}({\bf
u}-{\bf u}^{i+1}_{h})\cdot{\bf
v}\,d{\bf{x}}+\frac{\beta}{\rho}\int_{\kappa}(|{\bf u}|{\bf u}-|{\bf
u}^{i}_{h}|{\bf u}^{i+1}_{h})\cdot{\bf v}\,d{\bf{x}}\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt\hskip
56.9055pt\displaystyle+\displaystyle\int_{\kappa}\nabla(p-p_{h}^{i+1})\cdot\textbf{v}\,d{\bf{x}}-\displaystyle\int_{\kappa}({\bf
f}-{\bf f}_{h})\cdot{\bf v}\,d{\bf{x}}\,-\alpha\displaystyle\int_{K}({\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i})\cdot{\bf v}\,d{\bf{x}},\end{array}$

and then by using Lemma 3.9 and Properties 3.17 and 3.18, we get the following
bound:

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus 2.0pt\displaystyle\|-\nabla
p_{h}^{i+1}-\alpha({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})-\frac{\mu}{\rho}K^{-1}_{h}{\bf
u}^{i+1}_{h}-\frac{\beta}{\rho}|{\bf u}^{i}_{h}|{\bf u}^{i+1}_{h}+{\bf
f}_{h}\|_{L^{2}(\kappa)}\leq\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\hskip
56.9055pt\displaystyle C\Big{(}\|K^{-1}-K_{h}^{-1}\|_{L^{3}(\kappa)}+\|{\bf
u}-{\bf u}^{i+1}_{h}\|_{L^{2}(\kappa)}+\|\frac{\beta}{\rho}\big{(}|{\bf
u}|{\bf u}-|{\bf u}^{i}_{h}|{\bf
u}^{i+1}_{h}\big{)}+\nabla(p-p_{h}^{i+1})\|_{L^{2}(\kappa)}\\\ \hskip
85.35826pt\displaystyle+\displaystyle\|{\bf f}-{\bf
f}_{h}\|_{L^{2}(\kappa)}+\alpha\|{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}\|_{L^{2}(\kappa)}\Big{)}.\end{array}$ (3.67)

Thus we get the result by using the following triangle inequality and Lemma
3.9 :

$\begin{array}[]{ll}\vskip 6.0pt plus 2.0pt minus 2.0pt\displaystyle\|-\nabla
p_{h}^{i+1}-\alpha({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})-\frac{\mu}{\rho}K^{-1}{\bf u}^{i+1}_{h}-\frac{\beta}{\rho}|{\bf
u}^{i}_{h}|{\bf u}^{i+1}_{h}+{\bf f}_{h}\|_{L^{2}(\kappa)}\leq\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\hskip 14.22636pt\displaystyle\|-\nabla
p_{h}^{i+1}-\alpha({\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i})-\frac{\mu}{\rho}K^{-1}_{h}{\bf
u}^{i+1}_{h}-\frac{\beta}{\rho}|{\bf u}^{i}_{h}|{\bf u}^{i+1}_{h}+{\bf
f}_{h}\|_{L^{2}(\kappa)}+\frac{\mu}{\rho}\|K^{-1}-K^{-1}_{h}\|_{L^{3}(\kappa)}\|{\bf
u}^{i+1}_{h}\|_{L^{6}(\kappa)}.\end{array}$ (3.68)

∎

##  4\. Numerical simulation

We validate the theory developed here by showing numerical simulations using
Freefem++ (see [27]).
We consider the iterative scheme (3.9). For the stopping criterion given
later, we define the iterative error

$Err_{L}=\displaystyle\Big{(}\frac{||{\bf u}_{h}^{i+1}-{\bf
u}_{h}^{i}||_{L^{3}(\Omega)}+||\nabla(p_{h}^{i+1}-p_{h}^{i})||_{L^{3/2}(\Omega)}}{||{\bf
u}_{h}^{i+1}||_{L^{3}(\Omega)}+||\nabla
p_{h}^{i+1}||_{L^{3/2}(\Omega)}}\Big{)}.$

In the definition of $Err_{L}$, we consider the iterative error $||{\bf
u}_{h}^{i+1}-{\bf u}_{h}^{i}||_{L^{3}(\Omega)}$ in the natural space of the
velocity $L^{3}(\Omega)$ although in the definition of
$\eta_{\kappa,i}^{(L)}$, we used the error in $L^{2}(\kappa)$.

###  4.1. First test case

In this section, the domain ${\Omega}$ is the square $\Omega=]0,1[^{2}$ and
all computations start on a uniform initial triangular mesh obtained by
dividing $\Omega$ into $N^{2}$ equal squares, each one subdivided into $2$
triangles, so that the initial triangulation consists of $2N^{2}$ triangles.
The theory is tested by applying the numerical scheme (3.9) to the exact
solution $({\bf u},p,T)=({\bf curl}\,\psi,p,T)$ where $\psi$ and $p$ are given
by

$\displaystyle\psi(x,y)=e^{-\gamma((x-0.5)^{2}+(y-0.5)^{2})}$ (4.1)

and

$\displaystyle p(x,y)=x*(x-2./3.)*y*(y-2./3.),$ (4.2)

with the choice $K=I$, $\mu=\rho=1$ and $\gamma=50$. Here we have ${\bf
u}.{\bf n}=0$ and $b=\operatorname{div}({\bf u})=0$.

We begin by testing the dependency of the convergence of the iterative scheme
(3.9) with respect to $\alpha$. We consider $N=60$ and for each $\alpha$, we
stop the algorithm (3.9) when the error $Err_{L}<1e^{-5}$.
To discribe the convergence of Algorithm (3.9), we consider also the error

$Err=\displaystyle\Big{(}\frac{||{\bf u}_{h}^{i}-{\bf
u}||_{L^{3}(\Omega)}+||\nabla(p_{h}^{i}-p)||_{L^{3/2}(\Omega)}}{||{\bf
u}||_{L^{3}(\Omega)}+||\nabla p||_{L^{3/2}(\Omega)}}\Big{)}.$

We mention that in the definition of $Err$, we considered $\|{\bf
u}_{h}^{i}-{\bf u}\|_{L^{3}(\Omega)}$ despite that in Theorem 3.11 it figures
to the power $3/2$.

Tables 2 and 2 show, for $\beta=1$ and $\beta=10$, the error $Err$ and the
number of iterations $Nbr$ which describe the convergence of Algorithm (3.9)
with respect to $\alpha$ for ${\bf u}_{h}^{0}={\bf 0}$. We remark that the
best convergence is obtained for $\alpha_{min}=2.3$ when $\beta=1$ and for
$\alpha_{min}=12$ when $\beta=10$.

$\alpha$ | 0.001 | 0.01 | 0.1 | 1 | 1.4 | 1.9 | 2.1 | 2.3 | 2.5 | 2.7 | 3 | 3.7 | 5 | 10 | 100 | 1000
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Nbr | 75 | 74 | 63 | 25 | 20 | 16 | 15 | 14 | 15 | 16 | 17 | 20 | 24 | 41 | 231 | 1134

Table 1.  Number of iterations $Nbr$ for each $\alpha$. ($\beta=1$ and ${\bf
u}_{h}^{0}={\bf 0}$)). In all these cases, $Err=-0.939$ (in logarithmic
scale).

$\alpha$ | 0.001 | 0.01 | 0.1 | 1 | 6 | 8 | 10 | 11 | 12 | 13 | 14 | 15 | 18 | 21 | 35 | 100 | 1000
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Nbr | 693 | 681 | 579 | 232 | 55 | 42 | 34 | 32 | 30 | 31 | 32 | 34 | 38 | 42 | 59 | 116 | 465

Table 2. Number of iterations $Nbr$ for each $\alpha$. ($\beta=10$ and ${\bf
u}_{h}^{0}={\bf 0}$). In all these cases, $Err=-0.475$ (in logarithmic scale).

To go far with our numerical investigations, we test Algorithm (3.9) where the
initial guess ${\bf u}_{h}^{0}$ is calculated by using the Darcy’s problem
(which corresponds to $\beta=\alpha=0$). Tables 4 and 4 show, for $\beta=1$
and $\beta=10$, the error $Err$ and the number of iterations $Nbr$ with
respect to $\alpha$. We remark that here also the best convergence is obtained
for $\alpha_{min}=1.4$ when $\beta=1$ and for $\alpha_{min}=14$ when
$\beta=10$. In this case of the initial guess, we remark that the number of
the iterations is slightly smaller than that obtained for ${\bf
u}_{h}^{0}={\bf 0}$. Thus, in the following, all the numerical investigations
will be performed with the initial guess calculated by using the Darcy’s
problem.

$\alpha$ | 0.001 | 0.01 | 0.1 | 1 | 1.4 | 1.9 | 2.3 | 2.6 | 2.9 | 3.3 | 3.7 | 4 | 5 | 6 | 10 | 100 | 1000
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Nbr | 74 | 73 | 62 | 25 | 20 | 16 | 13 | 12 | 13 | 14 | 15 | 16 | 18 | 21 | 30 | 174 | 904

Table 3. Number of iterations $Nbr$ for each $\alpha$. ($\beta=1$). ${\bf
u}_{h}^{0}$ is calculated by using the Darcy’s problem. In all these cases,
$Err=-0.939$ (in logarithmic scale).

$\alpha$ | 0.001 | 0.01 | 0.1 | 1 | 6 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 18 | 21 | 28 | 100 | 1000
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Nbr | 692 | 680 | 578 | 232 | 55 | 34 | 31 | 29 | 27 | 26 | 27 | 28 | 30 | 33 | 40 | 92 | 458

Table 4. Number of iterations $Nbr$ for each $\alpha$. ($\beta=10$). ${\bf
u}_{h}^{0}$ is calculated by using the Darcy’s problem. In all these cases,
$Err=-0.475$ (in logarithmic scale).

Furthermore, Table 5 shows the dependancy of $\alpha_{min}$ whith respect to
$h$ for $\beta=100$. We remark that $\alpha_{min}$ increases when $h$
decreases, which is consistent with the results of the Theorems 3.2 and 3.4
(see Remarks 3.3 and 3.5).

$h$ | 0.1414 | 0.0708 | 0.0353 | 0.0283 | 0.0177 | 0.007889 | 0.00543
---|---|---|---|---|---|---|---
$\alpha_{min}$ | 66 | 68.5 | 80.5 | 82.25 | 82.75 | 85.5 | 87
Table 5. $\alpha_{min}$ with respect to $h$. ($\beta=100$ and ${\bf
u}_{h}^{0}=0$).

In the following of this section, we will show numerical investigations
corresponding to the a posteriori error estimate. We take $\beta=10$,
$\alpha=10$ and $N=10$ for the initial mesh.
On a given mesh and for the numerical calculation, it is convenient to compute
the following indicators:

$\eta_{i}^{(D)}=\big{(}\sum_{K\in\mathcal{T}_{h}}(\eta_{K,i}^{(D_{1})})^{2}+(\eta_{K,i}^{(D_{2})})^{2}\big{)}^{\frac{1}{2}}$

and

$\eta_{i}^{(L)}=\big{(}\sum_{K\in\mathcal{T}_{h}}((\eta_{K,i}^{(L)})^{2}\big{)}^{\frac{1}{2}}.$

The iterations are stopped following the criteria

$\eta_{i}^{(L)}\leq\tilde{\gamma}\eta_{i}^{(D)},$ (4.3)

where $\tilde{\gamma}=0.001$. For the study of the dependence of the stopping
criteria $\eqref{stoppingg}$ with $\tilde{\gamma}$, we refer to [21] and [23]
where the authors introduce this new stopping criterion.
For the adaptive mesh (refinement and coarsening), we use routines in
FreeFem++. The indicators (3.36) are used for mesh adaptation by the adapted
mesh algorithm introduced in [9].

In Figure 1, we present the evolution of the mesh during the iterations
(initial, second and fourth refinement levels). We notice that the mesh is
concentrated in the region where the solution needs to be well described.

Figure 1. Evolution of the mesh during the refinement levels (initial, second
and fourth).

Next, we plot and study the error curves between the exact and numerical
solutions corresponding to uniform and adaptive method.

Figure 2 plots a comparison of the global error curves $Err$ versus the total
number of vertices in logarithmic scales for the uniform and adapt methods;
global in the sense that they depict the sum of the velocity and pressure
errors. We notice that the errors of the adaptive mesh method are much smaller
than that obtained with the uniform method, hence the efficiency of this
method.

Figure 2. Comparison of the errors $Err$ with respect to the total number of
vertices in logarithmic scale.

In table 6, we present the effectivity index defined as

$EI=\displaystyle\frac{\displaystyle\eta_{i}^{(L)}+\eta_{i}^{(D)}}{\displaystyle\|\textbf{u}-\textbf{u}_{h}^{i+1}\|_{L^{3}(\Omega)}+||\nabla(p-p_{h}^{i+1|})||_{L^{3/2}(\Omega)}}$

with respect to the number of vertices during the refinement levels. This
effectivity index is calculated on each mesh level after the convergence on
the iterations i by using the stopping criteria (4.3). Table 6 shows that it
is between $38.51$ and $25.47$.

Refinement Level | initial | first | second | third | fourth | fifth | sixth
---|---|---|---|---|---|---|---
Number of vertices | 121 | 313 | 973 | 2638 | 6197 | 15358 | 37703
Effectivity index | 38.51 | 32.53 | 30.54 | 27.91 | 28.70 | 27.15 | 25.47
Table 6.  EI with respect to the refinement levels.

###  4.2. Second test case

In this case, we consider a more complicated geometry 3 presenting reentrant
corners to show the efficiency of the adaptive method proposed in this work.
Furthermore, we take $\mu=\rho=1,\alpha=10,\beta=10,{\bf u}.{\bf
n}=0,b=\operatorname{div}({\bf u})=0$, ${\bf f}=(f,0)$ where

$f=\left\\{\begin{array}[]{rcl}0&&\quad\mbox{if }y>1,\\\ -2&&\quad\mbox{if
}y<=1,\end{array}\right.$

and $K$ such that

$K^{-1}=\left(\begin{array}[]{lcl}2+\sin(\pi x)\sin(\pi y)&&\quad 0.2x\\\
0.2x&&\quad 3+\sin(\pi x)\sin(\pi y)\end{array}\right)$

for all the numerical simulations of this section.

Figure 3. Geometry.

We begin by showing comparisons between the uniform and the adaptive methods
corresponding to iterative system (3.9).
Figures 7-7 present the evolution of the mesh during the iterations. We remark
that, from an iteration to another, the concentration of the refinement is on
the complex vorticity regions, namely at the reentrant corner and some regions
of $\Omega$.

Figure 4. Initial mesh (912 triangles)

Figure 5.  Second level mesh (3956 triangles)

Figure 6.  Third level mesh (10389 triangles)

Figure 7.  Fourth level mesh (20361 triangles)

Figures 9 and 9 show color velocity and pressure at the fourth refinement
level. We can clearly see that the velocity in Figure 9 justifies the
concentration of the refinement showed in Figure 7.

Figure 8. numerical velocity at the Fourth refinement level.

Figure 9.  numerical pressure at the Fourth refinement level.

Next, we introduce the relative total error indicator given

$E_{tot}=\displaystyle\frac{\displaystyle\eta_{i}^{(D)}}{\displaystyle\|\textbf{u}_{h}^{i}\|_{L^{3}(\Omega)}+||\nabla(p_{h}^{i})||_{L^{3/2}(\Omega)}}$

where here also $\eta_{i}^{(D)}$ is computed after convergence on the
iterations $i$ (by using the stopping criteria (4.3)). Figure 10 shows and
compares the relative total error indicator given by $E_{tot}$ between the
uniform and adaptive methods.

Figure 10. Comparison of the errors $Err$ with respect to the total number of
vertices in logarithmic scale for the second case.

##  5\. Conclusion

In this article, we discretize a steady Darcy-Forchheimer problem. We
introduce error indicators and establish optimal a posteriori error estimates.
We perform several numerical simulations where the indicators are used for
mesh adaptation, and we show the efficiency of these adaptive methods.

## References

  * [1] Abboud H., El Chami F. & Sayah T., A priori and a posteriori estimates for three dimentional Stokes equations with non standard boundary conditions, Numer. Methods Partial Differential Equations, 28, pp. 1178-1193, (2012).
  * [2] Alonso A., Error estimators for a mixed method, Numerische Mathematik, 74(4), 385-395, (1996).
  * [3] Aziz K., Settari A., Petroleum Reservoir Simulation, Applied Science Publishers LTD, London, (1979).
  * [4] Ainsworth, M., Oden, J.T., A posteriori error estimation in finite element analysis.Computer methods in applied mechanics and engineering 142.1-2 (1997): 1-88.
  * [5] Babuska I. and Rheinboldt W.C., Error estimates for adaptive finite element computations, SIAM J. Numer. Anal., 15(4), 736-754, (1978).
  * [6] Babuska, I., The finite element method with Lagrangian multipliers.Numer. Math. 20, 179-192 (1973)

  * [7] Belhachmi, Z., Bernardi, C. and Deparis, Weighted Clément operator and application to the finite element discretization of the axisymmetric Stokes problem, Numer. Math., 105 (2), 105-217, (2006).
  * [8] Bao, G., Jiang, X., Li, P., and Yuan, X.  An adaptive finite element DtN method for the elastic wave scattering by biperiodic structures. ESAIM: Mathematical Modelling and Numerical Analysis, 55(6), 2921-2947. (2021)
  * [9] Bernardi C., Dakroub J., Mansour G., and Sayah T., A posteriori analysis of iterative algorithms for Navier-Stokes problem, ESAIM: Mathematical Modelling and Numerical Analysis, 50(4), 1035-1055, (2016).
  * [10] Bernardi C., Hecht F. & Verfürth R., Finite element discretization of the three-dimensional Navier-Stokes equations with mixed boundary conditions, Math. Model. and Numer. Anal., 3, 1185-1201, (2009).
  * [11] Braess D. and Verf$\ddot{u}$rth R., A posteriori error estimators for the Raviart-Thomas element, SIAM J. Numer. Anal., 33(6), 2431-2444, (1996).
  * [12] Bernardi C. & Sayah T., A posteriori error analysis of the time dependent Stokes equations with mixed boundary conditions, IMA J. Numer. Anal. (2014), doi: 10.1093/imanum/drt06.
  * [13] Bernardi C. & Sayah T., A posteriori error analysis of the time dependent Navier-Stokes equations with mixed boundary conditions, SeMA, 69, 1-23, (2015).
  * [14] Carstensen C., Aposteriori error estimate for the mixed finite element method, Mathematics of Computation, 66(218), 465-476, (1997).
  * [15] Chen W. and Wang Y., A posteriori estimate for the $H(div;\Omega)$ conforming mixed finite element for the coupled Darcy-Stokes system, Journal of Computational and Applied Mathematics, 255, 502-516, (2014).
  * [16] Ciarlet P. G., Basic error estimates for elliptic problems,In Handbook of Numerical Analysis, Vol. II, Handbook of Numerical Analysis, pages 17-351. North-Holland, Amsterdam, (1991).
  * [17] Girault V., Wheeler M.F., Numerical discretization of a Darcy-Forchheimer model, Numer. Math. 110(2), 161-198, (2008).
  * [18] Clément P., Approximation by finite element functions using local regularization, RAIRO Anal. Numér. 9, R2,77-84, (1975).
  * [19] Dakroub J., Faddoul J. & Sayah T., A posteriori analysis of the Newton method applied to the Navier-Stokes problem , Journal of Applied Mathematics and Computing, https://doi.org/10.1007/s12190-020-01323-w, (2020).
  * [20] Dib S., Girault V., Hecht F. and Sayah T., A posteriori error estimates for Darcy’s problem coupled with the heat equation, ESAIM Mathematical Modelling and Numerical Analysis, DOI:10.1051/m2an/2019049, (2019).
  * [21] El Alaoui L., Ern A.,Vohralík M., Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems, Computable Methods in Applied Mechanics and Engineering 200, 2782-2795, (2011).
  * [22] Ern A. and Guermond J.L., Finite element quasi-interpolation and best approximation, ESAIM: Mathematical Modelling and Numerical Analysis, 51 (4), 1367-1385, (2017).
  * [23] Ern A., Vohralík M., Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs, SIAMJ. Sci. Comput. 35, 4, A1761-A1791, (2013).
  * [24] Fabrie P., Regularity of the solution of Darcy-Forchheimer’s equation. Nonlinear Anal. Theory Methods, 13, 1025-1045, (1989).
  * [25] Forchheimer P., Wasserbewegung durch Boden, Z. Ver. Deutsh. Ing. 45, 1782-1788, (1901).
  * [26] Gatica G.N., Ruiz-Baier R., and Tierra G., A mixed finite element method for Darcy’s equations with pressure dependent porosity, Mathematics of Computation, 85, 1-33, (2016).
  * [27] Hecht F., New development in FreeFem++, Journal of Numerical Mathematics, 20, (2012), 251-266.
  * [28] Lopez H., Molina B., Jose J.S., Comparison between different numerical discretizations for a Darcy-Forchheimer model, Electron. Trans. Numer. Anal. 34, 187-203, (2009).
  * [29] Lovadina C. and Stenberg R., Energy norm a posteriori error estimates for mixed finite element methods, Mathematics of Computation, 75(256), 1659-1674, (2006).
  * [30] Monk P., A posterior error indicators for Maxwell’s equations,J. Comput. Appl. Math., 100 (1998), 173-190.
  * [31] Pan H., Rui H., Mixed Element Method for Two-Dimensional Darcy-Forchheimer Model, J. Sci. Comput. 52, 563-587, (2012).
  * [32] Ruth D., Ma H., On the derivation of the Forchheimer equation by means of the averaging theorem, Transp. Porous Media 7(3), 255-264, (1992).
  * [33] Salas J.J., Lopez H., Molina B., An analysis of a mixed finite element method for a Darcy-Forchheimer model, Mathematical and Computer Modelling 57, 2325-2338, (2013).
  * [34] Sayah T., Semaan G. and Triki F., Finite element methods for the Darcy-Forchheimer problem coupled with the convection-diffusion-reaction problem, ESAIM Mathematical Modelling and Numerical Analysis, DOI:10.1051/m2an/2021066, (2021).
  * [35] Temam, R., Navier-Stokes equations: theory and numerical analysis. Vol. 343. American Mathematical Soc., 2001.
  * [36] Triki F., and T. Yin Inverse Conductivity Equation with Internal Data, Journal of Computational Mathematics, 2022.
  * [37] Verf$\ddot{u}$rth R., A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques, Mathematics, Wiley and Teubner, New-York, 1996.
  * [38] Verfürth R. , A posteriori Error Estimation Techniques for finite Element Methods,,Numerical Mathematics And Scientific Computation, Oxford, (2013).