Convergence Acceleration of Iterative Algorithms for Solving Navier–Stokes Equations on Structured Grids

Basic tendency in computational fluid dynamics (CFD) consists in development of black box software for solving scientific and engineering problems. Numerical methods for solving nonlinear partial differential equations in black box manner should satisfy to the requirements: a) the least number of the problem-dependent components b) high computational efficiency c) high parallelism d) the least usage of the computer resources. We continue with the 2D (N = 2) Navier–Stokes equations governing flow of a Newtonian, incompressible viscous fluid. Let Ω ∈ RN be a bounded, connected domain with a piecewise smooth boundary ∂Ω. Given a boundary data, the problem is to find a nondimensional velocity field and nondimensional pressure such that: a) continuity equation ∂u ∂x + ∂v ∂y = 0 , (1)


Introduction
Basic tendency in computational fluid dynamics (CFD) consists in development of black box software for solving scientific and engineering problems.Numerical methods for solving nonlinear partial differential equations in black box manner should satisfy to the requirements: a) the least number of the problem-dependent components b) high computational efficiency c) high parallelism d) the least usage of the computer resources.We continue with the 2D (N = 2) Navier-Stokes equations governing flow of a Newtonian, incompressible viscous fluid.Let Ω ∈ R N be a bounded, connected domain with a piecewise smooth boundary ∂Ω.Given a boundary data, the problem is to find a nondimensional velocity field and nondimensional pressure such that: a) continuity equation ∂u ∂x b) X-momentum Reynold number Re is defined as where ρ and μ are density and viscosity, respectively.Choice of the velocity scale u s and geometric scale l s depends on the given problem.
Equations ( 1)-( 3) can be rewritten in the operator form where N is nonlinear convection-diffusion operator, F and G are source terms, V and P are velocity and pressure, respectively.It is assumed that the operator N accounts boundary conditions.Note that 2D and 3D Navier-Stokes equations can be written as equation ( 4), where first and second equations abbreviate momentum and continuity equations.Linearized discrete Navier-Stokes equations can be written in the matrix form in which α and β represent the discrete velocity and discrete pressure, respectively.
Here nonsymmetric A is a block diagonal matrix corresponding to the linearized discrete convection-diffusion operator N .The rectangular matrix B T represents the discrete gradient operator while B represents its adjoint, the divergence operator.
Large linear system of saddle point type ( 5) cannot be solved efficiently by standard methods of computational algebra.Due to their indefiniteness and poor spectral properties, such systems represent a significant challenge for solver developers Benzi et al. (2005).Preconditioned Uzawa algorithm enjoys considerable popularity in computational fluid dynamics.The iterations for solving the saddle point system (5) are given by ⎧ ⎨ ⎩ Aα (k+1) = −B T β (k) + f where the matrix Q is some preconditioner.Preconditioned Uzawa algorithm (6) defines the following way for improvement of the solvers for the Navier-Stokes equations: 1) development of numerical methods for solving the boundary value problems.Uzawa iterations require fast numerical inversion of the matrices A and Q.Now algebraic and geometric multigrid methods are often used for the given purpose Wesseling (1991).Multigrid methods give algorithms that solve sparse linear system of N unknowns with O(N) computational complexity for large classes of problems.Variant of geometric multigrid methods with the problem-independent transfer operators for black box or/and parallel implementation is proposed in Martynenko (2006;2010).
Error vector in Uzawa iterations satisfies to the condition

