Numerical Analysis and Comparison of Four Stabilized Finite Element Methods for the Steady Micropolar Equations

Abstract

In this paper, four stabilized methods based on the lowest equal-order finite element pair for the steady micropolar Navier–Stokes equations (MNSE) are presented, which are penalty, regular, multiscale enrichment, and local Gauss integration methods. A priori properties, existence, uniqueness, stability, and error estimation based on Fem approximation of all the methods are proven for the physical variables. Finally, some numerical examples are displayed to show the numerical characteristics of these methods.

1. Introduction

The MNSE is a coupling system of the conservation of mass, the conservation of linear momentum, and the conservation of angular momentum. Over the past few decades, micropolar fluids are frequently used in chemistry, physics, mechanical engineering, and medicine with the development of engineering applications. For instance, it can be used in modern lubrication theory, porous media theory, liquid crystals, suspensions, and animal blood. In this paper, let $\mathsf{\Omega }$ be a bounded domain in ${\mathbb{R}}^n,n=2,3,$ with a sufficiently smooth boundary $\partial \mathsf{\Omega }$, the following MNSE in a dimensionless form will be considered:

$$\left\{\begin{array}{cc}-{\nu }_1\mathsf{\Delta }u+\nabla p=2{\nu }_r\mathrm{rot}\omega +f,\hfill & \mathrm{in}\mathsf{\Omega },\hfill \\ \mathrm{div}u=0,\hfill & \mathrm{in}\mathsf{\Omega },\hfill \\ -{\nu }_2\mathsf{\Delta }\omega +4{\nu }_r\omega =2{\nu }_r\mathrm{rot}u+g,\hfill & \mathrm{in}\mathsf{\Omega },\hfill \end{array}\right.$$

(1)

where the fluid variables $u,\omega ,p$ are the linear velocity, angular velocity, and pressure, respectively. The symbols $\nu ,{\nu }_r,c_a,c_d$ are used to represent some definitely given physical parameters, where ${\nu }_1=\nu +{\nu }_r,{\nu }_2=c_a+c_d$ [1]. The external forces f and g are predefined. When appropriate boundary conditions are supplied to (1), the equations are well-posed. For the simplicity, we consider the following homogeneous boundary conditions:

$$u=0,\omega =0,\mathrm{on}\partial \mathsf{\Omega }.$$

(2)

There are many relevant results about mathematical analysis of the problem, such as the existence, uniqueness of the solution, regularity, and so on see reference [1]. In this paper, we mainly focus on numerical methods and numerical simulations of problem (1). Noting that the Galerkin variational problem of problem (1) is still a saddle-point problem, so from the viewpoint of theoretical analysis and numerical simulation, the variables, velocity, and pressure must satisfy the LBB condition either in discrete version or continuous version. The correspondingly mixed finite element spaces for these methods must be carefully chosen so that they satisfy the LBB condition. Although some stable pairs of finite elements have been studied and used widely for many years [2,3], the lowest-order finite element pairs with some supplied stabilized terms that do not meet the LBB condition also perform well. The relevant literature about stabilized mixed finite element methods for the Navier–Stokes equations are abundant, the reader can refer to [4,5,6,7] and the corresponding literature. However, the stabilized methods of the MNSE are not currently available in the literature, so it is of great significance to study these methods about the micropolar fluid.

The analysis of solutions of the MNSE has the following typical difficulties, such as incompressibility, nonlinearity, strong coupling, and multi-field coupling. Based on the above difficulties, direct numerical simulation of the MNSE will lead to a large-scale nonlinear discrete system. Therefore, it is necessary to design an efficient, accurate, and unconditionally stable numerical algorithm for the MNSE.

In this paper, the following stabilization methods, penalty method, regular method, multiscale enrichment method, and local Gauss integration method, are mainly considered [8]. Stabilized finite element methods usually have the desirable property of improving numerical stability of the standard Galerkin methods while maintaining accuracy. The stability of regular and local Gauss integration methods is achieved by introducing symmetric definite or semidefinite stabilized terms [9,10]. In [11,12,13,14], the multiscale enrichment method stabilizes the $p_1-p_1$ pair on both finite elements and their boundaries for the Stokes equations by using a multiscale approach. In [15], variational multiscale methods combined with artificial compressibility are presented for the nonstationary Navier–Stokes equations. In this paper, the above stable finite element methods will be employed for the steady micropolar flow based on the lowest order pair. Furthermore, a numerical comparison between these methods will be presented. In the case of nonlinear and nonhomogeneous boundaries, we also do some related examples, and the related theoretical results are discussed in another paper.

The rest of this paper is organized as follows: In Section 2, some notation and preliminary results for the stationary micropolar equations are introduced. Then several stabilized mixed finite element methods and their key stabilization techniques are presented. Stability and error estimates of these stabilized finite element solutions are derived in Section 3. Comparisons between these stabilized methods are performed numerically in Section 4. Finally, conclusions are stated in Section 5. Hereafter, $c,c_i(i=0,1,\dots )$ is used to indicate a generic constant, which may represent different values in different situations.

2. Problem Statement

In this section, we introduce some notations and the well-posedness of the weak solution for continuous and discrete variational formulations of problem (1) and (2). The norm in the standard Sobolev $H^m\left(\mathsf{\Omega }\right)$, $m=0,1,2,\dots ,$ is denoted by ${\parallel ·\parallel }_m$. Specially, if $m=0$, the space $H^0\left(\mathsf{\Omega }\right)$ is the general Hilbert $L^2\left(\mathsf{\Omega }\right)$, endowed with $L^2\left(\mathsf{\Omega }\right)$-scalar product $(·,·)$ and norm ${\parallel ·\parallel }_0$. The Sobolev space $H_0^1\left(\mathsf{\Omega }\right)$ is the subspace of $H^1\left(\mathsf{\Omega }\right)$ with homogeneous boundary condition. For the convenience of analysis, the following Sobolev spaces are introduced,

$$\begin{array}{cc}\hfill X& =H_0^1{\left(\mathsf{\Omega }\right)}^n=\{v\in H^1{\left(\mathsf{\Omega }\right)}^n:v=0\mathrm{on}\partial \mathsf{\Omega }\},\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill M& =L_0^2\left(\mathsf{\Omega }\right)=\{q\in L^2\left(\mathsf{\Omega }\right):{\int }_{\mathsf{\Omega }}q\mathrm{d}\mathrm{x}=0\}.\hfill \end{array}$$

(4)

From the Poincaré inequality, we know that $\forall v\in X$, ${\parallel \nabla v\parallel }_0$ and ${\parallel v\parallel }_1$ are equivalent norm in $H^1{\left(\mathsf{\Omega }\right)}^n$. And $H^{-1}{\left(\mathsf{\Omega }\right)}^n$ denotes the dual space of X with the norm:

$$\begin{array}{c}\hfill {\parallel f\parallel }_{-1}=\underset{0\ne v\in X}{sup}\frac{|\langle f,v\rangle |}{{\parallel v\parallel }_1},\end{array}$$

where $\langle ·,·\rangle$ denotes duality product. Next, let introduce the following integration by parts formula for the rot operator,

$$(\mathrm{rot}\omega ,u)=(\omega ,\mathrm{rot}u),\forall u,\omega \in X.$$

(5)

The weak formulation for the steady micropolar Equation (1) reads: Find $(u,\omega ,p)\in X\times X\times M$, such that for all $(v,s,q)\in X\times X\times M$,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{a}_1(u,v)-d(v,p)+d(u,q)=2{\nu }_r(\mathrm{rot}\omega ,v)+(f,v),\hfill \\ a_2(\omega ,s)+4{\nu }_r(\omega ,s)=2{\nu }_r(\mathrm{rot}u,s)+(g,s),\hfill \end{array}\right.\end{array}$$

(6)

where

$$\begin{array}{cc}\hfill a_1(u,v)& ={\nu }_1(\nabla u,\nabla v),\hfill \\ \hfill a_2(\omega ,s)& ={\nu }_2(\nabla \omega ,\nabla s),\hfill \\ \hfill d(v,p)& =(p,\mathrm{div}v).\hfill \end{array}$$

In addition, Young’s inequality will be frequently used in our analysis below:

$$\begin{array}{c}\hfill ab\le \frac{ϵ}{p}{a}^p+\frac{{ϵ}^{-\frac{p}{q}}}{q}{b}^q,a,b,p,q,ϵ\in {\mathbb{R}}^{+},\frac{1}{p}+\frac{1}{q}=1,p,q\in (1,\infty ).\end{array}$$

(7)

In order to obtain the existence and uniqueness of the weak solution of the problem (6), we should introduce the following LBB condition [2].

$$\underset{q\in M}{inf}\underset{v\in X}{sup}\frac{d(v,q)}{\parallel v\parallel \parallel q\parallel }\ge \beta ,$$

(8)

where $\beta >0$ is a constant. Based on the general theory of saddle-point problem, the variational problem (6) is well-posed, and the following theorem holds [1,2,16]:

Theorem 1.

Let $f\in L^2{\left(\mathsf{\Omega }\right)}^n,g\in L^2{\left(\mathsf{\Omega }\right)}^n$, then there exist a unique solution $(u,\omega ,p)\in X\times X\times M$ of the problem (6), which satisfies

$$\begin{array}{c}\hfill {\parallel u\parallel }_1+{\parallel \omega \parallel }_1+{\parallel p\parallel }_0\le {c(\parallel f\parallel }_{-1}+{\parallel g\parallel }_{-1}).\end{array}$$

(9)

Moreover, if Ω is regular of class $C^2$, then the following regularity result holds [17]:

$${\parallel u\parallel }_2+{\parallel \omega \parallel }_2+{\parallel p\parallel }_1\le {c(\parallel f\parallel }_0+{\parallel g\parallel }_0).$$

(10)

Next, we consider the discrete problem. Let ${\mathcal{T}}_h$ be a regular triangulation of $\mathsf{\Omega }$ made up of triangles K. Let $h_K=\mathrm{diam}\left\{K\right\},h=\max \{h_K:K\in {\mathcal{T}}_h\}$. Associated with ${\mathcal{T}}_h$, the finite element spaces $(X_h,M_h)$ are defined:

$$\begin{array}{cc}\hfill X_h& =\{v\in C^0{\left(\overline{\mathsf{\Omega }}\right)}^n\cap X:v{|}_K\in {\left(P_1\left(K\right)\right)}^n,\forall K\in {\mathcal{T}}_h\},\hfill \\ \hfill M_h& =\{q\in C^0\left(\overline{\mathsf{\Omega }}\right)\cap M:q{|}_K\in P_1\left(K\right),\forall K\in {\mathcal{T}}_h\}.\hfill \end{array}$$

Then, the Galerkin finite element formulation of the problem (6) is to seek $(u_h,{\omega }_h,p_h)\in (X_h,X_h,M_h)$, $\forall (v_h,s_h,q_h)\in (X_h,X_h,M_h)$ such that

$$\left\{\begin{array}{c}{a}_1(u_h,v_h)-d(v_h,p_h)+d(u_h,q_h)=2{\nu }_r(\mathrm{rot}{\omega }_h,v_h)+(f,v_h),\hfill \\ a_2({\omega }_h,s_h)+4{\nu }_r({\omega }_h,s_h)=2{\nu }_r(\mathrm{rot}{u}_h,s_h)+(g,s_h).\hfill \end{array}\right.$$

(11)

From the classical finite element theory [2,18,19], the following assumption is reasonable

Assumption 1.

There exist interpolation operators $I_{h}v\in X_h$ such that

$$\begin{array}{cc}\hfill & \parallel v-I_h{v\parallel }_0+h\parallel v-I_h{v\parallel }_1\le c{h}^2{\parallel v\parallel }_2,\hfill \end{array}$$

(12)

Furthermore, there exist the projection operator $r_h:M\to M_h$, for any given $p\in M$,

$$\begin{array}{c}\hfill (p-r_{h}p,q)=0,\forall q\in M_h,\end{array}$$

(13)

then the following inequality holds:

$$\begin{array}{c}\hfill \parallel p-r_h{p\parallel }_0\le ch{\parallel p\parallel }_1.\end{array}$$

(14)

Now, we will consider the local interpolation operator.

$$\begin{array}{c}\hfill \parallel v - I_h{v\parallel }_{0,K}+h_K\parallel v - I_h{v\parallel }_{1,K}+h_K^2\parallel v - I_h{v\parallel }_{2,K}\le c{h}_K^2{\parallel v\parallel }_{2,K},\forall v\in H^n\left(K\right),\\ \hfill \parallel v - I_h{v\parallel }_{0,Z}+h_Z\parallel v - I_h{v\parallel }_{1,Z}\le c{h}_Z^{\frac{3}{2}}{\parallel v\parallel }_{2,\widehat{K}},\forall v\in H^n(\widehat{K}),\end{array}$$

(15)

for $\forall K\in {\mathcal{T}}_h,Z\in {\mathsf{\Gamma }}_{kj}$, where $\widehat{K}=\cup \{K\in {\mathcal{T}}_h,Z\subset \partial K\}$. For pressure interpolation, we will use the Clément interpolation operator [2,20], ${\mathcal{C}}_h:H^1\left(\mathsf{\Omega }\right)\to M_h$ satisfying

$$\begin{array}{c}\hfill \parallel q - {\mathcal{C}}_h{\left(q\right)\parallel }_{t,\mathsf{\Omega }}\le c{h}^{1-t}{\parallel q\parallel }_{1,\mathsf{\Omega }},\forall q\in H^1\left(\mathsf{\Omega }\right),\end{array}$$

(16)

for $t=0,1$.

3. Stabilized Mixed Element Methods

Noticing that the variational problem (6) is still a saddle-point problem, u and p are restricted by the LBB condition. In order to avoid this restriction, we introduce four stabilized methods to convert the saddle-point problem into an elliptical problem.

3.1. Penalty Method

The penalty method for the micropolar problem (6) is defined as: find $(u_{\epsilon },{\omega }_{\epsilon },p_{\epsilon })\in (X,X,M)$, such that $\forall (v,s,q)\in (X,X,M)$,

$$\left\{\begin{array}{c}{a}_1(u_{\epsilon },v)-d(v,p_{\epsilon })=2{\nu }_r(\mathrm{rot}{\omega }_{\epsilon },v)+\langle f,v\rangle ,\hfill \\ d(u_{\epsilon },q)+\epsilon (p_{\epsilon },q)=0,\hfill \\ a_2({\omega }_{\epsilon },s)+4{\nu }_r({\omega }_{\epsilon },s)=2{\nu }_r(\mathrm{rot}{u}_{\epsilon },s)+\langle g,s\rangle .\hfill \end{array}\right.$$

(17)

or equivalent

$$\begin{array}{c}\hfill B((u_{\epsilon },{\omega }_{\epsilon },p_{\epsilon }),(v,s,q))=F(v,s,q),\end{array}$$

(18)

where

$$\begin{array}{cc}\hfill B((u_{\epsilon },{\omega }_{\epsilon },p_{\epsilon }),(v,s,q))& =a_1(u_{\epsilon },v)-d(v,p_{\epsilon })-2{\nu }_r(\mathrm{rot}{\omega }_{\epsilon },v)\hfill \\ \hfill +& d(u_{\epsilon },q)+\epsilon (p_{\epsilon },q)+a_2({\omega }_{\epsilon },s)+4{\nu }_r({\omega }_{\epsilon },s)-2{\nu }_r(\mathrm{rot}{u}_{\epsilon },s),\hfill \end{array}$$

(19)

$$\begin{array}{c}\hfill F(v,s,q)=\langle f,v\rangle +\langle g,s\rangle .\end{array}$$

(20)

Theorem 2.

Suppose $f,g\in H^{-1}{\left(\mathsf{\Omega }\right)}^n$, then there exists a unique solution $(u_{\epsilon },{\omega }_{\epsilon },p_{\epsilon })\in (X,X,M)$ of the problem (17) which satisfies

$$\begin{array}{c}\hfill \parallel u_{\epsilon }{\parallel }_1+ \parallel {\omega }_{\epsilon }{\parallel }_1+ \parallel p_{\epsilon }{\parallel }_0 \le {c(\parallel f\parallel }_{-1}+{\parallel g\parallel }_{-1}).\end{array}$$

(21)

Proof. 

Obviously, $B((u_{\epsilon },{\omega }_{\epsilon },p_{\epsilon }),(v,s,q))$ is continuous and coercive, then by using the Lax–Milgram theorem, we get the existence and uniqueness of solution. Next, taking $v=u_{\epsilon },s={\omega }_{\epsilon },q=p_{\epsilon }$ in problem (17) to get

$$\begin{array}{cc}\hfill & {\nu }_1\parallel u_{\epsilon }{\parallel }_1^2+ {\nu }_2\parallel {\omega }_{\epsilon }{\parallel }_1^2+ 4{\nu }_r\parallel {\omega }_{\epsilon }{\parallel }_0^2- 4{\nu }_r\parallel u_{\epsilon }{\parallel }_1\parallel {\omega }_{\epsilon }{\parallel }_0+ \epsilon {\parallel p_{\epsilon }\parallel }_0^2\hfill \\ \hfill & \ge \nu \parallel u_{\epsilon }{\parallel }_1^2+ {\nu }_2\parallel {\omega }_{\epsilon }{\parallel }_1^2 \ge \min \{\nu ,{\nu }_2\}(\parallel u_{\epsilon }{\parallel }_1^2+ \parallel {\omega }_{\epsilon }{\parallel }_1^2).\hfill \end{array}$$

$$\begin{array}{c}\hfill {\nu }_1\parallel u_{\epsilon }{\parallel }_1^2+ \epsilon \parallel p_{\epsilon }{\parallel }_0^2 \le 2{\nu }_r\parallel u_{\epsilon }{\parallel }_1\parallel {\omega }_{\epsilon }{\parallel }_0+ {\parallel f\parallel }_{-1}{\parallel u_{\epsilon }\parallel }_1,\end{array}$$

(22)

$$\begin{array}{c}\hfill {\nu }_2\parallel {\omega }_{\epsilon }{\parallel }_1^2+ 4{\nu }_r\parallel {\omega }_{\epsilon }{\parallel }_0^2 \le 2{\nu }_r\parallel u_{\epsilon }{\parallel }_1\parallel {\omega }_{\epsilon }{\parallel }_0+ {\parallel g\parallel }_{-1}{\parallel {\omega }_{\epsilon }\parallel }_1.\end{array}$$

(23)

Adding (22) and (23), we get

$$\begin{array}{cc}\hfill & \min \{\nu ,{\nu }_2\}(\parallel u_{\epsilon }{\parallel }_1^2+ \parallel {\omega }_{\epsilon }{\parallel }_1^2) + \epsilon \parallel p_{\epsilon }{\parallel }_0^2 \le {\parallel f\parallel }_{-1}\parallel u_{\epsilon }{\parallel }_1+ {\parallel g\parallel }_{-1}{\parallel {\omega }_{\epsilon }\parallel }_1\hfill \\ \hfill \le & \frac{1}{2}\min \{\nu ,{\nu }_2\}\parallel u_{\epsilon }{\parallel }_1^2+ \frac{1}{2}\min {\{\nu ,{\nu }_2\}}^{-1}{\parallel f\parallel }_{-1}^2+\frac{1}{2}\min \{\nu ,{\nu }_2\}\parallel {\omega }_{\epsilon }{\parallel }_1^2+ \frac{1}{2}\min {\{\nu ,{\nu }_2\}}^{-1}{\parallel g\parallel }_{-1}^2,\hfill \end{array}$$

thus, we obtain

$$\begin{array}{c}\hfill \min \{\nu ,{\nu }_2\}(\parallel u_{\epsilon }{\parallel }_1^2+ \parallel {\omega }_{\epsilon }{\parallel }_1^2) + 2\epsilon \parallel p_{\epsilon }{\parallel }_0^2 \le \min {\{\nu ,{\nu }_2\}}^{-1}{(\parallel f\parallel }_{-1}^2+{\parallel g\parallel }_{-1}^2).\end{array}$$

Then,

$$\begin{array}{c}\hfill \parallel u_{\epsilon }{\parallel }_1+ \parallel {\omega }_{\epsilon }{\parallel }_1 \le \min {\{\nu ,{\nu }_2\}}^{-\frac{1}{2}}{(\parallel f\parallel }_{-1}+{\parallel g\parallel }_{-1}).\end{array}$$

On the other hand, by the LBB condition, taking $q=0$ in (17),

$$\begin{array}{c}\hfill \beta \parallel p_{\epsilon }{\parallel }_0 \le \max \{{\nu }_1,2{\nu }_r\}(\parallel u_{\epsilon }{\parallel }_1+ \parallel {\omega }_{\epsilon }{\parallel }_1{) +\parallel f\parallel }_{-1}\le {c(\parallel f\parallel }_{-1}+{\parallel g\parallel }_{-1}).\end{array}$$

□

Let consider the relationship between the penalty problem and the solution of the variational problem.

Theorem 3.

Suppose $(u,\omega ,p)$, $(u_{\epsilon },{\omega }_{\epsilon },p_{\epsilon })$ are the solutions of (6) and (17) respectively, then there exist

$$\begin{array}{c}\hfill \parallel u - u_{\epsilon }{\parallel }_1+\parallel \omega - {\omega }_{\epsilon }{\parallel }_1+ {\parallel p-p_{\epsilon }\parallel }_0\le c\sqrt{\epsilon }.\end{array}$$

(24)

Proof. 

Subtracting (17) from (11), we get

$$\begin{array}{c}\hfill a_1(u-u_{\epsilon },v)-d(v,p-p_{\epsilon })+d(u-u_{\epsilon },q)-\epsilon (p_{\epsilon },q)=2{\nu }_r(\mathrm{rot}(\omega -{\omega }_{\epsilon }),v),\\ \hfill a_2(\omega -{\omega }_{\epsilon },s)+4{\nu }_r(\omega -{\omega }_{\epsilon },s)=2{\nu }_r(\mathrm{rot}(u-u_{\epsilon }),s).\end{array}$$

(25)

Taking $v=u-u_{\epsilon }$, $s=\omega -{\omega }_{\epsilon }$, $q=p-p_{\epsilon }$ in (25),

$$\begin{array}{cc}\hfill & a_1(u-u_{\epsilon },u-u_{\epsilon })+a_2(\omega -{\omega }_{\epsilon },\omega -{\omega }_{\epsilon })+4{\nu }_r(\omega -{\omega }_{\epsilon },\omega -{\omega }_{\epsilon })\hfill \\ \hfill & -4{\nu }_r(\omega -{\omega }_{\epsilon },\mathrm{rot}(u-u_{\epsilon }))-\epsilon (p_{\epsilon },p-p_{\epsilon })=0,\hfill \end{array}$$

(26)

Furthermore, we have

$$\begin{array}{cc}\hfill & a_1(u-u_{\epsilon },u-u_{\epsilon })+a_2(\omega -{\omega }_{\epsilon },\omega -{\omega }_{\epsilon })+4{\nu }_r(\omega -{\omega }_{\epsilon },\omega -{\omega }_{\epsilon })\hfill \\ \hfill & -4{\nu }_r(\omega -{\omega }_{\epsilon },\mathrm{rot}(u-u_{\epsilon }))+\epsilon (p_{\epsilon },p_{\epsilon })=\epsilon (p_{\epsilon },p),\hfill \end{array}$$

Thus

$$\begin{array}{c}\hfill \nu \parallel u - u_{\epsilon }{\parallel }_1^2+ {\nu }_2\parallel \omega - {\omega }_{\epsilon }{\parallel }_1^2 \le \frac{\epsilon }{2}{\parallel p\parallel }_0^2+\frac{\epsilon }{2}{\parallel p_{\epsilon }\parallel }_0^2.\end{array}$$

By using (9) and (21),

$$\begin{array}{c}\hfill \nu \parallel u - u_{\epsilon }{\parallel }_1^2+ {\nu }_2\parallel \omega - {\omega }_{\epsilon }{ \parallel }_1^2\le {\epsilon (\parallel f\parallel }_{-1}+{\parallel g\parallel }_{-1}{)}^2.\end{array}$$

On the other hand, taking $q=0$ in (25). For the penalty method, $d(v,q)$ still satisfies the LBB conditions, an estimate of the pressure can be obtained

$$\begin{array}{cc}\hfill \beta \parallel p - p_{ϵ}{\parallel }_0& \le \underset{v\in X}{sup}\frac{(\nabla ·v,\parallel p-p_{ϵ}\parallel )}{{\parallel v\parallel }_1}\hfill \\ \hfill & \le \underset{v\in X}{sup}\frac{|a_1(u-u_{\epsilon },v)| + |2{\nu }_r(\mathrm{rot}(\omega -{\omega }_{\epsilon }),v)|}{{\parallel v\parallel }_1}\hfill \\ \hfill & \le {\nu }_1\parallel u - u_{ϵ}{\parallel }_1+ 2\sqrt{2}{\nu }_r{\parallel \omega - {\omega }_{\epsilon }\parallel }_1\hfill \\ \hfill & \le \max \{{\nu }_1,2\sqrt{2}{\nu }_r\}(\parallel u-u_{ϵ}{\parallel }_1+\parallel \omega - {\omega }_{\epsilon }{\parallel }_1),\hfill \end{array}$$

which implies that

$$\begin{array}{c}\hfill \parallel p - p_{ϵ}{\parallel }_0 \le \frac{\max \{{\nu }_1,2\sqrt{2}{\nu }_r\}}{\beta }(\parallel u-u_{ϵ}{\parallel }_1+\parallel \omega - {\omega }_{\epsilon }{\parallel }_1).\end{array}$$

(27)

□

Next, let consider the corresponding discrete weak form of the penalty method (17): find $(u_{\epsilon h},{\omega }_{\epsilon h},p_{\epsilon h})\in (X_h,X_h,M_h)$ such that for all $(v,q,s)\in (X_h,X_h,M_h)$,

$$\left\{\begin{array}{c}{a}_1(u_{\epsilon h},v)-d(v,p_{\epsilon h})+d(u_{\epsilon h},q)+\epsilon (p_{\epsilon h},q)=2{\nu }_r(\mathrm{rot}{\omega }_{\epsilon h},v)+(f,v),\hfill \\ a_2({\omega }_{\epsilon h},s)+4{\nu }_r({\omega }_{\epsilon h},s)=2{\nu }_r(\mathrm{rot}{u}_{\epsilon h},s)+(g,s).\hfill \end{array}\right.$$

(28)

Theorem 4.

The problem (28) exists a unique solution $(u_{\epsilon h},{\omega }_{\epsilon h},p_{\epsilon h})\in (X_h,W_h,M_h)$ such that

$$\begin{array}{c}\hfill \parallel u_{\epsilon h}{\parallel }_1+ \parallel {\omega }_{\epsilon h}{\parallel }_1+ \sqrt{\epsilon }\parallel p_{\epsilon h}{\parallel }_0 \le {c(\parallel f\parallel }_{-1}+{\parallel g\parallel }_{-1}).\end{array}$$

(29)

Proof. 

Choosing $(v,s,q)=(u_{\epsilon h},{\omega }_{\epsilon h},p_{\epsilon h})$ in (28), we obtain

$$\begin{array}{cc}\hfill & {\nu }_1\parallel u_{\epsilon h}{\parallel }_1^2+ {\nu }_2\parallel {\omega }_{\epsilon h}{\parallel }_1^2+ \epsilon \parallel p_{\epsilon h}{\parallel }_0^2+ 4{\nu }_r{\parallel {\omega }_{\epsilon h}\parallel }_0^2\hfill \\ \hfill & \le 4{\nu }_r\parallel u_{\epsilon h}{\parallel }_1\parallel {\omega }_{\epsilon h}{\parallel }_0+ {\parallel f\parallel }_{-1}\parallel u_{\epsilon h}{\parallel }_1+ {\parallel g\parallel }_{-1}{\parallel {\omega }_{\epsilon h}\parallel }_1.\hfill \end{array}$$

By using Young’s inequality (7), we have

$$\begin{array}{cc}\hfill 4{\nu }_r\parallel u_{\epsilon h}{\parallel }_1{\parallel {\omega }_{\epsilon h}\parallel }_0\le & {\nu }_r\parallel u_{\epsilon h}{\parallel }_1^2+ 4{\nu }_r{\parallel {\omega }_{\epsilon h}\parallel }_0^2,\hfill \\ \hfill {\parallel f\parallel }_{-1}{\parallel u_{\epsilon h}\parallel }_1\le & \frac{1}{2}\min {\{\nu ,{\nu }_2\}}^{-1}{\parallel f\parallel }_{-1}^2+\frac{1}{2}\min \{\nu ,{\nu }_2\}{\parallel u_{\epsilon h}\parallel }_1^2,\hfill \\ \hfill {\parallel g\parallel }_{-1}{\parallel {\omega }_{\epsilon h}\parallel }_1\le & \frac{1}{2}\min {\{\nu ,{\nu }_2\}}^{-1}{\parallel g\parallel }_{-1}^2+\frac{1}{2}\min \{\nu ,{\nu }_2\}{\parallel {\omega }_{\epsilon h}\parallel }_1^2,\hfill \end{array}$$

Then, carrying the above inequalities to get

$$\begin{array}{c}\hfill \nu \parallel u_{\epsilon h}{\parallel }_1^2+ {\nu }_2\parallel {\omega }_{\epsilon h}{\parallel }_1^2+ \epsilon \parallel p_{\epsilon h}{\parallel }_0^2 \le {\parallel f\parallel }_{-1}\parallel u_{\epsilon h}{\parallel }_1+ {\parallel g\parallel }_{-1}{\parallel {\omega }_{\epsilon h}\parallel }_1,\end{array}$$

(30)

and hence,

$$\begin{array}{cc}\hfill & \min \{\nu ,{\nu }_2\}(\parallel u_{\epsilon h}{\parallel }_1^2+ \parallel {\omega }_{\epsilon h}{\parallel }_1^2) + \epsilon \parallel p_{\epsilon h}{\parallel }_0^2\hfill \\ \hfill \le & \frac{1}{2}\min {\{\nu ,{\nu }_2\}}^{-1}{\parallel f\parallel }_{-1}^2+\frac{1}{2}\min \{\nu ,{\nu }_2\}\parallel u_{\epsilon h}{\parallel }_1^2+ \frac{1}{2}\min {\{\nu ,{\nu }_2\}}^{-1}{\parallel g\parallel }_{-1}^2\hfill \\ \hfill & +\frac{1}{2}\min \{\nu ,{\nu }_2\}{\parallel {\omega }_{\epsilon h}\parallel }_1^2,\hfill \end{array}$$

thus,

$$\begin{array}{c}\hfill \min \{\nu ,{\nu }_2\}(\parallel u_{\epsilon h}{\parallel }_1^2+ \parallel {\omega }_{\epsilon h}{\parallel }_1^2) + 2\epsilon \parallel p_{\epsilon h}{\parallel }_0^2 \le \min {\{\nu ,{\nu }_2\}}^{-1}{(\parallel f\parallel }_{-1}^2+{\parallel g\parallel }_{-1}^2),\end{array}$$

the above inequality implies that

$$\begin{array}{cc}\hfill \min \{\nu ,{\nu }_2\}{\parallel u_{\epsilon h}\parallel }_1^2\le & \min {\{\nu ,{\nu }_2\}}^{-1}{(\parallel f\parallel }_{-1}^2+{\parallel g\parallel }_{-1}^2),\hfill \\ \hfill \parallel u_{\epsilon h}{\parallel }_1\le & \min {\{\nu ,{\nu }_2\}}^{-1}{(\parallel f\parallel }_{-1}+{\parallel g\parallel }_{-1});\hfill \\ \hfill \min \{\nu ,{\nu }_2\}{\parallel {\omega }_{\epsilon h}\parallel }_1^2\le & \min {\{\nu ,{\nu }_2\}}^{-1}{(\parallel f\parallel }_{-1}^2+{\parallel g\parallel }_{-1}^2),\hfill \\ \hfill \parallel {\omega }_{\epsilon h}{\parallel }_1\le & \min {\{\nu ,{\nu }_2\}}^{-1}{(\parallel f\parallel }_{-1}+{\parallel g\parallel }_{-1}).\hfill \end{array}$$

On the other hand, we obtain from (30)

$$\begin{array}{c}\hfill \sqrt{\epsilon }\parallel p_{\epsilon h}{\parallel }_0 \le {c(\parallel f\parallel }_{-1}+{\parallel g\parallel }_{-1}).\end{array}$$

□

Theorem 5.

Let $(u_{\epsilon },{\omega }_{\epsilon },p_{\epsilon })$ and $(u_{\epsilon h},p_{\epsilon h},{\omega }_{\epsilon h})$ be the solution of (17) and (28), respectively. Then the error satisfies

$$\parallel u_{\epsilon }-u_{\epsilon h}{\parallel }_1+ \parallel {\omega }_{\epsilon }-{\omega }_{\epsilon h}{\parallel }_1+ \sqrt{\epsilon }{\parallel p_{\epsilon }-p_{\epsilon h}\parallel }_0\le c(h+\frac{h}{\sqrt{\epsilon }}),$$

(31)

Proof. 

Subtracting (28) from (17) yields

$$\begin{array}{cc}\hfill a_1(u_{\epsilon }-u_{\epsilon h},v_h)-d(v_h,p_{\epsilon }-p_{\epsilon h})+d(u_{\epsilon }-u_{\epsilon h},q)+\epsilon (p_{\epsilon }-p_{\epsilon h},q)=& 2{\nu }_r(\mathrm{rot}({\omega }_{\epsilon }-{\omega }_{\epsilon h}),v_h),\hfill \\ \hfill a_2({\omega }_{\epsilon }-{\omega }_{\epsilon h},s_h)+4{\nu }_r(({\omega }_{\epsilon }-{\omega }_{\epsilon h}),s)=& 2{\nu }_r(\mathrm{rot}({\omega }_{\epsilon }-{\omega }_{\epsilon h}),s_h).\hfill \end{array}$$

(32)

Denoting $(e,\theta ,\eta )=(I_h{u}_{\epsilon }-u_{\epsilon h},I_h{\omega }_{\epsilon }-{\omega }_{\epsilon h},r_h{p}_{\epsilon }-p_{\epsilon h})$ with $(e,\theta ,\eta )=(v,s,q)$, we have

$$\begin{array}{cc}\hfill {\nu }_1{\parallel e\parallel }_1^2+\epsilon {\parallel \eta \parallel }_0^2\le & |{\nu }_1(\nabla (u_{\epsilon }-I_h{u}_{\epsilon }),\nabla e)| + 2{\nu }_r(\mathrm{rot}\theta ,e)+\left|2{\nu }_r(\mathrm{rot}({\omega }_{\epsilon }-I_h{\omega }_{\epsilon }),e)\right|\hfill \\ \hfill & + |d(e,p_{\epsilon }-r_h{p}_{\epsilon h})| + |d(u_{\epsilon }-I_h{u}_{\epsilon h},\eta )| + \epsilon |(p_{\epsilon }-r_h{p}_{\epsilon },\eta )|,\hfill \\ \hfill {\nu }_2{\parallel \theta \parallel }_1^2+4{\nu }_r{\parallel \theta \parallel }_0^2\le & |{\nu }_2(\nabla ({\omega }_{\epsilon }-I_h{\omega }_{\epsilon }),\theta )| + |4{\nu }_r({\omega }_{\epsilon }-I_h{\omega }_{\epsilon }),\theta \left)\right|+2{\nu }_r(\mathrm{rot}e,\theta )\hfill \\ \hfill & + 2{\nu }_r(\mathrm{rot}(u_{\epsilon }-I_h{u}_{\epsilon }),\theta ).\hfill \end{array}$$

Adding the above two formulas together and simplifying, we get

$$\begin{array}{cc}\hfill {\nu \parallel e\parallel }_1^2+{\nu }_2{\parallel \theta \parallel }_1^2+\epsilon {\parallel \eta \parallel }_0^2\le & |{\nu }_1(\nabla (u_{\epsilon }-I_h{u}_{\epsilon }),\nabla e)| + |2{\nu }_r(\mathrm{rot}({\omega }_{\epsilon }-I_h{\omega }_{\epsilon }),e)+d(e,p_{\epsilon }-r_h{p}_{\epsilon h})\hfill \\ \hfill & + |d(u_{\epsilon }-I_h{u}_{\epsilon h},\eta )| + \epsilon |(p_{\epsilon }-r_h{p}_{\epsilon },\eta )| + {\nu }_2\left|(\nabla ({\omega }_{\epsilon }-I_h{\omega }_{\epsilon }),\theta )\right|\hfill \\ \hfill & + 4{\nu }_r(\mathrm{rot}e,\theta )+2{\nu }_r(\mathrm{rot}(u_{\epsilon }-I_h{u}_{\epsilon }),\theta )\hfill \\ \hfill \le & {ch(\parallel e\parallel }_1+{\parallel \theta \parallel }_1+{ch\parallel \eta \parallel }_0\left)\right(\parallel u_{\epsilon }-I_h{u}_{\epsilon }{\parallel }_1+ \parallel {\omega }_{\epsilon }-I_h{\omega }_{\epsilon }{\parallel }_1+ \parallel p_{\epsilon }-r_h{p}_{\epsilon }{\parallel }_0),\hfill \end{array}$$

Now, by using Hölder and Young’s inequalities, we obtain

$$\begin{array}{c}\hfill {\parallel e\parallel }_1+{\parallel \theta \parallel }_1+\sqrt{\epsilon }{\parallel \eta \parallel }_0\le c\frac{h}{\sqrt{\epsilon }}.\end{array}$$

Applying the triangle inequality to gain

$$\begin{array}{c}\hfill \parallel u_{\epsilon }-u_{\epsilon h}{\parallel }_1+ \parallel {\omega }_{\epsilon }-{\omega }_{\epsilon h}{\parallel }_1+ \sqrt{\epsilon }{\parallel p_{\epsilon }-p_{\epsilon h}\parallel }_0\le c(h+\frac{h}{\sqrt{\epsilon }}).\end{array}$$

□

Theorem 6.

Under the Theorems 3 and 5, the errors $\parallel u - u_{\epsilon h}\parallel ,\parallel \omega -{\omega }_{\epsilon h}\parallel$ and $\parallel p - p_{\epsilon h}\parallel$ satisfy

$$\begin{array}{c}\hfill \parallel u - u_{\epsilon h}{\parallel }_1+\parallel \omega - {\omega }_{\epsilon h}{\parallel }_1+ \sqrt{\epsilon }{\parallel p-p_{\epsilon h}\parallel }_0\le c(\epsilon +h+\frac{h}{\sqrt{\epsilon }}).\end{array}$$

Proof. 

The result based on (24) and (31), we get the above inequality. □

3.2. Regular Method

The variational problem (6) of the regular method is: find $(u_{Rh},{\omega }_{Rh},p_{Rh})\in (X_h,X_h,M_h)$ such that

$$\begin{array}{c}\hfill B_2((u_{Rh},{\omega }_{Rh},p_{Rh}),(v,s,q))=F(v,s,q),\forall (v,s,q)\in (X_h,X_h,M_h),\end{array}$$

(33)

where

$$\begin{array}{cc}\hfill B_2((u_{Rh},{\omega }_{Rh},p_{Rh}),(v,s,q))=& a_1(u_{Rh},v)+a_2({\omega }_{Rh},s)+4{\nu }_r({\omega }_{Rh},s)-(\mathrm{div}v,p_{Rh})\hfill \\ \hfill & +(\mathrm{div}{u}_{Rh},q)-2{\nu }_r(\mathrm{rot}{\omega }_{Rh},v)-2{\nu }_r(\mathrm{rot}{u}_{Rh},s)\hfill \\ \hfill & +{\tau }_K\sum _{K\in {\mathcal{T}}_h}{(\nabla p_{Rh}-2{\nu }_r\mathrm{rot}{\omega }_{Rh},\nabla q)}_K.\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill F(v,s,q)& =(f,v)+(g,s)+\sum _{K\in {\mathcal{T}}_h}(f,\nabla q),\hfill \end{array}$$

where the stabilized parameter ${\tau }_K={\beta }_K\frac{h_K^2}{{\nu }_1}$. For the pressure variable $p\in M$, we define the mesh-dependent norm

$${\parallel p\parallel }_h={\left\{{\tau }_K\sum _{K\in {\mathcal{T}}_h}{\parallel p\parallel }_{1,K}^2\right\}}^{\frac{1}{2}}.$$

(35)

Next, let prove the following continuous and elliptical properties for bilinear form $B_2((u_{Gh},{\omega }_{Gh},p_{Gh}),(v,s,q))$.

Lemma 1.

The bilinear form $B_2((u_{Rh},{\omega }_{Rh},p_{Rh}),(v,s,q))$, $\forall (v,s,q)\in X_h\times W_h\times M_h$ satisfies the continuous property:

$$\begin{array}{c}\hfill B_2((u_{Rh},{\omega }_{Rh},p_{Rh}),(v,s,q))\le c(\parallel u_{Rh}{\parallel }_1+ \parallel {\omega }_{Rh}{\parallel }_1+ \parallel p_{Rh}{\parallel }_0{\left)\right(\parallel v\parallel }_1+{\parallel s\parallel }_1+{\parallel q\parallel }_0).\end{array}$$

(36)

and the elliptical property:

$$\begin{array}{cc}\hfill & B_2((u_{Rh},{\omega }_{Rh},p_{Rh}),(u_{Rh},{\omega }_{Rh},p_{Rh}))\ge c(\parallel u_{Rh}{\parallel }_{1,K}^2+ \parallel {\omega }_{Rh}{\parallel }_{1,K}^2+ \parallel p_{Rh}{\parallel }_h^2).\hfill \end{array}$$

(37)

$\forall (u_{Rh},{\omega }_{Rh},p_{Rh})\in X_h\times W_h\times M_h$.

Proof. 

The continuous property is obvious. Let first consider the restriction of $B_2(·,·)$ on each element K, by using Schwarz’s inequality,

$$\begin{array}{cc}\hfill & B_2{((u_{Rh},{\omega }_{Rh},p_{Rh}),(u_{Rh},{\omega }_{Rh},p_{Rh}))}_K\hfill \\ \hfill =& {\nu }_1\parallel u_{Rh}{\parallel }_{1,K}^2+ {\nu }_2\parallel {\omega }_{Rh}{\parallel }_{1,K}^2+ 4{\nu }_r\parallel {\omega }_{Rh}{\parallel }_{0,K}^2+ {\tau }_K{\parallel p_{Rh}\parallel }_{1,K}^2\hfill \\ \hfill & - 2{\nu }_r{(\mathrm{rot}{\omega }_{Rh},u_{Rh})}_K-2{\nu }_r{(\mathrm{rot}{u}_{Rh},{\omega }_{Rh})}_K\hfill \\ \hfill \ge & {\nu }_1\parallel u_{Rh}{\parallel }_{1,K}^2+ {\nu }_2\parallel {\omega }_{Rh}{\parallel }_{1,K}^2+ 4{\nu }_r\parallel {\omega }_{Rh}{\parallel }_{0,K}^2+ {\tau }_K{\parallel p_{Rh}\parallel }_{1,K}^2\hfill \\ \hfill & - 4{\nu }_r{({\omega }_{Rh},\mathrm{rot}{u}_{Rh})}_K.\hfill \end{array}$$

(38)

Applying Young’s inequality with $ϵ=\frac{1}{2}$, we can get the following estimate of (38)

$$\begin{array}{cc}\hfill & B_2{((u_{Rh},{\omega }_{Rh},p_{Rh}),(u_{Rh},{\omega }_{Rh},p_{Rh}))}_K\hfill \\ \hfill \ge & ({\nu }_1-2{\nu }_rϵ)\parallel u_{Rh}{\parallel }_{1,K}^2+ {\nu }_2\parallel {\omega }_{Rh}{\parallel }_{1,K}^2+(4{\nu }_r-\frac{2{\nu }_r}{ϵ})\parallel {\omega }_{Rh}{\parallel }_{0,K}^2+ {\tau }_K{\parallel p_{Rh}\parallel }_{1,K}^2\hfill \\ \hfill \ge & \min \{\nu ,{\nu }_2\}(\parallel u_{Rh}{\parallel }_{1,K}^2+ \parallel {\omega }_{Rh}{\parallel }_{1,K}^2) + {\tau }_K{\parallel p_{Rh}\parallel }_{1,K}^2,\hfill \end{array}$$

on each K, and the proof is finished by adding $K\in T_h$. □

We introduce the following approximation properties [4].

Lemma 2.

Let $(v,s,q)\in {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^n\times {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^n\times [H^1\left(\mathsf{\Omega }\right)\cap L_0^2\left(\mathsf{\Omega }\right)]$, there exists an interpolation $I_{h}v$ such that

$$\begin{array}{c}\hfill \left|\right|v - I_h{v\left|\right|}_1^2+\sum _{K\in {\mathcal{T}}_h}{\tau }_K^{-1}\parallel v - I_h{v\parallel }_{0,K}^2\le c{h}^2{\parallel v\parallel }_{2,\mathsf{\Omega }}^2.\end{array}$$

(39)

Theorem 7.

Let $(u,\omega ,p)\in {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^n\times {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^n\times [H^1\left(\mathsf{\Omega }\right)\cap L_0^2\left(\mathsf{\Omega }\right)]$ be the solution of (1), then the finite element solution $(u_{Rh},{\omega }_{Rh},p_{Rh})$ to (33) satisfies the following error estimate

$$\begin{array}{c}\hfill \parallel u - u_{Rh}{\parallel }_{1,\mathsf{\Omega }}+\parallel \omega - {\omega }_{Rh}{\parallel }_{1,\mathsf{\Omega }}+\parallel p - p_{Rh}{\parallel }_h \le {ch(\parallel u\parallel }_{2,\mathsf{\Omega }}+{\parallel \omega \parallel }_{2,\mathsf{\Omega }}+{\parallel p\parallel }_{1,\mathsf{\Omega }})\end{array}$$

(40)

Proof. 

Let

$$\begin{array}{cccccc}\hfill e& =u_{Rh}-I_{h}u,\hfill & \hfill \eta & =p_{Rh}-r_{h}p,\hfill & \hfill \theta & ={\omega }_{Rh}-I_h\omega ,\hfill \\ \hfill e_h& =u-I_{h}u,\hfill & \hfill {\eta }_h& =p-r_{h}p,\hfill & \hfill {\theta }_h& =\omega -I_h\omega ,\hfill \end{array}$$

(41)

Then, from Lemma 1, using the definition of $B_2$, we obtain

$$\begin{array}{cc}\hfill & B_2((e,\theta ,\eta ),(e,\theta ,\eta )=B_2((e_h,{\theta }_h,{\eta }_h),(e,\theta ,\eta ))\hfill \\ \hfill =& {\nu }_1{(\nabla e_h,\nabla e)}_{\mathsf{\Omega }}+{\nu }_2{(\nabla {\theta }_h,\nabla \theta )}_{\mathsf{\Omega }}+4{\nu }_r{({\theta }_h,\theta )}_{\mathsf{\Omega }}-{(\mathrm{div}e,{\eta }_h)}_{\mathsf{\Omega }}\hfill \\ \hfill & + {(\mathrm{div}{e}_h,\eta )}_{\mathsf{\Omega }}-2{\nu }_r{(\mathrm{rot}{\theta }_h,e)}_{\mathsf{\Omega }}-2{\nu }_r{(\mathrm{rot}{e}_h,\theta )}_{\mathsf{\Omega }}\hfill \\ \hfill & + {\tau }_K\sum _{K\in {\mathcal{T}}_h}{(\nabla {\eta }_h-2{\nu }_r\mathrm{rot}{\theta }_h,\nabla \eta )}_K,\hfill \end{array}$$

but, since $e_h$ vanishes on $\partial \mathsf{\Omega }$ (since u belongs to $H_0^1{\left(\mathsf{\Omega }\right)}^n$), after integration by parts

$$\begin{array}{cc}\hfill & B_2((e,\theta ,\eta ),(e,\theta ,\eta )\hfill \\ \hfill =& {\nu }_1{(\nabla e_h,\nabla e)}_{\mathsf{\Omega }}+{\nu }_2{(\nabla {\theta }_h,\nabla \theta )}_{\mathsf{\Omega }}+4{\nu }_r{({\theta }_h,\theta )}_{\mathsf{\Omega }}-2{\nu }_r{(\mathrm{rot}{\theta }_h,e)}_{\mathsf{\Omega }}-2{\nu }_r{(\mathrm{rot}{e}_h,\theta )}_{\mathsf{\Omega }}\hfill \\ \hfill & - {(\mathrm{div}e,{\eta }_h)}_{\mathsf{\Omega }}-{(e_h,\nabla \eta )}_{\mathsf{\Omega }}+{\tau }_K\sum _{K\in {\mathcal{T}}_h}{(\nabla {\eta }_h,\nabla \eta )}_K\hfill \\ \hfill \le & [\sum _{K\in {\mathcal{T}}_h}{\nu }_1\parallel e_h{\parallel }_{1,K}^2+ {\nu }_2\parallel {\theta }_h{\parallel }_{1,K}^2+ 4{\nu }_r\parallel {\theta }_h{\parallel }_{0,K}^2+ 2{\nu }_r\parallel {\theta }_h{\parallel }_{1,K}^2+ 2{\nu }_r{\parallel e_h\parallel }_{1,K}^2\hfill \\ \hfill & + \parallel {\eta }_h{\parallel }_{0,K}^2+ \frac{1}{{\tau }_K}\parallel e_h{\parallel }_{0,K}^2+ {\tau }_K{\parallel {\eta }_h\parallel }_{1,K}{]}^{\frac{1}{2}}\hfill \\ \hfill & · [\sum _{K\in {\mathcal{T}}_h}{\nu }_1{\parallel e\parallel }_{1,K}^2+ {\nu }_2{\parallel \theta \parallel }_{1,K}^2+ 4{\nu }_r{\parallel \theta \parallel }_{0,K}^2+2{\nu }_r{\parallel e\parallel }_{0,K}^2+2{\nu }_r{\parallel \theta \parallel }_{0,K}^2\hfill \\ \hfill & + {\parallel \mathrm{div}e\parallel }_{0,K}^2+{\tau }_K{\parallel \eta \parallel }_{1,K}^2+{\tau }_K{\parallel \eta \parallel }_{1,K}^2{]}^{\frac{1}{2}}\hfill \\ \hfill \le & c[\sum _{K\in {\mathcal{T}}_h}({\nu }_1+2{\nu }_r)\parallel e_h{\parallel }_{1,K}^2+ \frac{1}{{\tau }_K}\parallel e_h{\parallel }_{0,K}^2+({\nu }_2+2{\nu }_r)\parallel {\theta }_h{\parallel }_{1,K}^2+ 4{\nu }_r{\parallel {\theta }_h\parallel }_{0,K}^2\hfill \\ \hfill & + \parallel {\eta }_h{\parallel }_{0,K}^2+ {\tau }_K{\parallel {\eta }_h\parallel }_{1,K}^2{]}^{\frac{1}{2}}\hfill \\ \hfill & · {\left[\sum _{K\in {\mathcal{T}}_h}({\nu }_1+1){\parallel e\parallel }_{1,K}^2+2{\nu }_r{\parallel e\parallel }_{0,K}^2+{\nu }_2{\parallel \theta \parallel }_{1,K}^2+6{\nu }_r{\parallel \theta \parallel }_{0,K}^2+{\tau }_K{\parallel \eta \parallel }_{1,K}^2\right]}^{\frac{1}{2}},\hfill \end{array}$$

from the above inequalities,

$$\begin{array}{cc}\hfill & \min \{\nu ,{\nu }_2\}\left({\parallel e\parallel }_1^2+{\parallel \theta \parallel }_1^2+{\parallel \eta \parallel }_h^2\right)\hfill \\ \hfill \le & c[\sum _{K\in {\mathcal{T}}_h}({\nu }_1+2{\nu }_r)\parallel e_h{\parallel }_{1,K}^2+ \frac{1}{{\tau }_K}\parallel e_h{\parallel }_{0,K}^2+({\nu }_2+2{\nu }_r)\parallel {\theta }_h{\parallel }_{1,K}^2+ 4{\nu }_r{\parallel {\theta }_h\parallel }_{0,K}^2\hfill \\ \hfill & + \parallel {\eta }_h{\parallel }_{0,K}^2+ {\tau }_K{\parallel {\eta }_h\parallel }_{1,K}^2{]}^{\frac{1}{2}}\hfill \\ \hfill & · {\left[\sum _{K\in {\mathcal{T}}_h}({\nu }_1+1){\parallel e\parallel }_{1,K}^2+2{\nu }_r{\parallel e\parallel }_{0,K}^2+{\nu }_2{\parallel \theta \parallel }_{1,K}^2+6{\nu }_r{\parallel \theta \parallel }_{0,K}^2+{\tau }_K{\parallel \eta \parallel }_{1,K}^2\right]}^{\frac{1}{2}}\hfill \\ \hfill \le & c\sqrt{\max \{{\nu }_1+1,2{\nu }_r,{\nu }_2,6{\nu }_r\}}[\sum _{K\in {\mathcal{T}}_h}({\nu }_1+2{\nu }_r)\parallel e_h{\parallel }_{1,K}^2+ \frac{1}{{\tau }_K}{\parallel e_h\parallel }_{0,K}^2\hfill \\ \hfill & +({\nu }_2+2{\nu }_r)\parallel {\theta }_h{\parallel }_{1,K}^2+ 4{\nu }_r\parallel {\theta }_h{\parallel }_{0,K}^2+ \parallel {\eta }_h{\parallel }_{0,K}^2+ {\tau }_K{\parallel {\eta }_h\parallel }_{1,K}^2{]}^{\frac{1}{2}}\hfill \\ \hfill & · {\left[{\parallel e\parallel }_{1,\mathsf{\Omega }}^2+{\parallel \theta \parallel }_{1,\mathsf{\Omega }}^2+{\parallel \eta \parallel }_{h,\mathsf{\Omega }}^2\right]}^{\frac{1}{2}},\hfill \end{array}$$

and hence, dividing by the last term we get

$$\begin{array}{cc}\hfill & {\parallel e\parallel }_1+{\parallel \theta \parallel }_1+{\parallel \eta \parallel }_h\hfill \\ \hfill \le & c\frac{\sqrt{\max \{{\nu }_1+1,2{\nu }_r,c_1,6{\nu }_r\}}}{\min \{\nu ,{\nu }_2\}}[\sum _{K\in {\mathcal{T}}_h}({\nu }_1+2{\nu }_r)\parallel e_h{\parallel }_{1,K}^2+ \frac{1}{{\tau }_K}{\parallel e_h\parallel }_{0,K}^2\hfill \\ \hfill & +({\nu }_2+2{\nu }_r)\parallel {\theta }_h{\parallel }_{1,K}^2+ 4{\nu }_r\parallel {\theta }_h{\parallel }_{0,K}^2+ \parallel {\eta }_h{\parallel }_{0,K}^2+ {\tau }_K{\parallel {\eta }_h\parallel }_{1,K}^2{]}^{\frac{1}{2}}\hfill \\ \hfill \le & ch(\parallel u{\parallel }_{2,\mathsf{\Omega }}+\parallel \omega {\parallel }_{2,\mathsf{\Omega }}+\parallel p{\parallel }_{1,\mathsf{\Omega }})\hfill \end{array}$$

Finally, since $u-u_h=e_h-e,\omega -{\omega }_h={\theta }_h-\theta$ and $p-p_h={\eta }_h-\eta$, with the triangle inequality yields

$$\begin{array}{c}\hfill \parallel u - u_h{\parallel }_{1,\mathsf{\Omega }}+\parallel \omega - {\omega }_h{\parallel }_{1,\mathsf{\Omega }}+\parallel p - p_h{\parallel }_h \le {ch(\parallel u\parallel }_{2,\mathsf{\Omega }}+{\parallel \omega \parallel }_{2,\mathsf{\Omega }}+{\parallel p\parallel }_{1,\mathsf{\Omega }}).\end{array}$$

□

An Improved Error Estimate

Noting that the above estimation of pressure is still depend on the mesh parameter h, we can modified it into mesh independent. The main idea of proof can refer to the similar result of [9,21].

Theorem 8.

Let $(u,\omega ,p)\in {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^n\times {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^n\times [H^1\left(\mathsf{\Omega }\right)\cap L_0^2\left(\mathsf{\Omega }\right)]$ be the solution of (1). Then, the error satisfies

$$\begin{array}{c}\hfill \parallel u - u_{Rh}{\parallel }_{1,\mathsf{\Omega }}+\parallel \omega - {\omega }_{Rh}{\parallel }_{1,\mathsf{\Omega }}+\parallel p - p_{Rh}{\parallel }_{0,\mathsf{\Omega }} \le {ch(\parallel u\parallel }_{2,\mathsf{\Omega }}+{\parallel \omega \parallel }_{2,\mathsf{\Omega }}+{\parallel p\parallel }_{1,\mathsf{\Omega }}).\end{array}$$

Proof. 

Noting that, if we choose ${\tau }_Z=0$ in the multiscale enrichment method of the next subsection, then we can recover the regular method from the multiscale enrichment method. The proof is omitted. □

3.3. Multiscale Enrichment Method

Find $(u_{Mh},{\omega }_{Mh},p_{Mh})\in (X_h,W_h,M_h)$, such that

$$B_3((u_{Mh},{\omega }_{Mh},p_{Mh}),(v,s,q))=F(v,s,q),$$

(42)

where

$$\begin{array}{cc}\hfill & B_3((u_{Mh},{\omega }_{Mh},p_{Mh}),(v,s,q))\hfill \\ \hfill =& a_1(u_{Mh},v)+a_2({\omega }_{Mh},s)+4{\nu }_r({\omega }_{Mh},s)-(\mathrm{div}v,p_{Mh})+(\mathrm{div}{u}_{Mh},q)\hfill \\ \hfill & - 2{\nu }_r(\mathrm{rot}{\omega }_{Mh},v)-2{\nu }_r(\mathrm{rot}{u}_{Mh},s)\hfill \\ \hfill & +\sum _{K\in {\mathcal{T}}_h}{\tau }_K{(\nabla p_{Mh}-2{\nu }_r\mathrm{rot}{\omega }_{Mh},\nabla q)}_K\hfill \\ \hfill & +\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{\left(\left[\nu {\partial }_n{u}_{Mh}\right],\left[\nu {\partial }_{n}v\right]\right)}_Z,\hfill \end{array}$$

(43)

$$F(v,s,q)={(f,v)}_{\mathsf{\Omega }}+{(g,s)}_{\mathsf{\Omega }}+\sum _{K\in {\mathcal{T}}_h}{\tau }_K{(f,\nabla q)}_K,$$

(44)

$${\tau }_K={\beta }_K\frac{h_K^2}{{\nu }_1},{\tau }_Z={\beta }_e\frac{h_Z}{{\nu }_1},$$

(45)

where ${\tau }_K$ and ${\tau }_Z$ are the positive stabilization parameters, the quantity $h_Z=\left|Z\right|$ is the length of the edge $Z,Z\subset \partial K$, and $\left[v\right]$ denotes the jump of v across ${\mathsf{\Gamma }}_{kj}$. Define the mesh-dependent norms

$$\begin{array}{cc}\hfill {\left|\right|\left|v\right|\left|\right|}_h^2& ={\nu }_1{\parallel v\parallel }_{1,\mathsf{\Omega }}^2+\sum _{z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{\parallel \left[\nu {\partial }_{n}v\right]\parallel }_{0,Z}^2,\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill {\parallel q\parallel }_h^2& =\sum _{K\in {\mathcal{T}}_h}{\tau }_K{\parallel q\parallel }_{1,K}^2.\hfill \end{array}$$

(47)

Before analyzing the stability of the method given by (42), we introduce the following local trace theorems

$$\begin{array}{c}\hfill {\parallel v\parallel }_{0,\partial K}^2\le c(h_K^{-1}{\parallel v\parallel }_{0,K}^2+h_K{\parallel v\parallel }_{1,K}^2).\end{array}$$

(48)

Lemma 3.

The bilinear form $B_3((u_{Mh},{\omega }_{Mh},p_{Mh}),(v,s,q))$ satisfies the continuous property

$$\begin{array}{c}\hfill B_3((u_{Mh},{\omega }_{Mh},p_{Mh}),(v,s,q))\le c(\parallel u_{Mh}{\parallel }_1+ \parallel {\omega }_{Mh}{\parallel }_{+}\parallel p_{Mh}{\parallel }_0{\left)\right(\parallel v\parallel }_1+{\parallel s\parallel }_1+{\parallel q\parallel }_0).\end{array}$$

(49)

for all $(u_{Mh},{\omega }_{Mh},p_{Mh})\in (X_h,X_h,M_h)$, $(v,s,q)\in (X_h,X_h,M_h)$, and the elliptical property

$$\begin{array}{cc}\hfill & B_3((u_{Mh},{\omega }_{Mh},p_{Mh}),(u_{Mh},{\omega }_{Mh},p_{Mh}))\hfill \\ \hfill \ge & c(\parallel u_{Mh}{\parallel }_{1,\mathsf{\Omega }}^2+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z\parallel \left[{\nu }_1{\partial }_n{u}_{Mh}\right]{\parallel }_{0,Z}^2+\parallel {\omega }_{Mh}{\parallel }_{1,\mathsf{\Omega }}^2+\parallel p_{Mh}{\parallel }_h^2).\hfill \end{array}$$

(50)

Proof. 

The continuous property follows using (48) and Cauchy-Schwartz inequality.

$$\begin{array}{cc}\hfill & B_3{((u_{Mh},{\omega }_{Mh},p_{Mh}),(u_{Mh},{\omega }_{Mh},p_{Mh}))}_K\hfill \\ \hfill =& {\nu }_1\parallel u_{Mh}{\parallel }_{1,K}^2+ {\nu }_2\parallel {\omega }_{Mh}{\parallel }_{1,K}^2+4 {\nu }_r\parallel {\omega }_{Mh}{\parallel }_{0,K}^2+ {\tau }_K{\parallel p_{Mh}\parallel }_{0,K}^2\hfill \\ \hfill & - 4{\nu }_r{({\omega }_{Mh},\mathrm{rot}{u}_{Mh})}_K+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{\parallel \left[{\nu }_1{\partial }_n{u}_{Mh}\right]\parallel }_{0,Z}^2\hfill \\ \hfill \ge & {\nu }_1\parallel u_{Mh}{\parallel }_{1,K}^2+ {\nu }_2\parallel {\omega }_{Mh}{\parallel }_{1,K}^2+ 4{\nu }_r\parallel {\omega }_{Mh}{\parallel }_{0,K}^2+ {\tau }_K{\parallel p_{Mh}\parallel }_{0,K}^2\hfill \\ \hfill & - 4{\nu }_r\parallel {\omega }_{Mh}{\parallel }_{0,K}\parallel \nabla u_{Mh}{\parallel }_{0,K}+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{\parallel \left[{\nu }_1{\partial }_n{u}_{Mh}\right]\parallel }_{0,Z}^2.\hfill \end{array}$$

Since $2ab\le \frac{1}{ϵ}{a}^2+ϵb^2$ with $ϵ>0$, we see that

$$\begin{array}{cc}\hfill & B_3{((u_{Mh},{\omega }_{Mh},p_{Mh}),(u_{Mh},{\omega }_{Mh},p_{Mh}))}_K\hfill \\ \hfill \ge & ({\nu }_1-2{\nu }_rϵ)\parallel u_{Mh}{\parallel }_{1,K}^2+ {\nu }_2\parallel {\omega }_{Mh}{\parallel }_{1,K}^2+(4{\nu }_r-\frac{2{\nu }_r}{ϵ}){\parallel {\omega }_{Mh}\parallel }_{0,K}^2\hfill \\ \hfill & +\sum _{K\in {\mathcal{T}}_h}{\tau }_K\parallel p_{Mh}{\parallel }_{1,K}^2+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{\parallel \left[{\nu }_1{\partial }_n{u}_{Mh}\right]\parallel }_{0,Z}^2.\hfill \end{array}$$

Finally, taking $ϵ=\frac{1}{2}$ to obtain

$$\begin{array}{cc}\hfill & B_3{((u_{Mh},{\omega }_{Mh},p_{Mh}),(u_{Mh},{\omega }_{Mh},p_{Mh}))}_K\hfill \\ \hfill \ge & \nu \parallel u_{Mh}{\parallel }_{1,K}^2+ {\nu }_2\parallel {\omega }_{Mh}{\parallel }_{1,K}^2+\sum _{K\in {\mathcal{T}}_h}{\tau }_K\parallel p_{Mh}{\parallel }_{1,K}^2+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{\parallel \left[{\nu }_1{\partial }_n{u}_{Mh}\right]\parallel }_{0,Z}^2\hfill \\ \hfill \ge & c(\parallel u_{Mh}{\parallel }_{1,K}^2+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z\parallel \left[{\nu }_1{\partial }_n{u}_{Mh}\right]{\parallel }_{0,Z}^2+\parallel {\omega }_{Mh}{\parallel }_{1,K}^2+\parallel p_{Mh}{\parallel }_h^2).\hfill \end{array}$$

□

Next, we introduce the interpolation error about the Clément interpolation operator.

Lemma 4.

Let $(v,s,q)\in {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^2\times {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^2\times [H^1\left(\mathsf{\Omega }\right)\cap L_0^2\left(\mathsf{\Omega }\right)]$, and ${\tilde{q}}_h={\mathcal{C}}_h\left(q\right)-\frac{{({\mathcal{C}}_h\left(q\right),1)}_{\mathsf{\Omega }}}{\left|\mathsf{\Omega }\right|}$, such that

$$\begin{array}{c}\hfill \left|\right||v-I_h{v\left|\right||}_h^2 +\sum _{K\in {\mathcal{T}}_h}{\tau }_K^{-1}\parallel v - I_h{v\parallel }_{0,K}^2\le c{h}^2{\nu }_1{\parallel v\parallel }_{2,\mathsf{\Omega }}^2,\end{array}$$

(51)

$$\begin{array}{c}\hfill \parallel q - {\tilde{q}}_h{\parallel }_h+ \frac{1}{\sqrt{{\nu }_1}}\parallel q - {\tilde{q}}_h{\parallel }_{0,\mathsf{\Omega }} \le ch\frac{1}{\sqrt{{\nu }_1}}{\parallel q\parallel }_{1,\mathsf{\Omega }}.\end{array}$$

(52)

Proof. 

The result comes from the norm definition and uses $\parallel q - {\tilde{q}}_h{\parallel }_{0,\mathsf{\Omega }} \le {\parallel q-{\mathcal{C}}_h\left(q\right)\parallel }_{0,\mathsf{\Omega }}$ to combine with (15) and (16). □

Theorem 9.

Let $(u,\omega ,p)\in {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^n\times {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^n\times [H^1\left(\mathsf{\Omega }\right)\cap L_0^2\left(\mathsf{\Omega }\right)]$ be the solution of (1), and $(u_{Mh},{\omega }_{Mh},p_{Mh})$ the solution of (42). Then, the following error estimate holds

$$\begin{array}{cc}\hfill & \left|\right||u-u_{Mh}{\left|\right||}_h +\parallel \omega - {\omega }_{Mh}{\parallel }_{1,\mathsf{\Omega }}+\parallel p - p_{Mh}{\parallel }_h \le ch(\sqrt{{\nu }_1}{\parallel u\parallel }_{2,\mathsf{\Omega }}+{\parallel \omega \parallel }_{2,\mathsf{\Omega }}+\frac{1}{\sqrt{{\nu }_1}}{\parallel p\parallel }_{1,\mathsf{\Omega }}).\hfill \end{array}$$

Proof. 

For the sake of simplicity, let ${\tilde{u}}_h=I_{h}u,{\tilde{\omega }}_h=I_h\omega$, $({\eta }_h^u,{\eta }_h^{\omega },{\eta }_h^p)=(u-{\tilde{u}}_h,\omega -{\tilde{\omega }}_h,p-{\tilde{p}}_h)$. Then

$$\begin{array}{cc}\hfill & {\nu }_1\parallel u_{Mh}-{\tilde{u}}_h{\parallel }_1^2+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z\parallel {\nu }_1{\partial }_n(u_{Mh}-{\tilde{u}}_h){\parallel }_{0,Z}^2+ {\nu }_2\parallel {\omega }_{Mh}-{\tilde{\omega }}_h{\parallel }_1^2+ {\parallel p_{Mh}-{\tilde{p}}_h\parallel }_h^2\hfill \\ \hfill \le & B_3((u_{Mh}-{\tilde{u}}_h,{\omega }_{Mh}-{\tilde{\omega }}_h,p_{Mh}-{\tilde{p}}_h),(u_{Mh}-{\tilde{u}}_h,{\omega }_{Mh}-{\tilde{\omega }}_h,p_{Mh}-{\tilde{p}}_h))\hfill \\ \hfill =& B_3(({\eta }_h^u,{\eta }_h^{\omega },{\eta }_h^p),(u_{Mh}-{\tilde{u}}_h,{\omega }_{Mh}-{\tilde{\omega }}_h,p_{Mh}-{\tilde{p}}_h))\hfill \\ \hfill =& {\nu }_1{(\nabla {\eta }_h^u,\nabla (u_{Mh}-{\tilde{u}}_h))}_{\mathsf{\Omega }}+{\nu }_2{(\nabla {\eta }_h^{\omega },\nabla ({\omega }_{Mh}-{\tilde{\omega }}_h))}_{\mathsf{\Omega }}+4{\nu }_r{({\eta }_h^{\omega },{\omega }_{Mh}-{\tilde{\omega }}_h)}_{\mathsf{\Omega }}\hfill \\ \hfill & - {(\mathrm{div}(u_{Mh}-{\tilde{u}}_h),{\eta }_h^p)}_{\mathsf{\Omega }}-{({\eta }_h^u,\nabla (p_{Mh}-{\tilde{p}}_h))}_{\mathsf{\Omega }}-2{\nu }_r{(\mathrm{rot}{\eta }_h^{\omega },u_{Mh}-{\tilde{u}}_h)}_{\mathsf{\Omega }}\hfill \\ \hfill & - 2{\nu }_r{(\mathrm{rot}{\eta }_h^u,{\omega }_{Mh}-{\tilde{\omega }}_h)}_{\mathsf{\Omega }}+\sum _{K\in T_h}{\tau }_K{(\nabla {\eta }_h^p,\nabla (p_{Mh}-{\tilde{p}}_h))}_K\hfill \\ \hfill & +\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{(\left[{\nu }_1{\partial }_n{\eta }_h^u\right],\left[{\nu }_1{\partial }_n(u_{Mh}-{\tilde{u}}_h)\right])}_Z\hfill \\ \hfill \le & [{\nu }_1\parallel {\eta }_h^u{\parallel }_{1,\mathsf{\Omega }}^2+ {\nu }_2\parallel {\eta }_h^{\omega }{\parallel }_{1,\mathsf{\Omega }}^2+ 4{\nu }_r\parallel {\eta }_h^{\omega }{\parallel }_{0,\mathsf{\Omega }}^2+ 2{\nu }_r\parallel {\eta }_h^{\omega }{\parallel }_{0,\mathsf{\Omega }}^2+ 2{\nu }_r\parallel {\eta }_h^u{\parallel }_{1,\mathsf{\Omega }}^2+ {\parallel {\eta }_h^p\parallel }_{0,\mathsf{\Omega }}^2\hfill \\ \hfill & + \parallel {\eta }_h^u{\parallel }_{0,\mathsf{\Omega }}^2+ \parallel {\eta }_h^p{\parallel }_h^2+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{\parallel \left[{\nu }_1{\partial }_n{\eta }_h^u\right]\parallel }_{0,Z}^2{]}^{\frac{1}{2}}\hfill \\ \hfill & · [{\nu }_1\parallel u_{Mh}-{\tilde{u}}_h{\parallel }_{1,\mathsf{\Omega }}^2+ {\nu }_2\parallel {\omega }_{Mh}-{\tilde{\omega }}_h{\parallel }_{1,\mathsf{\Omega }}^2+ 6{\nu }_r\parallel {\omega }_{Mh}-{\tilde{\omega }}_h{\parallel }_{0,\mathsf{\Omega }}^2+ 2{\nu }_1{\parallel u_{Mh}-{\tilde{u}}_h\parallel }_{1,\mathsf{\Omega }}^2\hfill \\ \hfill & + 2\sum _{K\in T_h}{\tau }_K\parallel \nabla (p_{Mh}-{\tilde{p}}_h){\parallel }_{0,K}^2+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{\parallel \left[{\nu }_1{\partial }_n(u_{Mh}-{\tilde{u}}_h)\right]\parallel }_{0,Z}^2{]}^{\frac{1}{2}}\hfill \\ \hfill \le & c{\left[|\parallel {\eta }_h^u{\parallel |}_h^2+{\nu }_2\parallel {\eta }_h^{\omega }{\parallel }_{1,\mathsf{\Omega }}^2+ \parallel {\eta }_h^u{\parallel }_{0,\mathsf{\Omega }}^2+ 6{\nu }_r\parallel {\eta }_h^{\omega }{\parallel }_{0,\mathsf{\Omega }}^2+ \parallel {\eta }_h^p{\parallel }_{0,\mathsf{\Omega }}^2+ {\parallel {\eta }_h^p\parallel }_h^2\right]}^{\frac{1}{2}}\hfill \\ \hfill & · {\left[\left|\right||u_{Mh}-{\tilde{u}}_h{\left|\right||}_h^2+{\nu }_2\parallel {\omega }_{Mh}-{\tilde{\omega }}_h{\parallel }_{1,\mathsf{\Omega }}^2+ {\parallel p_{Mh}-{\tilde{p}}_h\parallel }_h^2\right]}^{\frac{1}{2}}.\hfill \end{array}$$

Hence, dividing by the last term we get

$$\begin{array}{cc}\hfill & \left|\right||u_{Mh}-{\tilde{u}}_h{\left|\right||}_h+{\nu }_2\parallel {\omega }_{Mh}-{\tilde{\omega }}_h{\parallel }_{1,\mathsf{\Omega }}+ {\parallel p_{Mh}-{\tilde{p}}_h\parallel }_h\hfill \\ \hfill \le & c{\left[\left|\right||{\eta }_h^u{\left|\right||}_h^2+{\nu }_2\parallel {\eta }_h^{\omega }{\parallel }_{1,\mathsf{\Omega }}^2+ \parallel {\eta }_h^u{\parallel }_{0,\mathsf{\Omega }}^2+ 6{\nu }_r\parallel {\eta }_h^{\omega }{\parallel }_{0,\mathsf{\Omega }}^2+ \parallel {\eta }_h^p{\parallel }_{0,\mathsf{\Omega }}^2+ {\parallel {\eta }_h^p\parallel }_h^2\right]}^{\frac{1}{2}}\hfill \\ \hfill \le & ch(\sqrt{{\nu }_1}\parallel u{\parallel }_{2,\mathsf{\Omega }}+\parallel \omega {\parallel }_{2,\mathsf{\Omega }}+ \frac{1}{\sqrt{{\nu }_1}}\parallel p{\parallel }_{1,\mathsf{\Omega }}).\hfill \end{array}$$

Finally, combing above inequalities with the triangular inequality gives the following result.

$$\begin{array}{c}\hfill \left|\right||u-u_{Mh}{\left|\right||}_h+{\nu }_2\parallel \omega - {\omega }_{Mh}{\parallel }_{1,\mathsf{\Omega }}+\parallel p - p_{Mh}{\parallel }_h \le ch(\sqrt{{\nu }_1}{\parallel u\parallel }_{2,\mathsf{\Omega }}+{\parallel \omega \parallel }_{2,\mathsf{\Omega }}+\frac{1}{\sqrt{{\nu }_1}}{\parallel p\parallel }_{1,\mathsf{\Omega }}).\end{array}$$

□

Remark 1.

Because of the norm definition, we cannot guarantee convergence of pressure. The next result shows that we have an independent optimal error estimate of h in the natural norm of pressure.

Theorem 10.

Let $(u,\omega ,p)\in {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^n\times {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^2\times [H^1\left(\mathsf{\Omega }\right)\cap L_0^n\left(\mathsf{\Omega }\right)]$ be the solution of (1), and $(u_{Mh},{\omega }_{Mh},p_{Mh})$ the solution of (42). Then, the following error estimate holds

$$\begin{array}{c}\hfill \parallel p - p_{Mh}{\parallel }_{0,\mathsf{\Omega }} \le ch({\nu }_1{\parallel u\parallel }_{2,\mathsf{\Omega }}+{\parallel \omega \parallel }_{2,\mathsf{\Omega }}+{\parallel p\parallel }_{1,\mathsf{\Omega }}).\end{array}$$

Proof. 

Known by the continuous inf-sup condition, there exists $\phi \in X$ such that $\nabla ·\phi =p-p_{Mh}$ and ${\parallel \phi \parallel }_{1,\mathsf{\Omega }}\le c{\parallel p-p_{Mh}\parallel }_{0,\mathsf{\Omega }}$. Let ${\phi }_h={\mathcal{C}}_h\left(\phi \right)\in X_h$, we obtain

$$\begin{array}{cc}\hfill \parallel p - p_{Mh}{\parallel }_{0,\mathsf{\Omega }}^2 =& {(\nabla ·\phi ,p-p_{Mh})}_{\mathsf{\Omega }}={(\nabla ·(\phi -{\phi }_h),p-p_{Mh})}_{\mathsf{\Omega }}+{(\nabla ·{\phi }_h,p-p_{Mh})}_{\mathsf{\Omega }}\hfill \\ \hfill \le & -\sum _{K\in T_h}{(\phi -{\phi }_h,\nabla (p-p_{Mh}))}_K+{\nu }_1{(\nabla (u-u_{Mh}),\nabla {\phi }_h)}_{\mathsf{\Omega }}\hfill \\ \hfill & - 2{\nu }_r{(\mathrm{rot}(\omega -{\omega }_{Mh}),{\phi }_h)}_{\mathsf{\Omega }}+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{(\left[{\nu }_1{\partial }_n(u-u_{Mh})\right],\left[{\nu }_1{\partial }_n{\phi }_h\right])}_{0,Z}\hfill \\ \hfill \le & \sum _{K\in T_h}\parallel \phi - {\phi }_h{\parallel }_{0,K}\parallel p - p_{Mh}{)\parallel }_{1,K}+{\nu }_1\parallel u - u_{Mh}{\parallel }_{1,\mathsf{\Omega }}{\parallel {\phi }_{Mh}\parallel }_{1,\mathsf{\Omega }}\hfill \\ \hfill & + 2{\nu }_r\parallel \omega - {\omega }_{Mh}{\parallel }_{1,\mathsf{\Omega }}\parallel {\phi }_{Mh}{\parallel }_{0,\mathsf{\Omega }}+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z\parallel \left[{\nu }_1{\partial }_n(u-u_{Mh})\right]{\parallel }_{0,Z}{\parallel \left[{\nu }_1{\partial }_n{\phi }_h\right]\parallel }_{0,Z}\hfill \\ \hfill \le & {\left[\sum _{K\in T_h}{\tau }_K^{-1}\parallel \phi - {\phi }_h{\parallel }_{0,K}^2+ {\nu }_1\parallel {\phi }_h{\parallel }_{1,\mathsf{\Omega }}^2+ 2{\nu }_r\parallel {\phi }_h{\parallel }_{1,\mathsf{\Omega }}^2+\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{\parallel \left[{\nu }_1{\partial }_n{\phi }_h\right]\parallel }_{0,Z}^2\right]}^{\frac{1}{2}}\hfill \\ \hfill & · [\sum _{K\in T_h}{\tau }_K\parallel p - p_{Mh}{)\parallel }_{1,K}^2+{\nu }_1\parallel u - u_{Mh}{\parallel }_{1,\mathsf{\Omega }}^2+ 2{\nu }_r{\parallel \omega -{\omega }_{Mh}\parallel }_{1,\mathsf{\Omega }}^2\hfill \\ \hfill & +\sum _{Z\in {\mathsf{\Gamma }}_{kj}}{\tau }_Z{\parallel \left[{\nu }_1{\partial }_n(u-u_{Mh})\right]\parallel }_{0,Z}^2{]}^{\frac{1}{2}}\hfill \\ \hfill \le & c\sqrt{{\nu }_1}\left(\right|\left|\right|u-u_{Mh}{\left|\right||}_h +\parallel \omega - {\omega }_{Mh}{\parallel }_{1,\mathsf{\Omega }}+\parallel p - p_{Mh}{)\parallel }_h{\left)\right(\parallel \phi \parallel }_{1,\mathsf{\Omega }}^2 + \parallel {\phi }_h{{\parallel }_{1,\mathsf{\Omega }}^2)}^{\frac{1}{2}}\hfill \\ \hfill \le & ch\sqrt{{\nu }_1}(\sqrt{{\nu }_1}{\parallel u\parallel }_{2,\mathsf{\Omega }}+{\parallel \omega \parallel }_{2,\mathsf{\Omega }}+\frac{1}{\sqrt{{\nu }_1}}{\parallel p\parallel }_{1,\mathsf{\Omega }})\parallel p-p_{Mh}{\parallel }_{0,\mathsf{\Omega }},\hfill \end{array}$$

Then divide by the last term to get the result. □

3.4. Local Gauss Integration Method

The center idea of this method is to add two local Gauss integrals to the original discrete formulation, seek $(u_{Gh},{\omega }_{Gh},p_{Gh})\in (X_h,X_h,M_h)$ such that

$$\left\{\begin{array}{c}{a}_1(u_{Gh},v)-d(v,p_{Gh})+d(u_{Gh},q)+G(p_{Gh},q)=2{\nu }_r(\mathrm{rot}{\omega }_{Gh},v)+(f,v),\hfill \\ a_2({\omega }_{Gh},s)+4{\nu }_r({\omega }_{Gh},s)=2{\nu }_r(\mathrm{rot}{u}_{Gh},s)+(g,s),\hfill \end{array}\right.$$

(53)

$$B_4((u_{Gh},{\omega }_{Gh},p_{Gh}),(v,s,q))=F(v,s),$$

(54)

where

$$\begin{array}{cc}\hfill & B_4((u_{Gh},{\omega }_{Gh},p_{Gh}),(v,s,q))\hfill \\ \hfill =& B((u_{Gh},{\omega }_{Gh},p_{Gh}),\left((v,s,q)\right)-G(p_{Gh},q)\hfill \\ \hfill =& {\nu }_1(\nabla u_{Gh},\nabla v)+{\nu }_2(\nabla {\omega }_{Gh},\nabla s)-(\mathrm{div}v,p_{Gh})+(\mathrm{div}{u}_{Gh},q)\hfill \\ \hfill & + 4{\nu }_r({\omega }_{Gh},s)-2{\nu }_r(\mathrm{rot}{\omega }_{Gh},v)-2{\nu }_r(\mathrm{rot}{u}_{Gh},s)-G(p_{Gh},q),\hfill \end{array}$$

(55)

$$F(v,s)=(f,v)+(g,s),$$

(56)

$G(p_{Gh},q)$ is defined by

$$G(p_{Gh},q)=(p_{Gh}-\mathsf{\Pi }{p}_{Gh},q-\mathsf{\Pi }q),$$

(57)

for all $(v,s,q)\in (X_h,X_h,M_h)$, and $\mathsf{\Pi }$ is a $L^2$ projection operator with the following properties:

$$(p,q)=(\mathsf{\Pi }p,q),\forall p\in L^2\left(\mathsf{\Omega }\right),q\in R_0,$$

(58)

$${\parallel \mathsf{\Pi }p\parallel }_0\le c{\parallel p\parallel }_0,\forall p\in L^2\left(\mathsf{\Omega }\right),$$

(59)

$${\parallel p-\mathsf{\Pi }p\parallel }_0\le ch{\parallel p\parallel }_1,\forall p\in H^1\left(\mathsf{\Omega }\right),$$

(60)

where $R_0\subset L^2\left(\mathsf{\Omega }\right)$ denotes the piecewise constant space associated with the triangulation ${\mathcal{T}}_h$.

Obviously, this method does not require a complex computation or a stabilization parameter. It is only necessary to compute the block $G_e$ with simple Gauss integrals. The stabilization term is defined as follows:

$$G(p_{Gh},q)=\sum _{K\in K_h}\left\{{\int }_{K,2}{p}_{Gh}q\mathrm{dx}-{\int }_{K,1}{p}_{Gh}q\mathrm{dx}\right\},$$

(61)

where ${\int }_{K,i}{p}_{Gh}q\mathrm{d}\mathrm{x}$ indicates a local Gauss integral over K that is exact for polynomials of degree $i,i=1,2$. Obviously, the bilinear form $G(p_{Gh},q)$ is symmetric and semi-definitely generated on each local set K.

Theorem 11.

The bilinear form $B_4((u_{Gh},{\omega }_{Gh},p_{Gh}),(v,s,q))$ satisfies the continuous property

$$\begin{array}{c}\hfill B_4((u_{Gh},{\omega }_{Gh},p_{Gh}),(v,s,q))\le c(\parallel u_{Gh}{\parallel }_1+ \parallel {\omega }_{Gh}{\parallel }_1+ \parallel p_{Gh}{\parallel }_0{\left)\right(\parallel v\parallel }_1+{\parallel s\parallel }_1+{\parallel q\parallel }_0)\end{array}$$

(62)

where $\forall (u_{Gh},{\omega }_{Gh},p_{Gh}),(v,s,q)\in (X_h,X_h,M_h)$, and the coercive property

$$\begin{array}{cc}\hfill \underset{(v,s,q)\in X_h,W_h,M_h}{sup}\frac{|B_4((u_{Gh},{\omega }_{Gh},p_{Gh}),(v,s,q))|}{{(\parallel v\parallel }_1^2+{\parallel s\parallel }_1^2+{\parallel q\parallel }_0^2{)}^{\frac{1}{2}}}\ge & \tilde{\beta }(\parallel u_{Gh}{\parallel }_1^2+ \parallel {\omega }_{Gh}{\parallel }_1^2+ \parallel p_{Gh}{{\parallel }_0^2)}^{\frac{1}{2}},\hfill \\ \hfill (u_{Gh},{\omega }_{Gh},p_{Gh})\in (X_h,X_h,M_h),\end{array}$$

(63)

where $\tilde{\beta }$ is a positive constant depending only on Ω.

Proof. 

$$\begin{array}{cc}\hfill & |B_4((u_{Gh},{\omega }_{Gh},p_{Gh}),(v,s,q))|\hfill \\ \hfill =& |{\nu }_1(\nabla u_{Gh},\nabla v)+{\nu }_2(\nabla {\omega }_{Gh},\nabla s)-d(v,p_{Gh})+d(u_{Gh},q)\hfill \\ \hfill & + 4{\nu }_r({\omega }_{Gh},s)-2{\nu }_r(\mathrm{rot}{\omega }_{Gh},v)-2{\nu }_r(\mathrm{rot}{u}_{Gh},s)-G(p_{Gh},q)|\hfill \\ \hfill \le & c(\parallel u_{Gh}{\parallel }_1\parallel v{\parallel }_1+\parallel {\omega }_{Gh}{\parallel }_1\parallel s{\parallel }_1+\parallel u_{Gh}{\parallel }_1\parallel s{\parallel }_1+\parallel {\omega }_{Gh}{\parallel }_1\parallel v{\parallel }_1\hfill \\ \hfill & + {\parallel v\parallel }_1\parallel p_{Gh}{\parallel }_0+ \parallel u_{Gh}{\parallel }_1{\parallel q\parallel }_0 + \parallel p_{Gh}{\parallel }_0{\parallel q\parallel }_0)\hfill \\ \hfill \le & c(\parallel u_{Gh}{\parallel }_1+ \parallel {\omega }_{Gh}{\parallel }_1+ \parallel p_{Gh}{\parallel }_0{\left)\right(\parallel v\parallel }_1+{\parallel s\parallel }_1+{\parallel q\parallel }_0).\hfill \end{array}$$

Thus it suffices to show the continuous property.

For the coercive property of $B_4$, $\forall p_{Gh}\in M_h$, there exists a positive constant $C_0$ and $z\in X$ such that [2]

$$\begin{array}{cc}\hfill (\mathrm{div}z,p_{Gh})& = \parallel p_{Gh}{\parallel }_0^2,\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill {\parallel z\parallel }_1& \le c_0{\parallel p_{Gh}\parallel }_0.\hfill \end{array}$$

(65)

Setting the finite element approximation $z_{Gh}\in X_h$ of z, we have

$$\begin{array}{cc}\hfill \parallel z_{Gh}{\parallel }_1& \le c_1{\parallel p_{Gh}\parallel }_0.\hfill \end{array}$$

(66)

Then, for any $p_{Gh}\in M_h$, we choose any $(v,s,q)=(u_{Gh}-\lambda z_{Gh},{\omega }_{Gh},-p_{Gh})$, where $0<\lambda <\frac{2(1-c_1)}{c_1{\nu }_r}$.

Obviously, it follows from (58)–(61), (65) and (66) and the Young inequality that

$$\begin{array}{cc}\hfill & B_4((u_{Gh},{\omega }_{Gh},p_{Gh}),(v,s,q))=B_4((u_{Gh},{\omega }_{Gh},p_{Gh}),(u_{Gh}-\lambda z_{Gh},{\omega }_{Gh},p_{Gh}))\hfill \\ \hfill =& a_1(u_{Gh},u_{Gh})-\alpha a_1(u_{Gh},z_{Gh})+a_2({\omega }_{Gh},{\omega }_{Gh})+\lambda d(z_{Gh},p_{Gh})\hfill \\ \hfill & + 4{\nu }_r({\omega }_{Gh},{\omega }_{Gh})-4{\nu }_r({\omega }_{Gh},\mathrm{rot}{u}_{Gh})+2\alpha {\nu }_r(\mathrm{rot}{\omega }_{Gh},z_{Gh})+G(p_{Gh},p_{Gh})\hfill \\ \hfill \ge & {\nu }_1\parallel u_{Gh}{\parallel }_1^2+ {\nu }_2\parallel {\omega }_{Gh}{\parallel }_1^2+ G(p_{Gh},p_{Gh})+4{\nu }_r\parallel {\omega }_{Gh}{\parallel }_0^2- \lambda {\nu }_1\parallel u_{Gh}{\parallel }_1{\parallel z_{Gh}\parallel }_1\hfill \\ \hfill & - 4{\nu }_r\parallel {\omega }_{Gh}{\parallel }_0\parallel u_{Gh}{\parallel }_1- 2{\nu }_r\parallel {\omega }_{Gh}{\parallel }_0{\parallel z_{Gh}\parallel }_1+\lambda d(z_{Gh}-z,p_{Gh})+\lambda d(z,p_{Gh})\hfill \\ \hfill \ge & \frac{\nu }{2}\parallel u_{Gh}{\parallel }_1^2+ {\nu }_2\parallel {\omega }_{Gh}{\parallel }_1^2+ G(p_{Gh},p_{Gh})+\lambda \parallel p_{Gh}{\parallel }_0^2- \frac{{\nu }_r{\alpha }^2}{2}{\parallel z_{Gh}\parallel }_1^2\hfill \\ \hfill & - c_2\lambda \parallel p_{Gh}-\mathsf{\Pi }{p}_{Gh}{\parallel }_1{\parallel p_{Gh}\parallel }_0\hfill \\ \hfill \ge & \frac{\nu }{2}\parallel u_{Gh}{\parallel }_1^2+ {\nu }_2\parallel {\omega }_{Gh}{\parallel }_1^2+(\lambda -\frac{c_3{\nu }_r{\lambda }^2}{2}-\frac{c_4\lambda }{2}){\parallel p_{Gh}\parallel }_0^2+\frac{1}{2}G(p_{Gh},p_{Gh})\hfill \\ \hfill \ge & c_5(\parallel u_{Gh}{\parallel }_1^2+ \parallel {\omega }_{Gh}{\parallel }_1^2+ \parallel p_{Gh}{\parallel }_0^2).\hfill \end{array}$$

(67)

Finally, we remark that

$$\begin{array}{cc}\hfill {\parallel v\parallel }_1+{\parallel s\parallel }_1+{\parallel q\parallel }_0& = \parallel u_{Gh}-\lambda z_{Gh}{\parallel }_1+ \parallel {\omega }_{Gh}{\parallel }_1+{\parallel p_{Gh}\parallel }_0\hfill \\ \hfill & \le c_6(\parallel u_{Gh}{\parallel }_1^2+ \parallel {\omega }_{Gh}{\parallel }_1^2+ \parallel p_{Gh}{{\parallel }_0^2)}^{\frac{1}{2}},\hfill \end{array}$$

setting $\tilde{\beta }=\frac{c_5}{c_6}$, which finishes the proof. □

Theorem 12.

Let $(u,\omega ,p)\in {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^n\times {[H^2\left(\mathsf{\Omega }\right)\cap H_0^1\left(\mathsf{\Omega }\right)]}^n\times [H^1\left(\mathsf{\Omega }\right)\cap L_0^2\left(\mathsf{\Omega }\right)]$ be the solution of (1), the stabilized finite element solution $(u_{Gh},{\omega }_{Gh},p_{Gh})$ satisfies the error estimates:

$$\parallel u - u_{Gh}{\parallel }_1+\parallel \omega - {\omega }_{Gh}{\parallel }_1+\parallel p - p_{Gh}{\parallel }_0 \le {ch(\parallel u\parallel }_2+{\parallel \omega \parallel }_2+{\parallel p\parallel }_1).$$

(68)

Proof. 

Subtracting (54) from (6), we have

$$\begin{array}{c}\hfill B_4((u-u_{Gh},\omega -{\omega }_{Gh},p-p_{Gh}),(v,s,q))=G(p,q).\end{array}$$

(69)

By setting $(e,\theta ,\eta )=(I_{h}u-u_{Gh},I_h\omega -{\omega }_{Gh},r_{h}p-p_{Gh})$, it follows from, we can see that

$$\begin{array}{cc}\hfill & B_4((e,\theta ,\eta ),(v,s,q))\hfill \\ \hfill =& G(p,q)-B_4((u-I_{h}u,\omega -I_h\omega ,p-r_{h}p),(v,s,q))\hfill \\ \hfill \le & \left|G(p,q)\right| + |B_4((u-I_{h}u,\omega -I_h\omega ,p-r_{h}p),(v,s,q))|\hfill \\ \hfill \le & {c\parallel p-\mathsf{\Pi }p\parallel }_0{\parallel q\parallel }_0+c(\parallel u-I_h{u\parallel }_1 +\parallel \omega - I_h{\omega \parallel }_1+\parallel p - r_h{p\parallel }_0{\left)\right(\parallel v\parallel }_1+{\parallel s\parallel }_1+{\parallel q\parallel }_0)\hfill \\ \hfill \le & {ch\parallel p\parallel }_1{(\parallel v\parallel }_1+{\parallel s\parallel }_1+{\parallel q\parallel }_0{) + ch(\parallel u\parallel }_2+{\parallel \omega \parallel }_2+{\parallel p\parallel }_1{\left)\right(\parallel v\parallel }_1+{\parallel s\parallel }_1+{\parallel q\parallel }_0)\hfill \\ \hfill \le & {ch(\parallel u\parallel }_2+{\parallel \omega \parallel }_2+{\parallel p\parallel }_1{\left)\right(\parallel v\parallel }_1+{\parallel s\parallel }_1+{\parallel q\parallel }_0).\hfill \end{array}$$

(70)

Obviously, it follows from (63) and (70) that

$$\begin{array}{cc}\hfill \tilde{\beta }{(\parallel e\parallel }_1^2+{\parallel \theta \parallel }_1^2+{\parallel \eta \parallel }_0^2{)}^{\frac{1}{2}}\le & \underset{0\ne (v,s,q)\in X_h\times W_h\times M_h}{sup}\frac{|B_4\left((e,\theta ,\eta )\right),(v,s,q)\left)\right|}{{(\parallel v\parallel }_1^2+{\parallel s\parallel }_1^2+{\parallel q\parallel }_0^2{)}^{\frac{1}{2}}}\hfill \\ \hfill \le & ch(\parallel u{\parallel }_2+\parallel \omega {\parallel }_2+\parallel p{\parallel }_1).\hfill \end{array}$$

(71)

Thanks to (71) and the triangle inequality

$$\begin{array}{cc}\hfill & \parallel u - u_{Gh}{\parallel }_1+\parallel \omega - {\omega }_{Gh}{\parallel }_1+ {\parallel p-p_{Gh}\parallel }_1\hfill \\ \hfill \le & \parallel u - I_h{u\parallel }_1 +\parallel \omega - I_h{\omega \parallel }_1 +\parallel p - r_h{p\parallel }_0+{\parallel e\parallel }_1+{\parallel \theta \parallel }_1+{\parallel \eta \parallel }_0\hfill \\ \hfill \le & ch(\parallel u{\parallel }_2+\parallel \omega {\parallel }_2+\parallel p{\parallel }_1).\hfill \end{array}$$

□

4. Numerical Experiments

In this section, we numerically compare the performance of the various stabilized mixed finite element methods discussed in the previous section for two examples. It is numerically solved by the stabilized mixed methods on uniform meshes (see the first picture in Figure 1). The second example solves the same problem on unstructured meshes (see the last picture in Figure 1). We recall that in our computations, the pressure and velocity are approximated by piecewise linear finite elements. All these computations have been based on the package FreeFem++, with some of our additional codes. The stabilized term of the regular must be controlled by carefully designed mesh-dependent parameters, whose optimal values are often unknown.

We now report the convergence rates for the penalty, regular, multiscale enrichment, and local Gauss integration methods for the steady micropolar equations solved by using the lowest equal-order element $p_1-p_1-p_1$ on the first uniform triangular mesh (see the first panel in Figure 1). The errors in the relative $L^2\left(\mathsf{\Omega }\right)$- and $H^1\left(\mathsf{\Omega }\right)$-norms for the velocity, angular velocity, and the relative $L^2\left(\mathsf{\Omega }\right)$-norm for the pressure and their corresponding convergence rates are given in

Table 1. Typical structured mesh with $\epsilon =O(h^{\frac{1}{2}})$ for the penalty method.

$\frac{1}{\mathit{h}}$	CPU	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{ϵ}\mathit{h}}{\parallel }_1}{{\parallel \mathit{u}\parallel }_1}$	${\mathit{u}}_{\mathit{H}1}$-Rate	$\frac{\parallel \mathit{\omega } - {\mathit{\omega }}_{\mathit{ϵ}\mathit{h}}{\parallel }_1}{{\parallel \mathit{\omega }\parallel }_1}$	${\mathit{\omega }}_{\mathit{H}1}$-Rate	$\frac{\parallel \mathit{p} - {\mathit{p}}_{\mathit{ϵ}\mathit{h}}{\parallel }_0}{{\parallel \mathit{p}\parallel }_0}$	${\mathit{p}}_{\mathit{L}2}$-Rate
12	0.071		$1.512\times {10}^0$		$8.828\times {10}^{-1}$		$9.606\times {10}^{-2}$
24	0.320	0.4089	$1.139\times {10}^0$	0.4948	$5.615\times {10}^{-1}$	0.6529	$6.817\times {10}^{-2}$
36	0.846	0.4332	$9.555\times {10}^{-1}$	0.4873	$4.405\times {10}^{-1}$	0.5985	$5.595\times {10}^{-2}$
48	1.755	0.4440	$8.409\times {10}^{-1}$	0.4841	$3.740\times {10}^{-1}$	0.5683	$4.867\times {10}^{-2}$
60	3.185	0.4505	$7.605\times {10}^{-1}$	0.4826	$3.309\times {10}^{-1}$	0.5497	$4.370\times {10}^{-2}$
72	5.332	0.4551	$6.999\times {10}^{-1}$	0.4819	$3.000\times {10}^{-1}$	0.5373	$4.003\times {10}^{-2}$

,

Table 2. Typical structured mesh with $\epsilon ={10}^{-6}$ for the penalty method.

$\frac{1}{\mathit{h}}$	CPU	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{ϵ}\mathit{h}}{\parallel }_0}{{\parallel \mathit{u}\parallel }_0}$	${\mathit{u}}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{ϵ}\mathit{h}}{\parallel }_1}{{\parallel \mathit{u}\parallel }_1}$	${\mathit{u}}_{\mathit{H}1}$-Rate	$\frac{\parallel \mathit{\omega } - {\mathit{\omega }}_{\mathit{ϵ}\mathit{h}}{\parallel }_0}{{\parallel \mathit{\omega }\parallel }_0}$	${\mathit{\omega }}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{p} - {\mathit{p}}_{\mathit{ϵ}\mathit{h}}{\parallel }_0}{{\parallel \mathit{p}\parallel }_0}$	${\mathit{p}}_{\mathit{L}2}$-Rate
12	0.046	$8.656\times {10}^{-2}$		$2.848\times {10}^{-1}$		$3.419\times {10}^{-2}$		$1.912\times {10}^{-1}$
24	0.225	$2.137\times {10}^{-2}$	2.018	$1.412\times {10}^{-1}$	1.012	$8.601\times {10}^{-3}$	1.991	$1.012\times {10}^{-1}$	0.917
36	0.668	$9.447\times {10}^{-3}$	2.014	$9.386\times {10}^{-2}$	1.008	$3.825\times {10}^{-3}$	1.999	$7.059\times {10}^{-2}$	0.890
48	1.47	$5.298\times {10}^{-3}$	2.010	$7.027\times {10}^{-2}$	1.006	$2.152\times {10}^{-3}$	2.000	$5.549\times {10}^{-2}$	0.836
60	2.806	$3.385\times {10}^{-3}$	2.008	$5.615\times {10}^{-2}$	1.005	$1.377\times {10}^{-3}$	2.001	$4.671\times {10}^{-2}$	0.772
72	4.965	$2.348\times {10}^{-3}$	2.005	$4.675\times {10}^{-2}$	1.005	$9.560\times {10}^{-4}$	2.001	$4.109\times {10}^{-2}$	0.704

,

Table 3. Typical structured mesh for the regular method.

$\frac{1}{\mathit{h}}$	CPU	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{Rh}}{\parallel }_0}{{\parallel \mathit{u}\parallel }_0}$	${\mathit{u}}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{R}\mathit{h}}{\parallel }_1}{{\parallel \mathit{u}\parallel }_1}$	${\mathit{u}}_{\mathit{H}1}$-Rate	$\frac{\parallel \mathit{\omega } - {\mathit{\omega }}_{\mathit{R}\mathit{h}}{\parallel }_0}{{\parallel \mathit{\omega }\parallel }_0}$	${\mathit{\omega }}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{p} - {\mathit{p}}_{\mathit{R}\mathit{h}}{\parallel }_0}{{\parallel \mathit{p}\parallel }_0}$	${\mathit{p}}_{\mathit{L}2}$-Rate
12	0.053	$1.641\times {10}^{-1}$		$4.454\times {10}^{-1}$		$3.225\times {10}^{-2}$		$3.045\times {10}^{-2}$
24	0.281	$4.069\times {10}^{-2}$	2.012	$1.800\times {10}^{-1}$	1.307	$8.121\times {10}^{-3}$	1.990	$9.207\times {10}^{-3}$	1.726
36	0.734	$1.797\times {10}^{-2}$	2.016	$1.077\times {10}^{-1}$	1.268	$3.613\times {10}^{-3}$	1.998	$4.675\times {10}^{-3}$	1.671
48	1.517	$1.001\times {10}^{-2}$	2.013	$7.565\times {10}^{-2}$	1.226	$2.033\times {10}^{-3}$	1.999	$2.916\times {10}^{-3}$	1.641
60	2.768	$6.430\times {10}^{-3}$	2.011	$5.794\times {10}^{-2}$	1.195	$1.301\times {10}^{-3}$	2.000	$2.032\times {10}^{-3}$	1.619
72	4.65	$4.457\times {10}^{-3}$	2.010	$4.680\times {10}^{-2}$	1.171	$9.034\times {10}^{-4}$	2.000	$1.517\times {10}^{-3}$	1.604

and

Table 4. Typical structured mesh for the multiscale enrichment method.

$\frac{1}{\mathit{h}}$	CPU	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{Mh}}{\parallel }_0}{{\parallel \mathit{u}\parallel }_0}$	${\mathit{u}}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{M}\mathit{h}}{\parallel }_1}{{\parallel \mathit{u}\parallel }_1}$	${\mathit{u}}_{\mathit{H}1}$-Rate	$\frac{\parallel \mathit{\omega } - {\mathit{\omega }}_{\mathit{Mh}}{\parallel }_0}{{\parallel \mathit{\omega }\parallel }_0}$	${\mathit{\omega }}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{p} - {\mathit{p}}_{\mathit{Mh}}{\parallel }_0}{{\parallel \mathit{p}\parallel }_0}$	${\mathit{p}}_{\mathit{L}2}$-Rate
12	0.1	$1.723\times {10}^{-1}$		$4.387\times {10}^{-1}$		$3.394\times {10}^{-2}$		$3.116\times {10}^{-2}$
24	0.529	$4.429\times {10}^{-2}$	1.960	$1.777\times {10}^{-1}$	1.304	$8.847\times {10}^{-3}$	1.940	$9.466\times {10}^{-3}$	1.719
36	1.446	$2.021\times {10}^{-2}$	1.936	$1.067\times {10}^{-1}$	1.259	$4.063\times {10}^{-3}$	1.919	$4.814\times {10}^{-3}$	1.668
48	3.112	$1.170\times {10}^{-2}$	1.900	$7.520\times {10}^{-2}$	1.216	$2.361\times {10}^{-3}$	1.888	$3.005\times {10}^{-3}$	1.638
60	5.755	$7.718\times {10}^{-3}$	1.863	$5.775\times {10}^{-2}$	1.183	$1.561\times {10}^{-3}$	1.853	$2.095\times {10}^{-3}$	1.617
72	9.652	$5.533\times {10}^{-3}$	1.825	$4.676\times {10}^{-2}$	1.158	$1.121\times {10}^{-3}$	1.817	$1.564\times {10}^{-3}$	1.602

