A Posteriori Error Estimators for the Quasi-Newtonian Stokes Problem with a General Boundary Condition

Abstract

In this paper, we approach two nonlinear differential equations applied in fluid mechanics by finite element methods (FEM). Our objective is to approach the solution to these problems; the first one is the “p-Laplacian” problem and the second one is the “Quasi-Newtonian Stokes” problem with a general boundary condition. To study and analyze our solutions, we introduce the a posteriori error indicator; this technique allows us to control the error, and each is shown the equivalent between the true and the a posterior errors estimators. The performance of the finite element method by this type of general boundary condition is presented via different numerical simulations.

1. Introduction

One of the important objectives of numerical studies of differential problems is: to have a “realistic simulation”; for this, many researchers have concentrated on controlling the error by using the adaptive finite element methods (FEM) meshes. The adaptive meshes provide effective means of optimizing calculation with reasonable results; “the meshes are automatically modified by enhancing the scope of their applications”. The error estimation technique provides an assessment of the accuracy of the solutions obtained by the finite element solvers, see [1,2] for more precise details. These techniques can efficiently offer certain flow features (stagnation, reattachment points, and recirculation eddies, with small velocity magnitudes). Generally, this technique is based on a posteriori local error estimation, see [3].

Nonlinear differential equations are used to model complex problems in the sciences and engineering. Many studies have been developed to simplify these complex models. Among these models, we mention the “p-Laplacian” equation and the “Quasi-Newtonian Stokes” system. These nonlinear equations have had more attention in recent years, and one of the important domains using these equations is the glaciology domain. These equations model the dynamics of ice sheets or glaciers, see [4], and the evolution of glacier geometry, see [5,6]. Another application in the biological domain is the common use for blood flows, see [7] for details.

One of the important domains that uses these equations is fluid mechanics. For example, to model “viscoelastic” fluids we use the “Quasi-Newtonian Stokes” method, see [8] for more details. This quasi-Newtonian law was first proposed by Carreau et al. in [9]; a popular extension to the Carreau law is studied by Yasuda et al. in [10] and a closely related model, the Cross law, is investigated by Cross in [11]. The well-posedness of these problems is established in these papers [12,13]. Much mathematical modeling for complex fluids and numerical algorithms are applied to solve linear/nonlinear equations, see [14]. The fully non-linear elliptic problems in divergence form by using a mixed finite volume scheme are studied in [15]; the discontinuous Galerkin approximation was considered in [5] and the hybrid high-order scheme was studied in [16]. The paper where the authors studied the comportment of solutions of some non–linear diffusion problems and in the boundary–layer flow of a pseudo-plastic fluid is given [17]. The a posteriori error estimator for elliptic partial differential equations by using finite elements is presented in [18] and proves some (elliptic) a posteriori error estimators. The p-Laplacian problem is studied by using a FEM method in many papers, see [19,20,21,22] for details. In [23,24], the authors analyzed a posteriori error estimators for quasi-Newtonian problems with a homogeneous Dirichlet boundary condition. There are at least as many relevant references over the past 30 years on the quasi-Newtonian Stokes problem concerning the posterior estimate and its use in adaptivity; there were important developments in the Chinese [25,26,27], German [28,29,30], and French [31,32] schools for flow fluid in complex porous media using new boundary conditions [33] and combined mixed finite element [34].

This paper is organized as follows: in Section 2 and in order to study the p-Laplacian problem with the new boundary condition, we introduce some useful notations, state our main assumptions regarding the modeling equations, and following a finite element discretization, we calculate error estimators with respect to the true error. Section 3 contains the description of our second problem (quasi-Newtonian system) which is discretized by a mixed finite element scheme. The a posteriori error estimator is developed in terms of the Residus of variational formulas with respect to the real error. Section 4 presents two numerical examples using a velocity angle error indicator.

The authors propose the study of two nonlinear differential equations that have applications in fluid mechanics, using FEM. These problems are studied in two steps, within the “p-Laplacian” problem, and then solving the “Quasi-Newtonian Stokes” problem imposing a general boundary condition. To study and analyze our solutions, an a posteriori error indicator will be introduced. In this way, there is the possibility to control the error at each step and compare the calculated error with the a posteriori error. The obtained results are supported by several numerical simulations that we considered significant for the presentation of the research.

2. Approximation of p-Laplacian Equation by Finite Element Method

The section aims to study the p-Laplacian equation with a nonhomogeneous Robin boundary condition; we use FEM to approximate this model. To analyze the error, we use the a postriori estimates for the Dirichlet boundary condition, see [17,24,35]. First of all, we recall some useful properties of generalized nonlinear diffusion problems and we investigate the existence/uniqueness of the solution. Let $\mathsf{\Omega }$ be an open-bounded (connected) subset of ${ℝ}^d$ ($d=2,3$) whose boundary Γ = ∂Ω, and let $\beta \in {\mathbb{N}}^{*}$ with conjugate ${\beta }^{\prime }=\frac{\beta }{\beta +1}$.

2.1. Results of the p-Laplacian Operator

We consider the generalized nonlinear diffusion problem, defined by

$$-\nabla \cdot \phi \left(x,\nabla u\left(x\right)\right)=f\left(x\right)\mathrm{in}\mathsf{\Omega },$$

(1)

where the flux $\phi :\mathsf{\Omega }\times {ℝ}^2\to {ℝ}^2$ is assumed to satisfy the following assumptions:

There exist two positive constants $C_1,C_2$, and two functions $b_1\in L^1\left(\mathsf{\Omega }\right)$ and $b_2\in L^{{\beta }^{\prime }}\left(\mathsf{\Omega }\right)$ such that