Remarks on solvers for simplified Navier-Stokes equations
Limited characteristics of the first computers and absence of efficient numerical methods put difficulties for simulation of fluid flows based on the full Navier-Stokes equations.As a result, computational fluid dynamics started from simulation of the simplest flows described by the simplified Navier-Stokes equations.
As an example, we consider 2D laminar flow between parallel plates.Figure 1 represents geometry of the problem.Assuming that the pressure is not changed across the flow (p ′ y = 0 in case of L ≫ 1), full Navier-Stokes equations can be reduced to the simplified form: a) X-momentum and mass conservation equations b) continuity equation (1).Since the mass conservation equation follows from the continuity equation (1), system (7) must be solved first.Solution of system (7) gives velocity components u and pressure p.A f t e r that the continuity equation ( 1) is used for determination of v.The computations are repeated until the convergent solution will be obtained.
Let us consider solution of system (7) in details.Assume that an uniform computational grid (h = h x = h y ) is generated.Linearized finite-differenced equations with block unknowns where is the given inlet mass flow rate and superscript n denotes time layer.Missing the superscript (n + 1), the system (8) can be rewritten in the matrix form Comparison of systems ( 9) and ( 5) shows that solution of the simplified Navier-Stokes equations ( 7) also is reduced to solution of the saddle point system.The principal difference between systems ( 9) and ( 5) consists in size of zero block in the coefficient matrix.Since the zero block in system ( 9) has the least size 1 × 1 because of the pressure is independent on y, efficient iterative algorithms for solution of system ( 9) have been proposed and developed.The most promising of them is secant method Briley (1974), where error of the mass conservation equation Hydrodynamics -Optimizing Methods and Tools www.intechopen.com is used for computation of pressure by the iterative method where superscript k denotes the secant method iterations.Note that the approach requires two starting guesses p (0) i and p (1) i .First starting guess can be obtained by extrapolation.For example, for uniform grid we obtain p (0) i = 2p i−1 − p i−2 and compute F (0) .Second starting guess can be given by perturbation of the first one, for example p i .I tgiv esF (1) .Function F depends almost linearly on p (n+1) i , but the secant method is direct solver for linear problems.Usually it is required several secant iterations to reduce error of the discrete mass conservation equation down to roundoff error.Note that in 2D case the system ( 9) can be solved by direct methods, i.e. without the secant iterations.However in 3D case the direct methods require unpractical computational efforts due to five-diagonal structure of the coefficient matrix.As contrasted to the Uzawa algorithm (6), the method does not require some preconditioner(s), relaxation parameter(s), extra computer memory and has high convergence rate.Unfortunately, basic assumption p = p(t, x) does not allow apply the method directly for solving full Navier-Stokes equations ( 1)-(3).Accounting the attractive properties, the algorithm for solving the simplified Navier-Stokes equations can be used for convergence acceleration of the iterative methods intended for full Navier-Stokes equations.Reduction of system (5) to the saddle point system with zero block of the least size is popular approach in CFD.For example, similar reduction based on special unknown ordering is used in Vanka smoother Vanka (1986).

Principle of formal decomposition of pressure
In order to apply the abovementioned approach for solving full Navier-Stokes equations, it is necessary artificially extract «one-dimensional parts of pressure» from the pressure field.For the given purpose, let add and subtract items p x (t, x), p y (t, y) and p z (t, z) depending only on one spatial variable, i.e. p(t, x, y, z)=p x (t, x)+p y (t, y)+p z (t, z) (t, x, y, z) , where superscripts x, y and z denote dependence of the functions on the spatial variables.Let us introduce a new function (t, x, y, z).
Finally the pressure can be represented as p (t, x, y, z) 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 component of pressure» p x (t, x), p y (t, y) or p z (t, z) is dominant.In this case impressive reduction of computational work is expected as compared with traditional algorithms (i.e.p x (t, x)=p y (t, y)=p z (t, z)=0).However for 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» in equation ( 10) are computed only in coupled manner.
Velocity components and «multidimensional component» p xyz (t, x, y, z) in equation ( 10) can be computed in decoupled (segregated) or coupled manner.Remark 7. Gradients of the «one-dimensional components» can be obtained in analytical form for explicit schemes.Implicit schemes require formulation of an auxiliary problem for determination of gradients of the «one-dimensional components».

Development of explicit schemes
First, consider modification of the explicit schemes using well-known benchmark problem about rotated flow in a driven cavity (Figure 2).Let a staggered grid with grid spacing h x and h y has been generated.Classical three-stage splitting scheme is represented as where h t is time semispacing, V (n+1/2) is intermediate velocity field and n is a time layer.Stage I consists in solution of the momentum equations without pressure gradients.For simplicity X-momentum can be written as where ψ ij is the given function defined as It is easy to see that intermediate velocity field V (n+1/2) is independent on pressure.This disadvantage can be compensated partially by the pressure decomposition (10).Application of the decomposition requires two mass conservation equations for 2D problems.Integration of the continuity equation ( 1) over the control volumes V 1 and V 2 shown on Figure 2 gives Approximation of the mass conservation equations on the staggered grid is given by where As contrasted with equation ( 11) in the classical approach, the velocity component u and «one-dimensional component of pressure» p x should satisfy to the system i.e. u and p x are computed in the coupled manner using the discrete mass conservation equation ( 14).
It is clear that the system ( 16) can be written in form of ( 5), where A is the diagonal matrix for explicit schemes.This fact allows obtain analytic solution of the saddle point system ( 16).Multiplication of the first equation in system ( 16) on h y and summation give Left-hand side of the equation equals zero due to equation ( 14).Furthermore Equation ( 17) is reduced to Substitution of the equation into system (16) gives a new form of the system Velocity component u and «one-dimensional component of pressure» p x are computed in the coupled manner saving explicit nature of the computation.Other velocity components are computed in the similar way.
Accounting decomposition (10), other stages of the algorithm are written as Stage II: In the stages only «multidimensional component» p xy is used for computation of the velocity field.
For the numerical experiment law of the lid motion is taken as Reynolds number Re = 1000 is based on the cavity height and the lid velocity max U

