#### Abstract

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 integration method 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.

#### 1. Introduction

In this paper, we will consider the following steady incompressible flow problem with damping: seek such that where is a convex polygonal domain with a Lipschitz continuous boundary and the symbols , , and denote the Laplacian, gradient, and divergence operators, respectively; , , and represent the velocity vector, the pressure, and the prescribed body force, respectively. Further, represents the damping term with two constants and . In addition, , and we linearize the nonlinear term by allowing it to lag one step behind.

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–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–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–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.

#### 2. Notations and Preliminaries

We will use the following Hilbert spaces: Here and in what follows, the spaces () are equipped with the -scalar product and -norm or , respectively. Further, we will consider the standard definitions for Sobolev spaces which, for any integer and any number , are equipped with the norm Notice that Then the variational formulation of problem (1) is to seek such that where

Now, for convenience, we introduce the finite subspaces , 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 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 and satisfy the following conditions [7].(i) There is a constant such that where .(ii) There is a constant independent of such that Then we have a unique solution of (5) satisfying where is a positive constant.

#### 3. 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 and be finite dimensional spaces. We consider the following saddle point problems: where and are unknown variables; , ; is symmetric positive definite operator, is linear map, is the transpose operator of , and is symmetric semidefinite operator. If the initial values and are given, then and are defined by where is a given real number.

Let be the iteration error generated by the above method. It is easy to show that
Let denote the largest eigenvalue of matrix . Then converges to if is chosen such that . In this case, and converge, respectively, to **X** and **Y** with a rate of convergence bounded by the absolute value of . For more details about the saddle point problems, please see [22–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 is the set of all polynomials on of degree less than 1. Note that is not a subspace of . However, in this nonconforming case, the pair of finite element spaces is ; 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. Seek such that Then we can get the following equations from (16): Next, let and be the array of the velocity and pressure, respectively. Then it is easy to see that (17) can be written in matrix form: where , , and , ; , , and , ; , , , , and , ; ; ; and .

*Remark 2. *Here, . Because we linearized the nonlinear term , then is a known-term in the process of solving (18). And for (16), the process of linearizing is as follows.

Given , seek such that

In the nonconforming case, the discrete nonconforming formulation for the steady incompressible flow problem with damping is to seek such that
where

Note that (17) is a saddle point problem and the lowest equal-order pair does not satisfy the discrete inf-sup condition
where the constant is independent of and , for all .

*Algorithm I* (Penalty method). The penalty method compensates for the inf-sup condition deficiency by adding the penalty term as follows.

Seek such that where 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 and are presented and is deduced from , using the base for in the usual manner; that is, Let ; 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 such that
for all , where is a stabilization parameter and . The matrix form of the above stabilized version can be expressed as
where additional blocks and correspond to
respectively; that is,
and , where
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 such that
where and are the positive stabilization parameters, is the normal outward vector, is normal derivative operator, and denotes the jump of across . Moreover, a direct algebraic manipulation leads to the matrix form
where the matrix 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 such that
where is defined by
, and indicates a local Gauss integral over that is exact for polynomials of degree . Then the corresponding matrix form of this stabilized method is
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 -norm and of order for the pressure in the -norm, respectively, by using all these stabilized methods.

#### 4. 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.

##### 4.1. Numerical Test 1

In this example, we consider the exact solution problem firstly. Let the domain be the unit square . The exact solution for the velocity and pressure is given as follows: and the right-hand side is determined by the original problem (1). Our goal in this test is to compare CPU-time, the -error of the pressure, and -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 , where and are the relative errors corresponding to the meshes of sizes and , respectively.

From this experiment, we can learn the following several points. 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. 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 , , , , and , , , because they can deal with the considered problem well. 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 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 CPU-time, 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.

##### 4.2. 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 . We assume the normal component of the velocity to be zero on and the tangential component to be zero except along , where it is set to one.

In this example, we simulate the referred physics phenomena. In Figures 1, 2, 3, 4, and 5, we present the velocity streamlines and pressure level lines for , , , and 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.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

##### 4.3. Numerical Test 3

In this example, we choose the exact solution for the velocity and pressure in the unit square as follows: and the right-hand side is determined just as Test 1. In this example, parameters we choose are , , and . For the penalty and local Gauss integration methods, CPU-time, the -error of the pressure, -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.

**(a)**

**(b)**

#### 5. 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.

#### Acknowledgments

This work is in part supported by the NSF of China (nos. 11271313 and 61163027), the China Postdoctoral Science Foundation (nos. 2012M512056 and 2013M530438), the Key Project of Chinese Ministry of Education (no. 212197), the NSF of Xinjiang Province (no. 2013211B01), and the Doctoral Foundation of Xinjiang University (no. BS120102).