$$\{\begin{array}{lll}\hfill & \left(\phi \left(x,y\right)-\phi \left(x,z\right),y-z\right)>0,\hfill & \left(H_1\right)\hfill \\ \hfill & \left(\phi \left(x,y\right),y\right)\ge C_1{\left|y\right|}^{\beta }-b_1\left(x\right),\hfill & \left(H_2\right)\hfill \\ \hfill & \phi \left(x,y\right)\le C_2{\left|y\right|}^{\beta -1}-b_2\left(x\right),\hfill & \left(H_3\right)\hfill \end{array}$$

for all $x,y,z\in \mathsf{\Omega }$ with $y\ne z$.

Under these assumptions the problem (1) has a unique solution $u\in W^{1,\beta }\left(\mathsf{\Omega }\right)$, and the functional $u\mapsto -\mathrm{div}\left(\phi \left(\cdot ,\nabla u(\cdot )\right)\right)$ is a Leray–Lions operator which satisfies

$$u\in {\left({\mathrm{L}}^{\mathsf{\beta }}\left(\mathsf{\Omega }\right)\right)}^2\mapsto \mathsf{\phi }\left(\cdot ,\nabla \mathrm{u}(\cdot )\right)\in {\left({\mathrm{L}}^{{\mathsf{\beta }}^{\prime }}\left(\mathsf{\Omega }\right)\right)}^2.$$

The p-Laplacian equation is defined by

$$-\nabla \cdot a\left(\nabla u\right)=f\mathrm{in}\mathsf{\Omega },$$

(2)

where the vector field $a\left(\zeta \right)={\left|\zeta \right|}^{\beta -2}\zeta$ and $\beta \in {\mathbb{N}}^{*}$.

The p-Laplacian equation is a nonlinear diffusion problem. Now, we recall some key lemmas useful to prove the monotonicity and continuity properties of such operators, see [12,20].

Let $z=\left(z_1,z_2\right)\in {ℝ}^2$ and let us define the following operator

$$A(\cdot ):u\in V\to A\left(u\right)=-\nabla \left({\left|\nabla u\right|}^{\beta -2}\nabla u\right)\in V^{\prime },$$

where space $V$ depends on the boundary condition and $V^{\prime }$ its dual. In order to prove the ellipticity of the problem, we need the following lemma.

Lemma 1. 

For all $y,z\in {ℝ}^2$, we have

$$\{\begin{array}{l}\left({\left|z\right|}^{\beta -2}z-{\left|y\right|}^{\beta -2}y,z-y\right)\ge \alpha {\left|z-y\right|}^{\beta -1}\mathrm{if}\beta \ge 2,\hfill \\ {\left(\left|z\right|+\left|y\right|\right)}^{2-\beta }\left({\left|z\right|}^{\beta -2}z-{\left|y\right|}^{\beta -2}y,z-y\right)\ge \alpha {\left|z-y\right|}^2\mathrm{if}1<\beta \le 2,\hfill \end{array}$$

where $\alpha >0$ is independent of $y$ and $z$.

Proof. 

See [12]. □

Proposition 2 (Ellipticity). 

For all $u,v\in V$, we have the following ellipticity properties

$$\{\begin{array}{lll}\left(A\left(u\right)-A\left(v\right),u-v\right)\ge \alpha \Vert u-v{\Vert }^{\beta }\hfill & \mathrm{if}& \hfill \beta \ge 2,\\ {\left(\Vert u\Vert +\Vert v\Vert \right)}^{2-\beta }\left(A\left(u\right)-A\left(v\right),u-v\right)\ge \alpha \Vert u-v{\Vert }^2\hfill & \mathrm{if}& \hfill 1<\beta \le 2,\end{array}$$

(3)

Proof. 

A direct consequence Lemma 2.1 from [12] for a detailed proof. □

To proof the continuity property, it needs the following lemma,

Lemma 3. 

For all $y,z\in {ℝ}^2$, we have

$$\{\begin{array}{lll}{\left|z\right|}^{\beta -2}z-{\left|y\right|}^{\beta -2}y\le \alpha \left|z-y\right|{\left(\left|z\right|+\left|y\right|\right)}^{\beta -2}\hfill & \mathrm{if}\hfill & \hfill \beta \ge 2,\\ {\left|z\right|}^{\beta -2}z-{\left|y\right|}^{\beta -2}y\le \alpha {\left|z-y\right|}^{\beta -1}\hfill & \mathrm{if}\hfill & \hfill 1<\beta \le 2,\end{array}$$

where $\alpha >0$ is independent of $y$ and $z$.

Proof. 

See [12]. □

Proposition 4 (Continuity). 

For all $u,v\in V$, we have the following continuity properties

$$\{\begin{array}{lll}\Vert A\left(u\right)-A{\left(v\right)\Vert }_{*}\le \alpha \Vert u-v\Vert {\left(\Vert u\Vert +\Vert v\Vert \right)}^{\beta -2}\hfill & \mathrm{if}& \hfill \beta \ge 2,\\ \Vert A\left(u\right)-A{\left(v\right)\Vert }_{*}\le \alpha \Vert u-v{\Vert }^{\beta -1}\hfill & \mathrm{if}& \hfill 1<\beta \le 2,\end{array}$$

(4)

here $\alpha >0$, independent to $u$ and $v$.

Proof. 

A (direct) consequence of the lemma 2.1, see also [12] for a detailed proof. □

2.2. Mathematical Problem

The practical applications of the p-Laplacian equations [4,5,36] have previously been mentioned; which, proved useful in a wide range of applications. For example, glacier dynamics are an important topic in engineering and hydrology. Thus, the ice flow is assumed to be an incompressible fluid with nonlinear viscosity. Next, let us consider the problem of the model

$$-\nabla \cdot a\left(\nabla u\right)=f \mathrm{in}\mathsf{\Omega },$$

(5)

$\mathrm{with}a\left(\zeta \right)={\left|\zeta \right|}^{\beta -2}\zeta$ and $f\in L^{{\beta }^{\prime }}\left(\mathsf{\Omega }\right)$. Note that if $\beta =2$, it coincides with the linear Laplacian operator, i.e., $A=-\mathsf{\Delta }$.

This equation can then be carried with the following new boundary conditions

$$\{\begin{array}{ll}\alpha \left(u\right)u+\left({\left|\nabla u\right|}^{\beta -2}\nabla u\right)\cdot n=g\hfill & {\mathrm{on}\mathsf{\Gamma }}_{ND},\hfill \\ u=0\hfill & {\mathrm{on}\mathsf{\Gamma }}_D,\hfill \end{array}$$

(6)

where $\alpha \in L^{\infty }\left(\mathsf{\Gamma }\right)$ ($\alpha >0$) and $g\in W^{1-\frac{1}{\beta },\beta }$(∂Ω).

Theorem 1. 

The problems (5) and (6) have a unique solution $u\in W^{1,\beta }\left(\Omega \right)$.

Proof. 

The proof is a consequence of ellipticity–continuity assumptions on the operator A. □

The systems (5) and (6) are equivalent to a minimization problem (see [12,13]) defined by: find $u\in V$ such that

$$J\left(v\right)\le J\left(u\right),\forall v\in V,$$

(7)

where

$$J\left(u\right)=\frac{1}{{}^{\beta }}{\int }_{\mathsf{\Omega }}{\left|\nabla u\right|}^{\beta }+\frac{1}{{}^2}{\int }_{{\mathsf{\Gamma }}_{ND}}\alpha {\left|u\right|}^2-{\int }_{\mathsf{\Omega }}fv-{\int }_{{\mathsf{\Gamma }}_{ND}}gv,$$

(8)

for all $v\in V$.

The function $J$ is continuous, strictly convex, and differentiable operator with $\underset{\Vert v\Vert \to \infty }{\lim }J\left(v\right)=+\infty$; it results, see [37], that $J$ is Gateaux differentiable and (7) admits a unique solution characterized by its variational formulation: find $u\in V$ such that

$${\int }_{\mathsf{\Omega }}\left({\left|\nabla u\right|}^{\beta -2}\nabla u\right)\nabla v+{\int }_{{\mathsf{\Gamma }}_{ND}}\alpha uv={\int }_{\mathsf{\Omega }}fv+{\int }_{{\mathsf{\Gamma }}_{ND}}gv,$$

(9)

for all $v\in V$.

Now, in order to analyze the finite element approximation of the problem (9), we consider a regular mesh $T_h$, $h>0$, of the domain $\mathsf{\Omega }$. For any element $T\in T_h$, we define: ${\omega }_T$ the set of elements share at least one edge with $T$; ${\tilde{\omega }}_T$ the set of all elements sharing at least one vertex with $T$; $\epsilon \left(T\right)$ the set of edges of $T$; $h_T$ the diameter of the simplex $T$; and $h=\underset{T\in T_h}{\max }{h}_T$. Respectively, for edge elements $E\in \partial T$ with $T\in T_h$: ${\omega }_E$ is the set of elements sharing at least one edge with $E$; ${\tilde{\omega }}_E$ the set of all elements sharing at least one vertex with $E$; and $h_E$ the diameter of a face $E$ of $T$. Let ${\epsilon }_h=\bigcup _{T\in T_h}\epsilon \left(T\right)$ designate the set of all edges; hence, one can divide it into interior and exterior edges such that ${\epsilon }_h={\epsilon }_{h,\mathsf{\Gamma }}{\cup }^{}{\epsilon }_{h,\mathsf{\Omega }}$ with

$${\epsilon }_{h,\mathsf{\Omega }}=\left\{E\in {\epsilon }_h:E\subset \mathsf{\Omega }\right\},{\epsilon }_{h,\mathsf{\Gamma }}=\left\{E\in {\epsilon }_h:E\subset \mathsf{\Gamma }\right\}.$$

Let $V_h\subset V=W^{1,\beta }\left(\mathsf{\Omega }\right)$, the finite dimensional spaces associated to regular partition of $\mathsf{\Omega }$, and a discrete weak formulation is defined using finite dimensional spaces as: find the vector $u_h\in V_h$ such that

$${\int }_{\mathsf{\Omega }}\left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\right)\nabla v_h+{\int }_{\mathsf{\Gamma }}\alpha u_h{v}_h={\int }_{\mathsf{\Omega }}f{v}_h+{\int }_{\mathsf{\Gamma }}g{v}_h,$$

(10)

for all $v_h\in V_h$.

2.3. A Posteriori Error Estimator

We illustrate the proposed technique with results of the a posteriori error estimation for the p-Laplacian problem. Denoting $u$ the solution of (5) and (6) and $u_h$ the approched solution of (10). Our aim is to estimate the velocity error $e=u-u_h\in V$ by using some important results.

Lemma 6. 

There is a constant $C>0$ for all the elements $K\in T_h$ and $v\in W^{1,\beta }\left(\Omega \right)$, such that

$$h_T{\Vert v\Vert }_{0,\beta ,\partial K}^{\beta }\le C\left({\Vert v\Vert }_{0,\beta ,K}^{\beta }+{\Vert v\Vert }_{1,\beta ,K}^{\beta }\right),$$

(11)

Proof. 

See [24]. □

Lemma 7 (Clement interpolation estimate). 

There is a constant $C>0$, for any $K\in T_h$ and for all $E\in \partial K$, let $v\in V$ and ${\pi }_h$ the operator of the interpolation of discontinuous functions defined by Clement satisfying

$${\Vert v-{\pi }_{h}v\Vert }_{0,\beta ,K}\le C{h}_K^{1-m}\sum _{K^{\prime }\in S_K}{\Vert v\Vert }_{0,\beta ,{\tilde{\omega }}_{K^{\prime }}}$$

(12)

for all $v\in H^1\left(S_K\right)$, and $m=0$ or $1$.

Where $S_K=\bigcup \left\{K^{\prime },K\bigcap K^{\prime }\ne \varnothing \right\}$. In particular, for $v_h$, is the quasi-interpolant of $v$defined by averaging as in

$$\Vert v-v_h\Vert {}_{0,\beta ,K}\le C{h}_k^{1-m}{\left|v\right|}_{1,\beta ,{\tilde{\omega }}_{K^{\prime }}}$$

(13)

with $m=0$ or $1$.

Proof. 

See [12,24]. □

The residual error estimator $R:V\mapsto ℝ$ is given by

$$⟨R,v⟩={\int }_{\mathsf{\Omega }}{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\nabla v+{\int }_{\mathsf{\Gamma }}\alpha u_{h}v-{\int }_{\mathsf{\Omega }}fv-{\int }_{\mathsf{\Gamma }}gv,$$

(14)

By applying the Green formula, we obtain

$$\begin{array}{ll}⟨R,v⟩& ={\displaystyle \sum _{K\in T_h}}\left\{{\int }_K{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\nabla v+{\int }_{\mathsf{\Gamma }\cap K}\alpha u_{h}v-{\int }_{K}fv-{\int }_{\mathsf{\Gamma }\cap K}gv\right\}\\ & ={\displaystyle \sum _{K\in T_h}}\{{\int }_K\left(-\nabla \cdot \left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\right)-f\right)v+{\int }_{\partial K}{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot nv+{\int }_{\mathsf{\Gamma }\cap K}\alpha u_{h}v+{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot nv-gv\}\\ & ={\displaystyle \sum _{K\in T_h}}\{{\int }_K\left(-\nabla \cdot \left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\right)-f\right)v+{\displaystyle \sum _{l\in \partial K}}{\int }_l{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot nv+{\displaystyle \sum _{l\in \mathsf{\Gamma }\cap \partial K}}{\int }_l\left(\alpha u_h+{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n-g\right)v\}.\end{array}$$

(15)

Because $R,v_h=0$ for all $v_h\in V_h$, we use the interpolation operator ${\pi }_h$ to obtain

$$\begin{array}{ll}⟨R,v⟩\hfill & ={\displaystyle \sum _{K\in T_h}}\{{\int }_K\left(-\nabla \cdot \left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\right)-f\right)\left(v-{\pi }_{h}v\right)+{\displaystyle \sum _{l\in \partial K}}{\int }_l{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\left(v-{\pi }_{h}v\right)\hfill \\ \hfill & +{\displaystyle \sum _{l\in \mathsf{\Gamma }\cap \partial K}}{\int }_l\left(\alpha u_h+{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n-g\right)\left(v-{\pi }_{h}v\right)\}\hfill \end{array}$$

(16)

and

$$\begin{array}{ll}⟨R,v⟩\hfill & \le {\displaystyle \sum _{K\in T_h}}\{{\int }_K{\left|-\nabla \cdot \left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\right)-f\right|}_{0,{\beta }^{\prime },K}{\left|v-{\pi }_{h}v\right|}_{0,\beta ,K}+{\displaystyle \sum _{l\in \partial K}}{\int }_l{\left|{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right|}_{0,{\beta }^{\prime },l}{\left|v-{\pi }_{h}v\right|}_{0,\beta ,l}\hfill \\ \hfill & +{\displaystyle \sum _{l\in \mathsf{\Gamma }\cap \partial K}}{\int }_l{\left|\alpha u_h+{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n-g\right|}_{0,{\beta }^{\prime },l}{\left|v-{\pi }_{h}v\right|}_{0,\beta ,l}\}.\hfill \end{array}$$

(17)

So, by (12), we get

$$\begin{array}{ll}⟨R,v⟩\le & C{\displaystyle \sum _{K\in T_h}}{h}_K{\left|-\nabla \cdot \left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\right)-f\right|}_{0,{\beta }^{\prime },K}\cdot {\Vert v\Vert }_{1,\beta ,K}+{\displaystyle \sum _{K\in T_h}}\{C{\displaystyle \sum _{l\in \partial K}}{\left|{\left[{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right]}_l\right|}_{0,{\beta }^{\prime },l}.{\left|v-{\pi }_{h}v\right|}_{0,\beta ,l}\\ & +C{\displaystyle \sum _{l\in \mathsf{\Gamma }\cap \partial K}}{\left|\alpha u_h+{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n-g\right|}_{0,{\beta }^{\prime },l}.{\left|v-{\pi }_{h}v\right|}_{0,\beta ,l}\}.\end{array}$$

(18)

For the first term of the second member, we have

$$I_1={\displaystyle \sum _{K\in T_h}}{h}_K{\left|-\nabla \cdot \left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\right)-f\right|}_{0,{\beta }^{\prime },K}.{\Vert v\Vert }_{1,\beta ,K}\le C{\left({\displaystyle \sum _{K\in T_h}}{h}_K^{{\beta }^{\prime }}{\left|-\nabla \cdot \left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\right)-f\right|}_{0,{\beta }^{\prime },K}^{{\beta }^{\prime }}\right)}^{\frac{1}{{\beta }^{\prime }}}{\Vert v\Vert }_{1,\beta ,K}.$$

(19)

Now, we increase the second term as follows

$$\begin{array}{ll}{I}_2& ={\displaystyle \sum _{l\in \partial K}}{\left|{\left[{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right]}_l\right|}_{0,{\beta }^{\prime },l}\cdot {\left|v-{\pi }_{h}v\right|}_{0,\beta ,l}\\ & ={\displaystyle \sum _{l\in \partial K}}{h}_l^{\frac{1}{{\beta }^{\prime }}}{\left|{\left[{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right]}_l\right|}_{0,{\beta }^{\prime },l}.h_l^{-\frac{1}{{\beta }^{\prime }}}{\left|v-{\pi }_{h}v\right|}_{0,\beta ,l}\\ & \le {\left({\displaystyle \sum _{l\in \partial K}}{h}_l{\left|{\left[{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right]}_l\right|}_{0,{\beta }^{\prime },l}^{{\beta }^{\prime }}\right)}^{\frac{1}{{\beta }^{\prime }}}\times {\left({\displaystyle \sum _{l\in \partial K}}{h}_l^{-\frac{\beta }{{\beta }^{\prime }}}{\left|v-{\pi }_{h}v\right|}_{0,\beta ,l}^{\beta }\right)}^{\frac{1}{\beta }}.\end{array}$$

(20)

Using (11), we have

$$\begin{array}{ll}{\displaystyle \sum _{l\in \partial K}}{h}_l^{-\frac{\beta }{{\beta }^{\prime }}}{\left|v-{\pi }_{h}v\right|}_{0,\beta ,l}^{\beta }& \le C{\displaystyle \sum _{K\in T_h}}{h}_K^{-\frac{\beta }{{\beta }^{\prime }}}{\left|v-{\pi }_{h}v\right|}_{0,\beta ,\partial K}^{\beta }\\ & \le C{\displaystyle \sum _{K\in T_h}}{h}_K^{-\left(1+\frac{\beta }{{\beta }^{\prime }}\right)}\left({\left|v-{\pi }_{h}v\right|}_{0,\beta ,\partial K}^{\beta }+h_K^{\beta }{\left|v-{\pi }_{h}v\right|}_{1,\beta ,K}^{\beta }\right)\\ & \le C{\displaystyle \sum _{K\in T_h}}{h}_K^{-\left(1+\frac{\beta }{{\beta }^{\prime }}\right)}\left(h_K^{\beta }{\left|v\right|}_{1,\beta ,K}^{\beta }+h_K^{\beta }{\left|v\right|}_{1,\beta ,K}^{\beta }\right).\end{array}$$

(21)

As $1+\frac{\beta }{{\beta }^{\prime }}=\beta$, it implies that

$$\begin{array}{ll}{I}_2\hfill & \le {\displaystyle \sum _{l\in \partial K}}{\left|{\left[{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right]}_l\right|}_{0,{\beta }^{\prime },l}.{\left|v-{\pi }_{h}v\right|}_{0,\beta ,l}\hfill \\ \hfill & \le C{\left({\displaystyle \sum _{l\in \partial K}}{h}_l{\left|{\left[{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right]}_l\right|}_{0,{\beta }^{\prime },l}^{{\beta }^{\prime }}\right)}^{\frac{1}{{\beta }^{\prime }}}\times {\left|v\right|}_{1,\beta ,\mathsf{\Omega }}.\hfill \end{array}$$

(22)

Now, it is simple to increase the third term as follows

$$\begin{array}{ll}{I}_3=& {\displaystyle \sum _{l\in \mathsf{\Gamma }\cap \partial K}}{h}_l^{\frac{1}{{\beta }^{\prime }}}{\left|\alpha u_{h}v+{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot nv-g\right|}_{0,{\beta }^{\prime },l}\cdot h_l^{-\frac{1}{{\beta }^{\prime }}}{\left|v-{\pi }_{h}v\right|}_{0,\beta ,l}\\ & \le {\left({\displaystyle \sum _{l\in \mathsf{\Gamma }\cap \partial K}}{h}_l{\left|\alpha u_h+{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n-g\right|}_{0,{\beta }^{\prime },l}^{{\beta }^{\prime }}\right)}^{\frac{1}{{\beta }^{\prime }}}\times {\left({\displaystyle \sum _{l\in \partial K}}{h}_l^{-\frac{\beta }{{\beta }^{\prime }}}{\left|v-{\pi }_{h}v\right|}_{0,\beta ,l}^{\beta }\right)}^{\frac{1}{\beta }}\\ & \le {\left({\displaystyle \sum _{l\in \mathsf{\Gamma }\cap \partial K}}{h}_l{\left|\alpha u_h+{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n-g\right|}_{0,{\beta }^{\prime },l}^{{\beta }^{\prime }}\right)}^{\frac{1}{{\beta }^{\prime }}}\times {\left|v\right|}_{1,\beta ,\mathsf{\Omega }}.\end{array}$$

(23)

Hence, by combining (19), (22), and (23) we get

$$\begin{array}{ll}⟨R,v⟩\le C& \{{\left({\displaystyle \sum _{K\in T_h}}{h}_K{\left|\nabla \cdot \left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\right)+f\right|}_{0,{\beta }^{\prime },K}^{{\beta }^{\prime }}\right)}^{\frac{1}{{\beta }^{\prime }}}+{\left({\displaystyle \sum _{l\in \partial K}}{h}_l{\left|{\left[{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right]}_l\right|}_{0,{\beta }^{\prime },l}^{\beta }\right)}^{\frac{1}{{\beta }^{\prime }}}\\ & +C{\left({\displaystyle \sum _{l\in \mathsf{\Gamma }\cap \partial K}}{h}_l{\left|\alpha u_h+{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n-g\right|}_{0,{\beta }^{\prime },l}^{{}^{{\beta }^{\prime }}}\right)}^{\frac{1}{{\beta }^{\prime }}}\}\Vert v\Vert .\end{array}$$

(24)

As $a^{\frac{1}{{\beta }^{\prime }}}+b^{\frac{1}{{\beta }^{\prime }}}\le 2^{\frac{1}{\beta }}{(a+b)}^{\frac{1}{{\beta }^{\prime }}}$ for $a$ and $b\ge 0$, we have the following estimate

$${\Vert \mathrm{R}\Vert }_{*}=\Vert \mathrm{A}\left(\mathrm{u}\right)-\mathrm{A}{\left({\mathrm{u}}_{\mathrm{h}}\right)\Vert }_{*}\le {\left(\mathrm{C}{\sum }_{\mathrm{K}\in {\mathrm{T}}_{\mathrm{h}}}{\mathsf{\eta }}_{\mathrm{K}}^{{\mathsf{\beta }}^{\prime }}\right)}^{\frac{1}{{\mathsf{\beta }}^{\prime }}},$$

with the contribution element of ${\eta }_K$, the residual error estimator is given by

$${\eta }_K^{{\beta }^{\prime }}=h_K^{{\beta }^{\prime }}\Vert R_K\Vert {}_{0,{\beta }^{\prime },K}^{{\beta }^{\prime }}+{\sum }_{l\in \partial K}{h}_l\Vert R_l\Vert {}_{0,{\beta }^{\prime },l}^{{\beta }^{\prime }},$$

(25)

his components are given by

$$R_K={\left\{\nabla \cdot \left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\right)+f\right\}}_K,$$

(26)

and

$$R_l=\{\begin{array}{lll}\frac{1}{2}{\left[{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right]}_l\hfill & \mathrm{if}\hfill & l\in {\epsilon }_{h,\mathsf{\Omega }},\hfill \\ \alpha \left(u_h\right)u_h+{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n-g\hfill & \mathrm{if}\hfill & l\in {\epsilon }_{h,\mathsf{\Gamma }}.\hfill \end{array}$$

(27)

where ${[\cdot ]}_l$ is the jump of the derivative of $u_h$ over the interior edge $l=T\cap S$, defined by

$${\left[{\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right]}_l=\left(\left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right){|}_T-\left({\left|\nabla u_h\right|}^{\beta -2}\nabla u_h\cdot n\right){|}_S\right){\overrightarrow{n}}_{E,T}.$$

(28)

It remains to connect ${\Vert R\Vert }_{*}$ and $\Vert u-u_h\Vert$, which uses a coercivity property of the operator $A:V\mapsto W^{-1,{\beta }^{\prime }}$ for $\beta \in \left]1,2\right[$, where we have

$$⟨A\left(u\right)-A\left(v\right),u-v⟩\ge C\frac{{\Vert u-v\Vert }^2}{{\left(\Vert u\Vert +\Vert v\Vert \right)}^{2-\beta }} ,$$

(29)

for all $u,v\in W^{1,\beta }\left(\mathsf{\Omega }\right)$, where $C$ is a constant that does not depend on either $u$ or $v$. By taking $v=0$, in (29), with $u$ as the solution of problem (5), we get (since $A\left(0\right)=0$)

$$⟨A\left(u\right),u⟩\ge \alpha {\Vert u\Vert }^{\beta }.$$

As $A\left(u\right)=f$, then $\alpha {\Vert u\Vert }^{\beta }\le f,u\le {\Vert f\Vert }_{W^{-1,{\beta }^{\prime }}\left(\mathsf{\Omega }\right)}.\Vert u\Vert$, which implies (since the injection of $L^{{\beta }^{\prime }}\left(\mathsf{\Omega }\right)$ in $W^{-1,{\beta }^{\prime }}\left(\mathsf{\Omega }\right)$ is continuous)

$$\Vert u\Vert \le {\left(\frac{1}{C}{\Vert f\Vert }_{W^{-1,{\beta }^{\prime }}\left(\mathsf{\Omega }\right)}\right)}^{\frac{1}{{}^{\beta }-1}}\le C{\left({\Vert f\Vert }_{L^{{\beta }^{\prime }}\left(\mathsf{\Omega }\right)}\right)}^{\frac{1}{{}^{\beta }-1}}.$$

Similarly, for $u_h$, we have

$${\left(\Vert u\Vert +\Vert u_h\Vert \right)}^{2-\beta }\le C{\left({\Vert f\Vert }_{L^{{\beta }^{\prime }}\left(\mathsf{\Omega }\right)}\right)}^{\frac{2-\beta }{{}^{\beta }-1}}.$$

(30)

Finally, from (19), (22), and (23), we obtain the following result.

Theorem 8. 

For any mixed finite element approximation defined on regular grids $T_h$, the residual estimator satisfies

$$\Vert e\Vert \le C{\left({\displaystyle \sum _{K\in T_h}}{\eta }_K^{{\beta }^{\prime }}\right)}^{\frac{1}{{\beta }^{\prime }}},$$

where C is independent of $\beta$, $\Omega$, $a$, and $f$.

Proof of Theorem 8. 

A direct consequence of these inequalities (15) to (29). □

Remark 9. 

Some remarks are in order: for $\beta \ge 2$, we obtain

$$\Vert e\Vert \le C{\left({\displaystyle \sum _{K\in T_h}}{\eta }_K^{{\beta }^{\prime }}\right)}^{\frac{1}{\beta }}$$

where $C$ is independent of $u$ and $f$, in this case and instead of (29), we have

$$⟨A\left(u\right)-A\left(v\right),u-v⟩\ge \alpha {\Vert u-v\Vert }^{\beta },$$

$\forall u,v\in W^{1,\beta }\left(\Omega \right).$

This result also holds for stable (and unstable) mixed approximations defined on a regular triangulation.

3. Approximation of Quasi-Newtonian Stokes Problem by Finite Element Method

In reality, the comportment of fluids is more complex; this type of fluid is modeled by nonlinear operators or tensors. Let us start this part with some definitions of fluid that we will treat. The type of fluid that is more important is modeled by a nonlinear operator; the problem is the quasi-Newtonian Stokes flow. For many details on the different models and algorithms solving these problems, see [14].

3.1. Definition of Quasi-Newtonian Stokes Problem

We denote by $u$ the velocity vector, $p$ the pressure, $\sigma$ the stress tensor, $D\left(u\right)$ the symmetric gradient of the velocity vector, and $f$ the external forces.

Definition 1 (Quasi-Newtonian fluid). 

The fluid is said to be quasi-Newtonian when there exists a positive function $\mu :{ℝ}^{+}\mapsto {ℝ}^{+}$ called the “viscosity” function, such that the stress deviator $\sigma$ is expressed as

$$\sigma =-p\mathbb{I}+2\mu \left({\left|2D\left(u\right)\right|}^2\right)D\left(u\right).$$

(31)

There are two classical laws of quasi-Newtonian Stokes problems, the first one is the “Carreau law viscosity” and the second one is the “power law viscosity”.

Definition 2 (Carreau law viscosity). 

The Carreau law expresses the viscosity as

$$\mu \left(s\right)={\mu }_{\infty }+\left({\mu }_0-{\mu }_{\infty }\right){\left(1+\lambda s\right)}^{\frac{-1+n}{2}}, \forall s\in {ℝ}^{+},$$

(32)

where ${\mu }_0$, ${\mu }_{\infty }$, $\lambda ,n\in {ℝ}^{+*}$ are given real constants satisfying ${\mu }_0\ge {\mu }_{\infty }$ when $n\le 1$ and ${\mu }_0\le {\mu }_{\infty }$ when $n\ge 1$.

The Figure 1 presents Carreau law viscosity for $n<1$ and for $n>1$. We can see that under simplification conditions, the Carreau law viscosity can be rewritten as power law viscosity defined by

$$\mu \left(\zeta \right)=K{\zeta }^{\frac{-1+n}{2}}, \forall \zeta \in ℝ,$$

(33)

where $K$ and $n\in {ℝ}^{+*}$.

Note that, when $n=1$, both the power law and Carreau law are reduced to a Newtonian fluid model with constant viscosity. When $n<1$, the viscosity is decreasing with the shear rate and the fluid is said to be “shear thinning” or “pseudoplastic”. The tensor norm and the shear rate are defined as follows (see [14], Definition 1.9 for more precise details).

Definition 3 (Tensor norm). 

The following tensor norm is defined as

$${\left|\tau \right|}^2=\frac{\tau :\tau }{2}=\frac{1}{2}{\displaystyle \sum _{i=1}^3}{\displaystyle \sum _{j=1}^3}{\tau }_{i,j}^2,\hspace{1em}\hspace{1em}\forall \tau \in {ℝ}^{3\times 3}.$$

(34)

Definition 4 (Shear rate). 

The shear rate, denoted by $\dot{\gamma }$, is defined by

$$\dot{\gamma }=\left|2D\left(u\right)\right|.$$

Note that, the stress tensor $\sigma$ in Newtonian fluid flow is “generally” defined by $\sigma =-p\mathbb{I}+2\mu D\left(u\right)$, where $\mu$ is the fluid viscosity (a bound function).

3.2. Mathematical Problem

To simplify this study, the stationary case where the flow is incompressible remains; hence, we neglect the inertia term which will allow us to focus on the nonlinearity of results (from the viscosity law). The governing system can be written as

$$\{\begin{array}{ll}-\mathrm{div}\left(2\mu (|2D\left(u\right)){|}^2)D\left(u\right)\right)+\nabla p\hfill & =f \mathrm{in}\mathsf{\Omega },\hfill \\ \mathrm{div}u\hfill & =0\mathrm{in}\mathsf{\Omega },\hfill \end{array}$$

(35)

where $f\in X^{\prime }$, we define the general boundary condition $C_{a,\mu }$ by

$$C_{a,\mu }:a\left(u\right)u+\left(2\mu \left({\left|2D\left(u\right)\right|}^2\right)D\left(u\right)-pI\right)\cdot n=g \mathrm{on}\mathsf{\Gamma },$$

(36)

where $g\in \mathsf{\Gamma }$ and $a\left(u\right)$ are a bounded function defined in the boundary.

The existence and unique solution $\left(u,p\right)$ of the problem (35) and (36) is defined in the product space $X\times M$ with $X={\left(H^1\left(\mathsf{\Omega }\right)\right)}^2$ and $M=L^2\left(\mathsf{\Omega }\right).$ The operator $A=2\mu \left({\left|2D\left(u\right)\right|}^2\right)D\left(u\right):X\mapsto X^{\prime }$, appearing in (35), satisfies the following two propositions.

Proposition 10. 

For all $u,v\in X$, we have

$$⟨\mathrm{A}\left(\mathrm{u}\right)-\mathrm{A}\left(\mathrm{v}\right),\mathrm{u}-\mathrm{v}⟩\ge \mathrm{C}\left(\mathrm{u},\mathrm{v}\right){\Vert \mathrm{u}-\mathrm{v}\Vert }^2,$$

(37)

where

$$C\left(u,v\right)={\left\{{\alpha }_1+{\alpha }_2{\left(\Vert u\Vert +\Vert v\Vert \right)}^{2-\beta }\right\}}^{-1},$$

with ${\alpha }_1, {\alpha }_2>0$.

Proposition 11. 

For all $u,v$ in $X$ and all $\beta >1$, there exists a positive constant $C>0$ such that

$${\Vert A\left(u\right)-A\left(v\right)\Vert }_{X^{\prime }}\le C\left(u,v\right){\Vert u-v\Vert }_X.$$

(38)

Proof of Proposition 10 and 11. 

See [38]. □

It should be noted that conditions (37) and (38) are satisfied for a fluid with a power law viscosity, see, e.g., [12,39]. For the existence and uniqueness of the solution, we refer to the works [17,40] by Baranger and Najib. The variational formulation to the problem (35) and (36) is equivalent to: find $\left(u,p\right)\in X\times M$ such that

$$\{\begin{array}{ll}\langle A\left(u\right),v\rangle +\langle Bv,p\rangle \hfill & =\langle f,v\rangle ,\hfill \\ \langle Bu,q\rangle \hfill & =\langle g,q\rangle ,\hfill \end{array}$$

(39)

for all $\left(v,q\right)\in X\times M$, where

$$\langle A\left(u\right),v\rangle ={\int }_{\mathsf{\Omega }}\left(2\eta \left({\left|2D\left(u\right)\right|}^2\right):D\left(u\right)\right)D\left(v\right)dx+{\int }_{\mathsf{\Gamma }}\alpha uvd\gamma \left(x\right),$$

(40)

and

$$\langle Bv,q\rangle ={\int }_{\mathsf{\Omega }}q \mathrm{div}\left(v\right).$$

(41)

Theorem 12. 

The problem (39) has a unique solution$\left(u,p\right)\in X\times M$ if and only if (40) has a unique solution and $B$ is surjective satisfying the inf-sup condition.

Proof. 

See [17,41]. □

The functional spaces X and $M$ are Hilbert spaces, then $\left(u,p\right)$ (resp. $\left(U,P\right)$) is a solution of problem (39) associated to the couple $\left(f,g\right)$ (resp. $\left(F,G\right)$) such that

$$\Vert u-U\Vert +\Vert p-P\Vert \le C\left({\Vert f-F\Vert }_{*}+{\left|g-G\right|}_{*}\right),$$

where $C$ is a positive constant, see [42] for more details.

We have to describe the discrete approximation problem corresponding to (39)–(41) by using notations of Section 2; to this aim, we define the finite element spaces $\left(X_h,M_h\right)\subset \left(X,M\right)$, where the discrete approximation of the quasi-Newtonian flow problem (with $C_{a,\mu }$ boundary condition) can be equivalently written as: find $\left(u_h,p_h\right)\in X_h,M_h$ such that

$$\{\begin{array}{ll}{\int }_{\mathsf{\Omega }}\left(2\mu {\left|2D\left(u_h\right)\right|}^2\right)D\left(u_h\right):D\left(v_h\right)+p_h{v}_{h}d+{\int }_{\mathsf{\Gamma }}a{u}_h{v}_{h}d\gamma ={\int }_{\mathsf{\Omega }}fvdx+{\int }_{\mathsf{\Gamma }}g{v}_{h}d\gamma ,x\hfill & \hfill \\ \hfill & \hfill \\ {\int }_{\mathsf{\Omega }}{q}_h\mathrm{div}{u}_{h}dx=0,\hfill & \hfill \end{array}$$

(42)

for all $\left(v_h,q_h\right)\in X_h\times M_h$.

Under properties (3.2) and (3.2), and assuming the inf-sup condition is satisfied in the discrete approximation problem (42), we obtain one solution, $\left(u_h,p_h\right)\in X_h\times M_h$.

3.3. A Posteriori Error Estimator

A posteriori error indicator and residual error for quasi-Newtonian problems with Dirichlet and Newman boundary conditions are developed in these papers [17,24,35]. Based on these results, one can adapt their estimates to the case of our problem (42); by using the notation ${\mu }_h=2\mu {\left|2D\left(u_h\right)\right|}^2$, we get

$$⟨R,v⟩={\displaystyle \sum _{K\in T_h}}{\int }_{K}div\left(2\mu {\left|2D\left(u_h\right)\right|}^2\right)D\left(u_h\right)D\left(v\right)-{\displaystyle \sum _{K\in T_h}}{\int }_K{p}_h\mathrm{div}v-{\displaystyle \sum _{K\in T_h}}{\int }_{K}fv+{\displaystyle \sum _{K\in T_h}}{\int }_{\partial K\cap \mathsf{\Gamma }}\alpha u_{h}v+{\displaystyle \sum _{K\in T_h}}{\int }_{\partial K\cap \mathsf{\Gamma }}gv.$$

(43)

From Green’s formula, and for any element $K$, we obtain

$$\begin{array}{ll}⟨R,v⟩& ={\displaystyle \sum _{K\in T_h}}\left({\int }_K\left(A\left(u_h\right)+\nabla p_h-f\right)v+{\int }_{\partial K}\left({\mu }_{h}D\left(u_h\right)n-p_{h}n\right)v+{\int }_{\partial K\cap \mathsf{\Gamma }}\alpha u_{h}v+{\int }_{\partial K\cap \mathsf{\Gamma }}gv\right)\\ & ={\displaystyle \sum _{K\in T_h}}\left({\int }_K\left(A\left(u_h\right)+\nabla p_h-f\right)v+{\int }_{\partial K}\left({\mu }_{h}D\left(u_h\right)n-p_{h}n\right)v+{\int }_{\partial K\cap \mathsf{\Gamma }}\alpha u_{h}v+\left({\mu }_{h}D\left(u_h\right)n-p_{h}n\right)v-gv\right)\\ & ={\displaystyle \sum _{K\in T_h}}\left({\int }_K\left(A\left(u_h\right)+\nabla p_h-f\right)v+{\displaystyle \sum _{l\in {\epsilon }_h}}{\int }_l\left({\mu }_{h}D\left(u_h\right)n-p_{h}n\right)v+{\displaystyle \sum _{l\in {\epsilon }_{\mathsf{\Gamma }}}}{\int }_l\alpha u_{h}v+\left({\mu }_{h}D\left(u_h\right)n-p_{h}n\right)v-gv\right)\end{array}$$

where ${\epsilon }_{\mathsf{\Gamma }}$ is the set of all edges of all elements $T_h$ divided into interior and exterior edges ${\epsilon }_h={\epsilon }_{h,\mathsf{\Gamma }}\cup {\epsilon }_{h,\mathsf{\Omega }}$ with ${\epsilon }_{h,\mathsf{\Omega }}=\left\{E\in {\epsilon }_h:E\subset \mathsf{\Omega }\right\}$ and ${\epsilon }_{h,\mathsf{\Gamma }}=\left\{E\in {\epsilon }_h:E\subset \mathsf{\Gamma }\right\}$. As $⟨R,v_h⟩=0$ for all $v_h\in X_h$, and for $v_h={\pi }_{h}v$, we get

$$\begin{array}{ll}⟨R,v⟩=⟨R,v-{\pi }_{h}v⟩& ={\displaystyle \sum _{K\in T_h}}{\int }_K\left(A\left(u_h\right)+\nabla p_h-f\right)\left(v-{\pi }_{h}v\right)\\ & +{\displaystyle \sum _{l\in {\epsilon }_h}}{\int }_l\left({\mu }_{h}D\left(u_h\right)n-p_{h}n\right)\left(v-{\pi }_{h}v\right)\\ & +{\displaystyle \sum _{l\in {\epsilon }_{\mathsf{\Gamma }}}}{\int }_l\left(\alpha u_h+\left({\mu }_{h}D\left(u_h\right)n-p_{h}n\right)-g\right)\left(v-{\pi }_{h}v\right).\end{array}$$

(44)

Hence

$$\begin{array}{ll}⟨R,v⟩\hfill & \le {\displaystyle \sum _{K\in T_h}}{\left|A\left(u_h\right)+\nabla p_h-f\right|}_{0,K}{\left|v-{\pi }_{h}v\right|}_{0,K}+{\displaystyle \sum _{l\in {\epsilon }_h}}{\left|{\mu }_{h}D\left(u_h\right)n-p_{h}n\right|}_{0,l}{\left|v-{\pi }_{h}v\right|}_{0,l}\hfill \\ \hfill & \hfill \\ \hfill & +{\displaystyle \sum _{l\in {\epsilon }_{\mathsf{\Gamma }}}}{\left|\alpha u_h+\left({\mu }_{h}D\left(u_h\right)n-p_{h}n\right)-g\right|}_{0,l}{\left|v-{\pi }_{h}v\right|}_{0,l}.\hfill \end{array}$$

(45)

In order to deal with the three terms of the right-hand side, we write

$$\begin{array}{ll}\hfill & {\displaystyle \sum _{K\in T_h}}{\left|A\left(u_h\right)+\nabla p_h-f\right|}_{0,K}{\left|v-{\pi }_{h}v\right|}_{0,K}\le {\left({\displaystyle \sum _{K\in T_h}}{h}_K^2{\left|A\left(u_h\right)+\nabla p_h-f\right|}_{0,K}^2\right)}^{\frac{1}{2}}{\left({\displaystyle \sum _{K\in T_h}}{h}_K^{-2}{\left|v-{\pi }_{h}v\right|}_{0,K}^2\right)}^{\frac{1}{2}}\hfill \\ \hfill & \hfill \end{array}$$

(46)

we use Lemma 3 (when $m=0$) to obtain

$${\displaystyle \sum _{K\in T_h}}{\left|A\left(u_h\right)+\nabla p_h-f\right|}_{0,K}{\left|v-{\pi }_{h}v\right|}_{0,K}\le {\left({\displaystyle \sum _{K\in T_h}}{h}_K^2{\left|A\left(u_h\right)+\nabla p_h-f\right|}_{0,K}^2\right)}^{\frac{1}{2}}\Vert v\Vert .$$

(47)

Defining the vector $s$

$$s={\left|{\mu }_{h}D\left(u_h\right)n-p_{h}n\right|}_{0,l}={\left|{\left[\sigma \left(u_{h,}{p}_h\right)n\right]}_l\right|}_{0,l}={\left|{\left[{\sigma }^{h}n\right]}_l\right|}_{0,l},$$

the second term of the right-hand side is rewritten as

$${\displaystyle \sum _{l\in {\epsilon }_h}}s{\left|v-{\pi }_{h}v\right|}_{0,l}\le {\left({\displaystyle \sum _{l\in {\epsilon }_h}}{h}_l{\left|s\right|}^2\right)}^{\frac{1}{2}}{\left({\displaystyle \sum _{l\in {\epsilon }_h}}{h}_l^{-1}{\left|v-{\pi }_{h}v\right|}_{0,l}^2\right)}^{\frac{1}{2}}.$$

(48)

Observe that

$$\begin{array}{ll}{\left({\displaystyle \sum _{K\in T_h}}{\displaystyle \sum _{l\in {\epsilon }_h}}{h}_l^{-1}{\left|v-{\pi }_{h}v\right|}_{0,l}^2\right)}^{\frac{1}{2}}\hfill & \le {\left({\displaystyle \sum _{K\in T_h}}C{h}_K^{-1}{\left|v-{\pi }_{h}v\right|}_{0,l}^2\right)}^{\frac{1}{2}}\hfill \\ \hfill & =C{\left({\displaystyle \sum _{K\in T_h}}{h}_K^{-1}{\left|v-{\pi }_{h}v\right|}_{0,l}^2\right)}^{\frac{1}{2}},\hfill \end{array}$$

(49)

then, one can use lemma 3 (with $\beta =2$) to get

$$\begin{array}{ll}{\left({\displaystyle \sum _{K\in T_h}}{\displaystyle \sum _{l\in {\epsilon }_h}}{h}_l^{-1}{\left|v-{\pi }_{h}v\right|}_{0,l}^2\right)}^{\frac{1}{2}}\hfill & \le C{\left({\displaystyle \sum _{K\in T_h}}{h}_K^{-2}{\left|v-{\pi }_{h}v\right|}_{0,K}^2+{\displaystyle \sum _{K\in T_h}}{h}_K^2{\left|v-{\pi }_{h}v\right|}_{1,K}^2\right)}^{\frac{1}{2}}\hfill \\ \hfill & \le C{\left({\displaystyle \sum _{K\in T_h}}{\left|v\right|}_{1,K}^2\right)}^{\frac{1}{2}}=C\Vert v\Vert .\hfill \end{array}$$

(50)

Thus, the second term of (45) is increased by

$$C{\left({\displaystyle \sum _{l\in {\epsilon }_h}}{\left|{\left[\sigma \left(u_{h,}{p}_h\right)n\right]}_l\right|}_{0,l}^2\right)}^{\frac{1}{2}}\Vert v\Vert .$$

To treat the last term of (45), we write §

$${\displaystyle \sum _{l\in {\epsilon }_h}}{\left|\alpha u_h+s-g\right|}_{0,l}{\left|v-{\pi }_{h}v\right|}_{0,l}\le {\left({\displaystyle \sum _{l\in {\epsilon }_h}}{h}_l{\left|\alpha u_h+s-g\right|}^2\right)}^{\frac{1}{2}}{\left({\displaystyle \sum _{l\in {\epsilon }_h}}{h}_l^{-1}{\left|v-{\pi }_{h}v\right|}_{0,l}^2\right)}^{\frac{1}{2}},$$

by using lemma 3 and (50), the last term of (45) is increased by

$$C{\left({\displaystyle \sum _{l\in {\epsilon }_h}}{\left|{\left[\left(\alpha u_h+\sigma \left(u_{h,}{p}_h\right)-g\right)n\right]}_l\right|}_{0,l}^2\right)}^{\frac{1}{2}}\Vert v\Vert .$$

(51)

Finally, we conclude

$$\begin{array}{ll}⟨R,v⟩\hfill & \le C{\left({\displaystyle \sum _{K\in T_h}}{h}_K^2{\left|A\left(u_h\right)+\nabla p_h-f\right|}_{0,K}^2\right)}^{\frac{1}{2}}\Vert v\Vert +C{\left({\displaystyle \sum _{l\in {\epsilon }_h}}{\left|{\left[\sigma \left(u_{h,}{p}_h\right)n\right]}_l\right|}_{0,l}^2\right)}^{\frac{1}{2}}\Vert v\Vert \hfill \\ \hfill & \hfill \\ \hfill & +C{\left({\displaystyle \sum _{l\in {\epsilon }_h}}{\left|{\left[\left(\alpha u_h+\sigma \left(u_{h,}{p}_h\right)-g\right)n\right]}_l\right|}_{0,l}^2\right)}^{\frac{1}{2}}\Vert v\Vert ,\hfill \end{array}$$

(52)

where

$$\begin{array}{ll}{\Vert R\Vert }_{*}\hfill & \le C{\left({\displaystyle \sum _{K\in T_h}}{h}_K^2{\left|A\left(u_h\right)+\nabla p_h-f\right|}_{0,K}^2\right)}^{\frac{1}{2}}+C{\left({\displaystyle \sum _{l\in {\epsilon }_h}}{\left|{\left[\sigma \left(u_{h,}{p}_h\right)n\right]}_l\right|}_{0,l}^2\right)}^{\frac{1}{2}}\hfill \\ \hfill & \hfill \\ \hfill & +C{\left({\displaystyle \sum _{l\in {\epsilon }_h}}{\left|{\left[\left(\alpha u_h+\sigma \left(u_{h,}{p}_h\right)-g\right)n\right]}_l\right|}_{0,l}^2\right)}^{\frac{1}{2}}.\hfill \end{array}$$

(53)

Finally, from (47), (50), and (51) we obtain the following result:

Theorem 13. 

There exists a constant $C$ (independent of $h$) such that

$$\Vert u-u_h\Vert +\Vert p-p_h\Vert \le C{\left({\displaystyle \sum _{K\in T_h}}\eta {\left(K\right)}^2\right)}^{\frac{1}{2}}$$

(54)

where

$$\begin{array}{ll}\eta \left(K\right)\hfill & =h_K^2{\left|A\left(u_h\right)+\nabla p_h-f\right|}_{0,K}^2+{\displaystyle \sum _{l\in {\epsilon }_h\cap K}}{\left|{\left[\sigma \left(u_{h,}{p}_h\right)n\right]}_l\right|}_{0,l}^2\hfill \\ \hfill & +{\displaystyle \sum _{l\in {\epsilon }_{\mathsf{\Gamma }}\cap K}}{\left|{\left[\left(\alpha u_h+\sigma \left(u_{h,}{p}_h\right)-g\right)n\right]}_l\right|}_{0,l}^2+{\left|\mathrm{div}{u}_h\right|}_{0,K}^2.\hfill \end{array}$$

(55)

Proof of Theorem 13. 

A direct consequence of these inequalities (43) to (53). □

Remark. 

For more complex models where the dependence of the viscosity with respect to the second invariant is more strongly nonlinear (for example Carreau law with ${\mu }_{\infty }=0$ or power law), we should use both methods of Section 2 and Section 3 and a nonlinear version of Theorem 3.3, in [43].

These systems are written as a big matrix (his component is nonlinear). In these simulations, we used the GMRES (GMRES is a generalized minimal residual algorithm applied to solve nonsymmetric linear systems, see [31,41,42] for more details.) algorithm to accelerate simulation.

4. Numerical Simulations

To conclude this paper, and in order to see the performance of the finite element method for the nonlinear equations, two different numerical simulations are represented: the first one uses the finite element software package “IFISS toolbox “ (IFISS software library is an algorithm executed under MATLAB for the interactive numerical study of differential equations for incompressible flow problems) to solve the Navier–Stokes equations with different rectangular discretization sizes of meshes ($16\times 16$, $32\times 32$, $64\times 64$, $128\times 128$, $256\times 256$) and multiple elements ($Q_1–Q_1,Q_1–P_0,Q_2–Q_1$ and $Q_2–P_1$); while the second one, uses the engineering simulation software “COMSOL Multiphysics software (COMSOL Multiphysics is a cross-platform finite element analysis, solver, and multi-physics simulation software. It allows conventional physics-based user interfaces and coupled systems of partial differential equations (PDEs)“. Inspired by the model defined in this paper [44] (applied in a porous media), we consider a “nonlinear Brinkman equation” with an inhomogeneous boundary condition, where the change of the parameter represents the change of the velocity field $u$ and the pressure $p$ in a different figure.

4.1. Fist Experience

We consider a Poiseuille (the Poiseuille flow problem is a steady horizontal flow in a channel driven by a pressure difference between the two ends) channel flow solution with an analytic solution of the Navier–Stokes equations as $u=\left(1-y^2,0\right)$ and $p=-2\nu x$, see [45,46] for more details. Where the boundary conditions are of Dirichlet or Neumann type on all the boundary—the inflow boundary is considered in the part $[x=-1,-1<y<1]$—a no-flow Dirichlet condition, $u=0$, is applied on the characteristic boundaries $y=\left\{-1,1\right\}$ and an outflow condition is considered in the rest of the boundary (i.e., $[x=1,-1<y<1]$). Figure 2, Figure 3, Figure 4 and Figure 5 represent the uniform streamline and the pressure associated with, respectively, $16\times 16$ and $256\times 256$ in cases $Q_1–Q_1,Q_1–P_0,Q_2–Q_1$ and $Q_2–P_1$ mixed approximation.

Figure 6, Figure 7, Figure 8 and Figure 9 represent the errors, associated with a $16\times 16$ and a $256\times 256$ uniform square grid with a different element, $Q_1–Q_1,Q_1–P_0,Q_2–Q_1$, and $Q_2–P_1$.

In the previous figures, we can notice that the results are well-presented if we take the mesh to be very small.

Table 1. Residual error estimator, estimated velocity divergence error, and Stokes solution residual for $16\times 16$, $32\times 32$, $64\times 64,$ $128\times 128$, and $256\times 256$.

Elements	Errors	Number of Girds
		$16\times 16$	$32\times 32$	$64\times 64$	$128\times 128$	$256\times 256$
$Q_1–Q_1$	$\eta$	5.86 × 100	1.34 × 100	9.99 × 10−1	8.38 × 10−1	2.37 × 10−1
	$\parallel \nabla \cdot {\overleftarrow{u}}_h{\parallel }_{0,\mathsf{\Omega }}$	8.77 × 10−2	2.22 × 10−2	5.65 × 10−3	1.47 × 10−3	3.92 × 10−4
	${\eta }_S$	4.22 × 100	1.75 × 100	9.43 × 10−1	5.58 × 10−1	3.26 × 10−1
$Q_1–P_0$	$\eta$	5.98 × 100	1.90 × 100	1.28 × 100	9.93 × 10−1	2.47 × 10−1
	$\parallel \nabla \cdot {\overleftarrow{u}}_h{\parallel }_{0,\mathsf{\Omega }}$	1.29 × 10−1	4.09 × 10−2	1.35 × 10−2	7.39 × 10−3	1.00 × 10−3
	${\eta }_S$	1.11 × 101	3.90 × 100	1.45 × 100	1.65 × 100	3.51 × 10−1
$Q_2–Q_1$	$\eta$	5.36 × 100	2.05 × 100	6.09 × 10−1	2.25 × 10−1	1.02 × 10−1
	$\parallel \nabla \cdot {\overleftarrow{u}}_h{\parallel }_{0,\mathsf{\Omega }}$	1.05 × 10−1	3.23 × 10−2	6.06 × 10−3	1.31 × 10−3	3.2 × 10−4
	${\eta }_S$	3.24 × 100	2.14 × 100	1.30 × 100	7.61 × 10−1	4.29 × 10−1
$Q_2–P_1$	$\eta$	5.36 × 100	1.67 × 100	4.92 × 10−1	1.55 × 10−1	6.14 × 10−2
	$\parallel \nabla \cdot {\overleftarrow{u}}_h{\parallel }_{0,\mathsf{\Omega }}$	1.07 × 10−1	1.28 × 100	5.25 × 10−3	1.23 × 10−3	3.03 × 10−4
	${\eta }_S$	2.70 × 100	6.09 × 10−1	1.21 × 100	7.18 × 10−1	4.09 × 10−1

represents the residual error estimator $\eta$, the estimated velocity divergence error $\parallel \nabla \cdot {\overleftarrow{u}}_h{\parallel }_{0,\mathsf{\Omega }}$, and the Stokes solution residual ${\eta }_S$ for difference discretization ($16\times 16$, $32\times 32$, $64\times 64$, $128\times 128$, and $256\times 256$) in these elements, cases $Q_1–Q_1,Q_1–P_0,Q_2–Q_1$ and $Q_2–P_1$ mixed approximations.

A good way to explore the capabilities of this method is to see the convergence time.

Table 2. Solve time of finite element methods for $16\times 16$, $32\times 32$, $64\times 64,$ $128\times 128$, and $256\times 256$.

Elements	$16\times 16$	$32\times 32$	$64\times 64$	$128\times 128$	$256\times 256$
$Q_1–Q_1$	1.71 × 10−1 s	2.08 × 10−1 s	8.75 × 10−1 s	3.78 × 100 s	1.99 × 101 s
$Q_1–P_0$	8.24 × 10−2 s	3.58 × 10−1 s	3.69 × 10−1 s	5.53 × 100 s	2.42 × 101 s
$Q_2–Q_1$	3.19 × 10−2 s	1.88 × 10−1 s	1.30 × 100 s	6.55 × 100 s	2.66 × 101 s
$Q_2–P_1$	1.69 × 10−1 s	1.88 × 100 s	5.33 × 100 s	3.25 × 101 s	5.66 × 101 s

presents the solve time of finite element methods with different discretization and elements.

4.2. Second Experience

The second experience comes from the second example (pore–scale flow experiments) defined by Sirivithayapakorn and Keller in [44] to model the transport of colloids in saturated porous media, and by Auset and Keller in [47] to control the dispersion of colloids. Inspired by the same model examples, and considering the new boundary conditions, the domain covers [0.640 µm × 0.320 µm], see Figure 10. The flow in the pores does not penetrate the solid grains and the inlet fluid pressure is known as $p=0.715Pa$, see Figure 10, and assumes that the boundary is changed at the outlet. Top and bottom boundaries are modeled by $\alpha \left(u\right)u+\left(-pI+K\right)\cdot n=f$, where the change of the parameter $\alpha \left(u\right)$ represents the change of the velocity field and the pressure in the figures. The primary zone of interest is the rectangular region with an upper-left corner at $\left(0,0\right)$ and lower-right coordinates at [581.6 µm × 265.0 µm].

Instead of solving for the creeping flow in the channels, the incompressible, stationary Brinkman equations are used, with the Stokes—Brinkman assumptions used next. Figure 11 shows the reservoir and localization of the boundary conditions. Figure 11 and Figure 12 represent the velocity and pressure of the model defined for these simulations. We use the finite element method to approach the unknown functions (pressure and velocity), change the value of $\mu$ in the boundary condition to see the comportment of the fluid in the reservoir, and take the two cases $\mu >>1$ and $\mu <<1$.

In this simulation, we can see the difference between Figure 11 and Figure 12: the distribution of the lines is modified when the parameter is changed. For this, in modeling the linear or nonlinear problems, we can count on this boundary condition, i.e., we can treat complex boundary conditions. This method, “the finite element method”, remains valid in our case to solve these types of problems.

5. Conclusions & Perspectives

In this paper, we studied two nonlinear differential equations, the “p-Laplacian” problem and the “Quasi-Newtonian Stokes” problem. This model is applied in many domains, for example in fluid mechanics, and we approach these models with a general boundary condition using the finite element method (FEM). For the theoretical study, we introduced the a posteriori error indicator to control the errors. The performance of this method is presented via different numerical simulations. From these perspectives, we can apply the $P1/P1$ Bubble method to approach this model for the study of the element, see [48] for a linear Brinkman model. Another problem that we can treat is to couple this model (quasi-Newtonian Stokes problem (35)) with Darcy’s law, see [49] for how the authors use the finite element methods to approach the coupled Darcy—Stokes problem. Another model which we can choose for permeability is a tensor in the equation and in the boundary condition for the quasi-Newtonian Stokes problem in (35) and (36) (i.e., we can choose a more complex model), see for example [50,51,52,53] applied to the linear model.