III.1 Solution Algorithms

Up until now, the we have dealt with simplified equations, such as the generalised transport equation. Our solution method is as follows :

  1. Define the geometry of the case
  2. Split the region of flow into cells
  3. Integrate the equations of interest over each cell (discretisation)
  4. Invert the resulting matrix
  5. (optional) Repeat for as many timesteps as necessary

Solving the NSE causes additional problems however. In vector notation, for an incompressible fluid the NSE are

 \~/ .u- = 0 (III.1)
@ux @t +  \~/ .(u-ux) = -1
--
r@p
---
@x + n \~/ 2u x (III.2)
@uy @t +  \~/ .(u-uy) = -1-
r@p-
@y + n \~/ 2u y (III.3)
@uz @t +  \~/ .(u-uz) = -1
--
r@p
---
@z + n \~/ 2u z (III.4)
We will ignore the first of these equations for the moment, and concentrate on the momentum equation(s). These are in transport-equation form, with a source term
  1 @p
- ----
  r@x
and a diffusion term on the r.h.s. However there are two main problems :
  1. The equation is non-linear - we need to know u- to evaluate the transport of ux into the domain
  2. The source term involves p, which is one of the variables we want to solve for

Both of these problems relate to the tangled nature of the NSE. Of the 4 equations making up the NSE (continuity and the 3 components of velocity) all components of velocity appear in all the equations, and the pressure appears in the 3 velocity equations. We cannot evaluate the velocity until we know the pressure, and we cannot find the pressure until we know the velocity. In order to find both, we have to guess one and solve for the other, then go back and correct the first. We might start by guessing the pressure, and use this to get a better estimate of the velocity, then correct the pressure, etc.

In CFD there are 2 basic algorithms for doing this :

We will start with PISO. We note the following :

The PISO algorithm is as follows :

  1. Guess p, flux (use values from the previous timestep)
  2. Use equation (III.2) to find ux, uy etc.
  3. Solve the pressure equation for p
  4. Correct the flux to satisfy continuity

Steps 3 and 4 can be iterated if necessary, but this is usually not necessary. This advances the solution one timestep - the whole procedure is then repeated from 1 to 4 for the next timestep.


PIC
Figure 1: Flow chart for the PISO algorithm

The above algorithm is the one usually used to calculate transient flows. To solve for steady-state problems, one can simply apply this algorithm until the solution is unchanging from one timestep to the next. Alternatively one can use the SIMPLE algorithm, as described below. The NSE for a steady flow can be written as

 \~/ .u- = 0
 \~/ .(u-ux) = -1-
r@p-
@x + n \~/ 2u x
 \~/ .(u-uy) = -1-
r@p-
@y + n \~/ 2u y (III.5)
 \~/ .(u-uz) = -1
--
r@p
---
@z + n \~/ 2u z
These still have the same problems as before, namely they are non-linear and the velocity and pressure are so closely entangled. We can adopt a similar procedure to that used in the PISO algorithm, starting by finding approximate values p* and u
--*, and then finding corrections p', u
--' to update these.
  1. Guess an initial pressure field p*
  2. Use equation (III.5) to create a velocity field u
--* from this pressure field
  3. Find a correction p' to the pressure field
    p** = p* + p'
  4. Correct the velocity (flux) to obey continuity.

In theory, p and u- should now be the desired solution. In practice, it is necessary to repeat this procedure as an iterative process, largely because of the nonlinear nature of the velocity. In addition, if the new solution p**, u-** is adopted at each step, it turns out that the algorithm becomes unstable. To correct this, a technique known as underrelaxation is used. The updated value pn+1 is taken as

pn+1 =  ap**+  (1 - a)pn
where a is an underrelaxation parameter. At each step of the iteration, the error in the solution (the residual) should decrease - residuals for all variables should be monitored to ensure this.


PIC
Figure 2: Flow chart for the SIMPLE algorithm