. In addition, the theoretical convergence rates should be of order $O\left(h^2\right)$ and $O\left(h\right)$ for the velocity and angular velocity in the $L^2\left(\mathsf{\Omega }\right)$- and $H^1\left(\mathsf{\Omega }\right)$-norms, respectively, and of order $O\left(h\right)$ for the pressure in the $L^2\left(\mathsf{\Omega }\right)$-norm by using all these stabilized methods. It follows from Table 1, Table 2, Table 3 and Table 4 that the theoretical results are confirmed for the velocity. Except for the penalty method, the speed and convergence order of the other methods have reached the results of theoretical analysis.

4.1. Problems with Smooth Solutions

Setting $\mathsf{\Omega }={[0,1]}^2$, equation parameters $\nu ={\nu }_r=c_a=c_d=0.1$, we take penalty parameter $\epsilon$ = ${10}^{-6}$ and $\epsilon =O\left(h^{\frac{1}{2}}\right)$ in Table 1 and Table 2 respectively, ${\beta }_1=100$, ${\beta }_2=100$, ${\beta }_3=150$, and the exact solution $(u,\omega ,p)$ is with the right-hand side function f generated by the exact solution:

$$\begin{array}{cc}\hfill u_1(x,y)=& 10x^2{(x-1)}^{2}y(y-1)(2y-1),\hfill \\ \hfill u_2(x,y)=& -10x(x-1)(2x-1)y^2{(y-1)}^2,\hfill \\ \hfill \omega (x,y)=& 10x^2{(x-1)}^{2}y(y-1)(2y-1)-10x(x-1)(2x-1)y^2{(y-1)}^2,\hfill \\ \hfill p(x,y)=& 10(2x-1)(2y-1).\hfill \end{array}$$

