Robust Multigrid Method for Solving the Navier-Stokes Equations on Structured Grids

The paper represents an efficient multigrid algorithm without the problem-dependent components for solving the Navier–Stokes equations in primitive variables formulation on structured grids. The algorithm consists of Vanka smoother, pressure decomposition and robust multigrid technique. Detailed description of the proposed approach and results of numerical experiment are given.


Introduction
Let N    be a bounded, connected domain with a piecewise smooth boundary  .Navier-Stokes equations governing flow of a Newtonian, incompressible viscous fluid can be written in the following operator form ( ) where  is a nonlinear convection-diffusion operator, P  is a pressure gradient, F and G are source terms.Given a boundary data, the problem is to find the velocity field V  and pressure P .In follows we assume that the boundary conditions are included in the operators and source terms.
Discretization of (1) using the finite differences or the finite elements and Picard or Newton linearization result in a generalized saddle point system 0 in which  and  represent the discrete velocities and discrete pressure, respectively.Here nonsymmetric A is a block diagonal matrix, where each block corresponds to a discrete convection-diffusion operator  with appropriate boundary conditions.The rectangular matrix T B represents the discrete gradient operator while B represents its adjoint, the divergence operator.
The saddle point problem (2) can be solved by the preconditioned Uzawa algorithm defined as , where matrix Q is a preconditioner.Substitution of ( 1) It is result to the following estimation where  is an exact solution.Choice of the preconditioner Q so 1 1 1 T I Q BA B q leads to the Uzawa iteration convergence .
Preconditioning has been and remains a most active area of research, some the most widely used and promising methods are described in Benzi (2005).Total efforts needed for solving the Navier-Stokes equations by Uzawa method are sum of efforts for the algorithm adaptation to the given problem (choice of the preconditioner Q , problem-dependent components for the matrices A and Q inversion, etc.) and efforts for solution of the saddle point problem (2).Optimal adaptation of the algorithm to the given problem is difficult question, but it guarantees high computational efficiency of the Uzawa solver.
On the other hand, the algorithms with the least number of the problem-dependent components seem to be more preferable for practical applications because of reduction of efforts for their adaptation.However in this case it is difficult to obtain high computational efficiency.
The goal of this paper is to present and test an iterative robust solver for the Navier-Stokes equations on structured grids.The method is based on the problem-independent components: Vanka smoother, pressure decomposition and robust multigrid technique.Special attention is paid to robustness and efficiency of the offered algorithm.

Vanka Smoother
Two-dimensional steady Navier-Stokes equations (1) can be written as: a) continuity equation where Re is a Reynolds number.
Discrete linearized continuity and momentum equations can be written as: a) discrete continuity equation b) discrete linearized X-momentum Coefficients in the momentum equations (excepting discrete pressure gradients) are functions of the velocity components.Stencils for the approximation of the continuity (4a) and momentum (4b)-(4c) equations on a staggered grid are shown on Figure 1.  (at the four cell faces) and the pressure ij p (at the center of the cell) simultaneously for each cell (Figure 1, right) Vanka (1986).The discrete X and Y-momentum equations (5b)-(5c) for the four cell faces together with the continuity equation (5a) can be written as where the components of the right-side vector are given by The system (6) can be solved by Gaussian elimination.Accounting nonlinear nature of the momentum equations, the velocity components and pressure can be updated using underrelaxation as .
In practice the underrelaxation parameters are taken as u     and 1 p   (Thompson & Ferziger, 1989).If a starting guess is not close to the solution, it is recommended values 0.2 0.8 in the first iterations.
Vanka method has no the problem-dependent components (preconditioners, artificial boundary conditions for pressure, relaxation parameters (if a starting guess is close to the solution) etc.), but slow convergence rate is observed and the rate depends strongly on the unknown ordering.As a result, the method requires different acceleration techniques.Since Vanka-like methods have very good smoothing properties, often the approach is used in multigrid algorithms as smoother.

