CFD-FVM-10: Solution of discretised equations

Introduction

Thomas (1949) developed a technique for rapidly solving tri-diagonal systems that is now called the Thomas algorithm or the tri-diagonal matrix algorithm (TDMA). It is computationally inexpensive and has the advantage that it requires a minimum amount of storage.

The TDMA

For a tri-diagonal system

and are known from boundary conditions. A general form can be

The tri-diagonal system can be written as

These equations can be solved by forward elimination and back-substitution.

Forward elimination

Removing by substitution from above equation to give

We use

then

The procedure can be repeated up to the last equation of the set.

Back-substitution

The general form of recurrence relationship

where

Since the value of is known at boundary location the values for can be obtained in reverse order by means of the recurrence formula.

Application of the TDMA to two-dimensional problems

To solve the system the TDMA is applied along chosen, e.g. north–south , lines. A general form like the above equation can be written as

The right hand side of is assumed to be temporarily known.

Application of the TDMA to three-dimensional problems

For three-dimensional problems the TDMA is applied line by line on a selected plane and then the calculation is moved to the next plane, scanning the domain plane by plane.

The values at W and E as well as those at B and T on the right hand side are considered to be temporarily known.

Closing remarks

TDMA consists of a forward elimination and a back-substitution stage:

• Forward elimination
  • arrange system of equations in the form of
    
  • calculate coefficients , , and
  • starting at calculate and using
    
  • repeat until
• Back-substitution
  • starting at obtain by evaluating
    
  • repeat until

For two- and three-dimensional problems, sweeping from upstream to downstream along the flow direction producing faster convergence than sweeping against the flow or parallel to the flow direction.

Convergence problems can be alleviated by alternating the sweep directions, which is particularly useful in complex three-dimensional recirculating flows, where the dominant flow direction is not known in advance

When overall stability considerations require coupling between the values over the whole calculation domain the TDMA can be unsatisfactory for the solution of discretised equations.

Point-iterative methods

A simple example can be

Rewrite the equations as

Starting from guessed , and , the equation can be solved by iteration. One condition for the iteration process to be convergent is that the matrix must be diagonally dominant.

When general systems of equations are solved it is sometimes necessary to rearrange the equations, but the finite volume method yields diagonally dominant systems as part of the discretisation process, so this aspect does not require special attention.

Jacobi iteration method

In the Jacobi method, the values , etc. on the left hand side at iteration – indicated here by the bracketed superscript – are evaluated by substituting in the right hand side the last known values , etc., which were obtained at iteration .

Naive substitution

The general form of the Jacobi iteration is

Gauss-Seidel iteration method

Use the newly calculated values of current iteration to update following values.

The general form of the Gauss-Seidel iteration is

Relaxation methods

The general form of the Jocabi method can be rewritten as

We try to modify the convergence rate of the iteration sequence by multiplying the second and third terms on the right hand side by the relaxation parameter :

When we choose the procedure is an **under-relaxation** method, whereas is called over-relaxation.

When apply the relaxation method to the Gauss-Seidel iteration, we obtain

Unfortunately, the optimum value of the relaxation parameter is problem and mesh dependent, and it is difficult to give precise guidance.

Multigrid techniques

The convergence rate of iterative methods, such as the Jacobi and Gauss–Seidel, rapidly reduces as the mesh is refined.

Multigrid concept

For the coarse mesh, the longest possible wavelengths of error components (i.e. those of the order of the domain size) are just within the short-wavelength range of the mesh and, hence, all error components reduce rapidly. On the finer meshes, however, the longest error wavelengths are progressively further outside the short-wavelength range for which decay is rapid.

Multigrid methods are designed to exploit these inherent differences of the error behaviour and use iterations on meshes of different size.

An outline of a multigrid procedure

Two-stage multigrid procedure:

1. Fine grid iterations:
   Perform iterations on the finest grid with spacing to generate an intermediate solution to system with true solution vector , we have residual vector . The error vector satisfies .

2. Restriction:
   The solution is transferred from the fine mesh with spacing onto a coarse mesh with spacing , where . Instead of solving for the solution vector we work with the error equation on the coarse mesh starting with an initial guess of

3. Prolongation:
   The error solution is transferred back to the fine mesh with spacing . The transfered error solution is denoted by .

4. Correction and final iteration:
   Correct the intermediate fine grid solution .

The above description is for the two-stage procedure (one fine mesh, one coarse mesh). In practice, however, the restriction is carried out into a number of increasingly coarse levels. Then prolongation procedures are also performed at each stage back to the starting mesh.

An illustrative example can be found in the book.

Multigrid cycles

In practical CFD calculations the multigrid transfer process is more sophisticated and different cycles of coarsening and refinement are used with special schedules of restriction and prolongation at different refinement levels. Common choices of multigrid cycles are the so-called V-, W- and F-cycles

In the technique known as the full multigrid (FMG) method the calculations do not start at the finest grid, but instead at the coarsest level. Solutions are transferred to successively finer grid levels and on the finest level the prolonged solution is used as the initial guess for start of the iterative process.

Grid generation for the multigrid method

For the coarse grid system matrix , if the matrix is computed using actual geometrical properties of the coarse mesh, this type of multigrid procedure is called a geometric multigrid procedure.

If the linear combinations of coefficients of the fine grid equations are used, such multigrid methods are called algebraic multigrid and are widely used in commercial CFD solvers.

Summary

Multigrid acceleration of the Gauss–Seidel point-iterative method is currently the solution algorithm of choice for commercial CFD codes. The rate of convergence of this procedure can be optimised by specific choices for:

• interpolation of the residual vector and coefficient matrix from fine to coarse meshes during restriction,
• interpolation of the error vector from coarse to fine meshes during prolongation,
• cycles of coarsening and refinement with special schedules of restriction and prolongation at different refinement levels.