Development of implicit schemes
Application of the pressure decomposition (10) for improvement of the implicit schemes requires solution of an auxiliary problem because of the «pressure» gradients can not be determined in explicit form such as equation ( 18).

Auxiliary problem
Auxiliary problem is intended for fast computation of the «one-dimensional components» p x (t, x), p y (t, y) and p z (t, z) in decomposition (10).It is assumed that the solution of the auxiliary problem will be close to the solution of the Navier-Stokes equations.Formulation of the auxiliary problem is based on replacement of the continuity equation ( 1) by the mass conservation equations.For example, for the driven cavity (Figure 2) the auxiliary problem with the mass conservation equations ( 13) instead of the continuity equation (1) takes the form: a) X-momentum and mass conservation equations b) Y-momentum and mass conservation equations where square brackets mean that the «pressure» gradients (p xy ) ′ x and (p xy ) ′ y are fixed (i.e. its values are taken from previous iteration).Braces mean that the momentum and mass conservation equations are solved only in coupled manner.Since the systems ( 20) and ( 21) are similar to the simplified Navier-Stokes equations ( 7), the systems can be solved by the same numerical methods.Main difference consists in stopping criterion: auxiliary problem can be solved approximately, i.e. it is necessary to perform several iterations of line Seidel method with the secant iterations.As a result, extra computational work for approximated solution of the auxiliary problem is negligible small as compared with the total efforts.Note that the equations of the auxiliary problem are not pressure-linked.To illustrate influence of the auxiliary problem on convergence rate, we use Uzawa algorithm (6) for simulation of stationary flow in the driven cavity starting the iterand zero: u (0) = 0, v (0) = 0andp (0) = 0. Accounting zero boundary conditions for v, first equation of system ( 6) is reduced to and v = 0.In the auxiliary problem the system (20) takes the form and v = 0. Finally both problem ( 22) and ( 23) are reduced to systems of linear algebraic equations Ax = b.For clearness these equations are solved until The computations are performed with Re = 100 on uniform staggered grid 101 × 101 (h x = h y = 1/100).Figure 4 represents solution of the Navier-Stokes equations in "stream function-vorticity" (+) Ghia et al. (1982), primitive variables formulations (-) and solutions of equations ( 22) and (23) in the middle section of the cavity (x = 0.5) at Re = 100.It is easy to see that use of the mass conservation equations in the auxiliary problem makes it possible to obtain more accurate approximation to solution of the full Navier-Stokes equations ( 1)-(3).

Main problem
Accounting the pressure decomposition (10), the momentum equations in the main problem are written as where square brackets mean that the «pressure» gradients (p x ) ′ and (p y ) ′ are fixed (i.e. the gradients have been computed in the auxiliary problem (such as equation ( 20) and ( 21) for the driven cavity)).Main problem consists of momentum ( 24), ( 25) and continuity equations (1).Algorithm for simulation of the flows with given mass flow rate can be represented as: Stage I: auxiliary problem: several iterations of line (2D) or plane (3D) Seidel method with the secant iterations Stage II: main problem: iterations of basic method (SIMPLE, Uzawa or Vanka iterations, etc.) Stage III: check convergence, continue (go to 1) if necessary