It can be seen from Table 2, Table 3, Table 4 and

Table 5. Typical structured mesh for the local Gauss integration method.

$\frac{1}{\mathit{h}}$	CPU	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{Gh}}{\parallel }_0}{{\parallel \mathit{u}\parallel }_0}$	${\mathit{u}}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{G}\mathit{h}}{\parallel }_1}{{\parallel \mathit{u}\parallel }_1}$	${\mathit{u}}_{\mathit{H}1}$-Rate	$\frac{\parallel \mathit{\omega } - {\mathit{\omega }}_{\mathit{Gh}}{\parallel }_0}{{\parallel \mathit{\omega }\parallel }_0}$	${\mathit{\omega }}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{p} - {\mathit{p}}_{\mathit{Gh}}{\parallel }_0}{{\parallel \mathit{p}\parallel }_0}$	${\mathit{p}}_{\mathit{L}2}$-Rate
12	0.054	$2.127\times {10}^{-1}$		$4.612\times {10}^{-1}$		$3.313\times {10}^{-2}$		$2.909\times {10}^{-2}$
24	0.256	$5.273\times {10}^{-2}$	2.012	$1.753\times {10}^{-1}$	1.396	$8.336\times {10}^{-3}$	1.991	$8.306\times {10}^{-3}$	1.809
36	0.682	$2.335\times {10}^{-2}$	2.009	$1.036\times {10}^{-1}$	1.297	$3.707\times {10}^{-3}$	1.999	$4.103\times {10}^{-3}$	1.739
48	1.357	$1.311\times {10}^{-2}$	2.006	$7.269\times {10}^{-2}$	1.232	$2.085\times {10}^{-3}$	2.000	$2.561\times {10}^{-3}$	1.700
60	2.403	$8.385\times {10}^{-3}$	2.003	$5.575\times {10}^{-2}$	1.189	$1.334\times {10}^{-3}$	2.001	$1.733\times {10}^{-3}$	1.672
72	3.938	$5.823\times {10}^{-3}$	2.000	$4.513\times {10}^{-2}$	1.159	$9.264\times {10}^{-4}$	2.001	$1.282\times {10}^{-3}$	1.652

