Numerical Study on Several Stabilized Finite Element Methods for the Steady Incompressible Flow Problem with Damping

We discuss several stabilized finite element methods, which are penalty, regular, multiscale enrichment, and local Gauss integration method, for the steady incompressible flow problem with damping based on the lowest equal-order finite element space pair. Then we give the numerical comparisons between them in three numerical examples which show that the local Gauss integrationmethod has good stability, efficiency, and accuracy properties and it is better than the others for the steady incompressible flow problem with damping on the whole. However, to our surprise, the regular method spends less CPU-time and has better accuracy properties by using Crout solver.

These equations describe various physical situations such as porous media flow, drag or friction effects, and some dissipative mechanisms from the resistance to the motion of the flow.If the damping system is different, energy dissipation is different.So the application of these equations is very extensive in the daily life (see [1][2][3][4][5][6] and references therein).
It is well known that it is very difficult to compute some PDEs directly while numerical method plays an important role in these problems, so developing an efficient and effective computational method for solving the incompressible flow problem has practical significance and has drawn the attention of many researchers (see [4,[7][8][9][10][11][12] and the references cited therein).During this time, mixed finite element methods [13] are a natural choice for solving fluid mechanics equations because these equations naturally appear in mixed form in terms of velocity and pressure.
In the analysis and practice of employing mixed finite element methods in solving the incompressible flow problems, the inf-sup condition has played an important role because it ensures stability and accuracy of the underlying numerical schemes.Pairs of finite element spaces that are used to approximate the velocity and the pressure unknown are said to be stable if they satisfy the inf-sup condition.Intuitively speaking, the inf-sup condition is something that enforces a certain correlation between two finite element spaces so that they both have the required properties when employed for the incompressible flow problems.However, due to computational convenience and efficiency in practice, some mixed finite element pairs which do not satisfy the inf-sup condition are also popular.Thus, much attention has been paid to the study of the stabilized method for the Stokes problem [4,7].
In the present years, studies have focused on stabilization techniques [14,15], which include penalty method [16,17], regular method [15], multiscale enrichment method [18], and local Gauss integration method [19].There exist a lot of theoretical results for the stabilized mixed finite element methods for the Stokes equations and the comparisons between them are also given (see [4,7,15,16,[18][19][20][21] and the references cited therein).It is pointed out that authors considered the performance of several stabilized methods for the Stokes equations based on the lowest equal-order pairs in [21].And authors studied several stabilized methods for the Stokes eigenvalue problem by using different conforming and nonconforming lower-order pairs in [7].In this report, we will adopt four kinds of stabilized finite element methods but will mainly focus on the incompressible flow problem with damping based on the lowest equal-order pairs.Moreover, we present the comparisons between these methods of the considered problem.
A brief outline of the rest of our paper is organized as follows: in Section 2, we introduce some notations and present some preliminary materials and some well-known results of the steady incompressible flow problem with damping to be used in our subsequent sections; then in Section 3, we review several stabilized mixed finite element methods and recall their key stabilization techniques; in Section 4, comparisons between these stabilized methods are performed numerically; finally, we end with some short conclusions in Section 5.

Notations and Preliminaries
We will use the following Hilbert spaces: Here and in what follows, the spaces  2 (Ω)  ( = 1, 2) are equipped with the  2 -scalar product (⋅, ⋅) and  2 -norm ‖ ⋅ ‖  2 or ‖ ⋅ ‖ 0 , respectively.Further, we will consider the standard definitions for Sobolev spaces  , (Ω) which, for any integer  > 0 and any number  ≥ 1, are equipped with the norm Notice that Then the variational formulation of problem ( 1) is to seek where Now, for convenience, we introduce the finite subspaces X ℎ ×  ℎ ⊂ X × , assumed to be uniformly regular in the usual sense.Suppose that  ℎ is a triangular decomposition of the domain Ω and ℎ is the maximum mesh size of the partition.Therefore, we define where  1 () represents the space of linear functions on .
And we assume the following basis functions: where  and  are the dimensions of  ℎ and  ℎ , respectively.And the bilinear forms (k ℎ , k ℎ ) and (k ℎ ,  ℎ ) satisfy the following conditions [7].
(i) There is a constant  1 > 0 such that where (ii) There is a constant  2 > 0 independent of ℎ such that sup Then we have a unique solution (u ℎ ,  ℎ ) of ( 5) satisfying where  3 is a positive constant.

Stabilized Mixed Finite Element Methods
In this section, we will give several stabilized mixed finite element algorithms to show their different aspects and several ways have been used to stabilize the lowest equal-order finite element space pair as follows (see [7,21]).First we introduce the following classical Uzawa iterative algorithm.
Let H 1 and H 2 be finite dimensional spaces.We consider the following saddle point problems: where where  is a given real number.
Let  Y  = Y − Y  be the iteration error generated by the above method.It is easy to show that Let  1 denote the largest eigenvalue of matrix In this case, X  and Y  converge, respectively, to X and Y with a rate of convergence bounded by the absolute value of 1− 1 .
For more details about the saddle point problems, please see [22][23][24] and the references cited therein.
Next we present four kinds of stabilized finite element method for the steady incompressible flow problem with damping.
Remark 1 (nonconforming finite element space).For convenience, we let ℎ be a positive parameter and  ℎ = {  } a regular triangulation of Ω. Denote by Γ  = Ω ∩   the boundary edge and by Γ  = Γ  =   ∩   the interior boundary.Set the centers of Γ  and Γ  by   and   , respectively.The nonconforming finite element space can be defined as where  1 (  ) is the set of all polynomials on   of degree less than 1.Note that NC ℎ is not a subspace of X.However, in this nonconforming case, the pair of finite element spaces is NC ℎ ×  ℎ ; that is, the conforming space is still used for pressure.
So the corresponding discrete variational formulation of (5) for the Stokes equations with damping reads as follows.
Then we can get the following equations from ( 16): Next, let U ℎ and P ℎ be the array of the velocity and pressure, respectively.Then it is easy to see that ( 17) can be written in matrix form: where . Because we linearized the nonlinear term |u| −2 u, then U is a knownterm in the process of solving (18).And for ( 16), the process of linearizing is as follows. Given In the nonconforming case, the discrete nonconforming formulation for the steady incompressible flow problem with damping is to seek where Note that ( 17) is a saddle point problem and the lowest equal-order pair does not satisfy the discrete inf-sup condition sup or sup where the constant  4 > 0 is independent of ℎ and ‖∇k ℎ ‖ 0,ℎ = (∑  |k| 2 1,  ) 2 , for all k ∈ NC ℎ .
Algorithm I (Penalty method).The penalty method compensates for the inf-sup condition deficiency by adding the penalty term as follows.
where  > 0 is a penalty parameter.The performance of this method obviously depends on the choice of the penalty parameter .Then the matrix form of ( 23) can be expressed as where the matrixes A 1 and B 1 are presented and D 1 is deduced from ( ℎ , ), using the base for  ℎ in the usual manner; that is, Let  P ℎ  = P ℎ − P ℎ  ( = 1, 2, . ..); we use the above Uzawa iterative algorithm and get Algorithm II (Regular method).This method uses a simple way to stabilize the mixed finite element approximation without a loss of accuracy, that is, to seek for all (k, ) ∈ X ℎ ×  ℎ , where  = ℎ 2 /(]) is a stabilization parameter and  > 0. The matrix form of the above stabilized version can be expressed as where additional blocks D 2 and C ℎ 1 correspond to respectively; that is, and Similar to (26), we also have Algorithm III (Multiscale enrichment method).Another stabilized way is the multiscale enrichment approach which includes the usual Galerkin least squares stabilized terms on each finite element and positive jump terms at interelement boundaries.Namely, seek where  1 = ℎ 2 /( 1 ]) and  2 = ℎ/( 2 ]) are the positive stabilization parameters,  is the normal outward vector,   is normal derivative operator, and [k] denotes the jump of k across .Moreover, a direct algebraic manipulation leads to the matrix form where the matrix D 3 is deduced from the term Similar to (26), we have Algorithm IV (Local Gauss integration method).The local Gauss integration method is to add two Gauss integrals rather than any stabilization parameter to the original discrete formulation (16) as follows.Seek (u ℎ ,  ℎ ) ∈ X ℎ ×  ℎ such that where ( ℎ , ) is defined by where As (26), we get It is well known that, if  is well chosen, then   and   converge, respectively, to  and  with a rate of convergence based on (26), (32), (35), and (40).From these equations, we can find that the penalty method converges faster than the local Gauss integration method.We can obtain that coefficient matrices of penalty algorithm, regular algorithm, and local Gauss integration algorithm are all symmetric; however, the multiscale enrichment method's coefficient matrices are not symmetric from ( 10), ( 28), (34), and (38).What is more, it is easy to see that the matrix calculations of multiscale enrichment algorithm are more complex and cumbersome, so this method maybe costs more time.
By using the regularity assumptions and well-established techniques for velocity and pressure [4,7], the theoretical convergence rates should be of order (ℎ) for the velocity in the  2 -norm and of order (ℎ 2 ) for the pressure in the  1norm, respectively, by using all these stabilized methods.

