Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volume 2012 (2012), Article ID 307939, 12 pages
Research Article

An Alternative HSS Preconditioner for the Unsteady Incompressible Navier-Stokes Equations in Rotation Form

Department of Mathematics and Statistics, University of West Florida, Pensacola, FL 32514, USA

Received 2 November 2011; Accepted 27 January 2012

Academic Editor: Kok Kwang Phoon

Copyright © 2012 Jia Liu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


We study the preconditioned iterative method for the unsteady Navier-Stokes equations. The rotation form of the Oseen system is considered. We apply an efficient preconditioner which is derived from the Hermitian/Skew-Hermitian preconditioner to the Krylov subspace-iterative method. Numerical experiments show the robustness of the preconditioned iterative methods with respect to the mesh size, Reynolds numbers, time step, and algorithm parameters. The preconditioner is efficient and easy to apply for the unsteady Oseen problems in rotation form.

1. Introduction

We study the numerical solution methods of the incompressible viscous fluid problems with the following form:

Equations (1.1) to (1.4) are also known as the Navier-Stokes equations. Here is an open set of , where , or , with boundary ; the variable is a vector-valued function representing the velocity of the fluid, and the scalar function represents the pressure. The pressure field, , is determined up to an additive constant. To uniquely determine , we may impose some additional condition, such as

The source function is given on . Here is a given constant called the kinematic viscosity, which is . is the Reynolds number: , where denotes the mean velocity, and is the diameter of , see [1]. Also, is the (vector) Laplacian operator in dimensions, is the gradient operator, and is the divergence operator. In (1.3) is some boundary operator; for example, the Dirichlet boundary condition ; Neumann boundary condition , where denotes the outward-pointing normal to the boundary, or a mixture of the two.

We use fully implicit time discretization and picard linearization to obtain a sequence of Oseen problems, that is, linear problems of the form

Here is a known divergence-free vector obtained from the previous linearized step (e.g., ). We call the vector the wind function. In addition, where is the time step. Equations (1.6)–(1.8) are referred to as the Oseen problem.

We can use either finite element or finite different methods to discretize (1.6)–(1.8). The resulting discrete system has the form

In this paper, we are interested in an alternative linearization of the steady-state Navier-Stokes equation. Based on the identity In order to linearize it, we replace in one place with a known divergence-free vector which can be the solution obtained from the previous Picard iteration. In this case we have

After substituting the right-hand side into (1.6), we find that the corresponding linearized equations have the following form: where is the so-called Bernoulli pressure. For the two-dimensional case where is a scalar function.

In the three-dimensional case, we have here , where denotes the th component of . Assume , then we have the formal expression of

Here the divergence-free vector field again denotes the approximate velocity from the previous Picard iteration. Note that when the “wind” function is irrotational (), (1.12)–(1.14) reduce to the Stokes problem. It is not difficult to see that the linearizations (1.6)–(1.8) and (1.12)–(1.14), although both conservative, are not mathematically equivalent. The momentum equation (1.12) is called the rotation form. We can see that no first-order terms in the velocities appear in (1.12); on the other hand, the velocities in the scalar equations comprising (1.12) are now coupled due to the presence of the term . The disappearance of the convective terms suggests that the rotation form (1.12) of the momentum equations may be advantageous over the standard form (1.6) from the linear solution point of view. This observation was first made by Olshanskii and his coworkers in [25]. In their papers, they showed the advantages of the rotation form over the standard convection form in several aspects.

After we discretize the Oseen problem in rotation form (1.12)–(1.14), we obtain the sparse linear system , where

Here , where is the discretization of the operator , and matrix is the discretization of the term , where . In the 2D case, The rectangular matrix is the discretization of the negative divergence, and is the discretization of the gradient.

If we use a finite difference method, like Mac and Cell (MAC), see [6], then is a diagonal matrix where its diagonal elements are the values of evaluated at the grid edges. Matrix is a weighted mass matrix if a finite element method is used. In the 3D case, we have

Again matrices , and are all diagonal matrices or weighted mass matrices. Typical sparsity patterns for in the 2D and 3D case are displayed in Figures 1(a) and 1(b).