that, except for the multiscale enrichment method, the velocity and angular velocity convergence order of the other methods have reached the result of theoretical analysis. The velocity and angular velocity $L^2\left(\mathsf{\Omega }\right)$-norm convergence order of the multiscale enrichment method do not reach $O\left(h^2\right)$ with the refinement of the grid. The convergence speed of the $L^2\left(\mathsf{\Omega }\right)$-norm of the pressure of the regular method, the local Gauss integration method, and the multiscale enrichment method has almost reached $O\left(h^{1.5}\right)$, which are better than the expected results of the theoretical analysis.

Additionally, note the CPU time listed in Table 2, Table 3, Table 4 and Table 5, all our calculations are performed under the same computing platform. From Table 2, Table 3, Table 4 and Table 5, as the grid is encrypted, the multiscale enrichment method consumes the most time, while the local Gauss integration method consumes the least time.

In addition, take the other exact solution as follows:

$$\begin{array}{cc}\hfill u_1(x,y)& =\pi {sin}^2\left(\pi x\right)sin\left(2\pi y\right),\hfill \\ \hfill u_2(x,y)& =-\pi sin\left(2\pi x\right){sin}^2\left(\pi y\right),\hfill \\ \hfill p(x,y)& =cos\left(\pi x\right)sin\left(\pi y\right),\hfill \\ \hfill \omega (x,y)& =\pi {sin}^2\left(\pi x\right){sin}^2\left(\pi y\right).\hfill \end{array}$$