Flow over a backward-facing step
The next benchmark problem about backward-facing step flow is used for illustration of the impressive convergence acceleration for the directed fluid flows.Consider the stationary laminar flow over a backward-facing step, which is another well studied test case.Figure 5 shows the geometry of the flow.The fact that the solution of the incompressible Navier-Stokes equations over a backward-facing step at Re = 800 is steady and stable has been confirmed in a number of recent works.
No-slip boundary conditions are imposed on the step and the upper and lower walls, a parabolic velocity u profile is specified at the channel inlet (v = 0), and zero natural boundary conditions (v = 0a n du ′ x = 0) are imposed at the channel outlet.The Reynolds number Re is based on the channel height (H = 1) and the average inlet velocity in the parabolic profile.The channel length is L = 14.Numerical experiments show that execution time can be reduced in ∼ 400 times for the given problem (staggered grid 101 × 1401, unpreconditioned Uzawa algorithm, Re = 800).Figure 6 explains the impressive reduction of the computational efforts.It is easy to see that pressure is changed mainly in x direction except small subdomain near attachment point of bottom eddy (i.e.p(x, y) ≈ p x (x)).Since the «one-dimensional component of the pressure» p x (x) is computed in the auxiliary problem, the proposed algorithm is very efficient for solving the problem.
Table 1 represents comparison of obtained results.

Authors l B l T w T x TL
x TR Nodes Barton (1997) 6.0150 5.6600 -4.8200 10.4800 Gartling (1990) 6.1000 5.6300 -4.8500 10.4800 129681 Gresho et al. (1993) 6.0820 5.6260 -4.8388 10.4648 245760 Gresho et al. (1993) 6.1000 5.6300 -4.8600 10.4900 8000 Keskar & Lin (1999)   Nonuniform staggered grid 385 × 3150 is used for the flow simulation (Re = 350).Figure 9 represents distribution of the stream function and pressure near first column of the needles.Chess order of the needle location results in eddy-free flow inside the microcatalyst.However intensive eddy formation after last column of the needles is observed (Figure 10).10) can be computed by coupled or segregated method using density-based or pressure-based approach.
Consider application of the pressure decomposition for simulation of compressible flow in flat Laval micronozzle.Width of subsonic part of the micronozzle is 1 mm.Grid generation is based on mapping of the non-dimensional physical domain with nonuniform grid onto computational domain (unit square) with uniform grid.Direct (ABCD → Ā B C D) and reverse (ABCD ← Ā B C D) mappings are shown on Figure 11, where the function ϕ(x) describes the micronozzle profile.The mappings can be given by , where (x, y) and ( x, ȳ) are spatial variables in physical and computational domains, respectively.Parameter β > 0 is intended for the grid refinement near solid wall.Jacobian (J) of the mapping is non-singular (J = 0).In addition, J → 1/ϕ( x) at β → 0 for uniform grid in y direction.

Fig. 11. Non-dimensional physical and computational domains
Finally, non-dimensional compressible Navier-Stokes equations in the computational domain are written as where Parameter ǫ is the micronozzle width-to-length ratio.
First mass conservation equation is obtained by integration of the continuity equation as follows In the auxiliary problem for incompressible flows, iterations of line (2D) or plane (3D) Seidel method are stopped then the velocity component satisfies to the mass conservation equation.Computation of compressible flows requires updating of thermophysical properties (density ρ, coefficient of viscosity in Re and heat conductivity coefficient in Pe) using updated pressure in the line or plane.In 3D case values of thermophysical properties of the fluid for X-momentum should be updated using pressure temperature T and equation of state.Here square brackets mean that the pressure components p y , p z and p xyz are fixed.
Figure 12 represents isobars in the Laval micronozzles.It is easy to see that the isobars are almost vertical lines near throat and in supersonic part of the micronozzle.It means that the pressure is changed mainly along the micronozzle axis.In other words, «one-dimensional component of the pressure» p x in decomposition ( 10) is dominant in this problem.For

Pressure decomposition in geometric multigrid methods
Proposed convergence acceleration technique based on the pressure decomposition (10) should be incorporated with well-known algorithms for solving Navier-Stokes equations.Multigrid methods having (almost) optimal convergence rate for many applications seem to be the most promising solvers for many CFD problems.Our purpose is development of multigrid method with the least number of problem-dependent components for using in black box software.