Figure 1: Sparsity patterns for different types of the Oseen problem in rotation form.

For some discretization methods, a stabilization matrix needs to be added to the (2,2) block of , namely, a matrix , where is a symmetric positive semidefinite diagonal or scaled mass matrix, or scaled Laplacian with small norm. Figures 2(a) and 2(b) show the sparsity pattern for the coefficient matrix with or without stabilization term in the 2D case. Such a stabilization is not necessary for the MAC discretization.

Figure 2: Sparsity patterns for different types of the 2D Oseen problem in convection form.

To solve the system , we can consider the Krylov subspace methods with the preconditioning. Many powerful preconditioning techniques have been explored for the generalized Oseen problems, for example, Uzawa-type preconditioner, block and approximate schur complement preconditioner, pressure preconditioner, and so forth, see [711] for more details. However, there is no “best” preconditioner for the saddle point system. To find the “best” preconditioner, we would like to find a preconditioner , such that the rate of convergence of the preconditioned Krylov subspace matrix is low and bounded independent of the mesh size, viscosity and time step . In addition, the cost of the preconditioning steps must be low. In this paper, we describe such a new preconditioner that satisfies the above requirements in most of cases and demonstrate its utility.

A summary of the paper is as follows. Section 2 demonstrates the Alternative Hermitian and Skew-Hermitian (AHSS) preconditioner; studies some of its convergence properties and the application of the HSS preconditioner for Krylov subspace methods; Section 3 shows the results of a series of numerical experiments. Finally, section 4 summarizes the approach and future work.

2. The Alternative HSS Preconditioner

The alternative HSS preconditioner is based on the nonsymmetric formulation We have analyzed the advantages of the nonsymmetric formulation in [12]. We gain positive (semi)-definiteness in this case. By changing the sign in front of the and blocks, we obtain an equivalent linear system with a matrix whose spectrum is entirely contained in the half-plane . (Here we use to denote the real part of ).The spectra of the nonsymmetric formulation is more friendly to the convergence of Krylov subspace iterations. For an example, GMRES methods, see [13, 14].

We have investigated the preconditioner based on the Hermitian and Skew-Hermitian splitting methods for the Navier-Stokes problem, see [12, 15, 16]. However, the HSS preconditioner still has some problems. When the time step is not small enough or the viscosity is relative larger, the iteration number increases a lot. Therefore, we are trying to find another preconditioner which works better than the HSS preconditioner. We find out that if we use a different splitting of the coefficient matrix, we can get a very good results. The following splitting is the new preconditioner we will introduce in the paper. Letting and , we have the following splitting of into two parts:

We denote

Therefore, we defined the preconditioner as the following: Here denotes the identity matrix of order , and is a parameter.

Similar in spirit to the classical ADI (alternating-direction implicit) method, we consider the following two splittings of : Here denotes the identity matrix of order . Note that is the shifted discretized Stokes problem, where denotes the identity matrix of order , and denotes the identity matrix of order . We obtian that is nonsingular and has positive definite symmetric part.

Alternating between these two splittings leads to the the following iteration: (). Here denotes the right-hand side of (1.9); the initial guess is chosen arbitrarily. Elimination of from (2.8) leads to a stationary (fixed-point) iteration of the form where is the iteration matrix and . The iteration converges for arbitrary initial guesses and right-hand sides to the solution if and only if , where denotes the spectral radius of .

Theorem 2.1. Consider the problem (2.1), that is, We assume that is positive real, and has full rank. Then the iteration (2.9) from the splitting (2.3) is unconditionally convergent; that is, for all and .

Proof. Consider the splitting (2.3). The iteration matrix is similar to where is symmetric and .
By Kellogg's lemma, since is positive semidefinite. is the unitary matrix so . Therefore,
We claim that .
Assume that is one eigenvalue of the preconditioned linear system . We have Therefore, We claim that , otherwise, . It turns out that
However, , and is orthogonal similar with the matrix which is a skew symmetric matrix with only pure imaginary eigenvalues. Thus, if is full rank, we claim that .
We define that . Since , is welldefined. Thus, with , we consider the following equation: If , then, . Since , we need to show , which means that is not a pure imaginary number.
Next we will prove that is not a pure imaginary number.
Consider the system We can obtain the following system of the equations:
We solve from the second equation . Plug in into the first equation, we have Applying to the both sides, we can obtain the following equation: Therefore, is a Hermitian matrix. Suppose that , that is, , where . We denote that . While , which leads to . A contradiction. Because is the discretization of the Oseen problem which is not Hermitian.
Thus, is not a pure imaginary number.