Given unstructured triangular mesh of $\mathsf{\Omega }$ with different h = 0.270282, 0.14633, 0.073319, 0.0368164, 0.0188245, the errors are listed in

Table 6. Unstructured mesh for the penalty method.

$\frac{1}{\mathit{h}}$	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{ϵ}\mathit{h}}{\parallel }_0}{{\parallel \mathit{u}\parallel }_0}$	${\mathit{u}}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{ϵ}\mathit{h}}{\parallel }_1}{{\parallel \mathit{u}\parallel }_1}$	${\mathit{u}}_{\mathit{H}1}$-Rate	$\frac{\parallel \mathit{\omega } - {\mathit{\omega }}_{\mathit{ϵ}\mathit{h}}{\parallel }_0}{{\parallel \mathit{\omega }\parallel }_0}$	${\mathit{\omega }}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{p} - {\mathit{p}}_{\mathit{ϵ}\mathit{h}}{\parallel }_0}{{\parallel \mathit{p}\parallel }_0}$	${\mathit{p}}_{\mathit{L}2}$-Rate
0.270282	$1.0559\times {10}^0$		$1.0587\times {10}^{-0}$		$1.1639\times {10}^{-1}$		$9.5639\times {10}^{-2}$
0.14633	$1.7658\times {10}^{-1}$	2.91	$6.6212\times {10}^{-1}$	0.77	$2.6803\times {10}^{-2}$	2.39	$9.9276\times {10}^{-3}$	3.69
0.073319	$2.1140\times {10}^{-2}$	3.07	$1.7191\times {10}^{-1}$	1.95	$6.6651\times {10}^{-3}$	2.01	$2.5280\times {10}^{-3}$	1.98
0.0368164	$3.3790\times {10}^{-3}$	2.66	$5.6879\times {10}^{-2}$	1.61	$1.6494\times {10}^{-3}$	2.03	$7.4832\times {10}^{-4}$	1.77
0.0188245	$7.3197\times {10}^{-4}$	2.28	$2.5374\times {10}^{-2}$	1.20	$4.1029\times {10}^{-4}$	2.07	$1.8817\times {10}^{-4}$	2.06