Pressure Decomposition
In partial cases the fluid flows can be described by the simplified Navier-Stokes equations.Basic assumption is that the pressure is not changed across the flow.Numerical solution of the simplified Navier-Stokes equations can be reduced to solution of series of the saddle point problem (2) with zero block of the least size in the coefficient matrix.This fact allows develop very efficient algorithms for solving the simplified Navier-Stokes equations.Unfortunately these algorithms cannot be applied for solving the full Navier-Stokes equations because the pressure is changed in all spatial directions.However algorithms developed for the simplified Navier-Stokes equations can be used for convergence acceleration of the iterative methods for the full Navier-Stokes equations.
For the given purpose, it is necessary to artificially extract «one-dimensional parts of pressure» from the pressure field by adding and subtracting items ( ) ) pressure is represented as sum of 1 N  «components», therefore the method requires N extra conditions for computation of «the one-dimensional components».The convergence acceleration technique uses N mass conservation equations as a priori information of physical nature.
Remark 3. In spite of representation of the pressure as sum of 1 N  «components», all momentum equations have only two «pressure gradients».For example, for X-momentum we obtain Remark 4. Efficiency of the acceleration technique depends strongly on the flow nature.For directed fluid flows (for example, flows in nozzles, pipes etc.) gradient of one of «one-dimensional components of pressure» ( ) However for strongly rotated flows (for example, flow in a driven cavity) the approach shows the least efficiency.
Remark 5.In 3D case the method will be more efficient than in 2D case.
Remark 6. Velocity components and corresponding «one-dimensional components of pressure» in ( 7) are computed only in coupled manner.Velocity components and «multidimensional component» ( ) xyz p t,x, y,z in (7) can be computed in decoupled or coupled manner.
Remark 7. Gradients of the «one-dimensional components» can be obtained in analytical form for the explicit schemes.Implicit schemes require an auxiliary problem for determination of the «one-dimensional components» in (7).
The pressure decomposition can be used for fast computation of a starting guess close to the solution of the full Navier-Stokes equations or built-in in some basic solver (SIMPLE, Uzawa or Vanka method, etc.).
Similarly to Vanka method, the convergence acceleration technique does not require the problem-dependent components.Detailed description of the approach and its applications for solving benchmark and applied problems are given in Martynenko (2009;2011).

Model Problem
For clearness of the multigrid algorithm representation, we consider problem about 2D steady fluid flow in unit cavity.Geometry of the problem is shown on Figure 2. Boundary conditions for the velocity components in (4) are given by 100 (0 2 ) 0 2 (0 ) 0 02 Convergence acceleration based on the pressure decomposition ( 7) requires two mass conservation equations.Integration of the continuity equation (4a) over the control volumes 1 V and 2 V shown on Figure 2 gives The mass conservation equations (8a) and (8b) will be used for fast computation of the «one-dimensional components of pressure» in (7).
Figure 2. Geometry of the model problems and the control volumes 1 V and 2 V

Robust Multigrid Technique
Robust multigrid technique (RMT) has been developed as a variant of the geometric multigrid methods with the problem-independent transfer operators Martynenko (2006).To overcome problem of robustness, RMT consists of two parts: analytical (adaptation of the boundary value problems to RMT) and computational (solution of the adapted problems by original multigrid algorithm).RMT is intended for solving a large class of nonlinear applied problems on the structured grids with almost optimal convergence rate.Really the smoothing procedure is single problem-dependent component of RMT.