3. Application of the Preconditioner

To solve the preconditioner , we first solve the system for , followed by

The first system requires solving systems with coefficient matrix , which is a system from discretized Stokes problem. We have many efficient solvers to solve this type of the system, see [1719].

The second system requires solving the sparse tridiagonal matrix . This can be done by sparse LU factorization, preconditioned GMRES method. Notice that since is a tridiagnoal matrix, it is very easy to solve this system. In practice, we can solve (3.1) and (3.2) with inexact solvers. Our experience is that the rate of convergence of the outer Krylov subspace iteration is scarcely affected by the use of inexact inner solves. We can only use 1 step of pcg method for (3.1) and 1 step of gmres method for (3.2).

4. Numerical Experiments

In this section we report on several numerical experiments meant to illustrate the behavior of the HSS preconditioner on a wide range of model problems. We consider both Stokes and Oseen-type problems, steady and unsteady, in 2D. The use of inexact solves is also discussed. Our results include iteration counts, as well as some timings and eigenvalue plots.

All results were computed in Matlab 7.6.0 on one processor of an Intel core i7 with 8 GB of memory.

In all experiments, a symmetric diagonal scaling was applied before forming the preconditioner, so as to have all the nonzeros diagonal entries of the saddle point matrix equal to 1. We found that this scaling is beneficial to convergence, and it makes finding (nearly) optimal values of the shift easier. Of course, the right-hand side and the solution vector were scaled accordingly. We found, however, that without further preconditioning Krylov subspace solvers converge extremely slowly on all the problems considered here. Right preconditioning was used in all cases.

Here we consider linear systems arising from the discretization of the linearized Navier-Stokes equations in rotation form. The computational domain is the unit square . Homogeneous Dirichlet boundary conditions are imposed on the velocities; experiments with different boundary conditions were also performed, with results similar to those reported below. We experimented with different forms of the divergence-free field appearing (via ) in the rotation form of the unsteady Oseen problem. Here we present results for the choice (2D case) and (3D case). The 2D case corresponds to the choice of . The equations were discretized with a Marker-and-Cell (MAC) scheme with a uniform mesh size . The outer iteration (full GMRES) was stopped when where denotes the residual vector at step . For the results presented in this section, we use the exact solver to solve the two systems.

In Tables 1, 2, and 3, we present results for unsteady problems with to for different values of . We can see from the table that AHSS preconditioning results in fast convergence in all cases, and that the rate of convergence is virtually -independent. Here as in all other unsteady (or quasi-steady) problems that we have tested, the rate of convergence is not overly sensitive to the choice of , especially for small . A good choice is for the two most grids.

Table 1: Results for 2D unsteady Oseen problem different values of (exact solves) and .
Table 2: Results for 2D unsteady Oseen problem different values of (exact solves) and .
Table 3: Results for 2D unsteady Oseen problem different values of (exact solves) and .

5. Conclusions

In this paper, we have considered preconditioned iterative methods applied to discretizations of the Navier-Stokes equations in rotation form. We focus on the unsteady case of the linearized Navier-Stokes problem. We have compared the performance of the alternative HSS (AHSS) preconditioners with regard to the mesh size, the Reynolds number, the time step, and other problem parameters. We find that the AHSS preconditioner has a robust behavior especially for unsteady Oseen problems. Although our computational experience has been limited to uniform MAC discretizations and simple geometries, the preconditioner should be applicable to more complicated problems and discretizations, including unstructured grids.

Compared with HSS (Hermitian and Skew-Hermitian preconditioner), the AHSS preconditioner works better for relative large viscosity. For an example, . For the smaller viscosity, , HSS preconditioner will be recommended.

In the future study, we will investigate the performance the AHSS preconditioner based using the inexact solvers for the inner iteration. Also the picard's iteration will be tested.


