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

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.


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 Ω be a bounded domain in R n , n = 2, 3, with a sufficiently smooth boundary ∂Ω, the following MNSE in a dimensionless form will be considered: in Ω, div u = 0, in Ω, −ν 2 ∆ω + 4ν r ω = 2ν r rot u + g, in Ω, where the fluid variables u, ω, p are the linear velocity, angular velocity, and pressure, respectively. The symbols ν, ν r , c a , c d are used to represent some definitely given physical parameters, where ν 1 = ν + ν r , ν 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: 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.

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 (Ω), m = 0, 1, 2, . . . , is denoted by · m . Specially, if m = 0, the space H 0 (Ω) is the general Hilbert L 2 (Ω), endowed with L 2 (Ω)-scalar product (·, ·) and norm · 0 . The Sobolev space H 1 0 (Ω) is the subspace of H 1 (Ω) with homogeneous boundary condition. For the convenience of analysis, the following Sobolev spaces are introduced, From the Poincaré inequality, we know that ∀v ∈ X, ∇v 0 and v 1 are equivalent norm in H 1 (Ω) n . And H −1 (Ω) n denotes the dual space of X with the norm: where ·, · denotes duality product. Next, let introduce the following integration by parts formula for the rot operator, (rot ω, u) = (ω, rot u), ∀u, ω ∈ X.
In order to obtain the existence and uniqueness of the weak solution of the problem (6), we should introduce the following LBB condition [2].
where β > 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 ∈ L 2 (Ω) n , g ∈ L 2 (Ω) n , then there exist a unique solution (u, ω, p) ∈ X × X × M of the problem (6), which satisfies u 1 + ω 1 + p 0 ≤ c( f −1 + g −1 ). (9) Moreover, if Ω is regular of class C 2 , then the following regularity result holds [17]: Next, we consider the discrete problem. Let T h be a regular triangulation of Ω made up of triangles K. Let h K = diam{K}, h = max {h K : K ∈ T h }. Associated with T h , the finite element spaces (X h , M h ) are defined: X h = {v ∈ C 0 (Ω) n ∩ X : v| K ∈ (P 1 (K)) n , ∀K ∈ T h }, Then, the Galerkin finite element formulation of the problem (6) From the classical finite element theory [2,18,19], the following assumption is reasonable Assumption 1. There exist interpolation operators I h v ∈ X h such that Furthermore, there exist the projection operator r h : M → M h , for any given p ∈ M, Adding (22) and (23), we get Then, On the other hand, by the LBB condition, taking q = 0 in (17), Let consider the relationship between the penalty problem and the solution of the variational problem. (6) and (17) respectively, then there exist Proof. Subtracting (17) from (11), we get Taking Furthermore, we have By using (9) and (21), 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 Next, let consider the corresponding discrete weak form of the penalty method (17): (28) Proof. Choosing (v, s, q) = (u εh , ω εh , p εh ) in (28), we obtain By using Young's inequality (7), we have 4ν r u εh 1 ω εh 0 ≤ν r u εh Then, carrying the above inequalities to get and hence, the above inequality implies that On the other hand, we obtain from (30) √ ε p εh 0 ≤ c( f −1 + g −1 ). Theorem 5. Let (u ε , ω ε , p ε ) and (u εh , p εh , ω εh ) be the solution of (17) and (28), respectively. Then the error satisfies Proof. Subtracting (28) from (17) yields Adding the above two formulas together and simplifying, we get Now, by using Hölder and Young's inequalities, we obtain Applying the triangle inequality to gain ).
Theorem 6. Under the Theorems 3 and 5, the errors u − u εh , ω − ω εh and p − p εh satisfy Proof. The result based on (24) and (31), we get the above inequality.

Regular Method
The variational problem (6) of the regular method is: where . For the pressure variable p ∈ M, we define the mesh-dependent norm (35) Next, let prove the following continuous and elliptical properties for bilinear form B 2 ((u Gh , ω Gh , p Gh ), (v, s, q)). Lemma 1. The bilinear form B 2 ((u Rh , ω Rh , p Rh ), (v, s, q)), ∀(v, s, q) ∈ X h × W h × M h satisfies the continuous property: and the elliptical property: Proof. The continuous property is obvious. Let first consider the restriction of B 2 (·, ·) on each element K, by using Schwarz's inequality, Applying Young's inequality with = 1 2 , we can get the following estimate of (38) on each K, and the proof is finished by adding K ∈ T h .
We introduce the following approximation properties [4].

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].
be the solution of (1). Then, the error satisfies Proof. Noting that, if we choose τ 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.