Analytical Part: Adaptation of the Navier-Stokes Equations to RMT
Adaptation of the Navier-Stokes equations (1) to RMT consists of representation of the solution as sum of two functions where the functions V C  and P C will be coarse grid corrections to velocity and pressure, and the functions V  and P will be approximations to the solution in the following multigrid iterations.
Representation ( 9) is called  -modification of the solution.Substitution of ( 9) into (1) yields the following  -modified form of the Navier-Stokes equations (1) Since the nonlinear convection-diffusion operator can be represented as we obtain the  -modified form of the Navier-Stokes equations ( ) ( ) where the source terms take the form Convergence of RMT means reduction of the corrections V C  and P C , therefore the approximations to the solution V  and P will satisfy to the Navier-Stokes equations (1) i.e. 0 ( )+ ( ) 0 0 0 ( ) Two-dimensional steady Navier-Stokes equations (4) can be rewritten in the  -modified form as: a)  -modified continuity equation where discrete analogues of the functions u c , c  and p c will be coarse grid corrections for velocity components u ,  and pressure p and discrete analogues of the functions û ,  and p will be approximations to the solutions in the following multigrid iterations.
The source terms in the  -modified Navier-Stokes equations coincide with non-modified ones: Additional convection terms in  -modified momentum equations ( 10b) and (10c) are result of nonlinear nature of the convection-diffusion operator  (i.e.     ).Approximation of the source terms (11) in Eq. ( 10) defines the accuracy, monotonicity and conservatism of the numerical solutions.Approximation of other terms in (10) affects only on the multigrid convergence rate.
Boundary conditions for the velocity components are  -modified in the same manner.Since the Dirichlet boundary conditions are given for the model problem, we obtain zero Dirichlet boundary conditions for the velocity corrections u c and c  .
Also the mass conservation equations ( 8) should be rewritten in the  -modified form: x u ĉ x, y dy u x, y dy u , y dy , d For clearness the pressure decomposition (7) will be used for fast computation of a starting guess.First of all it is necessary to formulate an auxiliary problem, we hope that the solution of the auxiliary problem will be close to the solution of the Navier-Stokes equations (4) as compared with zero starting guess.Pressure decomposition for the 2D steady Navier-Stokes equations in «corrections-residuals» variables formulation ( 10)-( 11) is written as Formulation of the auxiliary problem is based on replacement of the  -modified continuity equation (10a) by the  -modified mass conservation equations (12).Accounting (13), we obtain: a)  -modified X-momentum (10b) and  -modified mass conservation equations (12a) where the braces mean that the momentum and mass conservation equations are solved only in coupled manner.
Let us summarize the main features of the auxiliary problem ( 14): 1) due to (13), the momentum equations in the auxiliary problem ( 14) are not pressure-linked.Systems (14a) and (14b) can be solved by line (2D) or plane (3D) Seidel method coupled with secant iterations.Previously similar algorithm had been proposed by Briley (1967) for solving the simplified Navier-Stokes equations.
2) solution of the auxiliary problem satisfies to the mass conservation equations ( 12), but not satisfies to the continuity equation (10a).It gives us a reason to hope that solution of ( 14) will be an accurate starting guess for the Navier-Stokes equations (4).
To illustrate abovementioned deduction, Figures 3 and 4 represent comparison of the solutions of the auxiliary problem ( 14) and the Navier-Stokes equations for the model problem (Re=100, grid 101´101).Good agreement is observed in the stream function behavior, but agreement in pressure is not accurate because of only «part of pressure» ( x y p p  ) is computed in the auxiliary problem ( 14).

Computational Part
First stage of the computational part of RMT consists of the finest grid generation in domain for the following control volume approximation of the modified equations.The finest grid 0 1 G consists of two sets v (0;1) G and f (0;1) G of the grid points defined as The finest grid 0 1 G with 0 H 8 x  is shown on Figure 5. G have no common points, i.e.
3) all grids are similar to each other, but a mesh size on the coarse grids is three times as large as than the mesh size on the finest grid.
4) the discrete functions can be assigned to the grid points v x or to the grid points f x , but in both cases the control volume on the coarse grids is union of three control volumes on the finest grid (Figure 7).The finest grid forms the zero level and three coarse grids form the first level.The coarse grid generation is further recurrently repeated: each grid 1 2 3  of a current level l is considered to be the finest grid for the coarse grids  of the next level 1 l  .Nine coarse grids derived from the three grids of the first level form the second level, etc.The coarse grid generation is finished when no further coarsening can be performed.The grid hierarchy will be called a multigrid structure (Figure 8).

Second level
First level Third level ) can be represented as product of N one-dimensional grids, similar triple coarsening is performed independently in each spatial direction.Therefore th l level consists of 3 Nl grids in multidimensional case.
The number of levels can be computed in advance.Assume that majority of the coarsest grids has three grid points.Then the number of the finest grid points is 0 H 1 , where L  is number of the coarsest level.Therefore where square brackets mean integer part.
In 3D case the finest grid