,

Table 7. Unstructured mesh for the regular method.

$\frac{1}{\mathit{h}}$	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{Rh}}{\parallel }_0}{{\parallel \mathit{u}\parallel }_0}$	${\mathit{u}}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{R}\mathit{h}}{\parallel }_1}{{\parallel \mathit{u}\parallel }_1}$	${\mathit{u}}_{\mathit{H}1}$-Rate	$\frac{\parallel \mathit{\omega } - {\mathit{\omega }}_{\mathit{R}\mathit{h}}{\parallel }_0}{{\parallel \mathit{\omega }\parallel }_0}$	${\mathit{\omega }}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{p} - {\mathit{p}}_{\mathit{R}\mathit{h}}{\parallel }_0}{{\parallel \mathit{p}\parallel }_0}$	${\mathit{p}}_{\mathit{L}2}$-Rate
0.270282	$1.8999\times {10}^{-1}$		$4.0094\times {10}^{-1}$		$1.4119\times {10}^{-1}$		$6.1590\times {10}^0$
0.14633	$5.0548\times {10}^{-2}$	2.16	$1.9939\times {10}^{-1}$	1.14	$3.3420\times {10}^{-2}$	2.35	$2.0315\times {10}^0$	1.81
0.073319	$1.2644\times {10}^{-2}$	2.01	$9.8450\times {10}^{-2}$	1.02	$8.3367\times {10}^{-3}$	2.01	$5.4865\times {10}^{-1}$	1.89
0.0368164	$3.1384\times {10}^{-3}$	2.02	$4.8764\times {10}^{-2}$	1.02	$2.0596\times {10}^{-3}$	2.03	$1.5511\times {10}^{-1}$	1.83
0.0188245	$7.9530\times {10}^{-3}$	2.05	$2.4324\times {10}^{-2}$	1.04	$5.1178\times {10}^{-4}$	2.08	$4.8886\times {10}^{-1}$	1.72