Multiscale Enrichment Method
where where τ K and τ Z are the positive stabilization parameters, the quantity h Z = |Z| is the length of the edge Z, Z ⊂ ∂K, and [v] denotes the jump of v across Γ kj . Define the meshdependent norms Before analyzing the stability of the method given by (42), we introduce the following local trace theorems , and the elliptical property Proof. The continuous property follows using (48) and Cauchy-Schwartz inequality.
Finally, taking = 1 2 to obtain Next, we introduce the interpolation error about the Clément interpolation operator.
Proof. The result comes from the norm definition and uses q −q h 0,Ω ≤ q − C h (q) 0,Ω to combine with (15) and (16).
] be the solution of (1), and (u Mh , ω Mh , p Mh ) the solution of (42). Then, the following error estimate holds Proof. For the sake of simplicity, Hence, dividing by the last term we get Finally, combing above inequalities with the triangular inequality gives the following result.
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.
] be the solution of (1), and (u Mh , ω Mh , p Mh ) the solution of (42). Then, the following error estimate holds Proof. Known by the continuous inf-sup condition, there exists ϕ ∈ X such that ∇ · ϕ = p − p Mh and ϕ 1, Then divide by the last term to get the result.

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 , ω Gh , where for all (v, s, q) ∈ (X h , X h , M h ), and Π is a L 2 projection operator with the following properties: where R 0 ⊂ L 2 (Ω) denotes the piecewise constant space associated with the triangulation 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: where K,i p Gh qdx 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 , ω Gh , p Gh ), (v, s, q)) satisfies the continuous property whereβ is a positive constant depending only on Ω. Proof.
Thus it suffices to show the continuous property.
For the coercive property of B 4 , ∀p Gh ∈ M h , there exists a positive constant C 0 and z ∈ X such that [2] (div z, p Gh ) = p Gh Setting the finite element approximation z Gh ∈ X h of z, we have Then, for any p Gh ∈ M h , we choose any (v, s, q) = (u Gh − λz Gh , ω Gh , −p Gh ), where 0 < λ < 2(1−c 1 ) c 1 ν r . Obviously, it follows from (58)-(61), (65) and (66) and the Young inequality that Finally, we remark that settingβ = c 5 c 6 , which finishes the proof.
be the solution of (1), the stabilized finite element solution (u Gh , ω Gh , p Gh ) satisfies the error estimates: Proof. Subtracting (54) from (6), we have By setting (e, θ, η) = (I h u − u Gh , I h ω − ω Gh , r h p − p Gh ), it follows from, we can see that Obviously, it follows from (63) and (70) that Thanks to (71) and the triangle inequality

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 (Ω)and H 1 (Ω)-norms for the velocity, angular velocity, and the relative L 2 (Ω)-norm for the pressure and their corresponding convergence rates are given in Tables 1-4. In addition, the theoretical convergence rates should be of order O(h 2 ) and O(h) for the velocity and angular velocity in the L 2 (Ω)and H 1 (Ω)-norms, respectively, and of order O(h) for the pressure in the L 2 (Ω)-norm by using all these stabilized methods. It follows from Tables 1-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.

Problems with Smooth Solutions
Setting Ω = [0, 1] 2 , equation parameters ν = ν r = c a = c d = 0.1, we take penalty parameter ε = 10 −6 and ε = O(h 1 2 ) in Tables 1 and 2 respectively, β 1 = 100, β 2 = 100, β 3 = 150, and the exact solution (u, ω, p) is with the right-hand side function f generated by the exact solution: It can be seen from Tables 2-5 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 (Ω)-norm convergence order of the multiscale enrichment method do not reach O(h 2 ) with the refinement of the grid. The convergence speed of the L 2 (Ω)-norm of the pressure of the regular method, the local Gauss integration method, and the multiscale enrichment method has almost reached O(h 1.5 ), which are better than the expected results of the theoretical analysis. Additionally, note the CPU time listed in Tables 2-5, all our calculations are performed under the same computing platform. From Tables 2-5, as the grid is encrypted, the multiscale enrichment method consumes the most time, while the local Gauss integration method consumes the least time.

The Lid-Driven Flow
This gives the lid-driven flow on a square area, which has a boundary condition of noslip, and only at the upper boundary {(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 ∂Ω, and the tangential component is zero, except that y = 1 is set to 1. In Figures 2-6, we show the pressure levels and velocity streamlines with ε = 10 −6 , 1/h = 20 based on these four methods. It can be seen from Figures 2-6 that only the Gauss method can get the resolved pressure.

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 Γ 1 with radius r 1 and inner boundary Γ 2 with radius r 2 . The outer boundary is steady, and the inner body surrounded by Γ 2 is supposed to rotating along axis with rotating angular velocity ω r . So the homogeneous boundary condition u x = u y = ω = 0 on Γ 1 , shearing velocity boundary condition u x = −r 2 ω r sin(ω r t), u y = r 2 ω r cos(ω r t) and ω = 0 on Γ 2 are added. In this example, the parameters are selected as r 1 = 1, r 2 = 0.5, mesh size h = 1/200, j = ν = ν r = c a = c d = 0.1, equation parameters ν = ν r = c a = c d = 0.1, ε = 10 −8 , β 1 = 100, β 2 = 100, β 3 = 150. Some numerical results with rotating angular velocity ω r = 100, 500, 1000 are shown in Figures 8-23. We observe that the pressure is well obtained by the local Gauss integration method.  . ω r = 100, 500, 1000, horizontal velocity for the penalty method.   Figure 11. ω r = 100, 500, 1000, angular velocity for the penalty method.        Figure 19. ω r = 100, 500, 1000, angular velocity for the multiscale enrichment 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 ω r = 100, the pressure of the local Gauss integration method performs better. Taking different values of ω 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.

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.