Nonlinear multigrid iterations
To overcome problem of robustness, the Navier-Stokes equations (4) should be adapted for the multigrid algorithm Martynenko (2006).Adaptation of the Navier-Stokes equations (so called Σ-modification) consists in representation of the velocity V and pressure P as sum of where discrete analogues of the functions C V and C P will be coarse grid corrections and discrete analogues of the functions V and P will be approximations to the solutions in the following multigrid iterations.As a result, the Navier-Stokes equations ( 4) can be rewritten in the Σ-modified form where It is clear that main difference between Σ-modified and initial forms of the Navier-Stokes equations consists of the nonlinear convection-diffusion operator (N * in equation ( 26) instead of N in equation ( 4)) and source terms (F * and G * in equation ( 26) instead of F and G in equation ( 4)).Note that Σ-modification does not require some linearization of the Navier-Stokes equations.Therefore modified Navier-Stokes equations with other transport equations can be solved in coupled manner on all coarse grids.For example, 2D Σ-modified Navier-Stokes equations are written as: a) Σ-modified continuity equation where discrete analogues of the functions c u , c v and c p will be coarse grid corrections and discrete analogues of the functions û, v and p will be approximations to the solutions in the following multigrid iterations.
pressure» on the finest grid.In other words solution of the modified auxiliary problem is additional smoothing on the finest grid.Six-level multigrid structure with uniform finest staggered grid 501 × 501 (h t = h x = h y = 1/500) is used for the test.Correction of the «multidimensional component of pressure» is computed using Uzawa algorithm with diagonal preconditioning.
Figures 16-18 show evolution of the flow in the cavity.The main vortex is located near the lid after finish of the lid acceleration.Then the vortex moves to upper right corner under influence of the lid motion.After that the vortex moves along diagonal of the cavity to the center (Figure 16).Motion of the main vortex generates two additional vortices (Figure 17).The first vortex is formed in the lower corner, but the second vortex is formed on right vertical wall of the cavity.Location of the corner vortex is stable because of its motion is limited by the cavity walls.However the wall vortex can moving under influence of the main vortex.Corner and wall vortices agglomerate in common vortex.Agglomerated vortex tends to stable corner position, but a new corner vortex is generated in the left lower corner of the cavity.

Acknowledgements
Work supported by Russian Foundation for the Basic Research (project no.09-01-00151).I wish to express a great appreciation to professor M.P. Galanin (Keldysh Institute of Applied Mathematics of Russian Academy of Sciences), who have guided and supported the researches.

Conclusion
«Part of pressure» (i.e.sum of the «one-dimensional components» in decomposition (10)) can be computed using the simplified (pressure-unlinked) Navier-Stokes equations in primitive variables formulation and the mass conservation equations.«One-dimensional components of pressure» and corresponding velocity components are computed only in coupled manner.As a result, there are not pure segregated algorithms and pure density-based approach on structured grids.Proposed method does not require preconditioners and relaxation Fig. 1.Flow between parallel plates orderingshownonFigure1arewrittenas

Fig. 2 .
Fig. 2. Driven cavity and location of the control volumes V 1 and V 2

) 181 Convergence
Acceleration of Iterative Algorithms for Solving Navier-Stokes Equations on Structured Grids www.intechopen.com

9 Fig. 3 .
Fig. 3. Ratio of execution time at flow simulation in the driven cavity where (p x ) 5 is used for the flow simulation.Ratio 183 Convergence Acceleration of Iterative Algorithms for Solving Navier-Stokes Equations on Structured Grids www.intechopen.com of the execution time T(n) m /T (n) cis used as a criterion of the convergence acceleration, whereT (n) m and T (n) care execution time for abovementioned and classical approaches, respectively.Figure3shows result of the numerical test.Obtained result for n = acceleration efficiency arising at simulation of the rotated flows.

Fig. 4 .
Fig. 4. Distribution of the velocity component u in the middle section of the cavity

Fig. 5 .
Fig. 5. Geometry of problem about the backward-facing step flow

Fig. 13 .
Fig. 13.Starting location of the plunger the given problem, the auxiliary problem makes it possible to compute the most «part of pressure» (i.e.p x (t, x)+p y (t, y)) and corresponding change of thermophysical properties of the fluid based on simplified (pressure-unlinked) momentum equations in primitive variables formulation and mass conservation equations.

193Convergence
Acceleration of Iterative Algorithms for Solving Navier-Stokes Equations on Structured Grids www.intechopen.comtwo functions V = C V + V, P = C P + P , Fig. 15.Multigrid cycles

199Convergence
Acceleration of Iterative Algorithms for Solving Navier-Stokes Equations on Structured Grids www.intechopen.com Source terms in the Σ-modified equations coincide with the