Each level { }
x y z l l ,l ,l  consists of the computational grids .Control volume on the grid   control volumes on the finest grid.
An example of the coarse grid generation (two level structure) for the finest grid with 0 H 30 x  is shown on Figure 9.Let us define the mapping of indices to simplify operations with the multigrid structure.The notion «a grid of the th l level» means the one-to-one mapping of indices of the coarse grid points onto the indices of the finest grid points.In the rest of the paper, the braces {} will denote the mapping.The mapping of indices of the grid points from the sets v G and f G is written as v { } i x and f { } i x respectively, where i and {i} are the coarse and finest grid indices.For example, for grid 1 1 G shown on Figure 9 we obtain The mapping of indices gives a close-to-the-finest-grid notation.For example, the second order derivative of the discrete function assigned to the grid points v x on the multigrid structure is approximated as Figure 9 shows that the approximation in the grid points .

Multigrid Iterations
Multigrid cycle of RMT for solving linear problems is shown on Figure 10.Computations start on the coarsest level.When the coarsest level solution has been obtained, the transfer to the next finer level is performed.It should be emphasized that the transfer does not add any interpolation errors to the correction c as shown on Figure 11.It means that RMT has problem-independent prolongation operator.Smooth parts of the error are deleted on all grids of the next finer levels in the same manner (computation of the coefficient matrix and right-hand side vector and the smoothing iterations).The coarse grid correction to be added to û on the finest grid is c ( ˆû u c   ).The multigrid iterations repeatedly improve the approximation to the solution û until the current approximation becomes accurate enough.In particular applications the finest grid should be reconstructed after each multigrid iteration for their adaptation to the solution singularities.
In RMT more computational work must be spent on coarse grids in order to allow for the best approximation to the solution on the finest grid.Smoothing before coarse grid correction (pre-smoothing) is deleted in RMT to simplify solution of nonlinear problems (like a sawtooth cycle in classical multigrid).It is clear that RMT takes the intermediate place between classical and cascadic multigrid algorithms.
For nonlinear problems when no a priori information of the solution is available to assist in the choice of the initial guess on the finest grid, it is obviously wasteful to start the computation of the finest grid as shown on Figure 12.It is similarly to nested iterations in classical multigrid methods.Main difference consists of formulation of discrete problems of the coarse grids.Integration of the  -modified X-momentum (10b) over the control volume on some coarse grid gives where convective fluxes are given by and the source term is an averaged residual of the X-momentum equation.
Indices w,e,n and s denote values on western, eastern, northern and southern faces of the control volume.The convection terms with the exception of the coefficients , are average values of approximations to the velocity components on the control volume faces.To overcome problem of robustness, all integrals should be evaluated on the finest grid.As an example, we consider evaluation of the integral . The function ( ) q x can be redefined as 0, ( ) ( ), Further we define a characteristic function ( )  is the number of control volumes on the finest grid forming the real part of the control volume on the given grid.Computation of is computed by some method (trapezoids, rectangles etc.), and  are computed by the recursive equation Scheme of the integral evaluation on coarse grids 1 1 G and 2 1 G is shown on Figure 13.Virtual nodes and faces on each grid are intended only for the computational procedure.
-0.5 0 0.5 1 1.5 As an example, we consider computation of the average value on the control volume faces, where ˆ( , ) x y u x y e   .Figure 14 represents the finest uniform staggered grid and location of the control volume, the function û is assigned to the grid points  .The average value on the control volume face on the finest grid can be computed as  ).Note that the computational procedure is based on the additivity property of the integrals and, as a result, the procedure is problem-independent component of RMT.
Source term ( 17) is computed in the same manner.At first, the source term is computed on the finest grid as Values of the approximations to the velocity components û and  on the control volume faces can be approximated by high order schemes.Such different representation of the convection terms in both sides of the modified momentum equations is similar to the defect correction procedure Hackbusch (1981).Finally the source term on the coarse grids is evaluated as a double integral.The abovementioned averaging is a restriction operator of RMT.