,

Table 8. Unstructured mesh for the multiscale enrichment method.

$\frac{1}{\mathit{h}}$	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{Mh}}{\parallel }_0}{{\parallel \mathit{u}\parallel }_0}$	${\mathit{u}}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{M}\mathit{h}}{\parallel }_1}{{\parallel \mathit{u}\parallel }_1}$	${\mathit{u}}_{\mathit{H}1}$-Rate	$\frac{\parallel \mathit{\omega } - {\mathit{\omega }}_{\mathit{Mh}}{\parallel }_0}{{\parallel \mathit{\omega }\parallel }_0}$	${\mathit{\omega }}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{p} - {\mathit{p}}_{\mathit{Mh}}{\parallel }_0}{{\parallel \mathit{p}\parallel }_0}$	${\mathit{p}}_{\mathit{L}2}$-Rate
0.270282	$5.3045\times {10}^{-1}$		$5.7977\times {10}^0$		$1.9495\times {10}^{-1}$		$3.9490\times {10}^0$
0.14633	$1.6679\times {10}^{-1}$	1.89	$2.8543\times {10}^0$	1.15	$5.3267\times {10}^{-2}$	2.11	$1.6592\times {10}^0$	1.41
0.073319	$4.9313\times {10}^{-2}$	1.76	$1.3982\times {10}^0$	1.03	$1.4854\times {10}^{-2}$	1.85	$4.9800\times {10}^{-1}$	1.74
0.0368164	$1.5745\times {10}^{-2}$	1.66	$6.9553\times {10}^{-1}$	1.01	$4.3683\times {10}^{-3}$	1.78	$1.6378\times {10}^{-1}$	1.61
0.0188245	$5.8796\times {10}^{-3}$	1.47	$3.5345\times {10}^{-1}$	1.01	$1.4678\times {10}^{-3}$	1.63	$6.3343\times {10}^{-2}$	1.42