Numerical Experiments
In this section, we will give three numerical tests to confirm the numerical theory developed in the previous section.In the given experiments, the pressure and velocity are approximated by the lowest equal-order finite element pairs defined with respect to the same uniform triangulation; that is, the mesh consists of triangular elements that are obtained by dividing Ω into subsquares of equal size and then drawing the diagonal in each subsquare.

Numerical Test 1.
In this example, we consider the exact solution problem firstly.Let the domain Ω be the unit square Ω = (0, 1) × (0, 1) ⊂ R 2 .The exact solution for the velocity u = ( 1 ,  2 ) and pressure  is given as follows: and the right-hand side f = ( 1 (, ),  2 (, )) is determined by the original problem (1).Our goal in this test is to compare CPU-time, the  2 -error of the pressure, and  1 -error of the velocity; the experimental rates of convergence for these methods with different values of ℎ are tabulated in Tables 1, 2, 3, 4, and 5. What is more, the rates of convergence are calculated by the formula log(  / +1 )/ log(ℎ  /ℎ +1 ), where   and  +1 are the relative errors corresponding to the meshes of sizes ℎ  and ℎ +1 , respectively.
From this experiment, we can learn the following several points.(1) For penalty method, the result of this method is well in which parameter "" can be 2, 3, and 4 and we can use UMFPACK or default solver in the process of calculation.(2) Results got from the penalty, regular, multiscale enrichment, and local Gauss integration methods by using conforming and nonconforming elements are presented in Tables 1-5, respectively.Here, we choose  = 1.0 − 6,  = 1.0 − 4,  = 3, Re = 10000, and  = 160,  1 = 160,  2 = 100, For regular method, it works well only we choose Crout solver; however, for other methods, several solvers can be used and their difference is not big.What is more, the value of "" has little influence on the results within a certain range; in this paper, we give the results of  = 3 presented in tables.
From Tables 1-5, we can see that these methods work well and keep the convergence rates just like the theoretical analysis except for the multiscale enrichment method.Meanwhile, it can be seen that the penalty method requires the least CPUtime, which validates the analysis in Section 3. As expected, we have an interesting observation that the error of the nonconforming local Gauss integration method is better than that of the conforming version, which is not surprising since the degree of freedom of the nonconforming method is nearly three times than that of the conforming one on uniform mesh.Hence, it is natural that the nonconforming local Gauss integration method is more accurate and costs more CPU-time.The penalty method should use less time than the regular method and local Gauss integration method theoretically but in fact not; maybe this is caused by the different solver.
Besides, from the convergence results on this example, we can see that regular and multiscale enrichment methods are not better than other methods.And the nonconforming local Gauss integration method shows the best numerical stability.