Numerical Experiments
In the application of the algorithm to the cavity problem shown on Figure 2, a number of flow Reynolds numbers and computational grid nodes have been considered.Computations have been made for Reynolds numbers of 100 and 500 with computational grids consisting of 101´101 and 1001´1001 nodes.The residuals are defined as  and R  are given by Eq. ( 11).In addition, the residuals of the mass conservation equations (8) are defines as The auxiliary problem ( 14) has been used for fast formulation of the starting guess for the solution of the full Navier-Stokes equations.Since all grid of the same level have no common nodes and control volume faces, discrete equations on all grids of the same level can be written as ) where l A is a block diagonal matrix, the number of blocks is equal to the number of grids forming the level (i.e.   3 , 0,1, 2,..., Nl l L).
Transfer to the level consisting of the finer grids is abbreviated as Accordingly to Figure 11, the permutation matrix   1 l l P corresponds to the problem-independent prolongation operator.
Smoothing iteration for (19) can be written as where matrix l W defines the smoothing procedure and  l is number of the smoothing iterations.
Convergence analysis of RMT is based on the following assumptions: Assumption 1. Smoothing procedure is convergent iterative method, i.e.
where l S is a matrix of the smoothing iterations.
Assumption 2. Eq. ( 20) is solved exactly on the coarsest grids ( Eq. ( 20) can be written as where l c is an exact solution of Eq. ( 19) Since the starting guess  (0 l c is obtained by "prolongation" of approximation to the solution from coarse grids , where Finally the multigrid iterations of RMT is written as ( 1) ( ) 1 0 0 0 ( ) where the matrix of the multigrid iterations takes the form To estimate the matrix norm, it is used the assumption: Assumption 3 (approximation property) where the constant A C is independent on l .
Accounting the approximation property and Assumption 1, the matrix norm can be estimated as , where The estimation shows that convergence of the multigrid iterations  ( 1 ) M can be obtained by performing sufficient number of the smoothing iterations on the multigrid structure.In addition, smoothing on the finer grids has more influence on the convergence of RMT.
Note the estimation does not account the smoothing properties of used iterative methods.As a result, convergence rate of RMT cannot be studied theoretically using the estimation.Computational cost of the multigrid iterations is proportional to N lg N arithmetical operations, where N is the number of the finest grid nodes Martynenko (2006).The problem-independent prolongation operator of RMT leads to the multiplier lg N .In sense of the computational work RMT loses to the problem-adapted classical multigrid methods required N arithmetical operations for solving the boundary value problems.

Conclusions
Proposed algorithm based on the problem-independent components (Vanka smoother, pressure decomposition and robust multigrid technique) makes it possible to solve many boundary value problems for the Navier-Stokes equations to within truncation error at a cost of NlgN c  arithmetic operations, where N is the number of unknowns and c is a constant which depends on the problem.

Figure 1 .
Figure 1.Stencils and the control volumes for discretization of the Navier-Stokes equations on the staggered grids: X-momentum (left), Y-momentum (centre), continuity equation (right) Vanka-like procedures solve for four velocity components

Figure 5 .G
Figure 5. Staggered grid in X direction Figure 6.Coarsening in RMT Properties 1, 2, and 4 give efficient parallelism, problem-independent prolongation operator and problem-independent restriction operator, respectively.

Figure 7 .
Figure 7.Control volumes on the finest and coarse grids Figure 8. Multigrid structure

Figure 9 .
Figure 9. Coarse grids of the first level and the index mapping

Figure 12 .
Figure 10.Multigrid cycle for the linear problems

Figure 14 .
Figure 14.Finest staggered grid (left) and computation of the average value coarse grid points can be located outside the domain [0finest and coarse grids is performed as Eqs.(16) and (17), respectively.Error of the computation is given = 0 ErrMAX = .388E-06whereLevelX ( x l ) and LevelY ( y l ) are numbers of grid levels in X and Y directions, respectively.It easy to see that the error of the integral computation weakly depends on the mesh size on the coarse grid ( Figure 15.Convergence history of RMT l a restriction operator, i.e. the operator transfers the residual  it was shown that the operator is problem-independent component of RMT. results in the recursive form of Eq. (21) grid visiting order corresponds to order shown on Figure10.The recursive Eq. (22) should be written for each level of the multigrid structure.Accounting the Assumption 2, Eq (22) takes the form finer grids, we obtain the following form of Eq. (22)