and

Table 9. Unstructured mesh for the local Gauss integration method.

$\frac{1}{\mathit{h}}$	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{Gh}}{\parallel }_0}{{\parallel \mathit{u}\parallel }_0}$	${\mathit{u}}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{u} - {\mathit{u}}_{\mathit{G}\mathit{h}}{\parallel }_1}{{\parallel \mathit{u}\parallel }_1}$	${\mathit{u}}_{\mathit{H}1}$-Rate	$\frac{\parallel \mathit{\omega } - {\mathit{\omega }}_{\mathit{Gh}}{\parallel }_0}{{\parallel \mathit{\omega }\parallel }_0}$	${\mathit{\omega }}_{\mathit{L}2}$-Rate	$\frac{\parallel \mathit{p} - {\mathit{p}}_{\mathit{Gh}}{\parallel }_0}{{\parallel \mathit{p}\parallel }_0}$	${\mathit{p}}_{\mathit{L}2}$-Rate
0.270282	$1.8272\times {10}^{-1}$		$3.9585\times {10}^{-1}$		$1.4450\times {10}^{-1}$		$4.6246\times {10}^0$
0.14633	$4.5578\times {10}^{-2}$	2.26	$1.9863\times {10}^{-1}$	1.12	$3.4643\times {10}^{-2}$	2.33	$1.9980\times {10}^0$	1.37
0.073319	$1.1183\times {10}^{-2}$	2.03	$9.8443\times {10}^{-2}$	1.02	$8.6578\times {10}^{-3}$	2.01	$6.1863\times {10}^{-1}$	1.70
0.0368164	$2.7425\times {10}^{-3}$	2.04	$4.8775\times {10}^{-2}$	1.02	$2.1423\times {10}^{-2}$	2.03	$1.9223\times {10}^{-1}$	1.70
0.0188245	$6.8186\times {10}^{-4}$	2.07	$2.4329\times {10}^{-2}$	1.04	$5.3335\times {10}^{-4}$	2.07	$6.5361\times {10}^{-2}$	1.61

.

4.2. The Lid-Driven Flow

This gives the lid-driven flow on a square area, which has a boundary condition of no-slip, and only at the upper boundary $\left\{\right(x,1):0<x<1\}$ satisfies $u_1=1,u_2=0,w=1$. We assume that the normal component of velocity is zero on $\partial \mathsf{\Omega }$, and the tangential component is zero, except that $y=1$ is set to 1. In Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6, we show the pressure levels and velocity streamlines with $\epsilon ={10}^{-6},1/h=20$ based on these four methods. It can be seen from Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6 that only the Gauss method can get the resolved pressure.

4.3. The Co-Axis Bearing Lubrication

In this example, we extend these stabilized methods to the nonlinear co-axis bearing lubrication problem. In Figure 7, the typical fluid domain, boundary conditions, and structured grid are drawn. The fluid domain is an angular domain between outer boundary ${\mathsf{\Gamma }}_1$ with radius $r_1$ and inner boundary ${\mathsf{\Gamma }}_2$ with radius $r_2$. The outer boundary is steady, and the inner body surrounded by ${\mathsf{\Gamma }}_2$ is supposed to rotating along axis with rotating angular velocity ${\omega }_r$. So the homogeneous boundary condition $u_x=u_y=\omega =0$ on ${\mathsf{\Gamma }}_1$, shearing velocity boundary condition $u_x=-r_2{\omega }_r\sin \left({\omega }_{r}t\right),u_y=r_2{\omega }_r\cos \left({\omega }_{r}t\right)$ and $\omega =0$ on ${\mathsf{\Gamma }}_2$ are added.

In this example, the parameters are selected as $r_1=1,r_2=0.5$, mesh size $h=1/200$, $j=\nu ={\nu }_r=c_a=c_d=0.1$, equation parameters $\nu ={\nu }_r=c_a=c_d=0.1$, $\epsilon ={10}^{-8}$, ${\beta }_1=100$, ${\beta }_2=100$, ${\beta }_3=150$. Some numerical results with rotating angular velocity ${\omega }_r=100,500,1000$ are shown in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, Figure 19, Figure 20, Figure 21, Figure 22 and Figure 23. We observe that the pressure is well obtained by the local Gauss integration method.

From the above figures, as the rotation speed increases, the magnitude of fluid components horizontal velocity, vertical velocity, angular velocity, and pressure also increase. When ${\omega }_r=100$, the pressure of the local Gauss integration method performs better. Taking different values of ${\omega }_r$, the fluid field is still smooth, correspondingly, the non-linearity on the boundary also increases. It can also be applied to practical problems and has a wide range of applications.

5. Conclusions

In this paper, based on the lowest equal-order finite element space pair, a variety of stable mixed finite element methods are numerically studied for the stationary micropolar equations. The following conclusions are drawn through numerical comparison. The stability and efficiency of all these methods depend on their parameter values. As far as the penalty method is concerned, The smaller the parameter value, the more stable the method. However, for the regular and multiscale enrichment methods, their performance largely depends on the choice of the stabilization parameters. In fact, it is difficult to choose fine parameters. The local Gauss integration method has no stable parameters and shows the best performance among the considered methods on the numerical results.