Numerical Test 2.
In this test, we test a popular benchmark problem, the lid-driven flow.Let the computation be carried out in the region Ω = {(, ) | 0 < ,  < 1}.We assume the normal component of the velocity to be zero on Ω and the tangential component to be zero except along  = 1, where it is set to one.
In this example, we simulate the referred physics phenomena.In Figures 1, 2   and  = 3.For the penalty and local Gauss integration methods, CPU-time, the  2 -error of the pressure,  1 -error of the velocity, and the experimental rates of convergence for these methods with different values of ℎ are tabulated in Tables 6 and 7. What is more, the numerical solutions are given in Figure 6.From this test, we can get conclusion that is similar to Test 1.However, the result got from multiscale method is not better.

Conclusions
We have used several stabilized mixed finite element methods in solving the steady incompressible flow problem with damping based on the lowest equal-order pairs in this paper.We give some conclusions by comparing numerically as follows.
All of those methods' stability and efficiency depends on their parameter values.In terms of the penalty method,   the smaller its parameter value, the more stable the method.However, it can not be too small, otherwise, the condition number of the system matrix arising from this method will become too large to be solve.For the regular and multiscale enrichment methods whose performance heavily depends on the choice of the stabilization parameters, however, it is difficult to choose fine parameters in fact.What is more, a poor choice of these stabilization parameters can also lead to serious deterioration in the convergence rates.The local Gauss integration method is free of stabilization parameters and shows numerically the best performance among the methods considered for the given problem.
, 3, 4, and 5, we present the velocity streamlines and pressure level lines for  = 10 − 6,  = 100, ℎ = 1/30, and  = 4 based on these five methods.From Figures 1(b)-5(b), only Gauss methods can obtain resolved pressure.For the velocity, from Figures 1(a)-5(a), we can see that Gauss methods can capture this model better than the other methods.
3 > 0, and ∫ , () indicates a local Gauss integral over  that is exact for polynomials of degree  ( = 1, 2).Then the corresponding matrix form of this stabilized method is