The authors would like to thank Michele Benzi for helpful suggestions.


  1. H. C. Elman, D. J. Silvester, and A. J. Wathen, Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics, Numerical Mathematics and Scientific Computation, Oxford University Press, New York, NY, USA, 2005.
  2. M. A. Olshanskii, “An iterative solver for the Oseen problem and numerical solution of incompressible Navier-Stokes equations,” Numerical Linear Algebra with Applications, vol. 6, no. 5, pp. 353–378, 1999. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  3. M. A. Olshanskii, “A low order Galerkin finite element method for the Navier-Stokes equations of steady incompressible flow: a stabilization issue and iterative methods,” Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 47-48, pp. 5515–5536, 2002. View at Publisher · View at Google Scholar · View at MathSciNet
  4. M. A. Olshanskii, “Preconditioned iterations for the linearized Navier-Stokes system in the rotation form,” in Proceedings of the Second MIT Conference on Computational Fluid and Solid Mechanics, K. J. Bathe, Ed., pp. 1074–1077, 2003.
  5. M. A. Olshanskii and A. Reusken, “Navier-Stokes equations in rotation form: a robust multigrid solver for the velocity problem,” SIAM Journal on Scientific Computing, vol. 23, no. 5, pp. 1683–1706, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  6. F. H. Harlow and J. E. Welch, “Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface,” Physics of Fluids, vol. 8, no. 12, pp. 2182–2189, 1965. View at Google Scholar
  7. H. C. Elman, D. J. Silvester, and A. J. Wathen, “Performance and analysis of saddle point preconditioners for the discrete steady-state Navier-Stokes equations,” Numerische Mathematik, vol. 90, no. 4, pp. 665–688, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  8. M. Benzi, M. J. Gander, and G. H. Golub, “Optimization of the Hermitian and skew-Hermitian splitting iteration for saddle-point problems,” BIT. Numerical Mathematics, vol. 43, pp. 881–900, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  9. M. Benzi and G. H. Golub, “A preconditioner for generalized saddle point problems,” SIAM Journal on Matrix Analysis and Applications, vol. 26, no. 1, pp. 20–41, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  10. M. Benzi, G. H. Golub, and J. Liesen, “Numerical solution of saddle point problems,” Acta Numerica, vol. 14, pp. 1–137, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  11. M. Benzi and V. Simoncini, “On the eigenvalues of a class of saddle point matrices,” Numerische Mathematik, vol. 103, no. 2, pp. 173–196, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  12. J. Liu, Preconditioned Krylov subspace methods for incompressible flow problems, Ph.D. thesis, Emory University, Ann Arbor, Mich, USA, 2006.
  13. Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, Pa, USA, 2nd edition, 2003.
  14. Y. Saad and M. H. Schultz, “GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM Journal on Scientific and Statistical Computing, vol. 7, no. 3, pp. 856–869, 1986. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  15. Z.-Z. Bai, G. H. Golub, and M. K. Ng, “Hermitian and skew-Hermitian splitting methods for non-Hermitian positive definite linear systems,” SIAM Journal on Matrix Analysis and Applications, vol. 24, no. 3, pp. 603–626, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  16. V. Simoncini and M. Benzi, “Spectral properties of the Hermitian and skew-Hermitian splitting preconditioner for saddle point problems,” SIAM Journal on Matrix Analysis and Applications, vol. 26, no. 2, pp. 377–389, 2004/05. View at Publisher · View at Google Scholar · View at MathSciNet
  17. O. A. Karakashian, “On a Galerkin-Lagrange multiplier method for the stationary Navier-Stokes equations,” SIAM Journal on Numerical Analysis, vol. 19, no. 5, pp. 909–923, 1982. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  18. J. Peters, V. Reichelt, and A. Reusken, “Fast iterative solvers for discrete Stokes equations,” SIAM Journal on Scientific Computing, vol. 27, no. 2, pp. 646–666, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  19. G. M. Kobelkov and M. A. Olshanskii, “Effective preconditioning of Uzawa type schemes for a generalized Stokes problem,” Numerische Mathematik, vol. 86, no. 3, pp. 443–470, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet