How to solve incompressible and compressible flows in CFD

This article will take a deep dive into how to solve the Navier-Stokes equations numerically. We will look at the differences between compressible and incompressible flows and what solution algorithms exist to solve them. Along the way, we will also look at the difference between the conservative and primitive (or non-conservative) variable formulation. We will see (through an example and a derivation) why we really want to use conservative variables for compressible flows while primitive variables are preferred for incompressible flows.

And then, there is the scientific misconduct that no one is talking about; the SIMPLE algorithm. We will discuss the SIMPLE algoirithm, and its friends (SIMPLEC, PISO, PIMPLE), and why these really are just inferior and plagiarised description of the pressure projection method of Chorin. We will derive the equations from scratch, and by the end of the article, you will be able to see how all of these methods are one and the same. We have a wild journey ahead of us, did you bring a hot cup of tea? You will need it! Let’s get going.

In this series

• Part 1: How to Derive the Navier-Stokes Equations: From start to end
• Part 2: What are Hyperbolic, parabolic, and elliptic equations in CFD?
• Part 3: How to discretise the Navier-Stokes equations
• Part 4: Explicit vs. Implicit time integration and the CFL condition
• Part 5: Space and time integration schemes for CFD applications
• Part 6: How to implement boundary conditions in CFD
• Part 7: How to solve incompressible and compressible flows in CFD
• Part 8: The origin of turbulence and Direct Numerical Simulations (DNS)
• Part 9: The introduction to Large Eddy Simulations (LES) I wish I had
• Part 10: All you need to know about RANS turbulence modelling in one article
• Part 11: Advanced RANS and hybrid RANS/LES turbulence modelling in CFD

In this article

Introduction

In our journey through key concepts in CFD, we have looked at how to derive the Navier-Stokes equation, how to discretise it, how to use numerical schemes to obtain unknown quantities, and how to apply boundary conditions. You might reasonably assume that this is it; we have all the concepts together, and we can start implementing these equations and start solving them. There is just one small problem: the Navier-Stokes equations are under-constrained.

What do I mean by that? Well, let’s look at the equations. For the general Navier-Stokes equations, we have the following continuity or mass conservation equation:

\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{u})=0

In this equation, we have the density \rho, as well as the velocity vector \mathbf{u}, which contains three velocity components u, v, and w (or u_x, u_y, and u_z). Thus, we have 4 unknowns in this equation. If we now also consider the conservation of momentum equation, we have:

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \nu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\nu\nabla(\nabla\cdot\mathbf{u})

This equation is solved for \mathbf{u}, so this equation needs to be solved for each velocity component u, v, and w. Great, now we have 3 additional equations, and together with the continuity equation, we have now 4 equations. But, we also have introduced the pressure gradient, and thus the pressure p is another quantity we don’t know. So even if we now have 4 equations, we have a total of 5 unknowns at this point, that is, the density \rho, the velocity components u, v, and w, and the pressure p.

So, we need to figure out an additional equation that ties all of these variables together in some way. And this is what we will discuss in the article today. The Navier-Stokes equations can’t be solved as they are, and we need to augment them with some equations to get a solution. We also saw in the article on what are hyperbolic, parabolic, and elliptic flows, that the flow fundamentally changes whether we deal with subsonic or supersonic flows. That is, if the Mach number is below or above one, respectively.

If the character changes, then it stands to reason that the underlying equations also have to adapt accordingly. Thus, we have to find now two separate types of equations. The way that we do that in fluid dynamics is by considering the change in density. The following image shows how much the density changes for air as a function of the Mach number.

Thus, for flows which have a maximum Mach number below Ma<0.3, we can see that the density changes by at most 5%. We simplify this and say that the density does not change below a Mach number of Ma<0.3. This is our demarcation. Anything to the left of Ma<0.3 is treated as an incompressible flow, and anything to the right is treated as a compressible flow.

We can generalise this and say that an incompressible flow is one in which the density does not change, while the density is allowed to change for compressible flows. This definition is independent of the Mach number, but the previous figure nicely demonstrates this for air.

We can formalise this with a property called the compressibility \beta, and this is given as

\beta=-\frac{1}{V}\frac{\partial V}{\partial p}

Here, V is the volume and p is the pressure. Thus, we measure how much the volume changes as a result of increasing the pressure. If you take a softball and you apply pressure, for example, you place it between your hands and you squeeze it, it will get smaller (its volume will decrease). So, we see that the gradient \partial V/\partial p must be negative. We normalise this by the volume itself, and we multiply this by -1 so that our compressibility factor \beta is positive.

If we look at materials such as water, we have \beta\approx 4.6\cdot 10^{-10} at room temperature. Compare that to air, which has \beta\approx 1 at room temperature. A material with a low compressibility will not show any compressible effects, while a material with high compressibility will exhibit compressible effects.

As a result, water, for example, will essentially behave like an incompressible fluid. If we tried to simulate water with a compressible solver, we would see painfully slow convergence. This is because a compressible solver, as we will see, requires density changes to update the remaining flow variables. If the compressibility is low, the density changes are slow, and thus convergence is slow. An incompressible solver ignores density changes and thus is able to converge much faster at low Mach numbers (or flows with little to no density changes).

There is one main modelling difference between compressible and incompressible flows, and that is how variables are written in our derivatives. We distinguish between primitive and conserved variable formulations, and because this is another topic not getting the attention in the literature it deserves, we’ll need to talk about this topic first before we can look at the different techniques to solve compressible and incompressible flows.

Primitive vs. conserved variable formulations

To see the difference between the primitive and conservative variable formulation, let’s start by looking at the momentum equation of the Navier-Stokes equations. We have:

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \nu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\nu\nabla(\nabla\cdot\mathbf{u})

Now, we simplify this equation and say that we have an inviscid fluid (that is, \nu=0), and we also do not consider the pressure. This means the entire right-hand side goes to zero, and we are left with:

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = 0

If we now assume a one-dimensional flow, then we are left with only the velocity component in the x direction. This can be written as:

\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = 0

This equation is known as the Burgers equation, named after the Dutch physicist Johannes Martinus Burgers, and it is an important model equation in CFD. It allows us to study how numerical schemes behave to non-linearities, as the second term is fully non-linear. We could create a linear equation out of this by saying that the u in front of the space derivative is a constant and set u=a. If we were to do that, then we have obtained the linear advection equation, which is another important model equation in CFD.

The Burgers equation above is written in non-conservative, or primitive variable, form. The primitive variables are the ones we are solving for. In this case, the Burgers equation is solving for the variable u as indicated in the time derivative, but we also see that all other derivatives (in this case, just one, i.e. the second term with the spatial derivative) evaluating derivatives of u. When all derivatives solve for the same variable, we say the equation is written in primitive variable formulation.

The origin of why we call them primitive variables comes from various directions, but the first reported source is from a 1921 paper by Bjerknes (“The Meteorology of the Temperate Zone and the General Atmospheric CIRCULATION. 1.” Monthly Weather Review 49, no. 1 (1921): 1-3) where he was using primitive equations to solve meteorological equations. Back in 1921, CFD was not an established field and much of what would later form the body of CFD was developed by meteorologists.

What happens if we write the spatial derivative in the Burgers equation, i.e. u(\partial u/\partial x) as a fully non-linear derivative, i.e. \partial uu/\partial x=\partial u^2/\partial x? Well, let’s see. The product rule of differentiation tells us that the derivative of a product can be expressed as:

\frac{\partial ab}{\partial x}=a\frac{\partial b}{\partial x}+b\frac{\partial a}{\partial x}

In our case, both a and b are equal to u, so we would have:

\frac{\partial uu}{\partial x}=\frac{\partial u^2}{\partial x}=u\frac{\partial u}{\partial x}+u\frac{\partial u}{\partial x}=2u\frac{\partial u}{\partial x}

So we can’t simply write u(\partial u/\partial x)=\partial u^2/\partial x, this would be missing a factor of two. Instead, we have to divide the equation by this factor of two, which would result in the following equality:

\frac{1}{2}\frac{\partial u^2}{\partial x}=u\frac{\partial u}{\partial x}

We can now change the spatial derivative in our Burgers equation and arrive at:

\frac{\partial u}{\partial t}+\frac{1}{2}\frac{\partial u^2}{\partial x}=0

Now we have the velocity component u in the time derivative, that is the primitive variable we are solving for, and we have u^2 in the derivative. The derivative is not multiplied by any other variable. Instead, it contains all variables in the derivative. We could say it has conserved all variables. Therefore, we call this approach the conserved variable approach.

This becomes even more clear when we consider the Navier-Stokes equation again. This was given as:

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \nu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\nu\nabla(\nabla\cdot\mathbf{u})

This equation is given in primitive variable or non-conservative form. Let’s write this equation in conservative form. First, let’s multiply by the density. This results in

\rho\frac{\partial \mathbf{u}}{\partial t} + (\rho\mathbf{u}\cdot\nabla)\mathbf{u} = -\nabla p + \mu\nabla^2\mathbf{u} + \mu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\mu\nabla(\nabla\cdot\mathbf{u})

Here, we are using the fact that \nu=\mu/\rho. We can use the same product rule to write the time derivative in conservative form:

\frac{\partial \rho\mathbf{u}}{\partial t}=\rho\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\frac{\partial \rho}{\partial t}

The continuity equation is given as

\frac{\partial \rho}{\partial t}+\nabla\cdot (\rho\mathbf{u})=0

We can solve this for the time derivative and then insert in the above equation, which provides:

\frac{\partial \rho\mathbf{u}}{\partial t}=\rho\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\left[-\nabla\cdot (\rho\mathbf{u})\right]

If we solve this for the first derivative on the right-hand side, we obtain:

\rho\frac{\partial \mathbf{u}}{\partial t} =\frac{\partial \rho\mathbf{u}}{\partial t}+\mathbf{u}\nabla\cdot (\rho\mathbf{u})

Inserting this into our Navier-Stokes equations, we obtain:

\frac{\partial \rho\mathbf{u}}{\partial t}+\mathbf{u}\nabla\cdot (\rho\mathbf{u}) + (\rho\mathbf{u}\cdot\nabla)\mathbf{u} = -\nabla p + \mu\nabla^2\mathbf{u} + \mu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\mu\nabla(\nabla\cdot\mathbf{u})

Now, let’s investigate the second and third terms on the left-hand side. If we set a=\rho\mathbf{u} and b=\mathbf{u}, then we have

\mathbf{u}\nabla\cdot (\rho\mathbf{u}) + (\rho\mathbf{u}\cdot\nabla)\mathbf{u}=\mathbf{b}\nabla\cdot \mathbf{a} + (\mathbf{a}\cdot\nabla)\mathbf{b}=\nabla\cdot(\mathbf{a}\otimes \mathbf{b})=\nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u})

The additional term from the time derivative, after writing this term in conservative form, combines with the non-linear term. We can write this non-conservative form, i.e. (\rho\mathbf{u}\cdot\nabla)\mathbf{u} in conservative from, i.e. \nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u}) with the help of this additional term.

Thus, our final Navier-Stokes equation can be written as:

\frac{\partial \rho\mathbf{u}}{\partial t} + \nabla\cdot (\rho\mathbf{u}\otimes\mathbf{u}) = -\nabla p + \mu\nabla^2\mathbf{u} + \mu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\mu\nabla(\nabla\cdot\mathbf{u})

Now, let’s inspect the time derivative. We are no longer simply solving for the primitive variable \mathbf{u}, but instead we solve now for \rho\mathbf{u}. This is the momentum, and the equation above is the conservation law for momentum. As a result, we call the above formulation the conserved variable approach, where we now solve for the quantity that is being conserved (here, the momentum) through the conservation law.

The next question that you will be asking is, which of these forms is to be preferred? Well, let’s return to the Burgers equation, and let’s discretise this equation with a first-order upwind scheme using the finite difference method (we’ll assume the velocities are always positive, which makes our discussion easier to follow). Then, both conservative and non-conservative (primitive variable) formulations result in:

• Conservative form:

\frac{u^{n+1}_i-u^n_i}{\Delta t}+\left[\frac{u^2_i-u^2_{i-1}}{\Delta x}\right]^n=0

• Primitive variable form:

\frac{u^{n+1}_i-u^n_i}{\Delta t}+u^n_i\frac{u^n_i-u^n_{i-1}}{\Delta x}=0

Let’s see how both of them compare. We will solve them on a domain that is going from 0 to 1, and we will test this equation for two separate initial conditions. The first one is a smooth profile, where we use a sine wave that has exactly one period between 0 and 1. The animation of this simulation is shown below:

As you can see, there are some slight differences between both formulations, but overall, they give the same result. If anything, the conservative formulation has a slightly better resolution near the discontinuity once that forms. So, for smooth profiles, either the conservative or non-conservative approach works well.

Note that the above velocity profile contains negative velocities, so our discretisation will not work here. To ensure we have only positive velocities, we either have to shift the initial data so that all values of u are positive or use backward differencing if the velocities are negative. I have used the following (Python-based) discretisation to obtain the plot above:

# conservative
if u_old[i] >= 0.0:
  u[i] = old_u[i] - dt * 0.5 * (pow(old_u[i], 2) - pow(old_u[i-1], 2))/dx
else:
  u[i] = old_u[i] - dt * 0.5 * (pow(old_u[i+1], 2) - pow(old_u[i], 2))/dx
  
# non-conservative (primitive variable approach)
if old_u[i] >= 0.0:
  u[i] = old_u[i] - dt * old_u[i] * (old_u[i] - old_u[i-1])/dx
else:
  u[i] = old_u[i] - dt * old_u[i] * (old_u[i+1] - old_u[i])/dx

But let’s see what happens if we are dealing with discontinuities. We consider an initial profile that is 1 between 0 and 0.25 and 0 everywhere else. Let’s look at the comparison of the primitive and conserved variable formulation:

This is the correct solution and not a programming error! If you want to have a mathematical derivation, the book by LeVeque has some more detail, but let me give you an intuitive description. Let’s look at what happens around the x location of x=0.25 at the first time step. This is schematically shown below:

Let’s say we are currently at location i, that is, all values before are set to u=1 and all values to the right are set to u=0, including the node we are currently at. This means we have u_{i-1}=1, u_i=0, and u_{i+1}=0. Now let’s insert that into our spatial derivative for both the conservative and non-conservative (primitive variable) formulation. This results in:

• Conservative form:

\frac{\partial u^2}{\partial x}\approx\frac{u_i^2-u_{i-1}^2}{\Delta x}=\frac{0^2-1^2}{\Delta x}=-\frac{1}{\Delta x}

• Non-conservative (primitive variable) form:

u\frac{\partial u}{\partial x}\approx u_i\frac{u_i-u_{i-1}}{\Delta x}=0\cdot\frac{0-1}{\Delta x}=0\cdot \frac{-1}{\Delta x}=0

We are getting very different results near the discontinuity! In this case, the conservative formulation produces the correct physical behaviour, and, as a result, it is the preferred formulation for compressible flows. This ensures that when we deal with discontinuities, we are capturing them without problems. This doesn’t mean we can’t use the primitive variable formulation for compressible flows (we can, and people do), but we may need to take extra steps to ensure discontinuities are resolved.

So, from this discussion, it is clear that the conservative form is preferred for compressible flows. But what about incompressible flows? Let’s take a look. Let’s write the Navier-Stokes equations out again. We have:

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \nu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\nu\nabla(\nabla\cdot\mathbf{u})

We also have the continuity equation given as:

\nabla \cdot\mathbf{u}=0

We can use this to eliminate the last two terms in the momentum equation, which results in:

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}

Let’s multiply by the density first, this results in:

\rho\frac{\partial \mathbf{u}}{\partial t} + (\rho\mathbf{u}\cdot\nabla)\mathbf{u} = -\nabla p + \mu\nabla^2\mathbf{u}

We can now write the time derivative again in conservative form, this results in:

\frac{\partial \rho \mathbf{u}}{\partial t}=\rho\frac{\partial\mathbf{u}}{\partial t}+\mathbf{u}\frac{\partial \rho}{\partial t}

An incompressible flow is one in which the density does not change. In other words, the density is constant and thus, the derivative of the density must be zero, both in space and time. This means we have:

\frac{\partial \rho \mathbf{u}}{\partial t}=\rho\frac{\partial\mathbf{u}}{\partial t}

We can switch between both conserved and primitive variable formulation without incurring additional terms, thanks to the continuity equation. What about the non-linear term? Well, we saw before that we can write:

\mathbf{u}\nabla\cdot (\rho\mathbf{u}) + (\rho\mathbf{u}\cdot\nabla)\mathbf{u}=\nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u})

But we can use the continuity equation again here to eliminate the first term (\nabla\cdot (\rho\mathbf{u})=0). This results in:

(\rho\mathbf{u}\cdot\nabla)\mathbf{u}=\nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u})

Again, we see that both the conserved variable and primitive variable form are equivalent in a mathematical sense. If we further consider that our flow is typically smooth and does not have any discontinuities (they do exist in some cases, for example, at interfaces between two different fluids), then using the primitive variable formulation makes it easier to work with the equations.

While some researchers have successfully used the conservative variable formulation for incompressible flows, most solvers make use of the primitive variable formulation, i.e. the opposite of what compressible solvers are doing. For this reason, it seems to me that people generally split into either incompressible or compressible flows in terms of their specialisation (rarely both), as the mathematical and numerical treatment is slightly different for both.

Hopefully, you now have a much better idea of why we have primitive and conserved variables, what their mathematical difference are and how that manifests itself in the solution. I have always looked for a clear description on that when I was a student and could not find a good definition, hopefully this will have helped you in clearing this up.

Now that we have discussed the differences let us continue in our discussion and see how we solve the compressible Navier-Stokes equations next.

Compressible flows

As discussed in the previous section, the compressible form of the Navier-Stokes equation is most commonly solved in the conservative form, though some prefer to solve it with primitive variables instead. In this section, we will stick with the convention and look at the conservative form.

The conservation of mass, or continuity equation, for a compressible flow is given as:

\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{u})=0

We derived the compressible form of the Navier-Stokes equations in a previous article, which is given as:

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \nu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\nu\nabla(\nabla\cdot\mathbf{u})

We also saw in the previous section that we can multiply it by the density and then rewrite it in conservative form. This is then given as:

\frac{\partial \rho\mathbf{u}}{\partial t} + \nabla\cdot (\rho\mathbf{u}\otimes\mathbf{u}) = -\nabla p + \mu\nabla^2\mathbf{u} + \mu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\mu\nabla(\nabla\cdot\mathbf{u})

And finally, we need to look at the energy equation. We have not yet derived it in any previous article, and doing so would probably warrant its own article. We’ll use the end result here and likely return to its derivation in a future article. For now, we consider the following equation for the conservation of total energy:

\frac{\partial \rho E}{\partial t}+\nabla\cdot[(\rho E + p)\mathbf{u}]=\mathbf{u} \cdot \left[\mu\nabla^2\mathbf{u} + \mu\nabla(\nabla\cdot\mathbf{u})-\frac{2}{3}\mu\nabla(\nabla\cdot\mathbf{u})\right]-\nabla\cdot\kappa\nabla T

Both the momentum and energy equation contain the viscous stresses. These can be identified as:

\nabla\cdot\tau=\mu\nabla^2\mathbf{u} + \mu\nabla(\nabla\cdot\mathbf{u})-\frac{2}{3}\mu\nabla(\nabla\cdot\mathbf{u})

Here, \tau is the viscous stress tensor. With these equations now written out in full, we can now write out the scalar form of each equation, which we will need to solve in a CFD solver. These are given as follows:

• Conservation of mass (continuity equation):

\frac{\partial \rho}{\partial t} +\frac{\partial \rho u}{\partial x} +\frac{\partial \rho v}{\partial y} +\frac{\partial \rho w}{\partial z}=0

• Conservation of x-momentum:

\frac{\partial \rho u}{\partial t} + \frac{\partial (\rho u^2 + p)}{\partial x} + \frac{\partial \rho u v}{\partial y} + \frac{\partial \rho u w}{\partial z} = \frac{\partial \tau_{xx}}{\partial x} + \frac{\partial \tau_{yx}}{\partial y} + \frac{\partial \tau_{zx}}{\partial z}

• Conservation of y-momentum:

\frac{\partial\rho v}{\partial t} + \frac{\partial \rho u v}{\partial x}  + \frac{\partial (\rho v^2 + p)}{\partial y}  + \frac{\partial \rho v w}{\partial z} = \frac{\partial \tau_{xy}}{\partial x} + \frac{\partial \tau_{yy}}{\partial y} + \frac{\partial \tau_{zy}}{\partial z}

• Conservation of z-momentum:

\frac{\partial \rho w}{\partial t} + \frac{\partial \rho u w}{\partial x} + \frac{\partial \rho v w}{\partial y} + \frac{\partial (\rho w^2 + p)}{\partial z} = \frac{\partial \tau_{xz}}{\partial x} + \frac{\partial \tau_{yz}}{\partial y} + \frac{\partial \tau_{zz}}{\partial z}

• Conservation of total energy:

\frac{\partial \rho E}{\partial t} + \frac{\partial \left[ (\rho E + p) u \right]}{\partial x}  + \frac{\partial \left[ (\rho E + p) v \right]}{\partial y} + \frac{\partial \left[ (\rho E + p) w \right]}{\partial z} = \\[0.5em]
\frac{\partial (u \tau_{xx} + v \tau_{xy} + w \tau_{xz} - q_x)}{\partial x}  + \frac{\partial (u \tau_{yx} + v \tau_{yy} + w \tau_{yz} - q_y)}{\partial y}  + \frac{\partial (u \tau_{zx} + v \tau_{zy} + w \tau_{zz} - q_z)}{\partial z} 

Here, the components of the viscous stress tensor in scalar form are given as:

\tau_{xx} = 2\mu \frac{\partial u}{\partial x} + \lambda (\nabla \cdot \mathbf{v}), \quad
\tau_{yy} = 2\mu \frac{\partial v}{\partial y} + \lambda (\nabla \cdot \mathbf{v}), \quad
\tau_{zz} = 2\mu \frac{\partial w}{\partial z} + \lambda (\nabla \cdot \mathbf{v})\\[1.0em]
\tau_{xy} = \mu \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right), \quad
\tau_{xz} = \mu \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right), \quad
\tau_{yz} = \mu \left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right)

The heat flux \mathbf{q} is given as:

q_x = -\kappa \frac{\partial T}{\partial x}, \quad
q_y = -\kappa \frac{\partial T}{\partial y}, \quad
q_z = -\kappa \frac{\partial T}{\partial z}

And finally, the divergence of the velocity field is given as:

\nabla\cdot\mathbf{u}=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} + \frac{\partial w}{\partial z}

If we look through the equations, then we see that we have the following unknowns: the density \rho, the three velocity components u, v, and w in the x, y, and z direction, the pressure p, the total energy E and the temperature T. This means that we have 7 unknowns. However, we can obtain the temperature from the total energy.

The total energy (per unit mass) is given as:

E=e+\frac{1}{2}(u^2+v^2+w^2)

Here, e is the internal energy. This, in turn, can be expressed as:

e=c_vT

c_v is the specific heat capacity at constant volume, and this expression holds true for an ideal gas. Thus, we can insert this into our definition of the total energy and obtain:

E=c_vT+\frac{1}{2}(u^2+v^2+w^2)

The temperature can now be expressed as a function of the total energy as:

T=\frac{E-0.5(u^2+v^2+w^2)}{c_v}

So, instead of 7 unknowns, we really only have 6 unknowns, or 6 independent variables. If we count the equations, then we see that we have one conservation of mass equation, three conservation of momentum equations, and one conservation of energy equation. This gives a total of 5 equations for 6 unknowns. This means we need an additional equation, which is typically given in the form of the equation of state. For an ideal gas, we have:

p=\rho RT

Where R is the specific gas constant. This approximation works well for low pressures and high temperatures. At high pressures, atoms are forced closer together, showing stronger repulsive/attractive forces, which are not considered by the ideal gas law. At low temperatures (near the gas’ boiling point), a gas will condense into liquid, changing its properties, which isn’t captured by the ideal gas law.

If one of these real gas effects is important, we need to use a more appropriate equation of state. Examples are the Van der Walls, Redlick-Kwong, and Peng-Robinson equation of state.

Now we have the 6 equations for 6 independent variables, and we can solve the equations. The scalar equations derived above are particularly suitable to be solved with the finite difference approximation. We just have to approximate derivatives with a suitable scheme and then solve each equation for the unknown quantities. For finite volume approximations, however, it is more convenient to write the equations in vector form. The general vector form is given by:

\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}^{Inv}}{\partial \mathbf{x}}=\frac{\partial\mathbf{G}^{Visc}}{\partial \mathbf{x}}

Here, \mathbf{U} represents the vector of the conserved quantities, \mathbf{F}^{Inv} the inviscid fluxes, and \mathbf{G}^{Visc} the viscous fluxes. We can expand the vector in the x, y, and z direction as:

\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}^{Inv}_x}{\partial x}+\frac{\partial \mathbf{F}^{Inv}_y}{\partial y}+\frac{\partial \mathbf{F}^{Inv}_z}{\partial z}=\frac{\partial \mathbf{G}^{Visc}_x}{\partial x}+\frac{\partial \mathbf{G}^{Visc}_y}{\partial y}+\frac{\partial \mathbf{G}^{Visc}_z}{\partial z}

The conserved variables and inviscid fluxes and the left-hand side are given as:

\mathbf{U}=\begin{bmatrix}\rho \\[0.5em] \rho u \\[0.5em] \rho v \\[0.5em] \rho w \\[0.5em] \rho E \end{bmatrix},\quad

\mathbf{F}^{Inv}_x=\begin{bmatrix}\rho u \\[0.5em] \rho u^2 + p \\[0.5em] \rho u v \\[0.5em] \rho u w \\[0.5em] u(\rho E + p) \end{bmatrix},\quad

\mathbf{F}^{Inv}_y=\begin{bmatrix}\rho v \\[0.5em] \rho u v \\[0.5em] \rho v^2 + p \\[0.5em] \rho v w \\[0.5em] v(\rho E + p) \end{bmatrix},\quad

\mathbf{F}^{Inv}_z=\begin{bmatrix}\rho w \\[0.5em] \rho u w \\[0.5em] \rho v w \\[0.5em] \rho w^2 + p \\[0.5em] w(\rho E + p) \end{bmatrix}

The viscous fluxes on the right-hand side of the vector equation are given as:

\mathbf{G}^{Visc}_x=\begin{bmatrix} 0 \\[0.5em] \tau_{xx} \\[0.5em] \tau_{xy} \\[0.5em] \tau_{xz} \\[0.5em] u\tau_{xx} + v\tau_{xy} + w \tau_{xz} - q_x \end{bmatrix},\quad

\mathbf{G}^{Visc}_y=\begin{bmatrix} 0 \\[0.5em] \tau_{yx} \\[0.5em] \tau_{yy} \\[0.5em] \tau_{yz} \\[0.5em] u\tau_{yx} + v\tau_{yy} + w \tau_{yz} - q_y \end{bmatrix},\quad

\mathbf{G}^{Visc}_z=\begin{bmatrix} 0 \\[0.5em] \tau_{zx} \\[0.5em] \tau_{zy} \\[0.5em] \tau_{zz} \\[0.5em] u\tau_{zx} + v\tau_{zy} + w \tau_{zz} - q_z \end{bmatrix}

We still need to use the equation of state to close this system of equations, but working with the vector form is easier than with the scalar form. Why is that? Well, in the next step of the finite volume approximation, we need to interpolate values from the cell centroids to the face centres. We saw in the article on the finite volume method that a time derivative can be approximated as:

\int_V\int_t^{t+\Delta t}\frac{\partial \phi}{\partial t}\mathrm{d}t\mathrm{d}V=(\phi^{n+1}-\phi^n)V

On the other hand, a first-order derivative can be approximated as:

\int_t^{t+\Delta t}\int_An\cdot\phi\,\mathrm{d}A\mathrm{d}t=\int_An\cdot\phi\,\mathrm{d}A\Big|_t^{t+\Delta t}=\Delta t\int_An\cdot\phi\,\mathrm{d}A\approx \Delta t\sum_{i=1}^{n_{faces}}n_f\cdot \phi_i A_i

This means that for an arbitrary unstructured grid, we can combine these two approximations with our generalised vector equation form as:

(\mathbf{U}^{n+1}-\mathbf{U}^n)V+\Delta t\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_x\right)_i A_i+\Delta t\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_y\right)_i A_i+\Delta t\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_z\right)_i A_i=\\
\Delta t\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{G}^{Visc}_x\right)_i A_i+\Delta t\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{G}^{Visc}_y\right)_i A_i+\Delta t\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{G}^{Visc}_z\right)_i A_i

We can clean this up a bit by dividing each term by the volume of the cell V and the time step \Delta t. This results in:

\frac{\mathbf{U}^{n+1}-\mathbf{U}^n}{\Delta t}+\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_x\right)_i A_i+\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_y\right)_i A_i+\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_z\right)_i A_i\right]=\\
\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{G}^{Visc}_x\right)_i A_i+\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{G}^{Visc}_y\right)_i A_i+\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{G}^{Visc}_z\right)_i A_i\right]

If we are dealing with a structure grid with regular elements, as shown in the figure below:

where the cell’s centroid is assumed to be at location i,j,k, then we can write the locations on the faces at which the integrals are to be evaluated. We can further state that the volume can be written as V=\Delta x\Delta y\Delta z. Let’s look at the second term in more detail. We can simplify this term as follows, still assuming that we work with the cell shown above:

\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_x\right)_i A_i\right]=\frac{1}{\Delta x\Delta y\Delta z}\left[\begin{bmatrix}1\
\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_x\right)_i A_i\right]=\frac{1}{\Delta x\Delta y\Delta z}\left[\begin{bmatrix}1\\0\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k}\Delta A_{i+1/2,j,k}+\begin{bmatrix}-1\\0\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}\Delta A_{i-1/2,j,k}\right]=\\[1.0em]
\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k}\Delta y \Delta z -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}\Delta y\Delta z}{\Delta x\Delta y\Delta z}=\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k} -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}}{\Delta x}
\
\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_x\right)_i A_i\right]=\frac{1}{\Delta x\Delta y\Delta z}\left[\begin{bmatrix}1\\0\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k}\Delta A_{i+1/2,j,k}+\begin{bmatrix}-1\\0\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}\Delta A_{i-1/2,j,k}\right]=\\[1.0em]
\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k}\Delta y \Delta z -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}\Delta y\Delta z}{\Delta x\Delta y\Delta z}=\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k} -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}}{\Delta x}
\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k}\Delta A_{i+1/2,j,k}+\begin{bmatrix}-1\
\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_x\right)_i A_i\right]=\frac{1}{\Delta x\Delta y\Delta z}\left[\begin{bmatrix}1\\0\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k}\Delta A_{i+1/2,j,k}+\begin{bmatrix}-1\\0\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}\Delta A_{i-1/2,j,k}\right]=\\[1.0em]
\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k}\Delta y \Delta z -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}\Delta y\Delta z}{\Delta x\Delta y\Delta z}=\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k} -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}}{\Delta x}
\
\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_x\right)_i A_i\right]=\frac{1}{\Delta x\Delta y\Delta z}\left[\begin{bmatrix}1\\0\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k}\Delta A_{i+1/2,j,k}+\begin{bmatrix}-1\\0\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}\Delta A_{i-1/2,j,k}\right]=\\[1.0em]
\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k}\Delta y \Delta z -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}\Delta y\Delta z}{\Delta x\Delta y\Delta z}=\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k} -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}}{\Delta x}
\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}\Delta A_{i-1/2,j,k}\right]=\\[1.0em]
\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k}\Delta y \Delta z -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}\Delta y\Delta z}{\Delta x\Delta y\Delta z}=\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k} -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}}{\Delta x}

Here, we multiply the inviscid flux vector with the normal direction, and we can see the direction of the normal vector in the figure above. It is always pointed out of the control volume, wo on some faces it will have a positive sign, on others, it will have a negative sign (when it points against the coordinate direction). On the face at location i+1/2,j,k, the normal vector points in the positive direction (along the x-axis), while at location i-1/2,j,k it points against the x coordinate direction, so its component will be negative.

We also see that the face area can be approximated by the length of the faces multiplied by each other. This can then be simplified with the volume to arrive at the expression we see above. For the inviscid fluxes in the y direction, we obtain in a similar way the following expression:

\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_y\right)_i A_i\right]=\frac{1}{\Delta x\Delta y\Delta z}\left[\begin{bmatrix}0\\
\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_y\right)_i A_i\right]=\frac{1}{\Delta x\Delta y\Delta z}\left[\begin{bmatrix}0\\1\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_y\right)_{i,j+1/2,k}\Delta A_{i,j+1/2,k}+\begin{bmatrix}0\\-1\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_y\right)_{i,j-1/2,k}\Delta A_{i,j-1/2,k}\right]=\\[1.0em]
\frac{\left(\mathbf{F}^{Inv}_y\right)_{i,j+1/2,k}\Delta x \Delta z -\left(\mathbf{F}^{Inv}_y\right)_{i,j-1/2,k}\Delta x\Delta z}{\Delta x\Delta y\Delta z}=\frac{\left(\mathbf{F}^{Inv}_y\right)_{i,j+1/2,k} -\left(\mathbf{F}^{Inv}_y\right)_{i,j-1/2,k}}{\Delta y}
\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_y\right)_{i,j+1/2,k}\Delta A_{i,j+1/2,k}+\begin{bmatrix}0\\-1\
\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_y\right)_i A_i\right]=\frac{1}{\Delta x\Delta y\Delta z}\left[\begin{bmatrix}0\\1\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_y\right)_{i,j+1/2,k}\Delta A_{i,j+1/2,k}+\begin{bmatrix}0\\-1\\0\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_y\right)_{i,j-1/2,k}\Delta A_{i,j-1/2,k}\right]=\\[1.0em]
\frac{\left(\mathbf{F}^{Inv}_y\right)_{i,j+1/2,k}\Delta x \Delta z -\left(\mathbf{F}^{Inv}_y\right)_{i,j-1/2,k}\Delta x\Delta z}{\Delta x\Delta y\Delta z}=\frac{\left(\mathbf{F}^{Inv}_y\right)_{i,j+1/2,k} -\left(\mathbf{F}^{Inv}_y\right)_{i,j-1/2,k}}{\Delta y}
\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_y\right)_{i,j-1/2,k}\Delta A_{i,j-1/2,k}\right]=\\[1.0em]
\frac{\left(\mathbf{F}^{Inv}_y\right)_{i,j+1/2,k}\Delta x \Delta z -\left(\mathbf{F}^{Inv}_y\right)_{i,j-1/2,k}\Delta x\Delta z}{\Delta x\Delta y\Delta z}=\frac{\left(\mathbf{F}^{Inv}_y\right)_{i,j+1/2,k} -\left(\mathbf{F}^{Inv}_y\right)_{i,j-1/2,k}}{\Delta y}

And, in the z direction, we also obtain in a similar way the following expression:

\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_z\right)_i A_i\right]=\frac{1}{\Delta x\Delta y\Delta z}\left[\begin{bmatrix}0\
\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_z\right)_i A_i\right]=\frac{1}{\Delta x\Delta y\Delta z}\left[\begin{bmatrix}0\\0\\1\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_z\right)_{i,j,k+1/2}\Delta A_{i,j,k+1/2}+\begin{bmatrix}0\\0\\-1\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_z\right)_{i,j,k-1/2}\Delta A_{i,j,k-1/2}\right]=\\[1.0em]
\frac{\left(\mathbf{F}^{Inv}_z\right)_{i,j,k+1/2}\Delta x \Delta y -\left(\mathbf{F}^{Inv}_z\right)_{i,j,k-1/2}\Delta x\Delta y}{\Delta x\Delta y\Delta z}=\frac{\left(\mathbf{F}^{Inv}_z\right)_{i,j,k+1/2} -\left(\mathbf{F}^{Inv}_z\right)_{i,j,k-1/2}}{\Delta z}
\\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_z\right)_{i,j,k+1/2}\Delta A_{i,j,k+1/2}+\begin{bmatrix}0\
\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{F}^{Inv}_z\right)_i A_i\right]=\frac{1}{\Delta x\Delta y\Delta z}\left[\begin{bmatrix}0\\0\\1\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_z\right)_{i,j,k+1/2}\Delta A_{i,j,k+1/2}+\begin{bmatrix}0\\0\\-1\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_z\right)_{i,j,k-1/2}\Delta A_{i,j,k-1/2}\right]=\\[1.0em]
\frac{\left(\mathbf{F}^{Inv}_z\right)_{i,j,k+1/2}\Delta x \Delta y -\left(\mathbf{F}^{Inv}_z\right)_{i,j,k-1/2}\Delta x\Delta y}{\Delta x\Delta y\Delta z}=\frac{\left(\mathbf{F}^{Inv}_z\right)_{i,j,k+1/2} -\left(\mathbf{F}^{Inv}_z\right)_{i,j,k-1/2}}{\Delta z}
\\-1\end{bmatrix}\cdot\left(\mathbf{F}^{Inv}_z\right)_{i,j,k-1/2}\Delta A_{i,j,k-1/2}\right]=\\[1.0em]
\frac{\left(\mathbf{F}^{Inv}_z\right)_{i,j,k+1/2}\Delta x \Delta y -\left(\mathbf{F}^{Inv}_z\right)_{i,j,k-1/2}\Delta x\Delta y}{\Delta x\Delta y\Delta z}=\frac{\left(\mathbf{F}^{Inv}_z\right)_{i,j,k+1/2} -\left(\mathbf{F}^{Inv}_z\right)_{i,j,k-1/2}}{\Delta z}

We can repeat the exact same steps for the viscous fluxes, as these are also first-order derivatives, and we obtain:

\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{G}^{Visc}_x\right)_i A_i\right]=\frac{\left(\mathbf{G}^{Visc}_x\right)_{i+1/2,j,k} -\left(\mathbf{G}^{Visc}_x\right)_{i-1/2,j,k}}{\Delta x}\\[1em]

\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{G}^{Visc}_y\right)_i A_i\right]=\frac{\left(\mathbf{G}^{Visc}_y\right)_{i,j+1/2,k} -\left(\mathbf{G}^{Visc}_y\right)_{i,j-1/2,k}}{\Delta y}\\[1em]

\frac{1}{V}\left[\sum_{i=1}^{n_{faces}}n_f\cdot \left(\mathbf{G}^{Visc}_z\right)_i A_i\right]=\frac{\left(\mathbf{G}^{Visc}_z\right)_{i,j,k+1/2} -\left(\mathbf{G}^{Visc}_z\right)_{i,j,k-1/2}}{\Delta z}

Putting all of this together now, our generalised vector equation can be written for a structured, cartesian grid as:

\frac{\mathbf{U}^{n+1}_{i,j,k}-\mathbf{U}^n_{i,j,k}}{\Delta t}+
\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k} -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}}{\Delta x} + \frac{\left(\mathbf{F}^{Inv}_y\right)_{i,j+1/2,k} -\left(\mathbf{F}^{Inv}_y\right)_{i,j-1/2,k}}{\Delta y} + \frac{\left(\mathbf{F}^{Inv}_z\right)_{i,j,k+1/2} -\left(\mathbf{F}^{Inv}_z\right)_{i,j,k-1/2}}{\Delta z}=
\\[1em]
\frac{\left(\mathbf{G}^{Visc}_x\right)_{i+1/2,j,k} -\left(\mathbf{G}^{Visc}_x\right)_{i-1/2,j,k}}{\Delta x} +
\frac{\left(\mathbf{G}^{Visc}_y\right)_{i,j+1/2,k} -\left(\mathbf{G}^{Visc}_y\right)_{i,j-1/2,k}}{\Delta y} +
\frac{\left(\mathbf{G}^{Visc}_z\right)_{i,j,k+1/2} -\left(\mathbf{G}^{Visc}_z\right)_{i,j,k-1/2}}{\Delta z}

Now, remember that we have 5 variables within each vector of \mathbf{U}, \mathbf{F}^{Inv}, and \mathbf{G}^{Visc}. This means we have to solve this equation 5 times for each cell and for each time step. All we have to do is find approximations of these quantities at the face centroids, e.g. i+1/2,j,k, insert them into the equation and solve for \mathbf{U}^{n+1}.

We find the approximations for the vector quantities at the face locations through interpolation, and we looked at suitable schemes in my article on numerical schemes for the finite volume method. A good choice here would be the MUSCL or WENO scheme. Let’s assume that we are using an explicit time integration for simplicity. Then the above vector form can be written as:

\mathbf{U}^{n+1}_{i,j,k}=\mathbf{U}^n_{i,j,k}-\Delta t\left[
\frac{\left(\mathbf{F}^{Inv}_x\right)_{i+1/2,j,k} -\left(\mathbf{F}^{Inv}_x\right)_{i-1/2,j,k}}{\Delta x} + \frac{\left(\mathbf{F}^{Inv}_y\right)_{i,j+1/2,k} -\left(\mathbf{F}^{Inv}_y\right)_{i,j-1/2,k}}{\Delta y} + \frac{\left(\mathbf{F}^{Inv}_z\right)_{i,j,k+1/2} -\left(\mathbf{F}^{Inv}_z\right)_{i,j,k-1/2}}{\Delta z}\right]
\\[1em]
-\Delta t\left[\frac{\left(\mathbf{G}^{Visc}_x\right)_{i+1/2,j,k} -\left(\mathbf{G}^{Visc}_x\right)_{i-1/2,j,k}}{\Delta x} +
\frac{\left(\mathbf{G}^{Visc}_y\right)_{i,j+1/2,k} -\left(\mathbf{G}^{Visc}_y\right)_{i,j-1/2,k}}{\Delta y} +
\frac{\left(\mathbf{G}^{Visc}_z\right)_{i,j,k+1/2} -\left(\mathbf{G}^{Visc}_z\right)_{i,j,k-1/2}}{\Delta z}\right]

Add the correct boundary conditions, and you have all you need to write a compressible Navier-Stokes solver.

Incompressible flows

Let’s start our discussion with why we want to use different algorithms for incompressible flows. At the beginning of this article, we said that incompressible flows are those where the density does not change. We also saw that there are materials, such as water, that behave like an incompressible fluid due to their low compressibility factor \beta. This means that if we were to use any equation of state, for example, the ideal gas law, we have a functional relation between the pressure p and the density \rho as:

p=\rho RT

If the density changes slowly, then so will the pressure. OK, moving on, let’s consider the CFL condition next. When we looked at the CFL number, we said that this is a measure to determine a stable time step for explicit time integration methods. For our discussion, let us write the CFL number in the following way:

CFL=u_{max}\frac{\Delta t}{\Delta x},\quad u_{max}=\max(\lambda),\quad \lambda=[u, u+a, u-a]

Here, \lambda represents the eigenvalues of the equation we are solving. The eigenvalues will give us the characteristic speeds at which information is travelling. I have only shown the eigenvalues here for a one dimensional flow, in two and three dimensions, we also get v, v+a, v-a and w, w+a,w-a, respectively. Here, a is the speed of sound. Let’s quickly review what these eigenvalues physically represent.

The first eigenvalue u represents the speed at which a contact wave is travelling. A contact wave is defined as an interface at which the density is discontinuous, but the pressure is not. For example, if you consider a flame and you look at an interface between burnt and non-burned gas, there will be a difference in the density across the flame front since the flame will have a higher temperature, and so the burned gas will have a different density to the unburned gas, where the temperature is lower. The pressure is the same across the interface, though.

Interfaces, in general, are good examples of contact discontinuities, and we can see where the name contact is derived from. The water-air interface on a body of water is another good example. At the interface, air will impose atmospheric pressure on the water, so they will have the same pressure on either side of the interface, but the density of air and water is very different. In general, interfaces in multiphase flows will have these contact discontinuities, and the interfaces move with the characteristic speed u (and v,\, w in a three-dimensional flow).

The two other eigenvalues represent acoustic pressure waves travelling through your domain. Consider the following jet travelling at speed u:

We see that the wavefront travels with a velocity of u+a and u-a. If we are on the ground and u=0, then these wavefronts would propagate in all directions equally (with speed a, i.e. the speed of sound). But once we start moving, the sound wave ahead of the jet will slow down (u-a) while the propagation of waves behind the jet will increase (u+a).

Once we reach u=a, or a Mach number of Ma=u/a=1, the wave fronts can no longer propagate in front of us. We will have overtaken them. As a result, if we place a tangent on all of the wave fronts that emanate from the tip of the jet, we get the shock cone shown above, which represents the shock waves that will form.

Returning to the discussion of our incompressible flow, we said that an incompressible flow typically satisfies Ma<0. In many cases, we have speeds much smaller than Ma=0.3. Think, for example, of flow within a boundary layer. Close to the wall, we get very low velocities. Thus, our local velocity is very small in comparison to the speed of sound a.

Returning to our CFL condition, we said that it is given as CFL=u_{max}\Delta t/\Delta x. Solving this for the time step gives us \Delta t= CFL \Delta x /u_{max}. If we consider the flow around an airfoil of chord length 1 meter, a Reynolds number of around a million, standard air and an explicit time stepping, then we have

• CFL\le 1
• \Delta x<<1 (if we resolve this down to a y^+ value of 1, we have \Delta x\approx10^{-6} meters)
• u_{max}=\max (u, u+a, u-a)\approx 10^2\, m/s. For standard air, we have a=340\, m/s.

If we just consider orders of magnitude, then we can compute the timestep as

\Delta t=\frac{1\cdot 10^{-6}}{10^2}=10^{-8}\, s

This is a rather small time step. If we are interested in studying the vortex shedding behind the airfoil (well, we need a three-dimensional wing rather than a two-dimensional airfoil), then we can look at the relation between the Reynolds number and the Strouhal number, which is given by the following image:

The Strouhal number is given as

St=\frac{fL}{u}

Here, f is the frequency at which periodic flow patterns occur (such as vortex shedding), L is the characteristic length (which is the same definition of L in the Reynolds number), and u is the characteristic velocity (again, the same as in the Reynolds number definition.

Solving this equation for f, we get:

f=\frac{St\cdot u}{L}

In the example for the airfoil, we assumed standard air and a Reynolds number of around 1 million. With a characteristic length of L=1 meter, we have a velocity of around u=15 meters per second. From the plot above, we can see that a Reynolds number of a million produces a Strouhal number of about St=0.2. This means that our frequency can be computed as:

f=\frac{0.2\cdot 15}{1}=3\,\frac{1}{s}=3\, Hz

This means that one vortex shedding frequency takes T=1/f=1/3 seconds. If we divide this by the time step \Delta t=10^{-8} seconds, then we get the total number of time steps required to resolve a single vortex shedding period. We have: (1/3)/10^{-8}=33,333,333 time steps. This sounds like a lot, but it may not even be the upper limit. We said before that the pressure changes with changes in density. But, if the density changes very slowly, the pressure will change slowly, too.

Thus, we may even need more than the 33 million time steps for a single vortex shedding period. We can see that using a compressible method requires a lot of time steps to get the solution for an essentially incompressible fluid. In the literature, this is typically referred to as a stiff system of equations. The stiffness here means that the system does not change quickly and so we need a lot of iterations to see changes in the solution.

What are the solutions, then? Well, we no longer want to depend on density changes to drive the solution. Instead, we want to compute changes in pressure directly and see how that affects the solution. For this reason, we also differentiate between density-based and pressure-based solution approaches to solve compressible and incompressible flows, respectively. For example, in ANSYS Fluent, you have to specify if you want to use a density-based or a pressure-based approach, which selects either the compressible or incompressible solver.

If we already have a working compressible solver, we might be tempted to use it for incompressible flows if we can introduce a small correction instead of having to write a new solver from scratch. Approaches have been introduced in the literature which fall under the topic of preconditioning or low Mach number correction approaches. These are applied top density-based formulations, i.e. those where we link pressure and density through the equation of state.

Let’s look at preconditioning. One of the issues where a stiff system of equations comes from is the discrepancy between the local flow velocity u and the speed of sound a. If we can scale our system of equations so that the propagation speed of acoustic waves becomes similar to the local flow velocities, then we get larger time steps, and as a result, fewer time steps are required to advance the solution in time. We will look at one such method in just a second.

On the other hand, pressure-based solution algorithms accept that there is no coupling between the pressure and the density, so the equation of state can no longer be used. Instead, we say the pressure and velocity are completely decoupled and we have to come up with a clever way of coupling them back together. These algorithms are known as pressure velocity coupling algorithms, and you may be more familiar with their street names like SIMPLE, PISO, PIMPLE (great name!), and, in general, projection, or pressure projection, methods.

The last sentence has 8 commas, a new record. Whoop whoop 🥳🎉. Let’s continue.

We will review a subset of each of these methods so that you get an idea of how each method works. We will start with my personal favourite and the highly underrated artificial compressibility method first. It is my favourite because it is the simplest incompressible method to implement. It allows for a fully explicit implementation, making it a great method to implement from scratch, while all other project-based methods require the solution of an implicit system of equations. This can be confusing if you just want to start out.

However, while the artificial compressibility method is easy to implement, it is not as efficient as its projection-based counterparts (like SIMPLE, PISO, etc.). So you won’t find it used much these days. But for your own personal fun and solver development, I’ll give it an A+. Once you have that implemented and working, you can always switch to a different method and validate it against your artificial compressibility method implementation. Sounds good? Let’s see what this method is all about!

Artificial Compressibility

The artificial compressibility method was first published in 1967 by one of the great minds in CFD. Chorin, to me, is the father of incompressible algorithms. He gave us not just the artificial compressibility method but also introduced us to projection methods (covered later), which are the foundation for methods you have probably already used (SIMPLE, PISO, etc.). He is right up there with the greats such as Lax, Wendroff, Godunov, Harten, and van Leer.

So let’s see what Chorin did. He said that we should just use the equations that govern compressible flows. Now, we already know that this won’t work for the continuity equation, as the density is constant. However, let’s review this for a second. The equation is given as:

\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{u})=0

Chorin said that we know that the pressure and density are related through the equation of state, i.e. we have p=\rho R T. While we know that this equation no longer holds for an incompressible flow, we can generalise this as a functional relationship. We can express that as p=f(\rho, \beta)=\rho\beta, where \beta is some unknown scaling factor. This may not be a constant.

Chorin now said that if we solve this for the density, i.e. \rho=p/\beta, then we can insert this into the continuity equation. This would result in

\frac{1}{\beta}\frac{\partial p}{\partial \tau}+\nabla\cdot(\rho\mathbf{u})=0

Here, we have written the value of 1/\beta in front of the time derivative. Why? Well, we don’t know its value, so we treat it as a constant, even though it may not be one. But, take a closer look at the time derivative. It is no longer a time derivative. We have exchanged t (real time) for a new variable \tau (pseudo time). What does that mean?

Well, the equation above no longer represents the continuity equation. But, if we weaken the equation a bit, we can insist that the above equation is still valid. We can say that the equation is only valid in the limit of steady-state flows. Steady-state flows require, by definition, that there are no changes in time. If there are no more changes in time, then we have \partial p/\partial\tau = 0. So the time derivative of the pressure will vanish, and we recover the continuity equation of an incompressible flow, i.e. \nabla\cdot (\rho\mathbf{\rho}).

This doesn’t mean that the time history itself is completely wrong; if you look at an animation of the artificial compressibility method, it will appear as if it resolves the transient behaviour correctly. However, if you were to take measurements that depend on the timestep, like measuring frequencies, these would be all wrong now. Thus, we only resolve a physically meaningful flow field once the time derivative vanishes and the flow has settled for a steady-state solution.

Let us write the three-dimensional scalar form of this equation in Cartesian coordinates. Then we get

\frac{1}{\beta}\frac{\partial p}{\partial \tau}+\frac{\partial \rho u}{\partial x}+\frac{\partial \rho v}{\partial y}+\frac{\partial \rho w}{\partial z}=0

If we perform a quick discretisation using the finite difference method with a first-order Euler in time and a second-order central scheme in space, we get:

\frac{1}{\beta}\frac{p^{n+1}_{i,j,k}-p^n_{i,j,k}}{\Delta \tau}+\rho\left(\frac{u^n_{i+1,j,k}-u^n_{i-1,j,k}}{2\Delta x}+\frac{v^n_{i,j+1,k}-v^n_{i,j-1,k}}{2\Delta y}+\frac{w^n_{i,j,k+1}-w^n_{i,j,k-1}}{2\Delta z}\right)=0

Solving this for p^{n+1}_{i,j,k}, we get

p^{n+1}_{i,j,k}=p^n_{i,j,k} - \Delta\tau \rho\beta\left(\frac{u^n_{i+1,j,k}-u^n_{i-1,j,k}}{2\Delta x}+\frac{v^n_{i,j+1,k}-v^n_{i,j-1,k}}{2\Delta y}+\frac{w^n_{i,j,k+1}-w^n_{i,j,k-1}}{2\Delta z}\right)

Thus, we have found a cheap and easy way to compute an updated pressure field using the continuity equation. We can use this pressure field now in our momentum equation, which was given as:

\frac{\partial \mathbf{u}}{\partial \tau} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}

Let’s write the scalar form of this equation as well. In three dimensions, we get three equations for \mathbf{u}=[u,v,w], which are:

\frac{\partial u}{\partial \tau}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}=-\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+\frac{\partial^2 u}{\partial z^2}\right)\\[1em]

\frac{\partial v}{\partial \tau}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}=-\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}+\frac{\partial^2 v}{\partial z^2}\right)\\[1em]

\frac{\partial w}{\partial \tau}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}=-\frac{1}{\rho}\frac{\partial p}{\partial z}+\nu\left(\frac{\partial^2 w}{\partial x^2}+\frac{\partial^2 w}{\partial y^2}+\frac{\partial^2 w}{\partial z^2}\right)

To make the discretisation a bit more easy to manage, let us introduce the following operator:

\mathrm{upwind}(\phi)=\max(u_{i,j,k},\,0)\frac{\phi_{i,j,k}-\phi_{i-1,j,k}}{\Delta x}+\min(u_{i,j,k},\, 0)\frac{\phi_{i+1,j,k}-\phi_{i, j, k}}{\Delta x}+\\[1em]

\max(v_{i,j,k},\,0)\frac{\phi_{i,j,k}-\phi_{i,j-1,k}}{\Delta y}+\min(v_{i,j,k},\, 0)\frac{\phi_{i,j+1,k}-\phi_{i, j, k}}{\Delta y}+\\[1em]

\max(w_{i,j,k},\,0)\frac{\phi_{i,j,k}-\phi_{i,j,k-1}}{\Delta z}+\min(w_{i,j,k},\, 0)\frac{\phi_{i,j,k+1}-\phi_{i, j, k}}{\Delta z}

This is the upwind discretisation. Here, \phi is one of the velocity components, i.e. [u,v,w]. We can also introduce a diffusion operator as:

\mathrm{diffusion}(\phi)=\nu\left(\frac{\phi_{i+1,j,k}-2\phi_{i,j,k}+\phi_{i-1,j,k}}{\Delta x^2}+\frac{\phi_{i,j+1,k}-2\phi_{i,j,k}+\phi_{i,j-1,k}}{\Delta y^2}+\frac{\phi_{i,j,k+1}-2\phi_{i,j,k}+\phi_{i,j,k-1}}{\Delta z^2}\right)

With these operators defined, let us write the discretised form of the three-dimensional momentum equation with first-order Euler in time and a central operator for the pressure. This results in:

\frac{u^{n+1}_{i,j,k}-u^n_{i,j,k}}{\Delta \tau}+\mathrm{upwind}(u)=-\frac{1}{\rho}\frac{p^{n+1}_{i+1,j,k}-p^{n+1}_{i-1,j,k}}{2\Delta x}+\mathrm{diffusion}(u)

We can solve this now for u^{n+1}_{i,j,k} to arrive at:

u^{n+1}_{i,j,k}=u^n_{i,j,k}+\Delta \tau\left[-\mathrm{upwind}(u)-\frac{1}{\rho}\frac{p^{n+1}_{i+1,j,k}-p^{n+1}_{i-1,j,k}}{2\Delta x}+\mathrm{diffusion}(u)\right]

We can find similar discretisations for v^{n+1}_{i,j,k} and w^{n+1}_{i,j,k} as:

v^{n+1}_{i,j,k}=v^n_{i,j,k}+\Delta \tau\left[-\mathrm{upwind}(v)-\frac{1}{\rho}\frac{p^{n+1}_{i,j+1,k}-p^{n+1}_{i,j-1,k}}{2\Delta y}+\mathrm{diffusion}(v)\right]

And:

w^{n+1}_{i,j,k}=w^n_{i,j,k}+\Delta \tau\left[-\mathrm{upwind}(w)-\frac{1}{\rho}\frac{p^{n+1}_{i,j,k+1}-p^{n+1}_{i,j,k-1}}{2\Delta z}+\mathrm{diffusion}(w)\right]

Since we obtain the pressure already, it is no longer an unknown and so we can solve the momentum equation. We do that until all pseudo-time derivatives vanish. How can we check that? Well, if the pseudo time derivatives are discretised as \partial \phi/\partial\tau=(\phi^{n+1}-\phi^n)/\Delta \tau, then \partial \phi/\partial is zero for \phi^{n+1}=\phi^n, since (\phi^{n+1}-\phi^n)/\Delta \tau=0/\Delta\tau=0.

In other words, if we evaluate the continuity and momentum equation until they no longer change, then we have found a valid solution that satisfies the Navier-Stokes equation again. We can introduce a convergence condition that states that if the solution does not change by more than 0.1%, for example, then we have reached convergence. This may be expressed as \phi^{n+1}-\phi^n<10^{-3}.

Another convergence condition we can impose is the incompressible continuity equation itself. We know that for an incompressible flow, we have \nabla\cdot\mathbf{u}=0. We can impose a similar constraint by requiring the continuity equation to be satisfied, for example, by stating \nabla\cdot\mathbf{u}<10^{-3}. This could be written as

\frac{u_{i+1,j,k}-u_{i-1,j,k}}{2\Delta x}+\frac{v_{i,j+1,k}-v_{i,j-1,k}}{2\Delta y}+\frac{w_{i,j,k+1}-w_{i,j,k-1}}{2\Delta z}<10^{-3}

We can now check that this condition is satisfied at each node i,j,k, or we can compute the average through the L2 vector norm, for example, and compare that against 10^{-3} (or a smaller convergence condition).

Hopefully everything up to this point makes sense. The artificial compressibility method does not get more complicated than this. It is very straightforward to implement, and it does give results reasonably well. It suffers from a slower convergence compared to other methods, but there are ways to stabilise the procedure to potentially bring it closer to the projection methods discussed next. But before we go there, I wanted to close the circle and bring back our discussion on preconditioning.

For this, let us write all scalar equations we have obtained thus far below:

\frac{1}{\beta}\frac{\partial p}{\partial \tau}+\frac{\partial \rho u}{\partial x}+\frac{\partial \rho v}{\partial y}+\frac{\partial \rho w}{\partial z}=0\\[1em]

\frac{\partial u}{\partial \tau}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z} + \frac{1}{\rho}\frac{\partial p}{\partial x}=\nu\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+\frac{\partial^2 u}{\partial z^2}\right)\\[1em]

\frac{\partial v}{\partial \tau}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z} + \frac{1}{\rho}\frac{\partial p}{\partial y}=\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}+\frac{\partial^2 v}{\partial z^2}\right)\\[1em]

\frac{\partial w}{\partial \tau}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z} + \frac{1}{\rho}\frac{\partial p}{\partial z}=\nu\left(\frac{\partial^2 w}{\partial x^2}+\frac{\partial^2 w}{\partial y^2}+\frac{\partial^2 w}{\partial z^2}\right)

We can also write this system of equations as:

\mathbf{P}^{-1}\frac{\partial \mathbf{U}}{\partial \tau}+\mathbf{A}^{Inv}\frac{\partial \mathbf{U}}{\partial x}+\mathbf{B}^{Inv}\frac{\partial \mathbf{U}}{\partial y}+\mathbf{C}^{Inv}\frac{\partial \mathbf{U}}{\partial z}=\mathbf{A}^{Visc}\frac{\partial^2 \mathbf{U}}{\partial x^2}+\mathbf{B}^{Visc}\frac{\partial^2 \mathbf{U}}{\partial y^2}+\mathbf{C}^{Visc}\frac{\partial^2 \mathbf{U}}{\partial z^2}

Here, we have introduced the primitive variable vector \mathbf{U}=[p,u,v,w]^T and matrices:

\mathbf{A}^{Inv}=\begin{bmatrix}0 & \rho & 0 & 0 \\ 1/\rho & u & 0 & 0 \\ 0 & 0 & u & 0 \\ 0 & 0 & 0 & u\end{bmatrix},\quad
\mathbf{B}^{Inv}=\begin{bmatrix}0 & 0 & \rho & 0 \\ 0 & v & 0 & 0 \\ 1/\rho & 0 & v & 0 \\ 0 & 0 & 0 & v\end{bmatrix},\quad
\mathbf{C}^{Inv}=\begin{bmatrix}0 & 0 & 0 & \rho \\ 0 & w & 0 & 0 \\ 0 & 0 & w & 0 \\ 1/\rho & 0 & 0 & w\end{bmatrix},\quad\\[1em]
\mathbf{A}^{Visc}=\mathbf{B}^{Visc}=\mathbf{C}^{Visc}=\begin{bmatrix}0 & 0 & 0 & 0 \\ 0 & \nu & 0 & 0 \\ 0 & 0 & \nu & 0 \\ 0 & 0 & 0 & \nu\end{bmatrix}

Does this vector equation look familiar? Yes, it is almost identical to the vector equation we saw in the discussion on compressible flows. The difference here is that we use primitive or non-conserved variables. So we no longer deal with the inviscid and viscous fluxes \mathbf{F}^{Inv} and \mathbf{G}^{Visc} but rather with the primitive variables \mathbf{U} directly. To see that these matrixes give us the right scalar equation again, we can write out the vector equation as:

\mathbf{P}^{-1}\begin{bmatrix}\partial p/\partial \tau \\ \partial u/\partial \tau \\ \partial v/\partial \tau \\ \partial w/\partial \tau\end{bmatrix}+
\begin{bmatrix}0 & \rho & 0 & 0 \\ 1/\rho & u & 0 & 0 \\ 0 & 0 & u & 0 \\ 0 & 0 & 0 & u\end{bmatrix}
\begin{bmatrix}\partial p/\partial x \\ \partial u/\partial x \\ \partial v/\partial x \\ \partial w/\partial x\end{bmatrix}+
\begin{bmatrix}0 & 0 & \rho & 0 \\ 0 & v & 0 & 0 \\ 1/\rho & 0 & v & 0 \\ 0 & 0 & 0 & v\end{bmatrix}
\begin{bmatrix}\partial p/\partial y \\ \partial u/\partial y \\ \partial v/\partial y \\ \partial w/\partial y\end{bmatrix}+
\begin{bmatrix}0 & 0 & 0 & \rho \\ 0 & w & 0 & 0 \\ 0 & 0 & w & 0 \\ 1/\rho & 0 & 0 & w\end{bmatrix}
\begin{bmatrix}\partial p/\partial z \\ \partial u/\partial z \\ \partial v/\partial z \\ \partial w/\partial z\end{bmatrix}
=\\[1em]
\begin{bmatrix}0 & 0 & 0 & 0 \\ 0 & \nu & 0 & 0 \\ 0 & 0 & \nu & 0 \\ 0 & 0 & 0 & \nu\end{bmatrix}
\begin{bmatrix}\partial^2 p/\partial x^2 \\ \partial^2 u/\partial x^2 \\ \partial^2 v/\partial x^2 \\ \partial^2 w/\partial x^2\end{bmatrix}+
\begin{bmatrix}0 & 0 & 0 & 0 \\ 0 & \nu & 0 & 0 \\ 0 & 0 & \nu & 0 \\ 0 & 0 & 0 & \nu\end{bmatrix}
\begin{bmatrix}\partial^2 p/\partial y^2 \\ \partial^2 u/\partial y^2 \\ \partial^2 v/\partial y^2 \\ \partial^2 w/\partial y^2\end{bmatrix}+
\begin{bmatrix}0 & 0 & 0 & 0 \\ 0 & \nu & 0 & 0 \\ 0 & 0 & \nu & 0 \\ 0 & 0 & 0 & \nu\end{bmatrix}
\begin{bmatrix}\partial^2 p/\partial z^2 \\ \partial^2 u/\partial z^2 \\ \partial^2 v/\partial z^2 \\ \partial^2 w/\partial z^2\end{bmatrix}

The advantage of doing this is that we can create an incompressible solver with very little effort if we already have a compressible solver implemented. We can reuse the same numerical schemes and everything else stays the same, we just have to change the the vector equation from conserved to non-conserved variable formulation and ensure that we have subsonic boundary conditions implemented.

The downside is that the artificial compressibility method retains the hyperbolic character of the compressible Navier-Stokes equations. As we discussed in the article on hyperbolic, parabolic, and elliptic partial differential equations, hyperbolic flows behave very differently from elliptic flows.

Is that a problem? Not necessarily, but consider the following analogy: Imagine you are extremely active, and you exercise several times per day. You do some running and some cycling, and you visit your local gym in the morning before heading to university or work. Now, your friend approaches you and asks you: “Do you want to run a marathon with me?” You probably don’t even have to think twice and accept the invite on the spot because this activity aligns with your character.

Now, imagine you are a couch potato. While you do occasionally enjoy outdoor walking and perhaps a bit of recreational cycling, you feel most comfortable in your home sweet home. You work or study from home when you can, and you are happy to do your groceries once or twice a week; you don’t need more activity in your life than that. Now your friend comes around again and asks: “Do you want to run a marathon with me?”. This does not align with your character, so your friend would have to spend a lot of time convincing you to join in.

In both cases, you can be convinced to run the marathon, but in once case, it is very fast, in the other case, it takes time. The same is true for the Navier-Stokes equation. If the equations you are solving are fully hyperbolic (artificial compressibility method) but the underlying equations they are solving are mixed elliptic/hyperbolic (incompressible Navier-Stokes equation), then it will take time to get convergence.

However, if you solve the compressible Navier-Stokes equations, which are fully hyperbolic, with a system of fully hyperbolic equations, then you get better and much faster convergence. Since we do not modify the equations for compressible flows but rather solve them as they are (while closing them with the equation of state), we retain the mathematical character and, therefore, get better convergence.

We haven’t looked yet at \mathbf{P}. This can be written as:

\mathbf{P}^{-1}=\begin{bmatrix}1/\beta & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0\end{bmatrix}

We see that this matrix only multiplies the time derivative in the continuity equation with the factor of 1/\beta. We said earlier that preconditioning methods achieve faster convergence for compressible methods by scaling the equations in some form. Here, we can think of \beta as a scaling parameter for the pressure. Its value can be anything we want since it multiplies the time derivative of the pressure, which will vanish for steady-state flows.

Think of the physical role of pressure; it dictates the propagation of pressure waves, which in turn is governed by the speed of sound. We said that incompressible flows suffer from slow convergence because local flow velocities are much smaller than the local speed of sound. So, by tampering (scaling) the pressure, we could reasonably argue that we would be tampering with the local speed of sound. If we wanted to see this in equation form, we would need to calculate the eigenvalues of the system. To do that, let’s write out our continuity equation for a moment:

\frac{1}{\beta}\frac{\partial p}{\partial \tau}+\frac{\partial \rho u}{\partial x}+\frac{\partial \rho v}{\partial y}+\frac{\partial \rho w}{\partial z}=0

Let’s modify it by multiplying this equation by \beta. This will result in:

\frac{\partial p}{\partial \tau}+\frac{\partial \beta\rho u}{\partial x}+\frac{\partial \beta\rho v}{\partial y}+\frac{\partial \beta\rho w}{\partial z}=0

Since \beta is a constant (or at least, we treat it as one), we can just place it within the derivative. Now, the preconditioning matrix \mathbf{P}^{-1} is a zero matrix, so we can remove it, but we need to modify our \mathbf{A}^{Inv}, \mathbf{B}^{Inv}, and \mathbf{C}^{Inv} matrix, which now become:

\mathbf{A}^{Inv}=\begin{bmatrix}0 & \beta\rho & 0 & 0 \\ 1/\rho & u & 0 & 0 \\ 0 & 0 & u & 0 \\ 0 & 0 & 0 & u\end{bmatrix},\quad
\mathbf{B}^{Inv}=\begin{bmatrix}0 & 0 & \beta\rho & 0 \\ 0 & v & 0 & 0 \\ 1/\rho & 0 & v & 0 \\ 0 & 0 & 0 & v\end{bmatrix},\quad
\mathbf{C}^{Inv}=\begin{bmatrix}0 & 0 & 0 & \beta\rho \\ 0 & w & 0 & 0 \\ 0 & 0 & w & 0 \\ 1/\rho & 0 & 0 & w\end{bmatrix}

Now we have to look at the eigenvalues of \mathbf{A}^{Inv}, \mathbf{B}^{Inv}, and \mathbf{C}^{Inv}. We compute the eigenvalues for each matrix as

\lambda^\mathbf{A}=\det(\mathbf{A}^{Inv}-\lambda\mathbf{I}),\quad \lambda^\mathbf{B}=\det(\mathbf{B}^{Inv}-\lambda\mathbf{I}),\quad \lambda^\mathbf{C}=\det(\mathbf{C}^{Inv}-\lambda\mathbf{I})

We can do this by hand, but for anything more than a 2 by 2 matrix, I can’t be bothered doing this analytically and rather write a piece of code that can do this for me. I trust that most know how to do it by hand, but perhaps not how to do this in code, so let’s see how we can achieve this. The following Python code uses the sympy package for symbolic math manipulation. Let’s look through it first and then discuss it afterwards.

import sympy as sp

# define symbols required to construct matrices
rho = sp.symbols('rho')
beta = sp.symbols('beta')
u = sp.symbols('u')
v = sp.symbols('v')
w = sp.symbols('w')

# define inviscid coefficient matrices
A = sp.Matrix([
    [0, beta * rho, 0, 0],
    [1/rho, u, 0, 0],
    [0, 0, u, 0],
    [0, 0, 0, u]
])

B = sp.Matrix([
    [0, 0, beta * rho, 0],
    [0, v, 0, 0],
    [1/rho, 0, v, 0],
    [0, 0, 0, v]
])

C = sp.Matrix([
    [0, 0, 0, beta * rho],
    [0, w, 0, 0],
    [0, 0, w, 0],
    [1/rho, 0, 0, w]
])

# compute eigenvalues and roots of eigenvalues
lam = sp.symbols('lambda')

eigenvalues_of_A = sp.det(A - lam * sp.eye(4))
eigenvalues_of_B = sp.det(B - lam * sp.eye(4))
eigenvalues_of_C = sp.det(C - lam * sp.eye(4))

roots_of_A = sp.roots(sp.Poly(eigenvalues_of_A, lam))
roots_of_B = sp.roots(sp.Poly(eigenvalues_of_B, lam))
roots_of_C = sp.roots(sp.Poly(eigenvalues_of_C, lam))

# print eigenvalues
print(roots_of_A)
print(roots_of_B)
print(roots_of_C)

After we have imported the module on line 1, we introduce our symbols on lines 4-8. Sympy needs these to know which variables to treat as symbols, in our case, we have the density and velocity components, as well as the \betaparameter. We then construct the coefficient matrices on lines 11-30 which we can now use for the eigenvalue calculation.

We define the additional symbol lambda on line 33, which we then use to compute the eigenvalues on lines 35-37. Here, sp.eye(4) produces a 4 by 4 identity matrix. With those eigenvalues found, we want to extract the roots of the polynomial that is generated by lines 35-37, which is done on lines 39-41. Printing the results for each matrix on lines 44-46 results in the following values:

A:\quad u,\,u,\,\frac{1}{2}\left(u+\sqrt{u^2+4\beta}\right),\,\frac{1}{2}\left(u-\sqrt{u^2+4\beta}\right)\\[1em]
B:\quad v,\,v,\,\frac{1}{2}\left(v+\sqrt{v^2+4\beta}\right),\,\frac{1}{2}\left(v-\sqrt{v^2+4\beta}\right)\\[1em]
C:\quad w,\,w,\,\frac{1}{2}\left(w+\sqrt{w^2+4\beta}\right),\,\frac{1}{2}\left(w-\sqrt{w^2+4\beta}\right)

We have again a few contact discontinuities in each direction, i.e. u, v, and w. And then, we have two additional eigenvalues. Let’s see what happens if we set \beta=0. In this case, the square root of, for example, \sqrt{u^2+4\beta} becomes \sqrt{u^2+4\cdot 0}=\sqrt{u^2}=u.

Let’s abbreviate this square root with the letter a. Then, we can rewrite these eigenvalues as 0.5(u+\sqrt{u^2+4\beta})=0.5(u+a) and 0.5(u-\sqrt{u^2+4\beta})=0.5(u-a). So we see, the content of the square root resembles the speed of sound in compressible flows, or, in other words, the propagation speed of characteristics (e.g. pressure waves, disturbances).

Now, we are in a better position to judge the real importance of the \beta parameter. Not only is it used now to represent the functional relationship between the pressure p and the density \rho, but also to scale the characteristic speed at which disturbances (pressure waves) travel. In other words, we are able to scale a=\sqrt{u^2+4\beta} by modifying \beta so that u and a become closer together.

Let’s return to our airfoil example above, where we calculated the number of timesteps required to resolve one vortex shedding cycle or period. We said that our time step is calculated as:

\Delta t=\frac{1\cdot 10^{-6}}{10^2}=10^{-8}\, s

Here, the denominator of 10^2 represents the max eigenvalue, e.g. u\pm a. If we are able to scale this value now by \beta, so that a becomes either the same order of magnitude or insignificant compared to u, then the denominator simply becomes u. Indeed, for incompressible flows, the CFL number is calculated based on the maximum local velocity rather than the maximum eigenvalue.

Thus, let’s say that our velocities are of the order of 1, then our time step becomes:

\Delta t=\frac{1\cdot 10^{-6}}{1}=10^{-6}\, s

Before we had 33,333,333 time steps required to resolve a single vortex shedding period, this reduces no by two orders of magnitude, i.e. to (1/f)/(\Delta t)=(1/3)/10^{-6}=333,333 time steps required per vortex shedding period. Even though this change seems small, with a good high-performance cluster, even a large and complex case, this small change can be the difference between 1 day of computation vs 100 days of computation.

Thus, by optimising the \beta value, we can (in theory) reduce the number of iterations required to reach convergence. But as mentioned before, a value of \beta=1 is typically best.

If you want to play around with this parameter, you can try the range of 0.1\le\beta\le 10, these tend to work well. A while ago, I did test it for a wide range of values. This is what I obtained:

On the left y-axis, you see the number of iterations required, and on the right y-axis, the time taken (in seconds) to get a solution. On the x-axis, we see different values for the \beta parameter. I have tested this for one test case and different numerical schemes (where RS stands for Riemann solver). For other test cases, this will likely look slightly different, but you will probably always get an optimum around \beta=1.

And this is all there is to the artificial compressibility method. It is a quick and easy method to implement for incompressible flows, and I would strongly recommend starting with this formulation first when you are writing an incompressible Navier-Stokes equation solver. It is very quick to get an implementation with this equation, and we don’t have to mess with linear systems here, which is a bonus. This cannot be said for the remaining equations we will look at, so let’s have a look at them next.

Exact projection methods: The fractional step, pressure projection

Chorin gave us the artificial compressibility method (ACM), and if you ever walk the Monaco Formula 1 track, you will lots of appreciation for the ACM. They have placed road signs everywhere and even created a medal of honour that any citizen receives who understands the ACM. Don’t believe me? Here is proof (go on, you have my permission to have a virtual walk around Monaco on google earth, I’ll wait here …)

Look, it’s either that or ACM stands for Auto Club de Monaco; you choose which conspiracy theory you believe most …

But Chorin didn’t stop there. He thought, why not introduce a second, completely different method? Why not create some competition for myself? Why not fight fire with fire? His theory was so good that it was not fully understood (in my personal view). In fact, many copy cats sprung up (long before we even used the term copy cats) and feasted on his research without either knowing it (i.e. incompetency) or acknowledging it (i.e. scientific misconduct (fancy word for fraud)). We’ll get to that, don’t you worry. The SIMPLE method will get some beating today …

But before we collectively circle the SIMPLE algorithm and stone it as it should have been a long, long time ago, let’s first indulge ourselves in the pressure projection method of Chorin. It is pretty ingenious. If it is not obvious by now, I do have a little man crush on Chorin, but don’t tell him.

The main problem with the momentum equation for incompressible flows is that it contains the pressure gradient. If the equation only contained the velocity, well, then there wouldn’t be any dependence on the pressure and we could solve the momentum equation directly. So Chorin said, “bin it” (I did not check but I’m sure these were his words!) and threw the pressure out. Problem solved. It must be this famous Polish mathematicians’ humour. Well, this is the momentum equation that we are getting:

\frac{\mathbf{u}^*-\mathbf{u}^n}{\Delta t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = \nu\nabla^2\mathbf{u}

The pressure gradient is gone, but wait, he did it again! Look at the time derivative, it is all messed up. Chorin appears to have a troublesome relation with time derivatives, but it is not as bad as it first seems.

In the artificial compressibility method, we scaled the time derivative in the continuity equation. This directly destroyed the time derivative and we no longer had real time but rather pseudo time. But in the equation above, we left the time derivative as it is, so we are still operating in real time, not pseudo time. However, since we removed the pressure, we can no longer claim that the velocity we obtain from this equation represents the velocity at the next time step n+1.

Instead, we solve the velocity for some intermediate velocity field \mathbf{u}^*, which is somewhere between time level n and n+1. For this reason, we have to write the time derivative in discretised form, and we call the equation above a semi-discretised equation as a result (you may see this term used in the literature).

We can’t ignore the pressure indefinitely, though, and eventually, we have to face it. What Chorin did was to use the so called fractional step procedure. This procedure allows you to split partial differential equation (or ordinary ones) into several equations and solve them individually. The way that it works is to write each dropped term (here the pressure gradient) as its own fractional step. Each fractional step will then advance the time derivative from the intermediate velocity field to the next time level.

We can introduce as many fractional steps as we want, we would just have more intermediate velocity fields. In our case, we have exactly two fractional step, and the second becomes:

\frac{\mathbf{u}^{n+1}-\mathbf{u}^*}{\Delta t} = -\frac{1}{\rho}\nabla p

Here, the time derivative advances the solution in time from the intermediate velocity field \mathbf{u}^*, not the old velocity field \mathbf{u}^n, and goes to the next time level at n+1. We have to collect terms on the right-hand side that we have dropped in the previous fractional step, i.e. here the pressure gradient.

So why does that work? Well, it’s easy to prove that both fractional steps are consistent with the original momentum equation, yet I have never seen this proof written down (perhaps I haven’t looked very far?!). But once you have seen it, this fractional step procedure feels very natural. So let’s write down the proof (well, proof is perhaps a strong word, I’ll show you how it works, I think this is more useful). To do that, let’s solve the above equation for the intermediate velocity \mathbf{u}^*. This results in:

\mathbf{u}^*= \frac{\Delta t}{\rho}\nabla p+\mathbf{u}^{n+1}

Let’s insert this expression into the first fractional step:

\frac{\left(\frac{\Delta t}{\rho}\nabla p+\mathbf{u}^{n+1}\right)-\mathbf{u}^n}{\Delta t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = \nu\nabla^2\mathbf{u}

We split the addition within the parenthesis into two fractions, which yields:

\frac{\frac{\Delta t}{\rho}\nabla p}{\Delta t}+\frac{\mathbf{u}^{n+1}-\mathbf{u}^n}{\Delta t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = \nu\nabla^2\mathbf{u}

We simplify the pressure gradient term, as \Delta t appears both in the nominator and denominator. It will vanish and so we are left with:

\frac{1}{\rho}\nabla p+\frac{\mathbf{u}^{n+1}-\mathbf{u}^n}{\Delta t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = \nu\nabla^2\mathbf{u}

Subtracting the pressure gradient, and writing the discretised time derivative as a partial derivative again results in:

\frac{\partial \mathbf{u}}{\partial  t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = - \frac{1}{\rho}\nabla p +\nu\nabla^2\mathbf{u}

As you can see, we are still solving the same equation, we have just split it into two separate equations. If you wanted, you could write as many fractional steps as you have terms. For example, the following fractional step procedure will also solve the same momentum equation:

\frac{\mathbf{u}^{*}-\mathbf{u}^{n}}{\Delta t}=-(\mathbf{u}\cdot\nabla)\mathbf{u} \\[1em]
\frac{\mathbf{u}^{**}-\mathbf{u}^{*}}{\Delta t}=- \frac{1}{\rho}\nabla p \\[1em]
\frac{\mathbf{u}^{n+1}-\mathbf{u}^{**}}{\Delta t}=\nu\nabla^2\mathbf{u} \\[1em]

Go ahead and solve for \mathbf{u}^* and \mathbf{u}^{**} as we did before, and insert into modify the equation until you recover the momentum equation. You will always get back to that equation.

This three stage algorithm looks, in equations, very different from Chorin’s pressure projection method, so why don’t we call it the tom algorithm? Obviously, I am taking the piss here, but this is, in essence, what the SIMPLE algorithm does, and this is, why I have some beef with it. It uses a different mathematical representation, but it solves, in the end, the same equations. Under specific conditions, the SIMPLE algorithm and the pressure projection method become almost identical!

Anyways, so why does this fractional step procedure help us? Well, as we have so any idiot (me) can split the equation in as many fractional steps as they want. The genius idea of Chorin was to isolate the pressure, and the pressure only (as we saw, we only have two fractional steps). This means that the first fractional step depends only on the velocity field, which we can solve. However, the second fractional step consisted of the pressure that we did not know. But let’s have another look at Chorin’s second fractional step. We saw it was given as:

\frac{\mathbf{u}^{n+1}-\mathbf{u}^*}{\Delta t} = -\frac{1}{\rho}\nabla p

Yes, we do not know the pressure, but, since we have a mechanism to compute the intermediate velocity field, \mathbf{u}^* is no longer unknown. So if the pressure is unknown, we have to turn our attention to the other unknown in this equation, i.e. \mathbf{u}^{n+1}. Do we know this value? Not precisely, but we know that the incompressible continuity equation is given as \nabla\cdot\mathbf{u}=0. So the idea now is to use the continuity equation to get rid of \mathbf{u}^{n+1}.

The idea here is that the continuity equation has to be satisfied for each time step. So the velocity field at n and n+1 has to satisfy \nabla\cdot\mathbf{u}=0. This means \nabla\cdot\mathbf{u}^n=\nabla\cdot\mathbf{u}^{n+1}=0. But what about \nabla\cdot\mathbf{u}^*?

The continuity equation is only valid for a fluid that follows the incompressible Navier-Stokes equation. But, we obtained this intermediate velocity field \mathbf{u}^* from the modified momentum equation and therefore the velocity field will not satisfy the Navier-Stokes equation. As a result, \nabla\cdot\mathbf{u}^*\ne 0. This is good news, because we can now multiply each term with the divergence operator as:

\frac{\overbrace{\nabla\cdot\mathbf{u}^{n+1}}^{=0}-\overbrace{\nabla\cdot\mathbf{u}^*}^{\ne 0}}{\Delta t} = -\frac{1}{\rho}\nabla\cdot\nabla p

If you look up the pressure projection method in the literature, you will also hear that we have used the Helmholtz-Hodge decomposition, where we have separated the velocity field into an irrotational and sinusoidal part. Instead of saying “we can decompose the velocity field”, people also say that “the velocity field is projected onto an irrotational and sinusodial space”.

Engineers tend to use the word decomposition, while Mathematicians think in terms of projections. They refer to one and the same thing (at least in this context), and this is where the term pressure projection is coming from. If you buy a book for £50 and there is a 20% VAT on the sales price, then we can decompose the book price into £40 profit for the publisher or author and £10 sales tax for the government. Or we could say that the sales price can be projected onto a publisher profit space and a sales tax space.

So if you hear the terminology Helmholtz-Hodge decomposition, or projection operator, think: \nabla\cdot\mathbf{u}^{n+1}=0. This is why we use it. Fancy words to describe one and the same thing.

As we can see from the equation above, we can remove the term involving \nabla\cdot\mathbf{u}^{n+1} and the right-hand side becomes the laplacian operator for the pressure.

\frac{-\nabla\cdot\mathbf{u}^*}{\Delta t} = -\frac{1}{\rho}\nabla^2 p

Solving this for the pressure results in:

\nabla^2 p=\frac{\rho}{\Delta t}\nabla\cdot\mathbf{u}^*

So, now that we have the intermediate velocity field \mathbf{u}^*, we can compute the pressure field p as a result. Let’s look at the second fractional step again:

\frac{\mathbf{u}^{n+1}-\mathbf{u}^*}{\Delta t} = -\frac{1}{\rho}\nabla p

We haven’t actually used this equation yet, we have only used a derived form of it. So, we do know what the pressure is now, as well as the intermediate velocity field. The only unknown is the velocity \mathbf{u}^{n+1}. So let’s solve this equation for it. This produces:

\mathbf{u}^{n+1}=\mathbf{u}^* -\frac{\Delta t}{\rho}\nabla p

And with that, we have found a way to compute both the pressure and the velocity at the next time step. Cool, let’s write the equation in a form that we can implement, i.e. a discretised form.

The first equation we need to solve is the momentum equation with the pressure removed. We had:

\frac{\mathbf{u}^*-\mathbf{u}^n}{\Delta t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = \nu\nabla^2\mathbf{u}

Using the \mathrm{upwind}(\phi) and \mathrm{diffusion}(\phi) operator from the previous section, we can write this equation as

\frac{\mathbf{u}^*_{i,j,k}-\mathbf{u}^n_{i,j,k}}{\Delta t}+\mathrm{upwind}(\mathbf{u})=\mathrm{diffusion}(\mathbf{u})\\[1em]
\mathbf{u}^*_{i,j,k} = \mathbf{u}^n_{i,j,k} + \Delta t\left[ -\mathrm{upwind}(\mathbf{u}) + \mathrm{diffusion}(\mathbf{u}) \right]

The scalar equations for \mathbf{u}=[u,\, v,\, w can then be written as:

u^*_{i,j,k} = u^n_{i,j,k} + \Delta t\left[ -\mathrm{upwind}(u) + \mathrm{diffusion}(u) \right]\\[1em]
v^*_{i,j,k} = v^n_{i,j,k} + \Delta t\left[ -\mathrm{upwind}(v) + \mathrm{diffusion}(v) \right]\\[1em]
w^*_{i,j,k} = w^n_{i,j,k} + \Delta t\left[ -\mathrm{upwind}(w) + \mathrm{diffusion}(w) \right]

With the intermediate velocity field available, we can solve the pressure Poisson equation, which we obtained as:

\nabla^2 p=\frac{\rho}{\Delta t}\nabla\cdot\mathbf{u}^*

There is one small problem here, can you spot it? We no longer have a time derivative, so we are unable to use a simple explicit time integration. As we have discussed in my article on how to discretise the steady state heat diffusion equation, i.e. \alpha(\partial^2 T/\partial x^2)=0, whenever we want to obtain the solution to an equation where there is no time derivative, we have to use an implicit discretisation. This means we are now evaluating the pressure in the equation above at p^{n+1}.

Using a second-order central scheme for the pressure and a central scheme for the intermediate velocity field, and bringing the term \rho/\Delta t on the left-hand side, the Poisson equation is written in discretised form as:

\frac{\Delta t}{\rho}\frac{p^{n+1}_{i+1,j,k}-2p^{n+1}_{i,j,k}+p^{n+1}_{i-1,j,k}}{\Delta x^2}+\frac{\Delta t}{\rho}\frac{p^{n+1}_{i,j+1,k}-2p^{n+1}_{i,j,k}+p^{n+1}_{i,j-1,k}}{\Delta y^2}+\frac{\Delta t}{\rho}\frac{p^{n+1}_{i,j,k+1}-2p^{n+1}_{i,j,k}+p^{n+1}_{i,j,k-1}}{\Delta z^2}=\\[1em]

\left(\frac{u^*_{i+1,j,k}-u^*_{i-1,j,k}}{2\Delta x}+\frac{u^*_{i,j+1,k}-u^*_{i,j-1,k}}{2\Delta y}+\frac{u^*_{i,j,k+1}-u^*_{i,j,k-1}}{2\Delta z}\right)

Now we factor out all the terms that go in front our pressure. This results in:

p^{n+1}_{i+1,j,k}\left(\frac{\Delta t}{\rho\Delta x^2}\right) + p^{n+1}_{i,j,k}\left(\frac{-2\Delta t}{\rho\Delta x^2}\right) + p^{n+1}_{i-1,j,k}\left(\frac{\Delta t}{\rho\Delta x^2}\right) + \\[1em]
p^{n+1}_{i,j+1,k}\left(\frac{\Delta t}{\rho\Delta y^2}\right) + p^{n+1}_{i,j,k}\left(\frac{-2\Delta t}{\rho\Delta y^2}\right) + p^{n+1}_{i,j-1,k}\left(\frac{\Delta t}{\rho\Delta y^2}\right) +\\[1em]
p^{n+1}_{i,j,k+1}\left(\frac{\Delta t}{\rho\Delta z^2}\right) + p^{n+1}_{i,j,k}\left(\frac{-2\Delta t}{\rho\Delta z^2}\right) + 
p^{n+1}_{i,j,k-1}\left(\frac{\Delta t}{\rho\Delta z^2}\right)
=\\[1em]
\frac{u^*_{i+1,j,k}-u^*_{i-1,j,k}}{2\Delta x}+\frac{u^*_{i,j+1,k}-u^*_{i,j-1,k}}{2\Delta y}+\frac{u^*_{i,j,k+1}-u^*_{i,j,k-1}}{2\Delta z}

Combining presure terms at location i,\, j,\, k, we get:

p^{n+1}_{i,j,k}\left(\frac{-2\Delta t}{\rho\Delta x^2}+\frac{-2\Delta t}{\rho\Delta y^2}+\frac{-2\Delta t}{\rho\Delta z^2}\right) +\\[1em]
p^{n+1}_{i+1,j,k}\left(\frac{\Delta t}{\rho\Delta x^2}\right) + p^{n+1}_{i-1,j,k}\left(\frac{\Delta t}{\rho\Delta x^2}\right) + \\[1em]
p^{n+1}_{i,j+1,k}\left(\frac{\Delta t}{\rho\Delta y^2}\right) + p^{n+1}_{i,j-1,k}\left(\frac{\Delta t}{\rho\Delta y^2}\right) +\\[1em]
p^{n+1}_{i,j,k+1}\left(\frac{\Delta t}{\rho\Delta z^2}\right) + p^{n+1}_{i,j,k-1}\left(\frac{\Delta t}{\rho\Delta z^2}\right)
=\\[1em]
\frac{u^*_{i+1,j,k}-u^*_{i-1,j,k}}{2\Delta x}+\frac{u^*_{i,j+1,k}-u^*_{i,j-1,k}}{2\Delta y}+\frac{u^*_{i,j,k+1}-u^*_{i,j,k-1}}{2\Delta z}

By setting the following coefficients:

a_p = \left(\frac{-2\Delta t}{\rho\Delta x^2}+\frac{-2\Delta t}{\rho\Delta y^2}+\frac{-2\Delta t}{\rho\Delta z^2}\right)\\[1em]
a_w = \left(\frac{\Delta t}{\rho\Delta x^2}\right)\\[1em]
a_e = \left(\frac{\Delta t}{\rho\Delta x^2}\right)\\[1em]
a_n = \left(\frac{\Delta t}{\rho\Delta y^2}\right)\\[1em]
a_s = \left(\frac{\Delta t}{\rho\Delta y^2}\right)\\[1em]
a_f = \left(\frac{\Delta t}{\rho\Delta z^2}\right)\\[1em]
a_b = \left(\frac{\Delta t}{\rho\Delta z^2}\right)\\[1em]
b = \frac{u^*_{i+1,j,k}-u^*_{i-1,j,k}}{2\Delta x}+\frac{u^*_{i,j+1,k}-u^*_{i,j-1,k}}{2\Delta y}+\frac{u^*_{i,j,k+1}-u^*_{i,j,k-1}}{2\Delta z}

We can write our pressure Poisson equation as:

a_p\cdot p^{n+1}_{i,j,k} + a_e\cdot p^{n+1}_{i+1,j,k} + a_w\cdot p^{n+1}_{i,j-1,k} + a_n\cdot p^{n+1}_{i,j+1,k} + a_s\cdot p^{n+1}_{i,j-1,k} + a_f\cdot p^{n+1}_{i,j,k+1} + a_b\cdot p^{n+1}_{i,j,k-1} = b

Or, in compact form:

a_p\cdot p^{n+1}_{i,j,k} + \sum_{f}^{nNeighbours} a_f \cdot p^{n+1}_f = b

Solving this for p^{n+1}_{i,j,k}, i.e. the varible we want to solve for, we get:

p^{n+1}_{i,j,k} = \frac{1}{a_p}\left(b - \sum_{f}^{nNeighbours} a_f \cdot p^{n+1}_f \right)

Now we have the problem that we have p^{n+1} on both the left and right-hand side of the equation. To find a solution to this problem, we modify this equation, which now reads:

p^{n+1,m+1}_{i,j,k} = \frac{1}{a_p}\left(b - \sum_{f}^{nNeighbours} a_f \cdot p^{n+1,m}_f \right)

Here, we have introduced another time index m for the pressure, which indicates that we have to solve the above equation iteratively until we have reached a converged solution, i.e. until p^{n+1,m+1}\approx p^{n+1, m}. We can check that by taking the difference between the two, and if the the difference is below a predefined convergence threshold, for example, \epsilon = 10^{-4}, then we can stop the iterative loop.

The algorithm described above is known as the Jacobi method. It is the slowest method to reach a converged solution for an implicit discretisation. There are much better solution algorithms available, but these require us to provide the coefficient matrix \mathbf{A}, i.e. we need to construct the following system of equations:

\mathbf{A}p=\mathbf{b}

We can write out this equation as:

\begin{bmatrix}
a_p & a_e &  \dots & a_n & \dots & a_f & \dots & 0\\
a_w & a_p & a_e & & \ddots & & \ddots & \vdots \\
\vdots & \ddots & \ddots & \ddots & & \ddots & & a_f\\
a_s & & \ddots & \ddots  & \ddots & & \ddots & \vdots \\
\vdots & \ddots & & \ddots & \ddots & \ddots & & a_n \\
a_b & & \ddots &  & \ddots & \ddots & \ddots & \vdots \\
\vdots & \ddots & & \ddots & & a_w & a_p  & a_e \\
0 & \dots & a_b & \dots &  a_s & \dots & a_w& a_p \\
\end{bmatrix}
p=
\mathbf{b}

The right-hand side vector will contain the divergence of the velocity, as well as boundary conditions close to the boundaries. Once we have constructed both the coefficient matrix \mathbf{A} and the right-hand side vector \mathbf{b}, we can use an algorithm such as the Conjugate Gradient algorithm to solve for the pressure p.

If you have missed it, I have written an entire series on how to implement the Conjugate Gradient algorithm and solve equations like the one above. If you want to speed up your convergence, then this is likely what you want to implement. The Conjugate Gradient algorithm is what you would find in commercial and more serious open-source CFD solvers to solve the implicit equation, in this case, the pressure Poisson equation.

Whichever way to choose for the pressure p, eventually, you will end up with values for p^{n+1}. Once we have these values, we saw that we can obtain the updated velocity field as

\mathbf{u}^{n+1}=\mathbf{u}^* -\frac{\Delta t}{\rho}\nabla p

Using a central difference approximation for the pressure gradient term, this can be written in discretised form as:

u^{n+1}_{i,j,k}=u^*_{i,j,k} -\frac{\Delta t}{\rho}\frac{p^{n+1}_{i+1,j,k}-p^{n+1}_{i-1,j,k}}{2\Delta x}\\[1em]
v^{n+1}_{i,j,k}=v^*_{i,j,k} -\frac{\Delta t}{\rho}\frac{p^{n+1}_{i,j+1,k}-p^{n+1}_{i,j-1,k}}{2\Delta x}\\[1em]
w^{n+1}_{i,j,k}=w^*_{i,j,k} -\frac{\Delta t}{\rho}\frac{p^{n+1}_{i,j,k+1}-p^{n+1}_{i,j,k-1}}{2\Delta x}

And now you know how to implement the pressure projection method as well. But, there is one more family of algorithms to explore. It is time to face the villain; the SIMPLE algorithm and its thug friends (SIMPLER, SIMPLEC, PISO, PIMPLE).

Approximate projection methods: The SIMPLE method and its derivatives

Chorin introduced the notion of pressure projection methods, and we said it was a projection method because we decomposed (or projected) the velocity into an irrotational and sinusoidal part (or irrotational and sinusoidal space). We can write this as \mathbf{u}=\mathbf{u}^\mathrm{irrotational} + \mathbf{u}^\mathrm{sinusodial}. Keep this decomposition/projection in mind for later.

Chorin’s pressure projection algorithm was introduced in 1968, and 4 years later, Patankar and Spalding gave us the SIMPLE algorithm in 1972. SIMPLE is an acronym standing for Semi-Implicit Method for Pressure Linked Equations (actually, that should be SIMfPLE … anyways). We will see later where this semi-implicit nature is coming from. As I have alluded to in the previous sections, the SIMPLE algorithm has so much resemblance to Chorin’s pressure projection method that, in my view, it is, at best, a modification but not a new method in its own right!

For crying out loud, SIMPLE, SIMPLER, SIMPLEC, PISO, and all of these variants all use the same fundamental building blocks, and all derive from Chorin’s pressure projection. We call these methods approximate projection method, and this approximation actually makes SIMPLE and co less rigorous compared to Chorin’s method! The reason we have so many variants of SIMPLE (i.e. SIMPER and SIMPLEC) is because of the approximations made in SIMPLE, and each new variant tries to bring back some rigour into the algorithm.

Well, there are still small differences, but we will explore those differences in detail, and you will see that these differences do not warrant new methods being released to the world.

The other issue I have with the SIMPLE algorithm, and this is not (entirely) Patankar’s fault, is the way the SIMPLE algorithm is introduced in the literature. Sure, in 1972, I understand why people were using staggered grids. And no, I’m not going to introduce staggered grids here as they are outdated. No one (serious about CFD) is still using staggered grids, and we have alternative cures for the issue stemming from co-located grids (storing and solving all variables at the same location, e.g. cell centroid or vertices), typically by using the Rhie-Chow interpolation.

No, we will look at how the SIMPLE algorithm works on a normal grid. Patankar did introduce the SIMPLE algorithm for staggered grids, and most textbooks that talk about the SIMPLE algorithm have pretty much just copied and pasted that description. Yet, any textbook talking about the pressure projection method (which suffers from exactly the same issue) does not talk about the need for staggered grids. So the presentation of the SIMPLE algorithm is just needlessly complex in textbooks, and I remember being very confused about it when I was a student.

To follow the textbook derivation as given by Patankar in his 1980 book Numerical Heat Transfer and Fluid Flow (which is also repeated in Versteeg and Malalasekera’s book, well, their book has pretty much copied and pasted chapter 4-6 from Patankar’s book, which appears in chapter 4-6 in their book), we need to start by discretising the momentum equation. I’ll do a one-dimension consideration here, which we will later generalise to three dimensions.

The momentum equation was given as:

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}

In the SIMPLE algorithm, we assume an implicit discretisation for the velocity. Thus, all velocity components are evaluated at time level n+1. Since we can only solve for one variable per equation, we can’t have an implicit treatment for the pressure, and it is evaluated at time level n (in other words, we assume the pressure to be known, either from the previous timestep or the initial solution). This is what SI in SIMPLE stands for: semi-implicit. It is implicit for velocity but explicit for pressure. The equation can be written as:

\frac{\mathbf{u}^{n+1}-\mathbf{u}^n}{\Delta t} + \left[(\mathbf{u}\cdot\nabla)\mathbf{u}\right]^{n+1} = -\frac{1}{\rho}\nabla p^n + \left[\nu\nabla^2\mathbf{u}\right]^{n+1}

For a one-dimensional flow (which doesn’t really solve anything of interest for incompressible flows, but it does make the derivation easier to follow), we can write the scalar equation as:

\frac{u^{n+1}-u^n}{\Delta t}+\phi\frac{\partial u^{n+1}}{\partial x}=-\frac{1}{\rho}\frac{\partial p^n}{\partial x} + \nu\frac{\partial^2 u^{n+1}}{\partial x^2}

Here, I have introduced the linearisation \phi=u^n. We have to do this because implicit equations are solved with linear systems of equations. While these can be extended to solve non-linear systems as well, they are too expensive to solve for very little additional benefit.

Let’s use a first-order upwind discretisation here for the convective flux, a central scheme for the pressure, and the second-order diffusive derivative. This results in:

\frac{u_i^{n+1}-u_i^n}{\Delta t}+\max(\phi_i,0)\frac{u^{n+1}_i-u^{n+1}_{i-1}}{\Delta x} + \min(\phi_i,0)\frac{u^{n+1}_{i+1}-u^{n+1}_{i}}{\Delta x}=-\frac{1}{\rho}\frac{p^n_{i+1}-p^n_{i-1}}{2\Delta x} + \nu\frac{u^{n+1}_{i+1}-2u^{n+1}_{i}+u^{n+1}_{i-1}}{\Delta x^2}

Next, we rearrange the equation so that all contributions at time level n+1 (in this case, all velocities) are on the left-hand side, while we collect all explicit contributions on the right-hand side. This results in:

\frac{u_i^{n+1}}{\Delta t}+\max(\phi_i,0)\frac{u^{n+1}_i-u^{n+1}_{i-1}}{\Delta x} + \min(\phi_i,0)\frac{u^{n+1}_{i+1}-u^{n+1}_{i}}{\Delta x} - \nu\frac{u^{n+1}_{i+1}-2u^{n+1}_{i}+u^{n+1}_{i-1}}{\Delta x^2}=\frac{u_i^n}{\Delta t}-\frac{1}{\rho}\frac{p^n_{i+1}-p^n_{i-1}}{2\Delta x}

Now we collect terms that multiply with u_{i+1}^{n+1}, u_{i}^{n+1}, and u_{i-1}^{n+1}, which can be written as:

u_{i+1}^{n+1}\left(\frac{\min(\phi_i, 0)}{\Delta x} - \frac{\nu}{\Delta x^2} \right) +
u_{i}^{n+1}\left(\frac{1}{\Delta t} + \frac{\max(\phi_i,0)}{\Delta x} - \frac{\min(\phi_i, 0)}{\Delta x} + \frac{2\nu}{\Delta x^2} \right) +
u_{i-1}^{n+1}\left(-\frac{\max(\phi_i,0)}{\Delta x} - \frac{\nu}{\Delta x^2} \right) = \\[1em]
\frac{u_i^n}{\Delta t}-\frac{1}{\rho}\frac{p^n_{i+1}-p^n_{i-1}}{2\Delta x}

We can now introduce coefficients to make the equation above a bit easier to write:

a_e = \left(\frac{\min(\phi_i, 0)}{\Delta x} - \frac{\nu}{\Delta x^2} \right)\\[1em]
a_w = \left(-\frac{\max(\phi_i,0)}{\Delta x} - \frac{\nu}{\Delta x^2} \right)\\[1em]
a_p = \left(\frac{1}{\Delta t} + \frac{\max(\phi_i,0)}{\Delta x} - \frac{\min(\phi_i, 0)}{\Delta x} + \frac{2\nu}{\Delta x^2} \right)\\[1em]
b_i = \frac{u_i^n}{\Delta t}

This will reduce the equation above to:

a_e\cdot u_{i+1}^{n+1} +
a_p \cdot u_{i}^{n+1} +
a_w \cdot u_{i-1}^{n+1} =
b_i-\frac{1}{\rho}\frac{p^n_{i+1}-p^n_{i-1}}{2\Delta x}

Now we say that velocities and coefficients at the east and west location, i.e. at i+1 and i-1, respectively, are neighbour coefficients to the central coefficient a_p at location i. We therefore write these neighbouring coefficients as:

- a_e\cdot u_{i+1}^{n+1}
- a_w \cdot u_{i-1}^{n+1} =
\sum^{numNeighbours}_{nb=0}a_{nb} \cdot u^{n+1}_{nb}=\sum a_{nb} \cdot u^{n+1}_{nb}

In this example, we have numNeighbours = 2. Inserting this into our equation, we arrive at:

a_p \cdot u_{i}^{n+1} -
\sum a_{nb} \cdot u^{n+1}_{nb}
=
b_i-\frac{1}{\rho}\frac{p^n_{i+1}-p^n_{i-1}}{2\Delta x}

Bringing the summation on the right-hand side, we obtain:

a_p \cdot u_{i}^{n+1}
=
\sum a_{nb} \cdot u^{n+1}_{nb} +
b_i - 
\frac{1}{\rho}\frac{p^n_{i+1}-p^n_{i-1}}{2\Delta x}

This is the starting point for the SIMPLE algorithm in many textbooks, and I just wanted to make sure we get to this equation without missing any steps in the derivation so you can follow this discussion. We solve this equation to obtain a velocity for the next time level n+1 by using the pressure from the previous time step. But we also need to find a way to solve for the pressure and to do so requires us to rewrite the equation a bit.

The main idea behind the SIMPLE algorithm is that we replace the velocity and pressure with guessed values. If we guess values for the pressure that satisfy the above-shown momentum equation, then we have a solution. To verify that the solution is physical, the velocity field also needs to satisfy the continuity equation. This is an iterative procedure so we have to loop over the momentum equation and iteratively update the pressure until both the continuity and momentum equation are satisfied.

Therefore, the first step in the SIMPLE algorithm to construct an equation for the pressure is to rewrite the above-shown momentum equation in terms of guessed velocity and pressure fields. This is achieved by replacing the time indices n and n+1 with *, which now indicate that we are using some guessed values. The equation can thus be written as:

a_p \cdot u_{i}^{*}
=
\sum a_{nb} \cdot u^{*}_{nb} +
b_i - 
\frac{1}{\rho}\frac{p^*_{i+1}-p^*_{i-1}}{2\Delta x}

A guessed velocity and pressure field, or more generally, an arbitrary value \phi^*, can be corrected by some value \phi' to obtain the correct value of \phi^c, This can be expressed as:

\underbrace{\phi^c}_\mathrm{correct}=\underbrace{\phi^*}_\text{guess}+\underbrace{\phi'}_\text{correction}

Therefore, we can also rewrite our momentum equation with velocities at u^{n+1} as the correct velocities, now expressed as u using the above notation. The correct momentum equation becomes:

a_p \cdot u_{i}^c
=
\sum a_{nb} \cdot u_{nb}^c +
b_i - 
\frac{1}{\rho}\frac{p^c_{i+1}-p^c_{i-1}}{2\Delta x}

If we now subtract the momentum equation containing \phi^{*} (i.e. the guessed velocity and pressure values) from the correct momentum equation (i.e. the one containing \phi^c), then we obtain the following equation:

a_p \cdot (u_{i}^c-u_i^*)
=
\sum a_{nb} \cdot (u_{nb}^c - u_i^* ) - 
\frac{1}{\rho}\frac{(p^c_{i+1} - p^*_{i+1})-(p^c_{i-1} - p^*_{i-1})}{2\Delta x}

Since both equations contain the term b_i, and we have subtracted one equation from the other, the b_i coefficient has vanished. Three equations earlier, we saw that we can express the correction of a quantity \phi as \phi'=\phi^c-\phi^*. Therefore, we can rewrite the above momentum equation as:

a_p \cdot u_{i}'
=
\sum a_{nb} \cdot u_{nb}' - 
\frac{1}{\rho}\frac{p'_{i+1}-p'_{i-1}}{2\Delta x}

This equation now expresses a correction equation. We still can’t solve it, but it shows us how the pressure field would have to change in response to a change in velocity. We will use this equation in a second to derive a correction equation for the pressure. But before we do that, let’s take a look at how the SIMPLE algorithm treats this equation. It drops the first term on the right-hand side, i.e. the summation over neighbouring velocities and their coefficients. This means the above equation simplifies to:

a_p \cdot u_{i}'
=
- \frac{1}{\rho}\frac{p'_{i+1}-p'_{i-1}}{2\Delta x}

This is the main approximation in the SIMPLE algorithm, and this term is reintroduced (in one way or another) in the SIMPLER (SIMPLE Revised) and SIMPLEC (SIMPLE Consistent) algorithms. In order to make use of this equation, we insert the definition for the correct velocity into this equation. That is, u'=u^c-u^*. This results in:

a_p \cdot (u_i^c - u_i^* )
=
- \frac{1}{\rho}\frac{p'_{i+1}-p'_{i-1}}{2\Delta x}

Solving this equation for the correct velocity u^* gives us:

u_i^c=u^*_i-\frac{1}{\rho a_p}\frac{p'_{i+1}-p'_{i-1}}{2\Delta x}

Let us write out the coefficient for a_p. As a reminder, a_p was given as:

a_p = \left(\frac{1}{\Delta t} + \frac{\max(\phi_i,0)}{\Delta x} - \frac{\min(\phi_i, 0)}{\Delta x} + \frac{2\nu}{\Delta x^2} \right)

This results in the following expanded velocity correction equation:

u_i^c=u^*_i-\left(\frac{\Delta t}{\rho} + \frac{\Delta x}{\max(\phi_i,0)\rho} - \frac{\Delta x}{\min(\phi_i,0)\rho} +\frac{\Delta x^2}{2\nu\rho}\right)\frac{p'_{i+1}-p'_{i-1}}{2\Delta x}

Let us now introduce the following coefficient:

a_{\Delta x}=\frac{\Delta x}{\max(\phi_i,0)\rho} - \frac{\Delta x}{\min(\phi_i,0)\rho} +\frac{\Delta x^2}{2\nu\rho}

Let us also write the velocity correction equation just once using u^c=u^{n+1}, as we assume that the correct velocity must be the velocity at the next time level n+1. Then, we can write the equation above as:

u_i^{n+1}=u^*_i-\left(\frac{\Delta t}{\rho} + a_{\Delta x} \right)\frac{p'_{i+1}-p'_{i-1}}{2\Delta x}

Let’s compare that to the velocity correction equation in Chorin’s pressure projection algorithm:

u^{n+1}_{i}=u^*_{i} -\frac{\Delta t}{\rho}\frac{p^{n+1}_{i+1}-p^{n+1}_{i-1}}{2\Delta x}

I don’t know about you, but to me, they look pretty similar. OK, but I hear the Pantankar fanboys (and girls) screaming that there is a coefficient of a_{\Delta x}, which depends on \Delta x. Good observation. Let’s look at the CFL number for a second, it is given as:

CFL=u\frac{\Delta t}{\Delta x}\quad \rightarrow\quad \Delta x = u\frac{\Delta t}{CFL}

Thus, the spatial distance \Delta x is proportional to \Delta t/CFL. In other words, as the CFL number is increasing, we have \Delta x << \Delta t. Therefore, only for small CFL numbers (of the order of one or less) does the coefficient a_{\Delta x} really contribute to the above velocity correction equation. For large CFL values (say, of the order of 10 and larger), both SIMPLE and the pressure projection method of Chorin behave essentially the same.

OK, let’s finish this little detour and come back to the SIMPLE algorithm. We have stolen from the pressure projection, apologies, I meant derived, the velocity correction equation, which was given as:

u_i^c=u^*_i-\frac{1}{\rho a_p}\frac{p'_{i+1}-p'_{i-1}}{2\Delta x}

Now, we take the continuity equation of an incompressible fluid, i.e. \nabla\cdot \mathbf{u} = 0, and we want to derive a correction equation for the pressure. The continuity equation will only hold for a physically correct velocity field. In the pressure projection algorithm, we said that using the Helmholtz-Hodge decomposition, we can project (decompose) the velocity field onto an irrotational and sinusoidal part. That is, we have:

\mathbf{u}=\mathbf{u}^\text{irrotational}+\mathbf{u}^\text{sinusodial}

This is exact, i.e. this holds true for any vector field in a mathematical sense. Therefore, we say that Chorin’s projection method is an exact pressure projection. How about the simple method? Well, we already introduced the correction equation, which was given as:

\phi^c=\phi^*+\phi'

We use different terms here, like guessed or initial value and corrected value, but in the end, both concepts represent the same thing. However, since in a strict mathematical sense, we need to iteratively satisfy the equation for \\phi^c, we say that this is an approximate pressure projection method.

For the pressure projection, we said that the velocity at time level n+1 satisfies the continuity equation, i.e. \nabla\cdot\mathbf{u}^{n+1}=0, but the intermediate velocity \mathbf{u}^* does not, i.e. we have \nabla \cdot \mathbf{u}^*\ne 0.

By analogy, the same is true for the SIMPLE algorithm. We have \nabla \cdot \mathbf{u}^c = 0, but \nabla \cdot \mathbf{u}^* \ne 0. Therefore, if we want to use the continuity equation, we have to use it with respect to the correct velocity u^c, for which we have derived the correction equation above.

Let us first write the scalar, discretised continuity equation using a central differencing scheme, still in one dimension. This results in:

\nabla\cdot\mathbf{u}=\nabla\cdot\mathbf{u}^c=\nabla \cdot u^c=\frac{\partial u^c}{\partial x}=\frac{u^c_{i+1}-u^c_{i-1}}{2\Delta x}=0

Now we insert the velocity correction equation, which results in:

\frac{1}{2\Delta x}\left[\left(u^*_{i+1}-\frac{1}{\rho a_p}\frac{p'_{i+2}-p'_{i}}{2\Delta x}\right)-\left(u^*_{i-1}-\frac{1}{\rho a_p}\frac{p'_{i}-p'_{i-2}}{2\Delta x}\right)\right]=0

In this equation, everything is known except for the pressure. The velocity u^* is the guessed velocity (i.e. from the previous time step or initial solution) while a_p contains only known coefficients. Therefore, we can solve this equation for the pressure now. We start by expanding the brackets. This leads to:

\frac{u^*_{i+1}}{2\Delta x}-\frac{1}{\rho a_p}\frac{p'_{i+2}-p'_{i}}{4\Delta x^2}-\frac{u^*_{i-1}}{2\Delta x}+\frac{1}{\rho a_p}\frac{p'_{i}-p'_{i-2}}{4\Delta x^2}=0

Now we multiply by \rho a_p and we merge fractions, this results in:

\rho a_p\frac{u^*_{i+1}-u^*_{i+1}}{2\Delta x}-\frac{p'_{i+2}-2p'_{i}+p'_{i-2}}{4\Delta x^2}=0

Now, we place the velocity on the right-hand side, which results in:

\frac{p'_{i+2}-2p'_{i}+p'_{i-2}}{4\Delta x^2}=\rho a_p\frac{u^*_{i+1}-u^*_{i+1}}{2\Delta x}

Can we express the pressure term on the left-hand side differently? Well, let’s see what happens if we expand a Taylor series around i+2 and i-2. Then we get:

p(x+2\Delta x)=p(x)+\frac{\partial p}{\partial x}(2\Delta x)+\frac{\partial^2 p}{\partial x^2}\frac{(2\Delta x)^2}{2!}+\mathcal{O}(\Delta x^3)\\[1em]
p(x-2\Delta x)=p(x)-\frac{\partial p}{\partial x}(2\Delta x)+\frac{\partial^2 p}{\partial x^2}\frac{(2\Delta x)^2}{2!}+\mathcal{O}(\Delta x^3)

Let’s add these two equation together:

p(x+2\Delta x) + p(x-2\Delta x) = 2p(x) + \frac{2}{2!}\frac{\partial^2 p}{\partial x^2}(2\Delta x)^2

If we solve this for the second-order derivative, we obtain (with (2\Delta x)^2=2^2\Delta x^2=4\Delta x^2):

\frac{\partial^2 p}{\partial x^2} = \frac{p(x+2\Delta x) - 2p(x) + p(x-2\Delta x)}{4\Delta x^2}

If we now introduce discrete locations for x, where we have location i at x, i+2 at x+2\Delta x and i-2 at x-2\Delta x, then we have:

\frac{\partial^2 p}{\partial x^2} = \frac{p_{i+2} - 2p_i + p_{i-2}}{4\Delta x^2}

Therefore, our pressure correction equation, i.e.:

\frac{p'_{i+2}-2p'_{i}+p'_{i-2}}{4\Delta x^2}=\rho a_p\frac{u^*_{i+1}-u^*_{i+1}}{2\Delta x}

can now be written as:

\frac{\partial^2 p'}{\partial x^2}=\rho a_p\frac{\partial u^*}{\partial x} \quad\rightarrow\quad \nabla^2 p'=\rho a_p\nabla\cdot \mathbf{u}^*

Thus, we have to solve this pressure Poisson equation to get values for the corrected pressure values p'. Looks familiar? Well, let’s see how we solved the pressure in Chorin’s pressure projection method, we had:

\nabla^2 p=\frac{\rho}{\Delta t}\nabla\cdot\mathbf{u}^*

Let us expand a_p again in the pressure correction equation. Then we have:

\nabla^2 p'=\left(\frac{\rho}{\Delta t}+\frac{1}{a_{\Delta x}}\right)\nabla\cdot \mathbf{u}^*

We see again the same coefficient as in the pressure projection method of Chorin. In addition, we get a_{\Delta x} again, though now we have to take the inverse, i.e. 1/a_{\Delta x}.

This is the SIMPLE algorithm of Patankar and Spalding. It appeared 4 years after the pressure projection method of Chorin was published, and people got accused and found guilty of far stronger cases of academic misconduct. I will say this, though: Chorin published his work in predominantly mathematical circles, while Patankar’s work was published in more engineering-orientated circles. If you had to conduct a literature review back in the 1970s (without the internet and Google), your literature review could only be as good as the printed journals available in your library.

While I’m sure the work of Chorin would have been available to Patankar and Spalding, they may not have thought of looking through applied mathematical journals. Thus, even though I was critical of the work of Patankar, I suppose he has plausible deniability. This doesn’t mean, though, that we have to keep using the SIMPLE algorithm or any of its siblings when the pressure projection algorithm is essentially the same, with a more rigorous mathematical foundation to support its derivation.

There is one main difference between these two methods that we have not yet looked at. And that is the momentum equation. In the SIMPLE algorithm, we retain the pressure gradient, while in the pressure projection algorithm of Chorin, we do not. We saw that in the SIMPLE algorithm, this leads to a semi-implicit behaviour in the pressure, where the velocity is solved at time level n+1 while the pressure is solved explicitly at time level n.

The pressure projection method of Chorin is fully implicit (if we want, we can also solve the equation fully explicit). Since only the velocity appears in the first fractional step of the momentum equation, we can solve the equation fully implicitly, which means that we can choose a CFL number that is arbitrarily large. In the SIMPLE algorithm, we have to solve the equations implicitly, but we are restricted in the CFL number due to the explicit contribution from the pressure.

Furthermore, since we have made quite a few approximations with the simple algorithm (e.g. dropping the neighbouring contributions of velocity and its coefficients), the SIMPLE algorithm is prone to slow convergence. Thus, we can summarise the SIMPLE algorithm as:

Slow convergence + limited CFL number (despite implicit treatment) + less mathematical rigour + plagiarism = most used pressure velocity coupling algorithm for incompressible flows.

CFD is fascinating if you look into the details, isn’t it?

There is one more thing I want to do before I let you go. I want to take the equations we have derived here and put them into a matrix form. This uses a different notation, but it is one you will find in the literature as well, especially when you are dealing with OpenFOAM. Seeing them in one place will help to connect the equations when you find them.

We start with our correction momentum equation. This was given as:

a_p \cdot u_{i}^{n+1}
=
\sum a_{nb} \cdot u^{n+1}_{nb} +
b_i - 
\frac{1}{\rho}\frac{p^n_{i+1}-p^n_{i-1}}{2\Delta x}

If we collect velocities at time level n+1 on the left-hand side, we get:

a_p \cdot u_{i}^{n+1} -
\sum a_{nb} \cdot u^{n+1}_{nb} 
=
b_i - 
\frac{1}{\rho}\frac{p^n_{i+1}-p^n_{i-1}}{2\Delta x}

Let us now expand the summation again so we see all velocities explicitly. This is given as:

a_p \cdot u_{i}^{n+1} -
a_{e} \cdot u^{n+1}_{e} - 
a_{w} \cdot u^{n+1}_{w}
=
b_i - 
\frac{1}{\rho}\frac{p^n_{i+1}-p^n_{i-1}}{2\Delta x}

Instead of using indices for the east and west nodes, we write this equation index i instead. Thus, the node at location e is located at i+1 and w is located at i-1. The coefficient a_p is now also evaluated at location i, just like the velocity that it multiplies. This is then given as:

-a_{i-1} \cdot u^{n+1}_{i-1} +
a_i \cdot u_{i}^{n+1} -
a_{i+1} \cdot u^{n+1}_{i+1}
=
b_i - 
\frac{1}{\rho}\frac{p^n_{i+1}-p^n_{i-1}}{2\Delta x}

This equation can then be written in matrix form:

\begin{bmatrix}
a_1 & -a_2 & 0 & \dots & & 0 & 0 & 0\\[1em]
-a_1 & a_2 & -a_3 &0 & \dots & & & 0\\[1em]
0 & -a_2 & a_3 & -a_4 & 0 & \dots & & 0\\[1em]
\vdots & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots\\[1em]
\vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots\\[1em]
0 &  & \dots & 0 & -a_{n-3} & a_{n-2} & -a_{n-1} & 0\\[1em]
0 &  & & \dots & 0 & -a_{n-2} & a_{n-1} & -a_{n}\\[1em]
0 & 0 & 0 & & \dots & 0 & -a_{n-1} & a_{n}\\[1em]
\end{bmatrix}

\begin{bmatrix}
u_1\\[1em]
u_2\\[1em]
u_3\\[1em]
\vdots\\[1em]
\vdots\\[1em]
u_{n-2}\\[1em]
u_{n-1}\\[1em]
u_n\\[1em]
\end{bmatrix}
=
\begin{bmatrix}
b_1-\frac{1}{\rho}\frac{p_2-p_\text{left boundary}}{2\Delta x}\\[1em]
b_2-\frac{1}{\rho}\frac{p_3-p_1}{2\Delta x}\\[1em]
b_3-\frac{1}{\rho}\frac{p_4-p_2}{2\Delta x}\\[1em]
\vdots\\[1em]
\vdots\\[1em]
b_{n-2}-\frac{1}{\rho}\frac{p_{n-1}-p_{n-3}}{2\Delta x}\\[1em]
b_{n-1}-\frac{1}{\rho}\frac{p_n-p_{n-2}}{2\Delta x}\\[1em]
b_n-\frac{1}{\rho}\frac{p_{\text{right boundary}}-p_{n-1}}{2\Delta x}\\[1em]
\end{bmatrix}

Now let’s say that we are at the second node, i.e. i=2. Looking at the matrix equation above, we could write the equation for node 2 as

-a_1\cdot u_1^{n+1} + a_2\cdot u_2^{n+1} - a_3 \cdot u_3^{n+1} = b_2 - \frac{1}{\rho}\frac{p_3^n - p_1^n}{2\Delta x}

We can see that this is equivalent to the equation we expressed for the general location i before. We can now write this matrix equation in a more compact form:

\mathbf{\mathcal{M}}\mathbf{U}^{n+1}=\mathbf{RHS}(\mathbf{b}(\mathbf{U}^n), \nabla p^n)

This is the equation we are solving to get an initial velocity field. But we also use this equation to derive the correction equation for the pressure. We saw that the first step requires us to derive a correction equation, which we saw was given as:

a_p \cdot u_{i}'
=
\sum a_{nb} \cdot u_{nb}' - 
\frac{1}{\rho}\frac{p'_{i+1}-p'_{i-1}}{2\Delta x}

We see that the term \mathbf{b} has vanished. Therefore, our corrected momentum equation can be written in matrix form as:

\mathbf{\mathcal{M}}\mathbf{U}'=-\frac{1}{\rho}\nabla p'

Now we separate the matrix \mathbf{\mathcal{M}} into two matrices. The first matrix will only contain diagonal coefficients, and we call this matrix \mathbf{\mathcal{A}}. The second matrix will contain all the off-diagonal components, and we shall name this matrix \mathbf{\mathcal{H}} (I am using the naming convention here that you will find in the literature). This can be written in matrix form as:

\mathbf{\mathcal{M}}\mathbf{U}'=\mathbf{\mathcal{A}}\mathbf{U}'-\mathbf{\mathcal{H}}

If we expand these matrixes, for example, for a simple 1D grid with 5 nodes, then we have:

\underbrace{\begin{bmatrix}
a_1 & -a_2 & 0 & 0 & 0\\[1em]
-a_1 & a_2 & -a_3 & 0 & 0\\[1em]
0 & -a_2 & a_3 & -a_4 & 0\\[1em]
0 & 0 & -a_3 & a_4 & -a_5\\[1em]
0 & 0 & 0 & -a_4 & a_5
\end{bmatrix}}_{\mathbf{\mathcal{M}}}
\underbrace{\begin{bmatrix}
u_1 \\[1em] u_2 \\[1em] u_3 \\[1em] u_4 \\[1em] u_5
\end{bmatrix}}_{\mathbf{U}'}
=
\underbrace{\begin{bmatrix}
a_1 & 0 & 0 & 0 & 0\\[1em]
0 & a_2 & 0 & 0 & 0\\[1em]
0 & 0 & a_3 & 0 & 0\\[1em]
0 & 0 & 0 & a_4 & 0\\[1em]
0 & 0 & 0 & 0 & a_5
\end{bmatrix}}_{\mathbf{\mathcal{A}}}
\underbrace{\begin{bmatrix}
u_1 \\[1em] u_2 \\[1em] u_3 \\[1em] u_4 \\[1em] u_5
\end{bmatrix}}_{\mathbf{U}'}
-
\underbrace{\begin{bmatrix}
0 & -a_2\cdot u_2' & 0 & 0 & 0\\[1em]
-a_1\cdot u_1' & 0 & -a_3\cdot u_3' & 0 & 0\\[1em]
0 & -a_2\cdot u_2' & 0 & -a_4\cdot u_4' & 0\\[1em]
0 & 0 & -a_3\cdot u_3' & 0 & -a_5\cdot u_5'\\[1em]
0 & 0 & 0 & -a_4\cdot u_4' & 0
\end{bmatrix}}_{\mathbf{\mathcal{H}}}

Inserting this definition, we have the momentum equation of the form:

\mathbf{\mathcal{A}}\mathbf{U}'-\mathbf{\mathcal{H}}=-\frac{1}{\rho}\nabla p'

We can also write this equation with the off-diagonal matrix on the right-hand side, which results in:

\mathbf{\mathcal{A}}\mathbf{U}'=\mathbf{\mathcal{H}}-\frac{1}{\rho}\nabla p'

Compare that to the discrete form we had obtained previously:

\underbrace{a_p}_\mathbf{\mathcal{A}} \cdot \underbrace{u_{i}'}_{\mathbf{U}'}
=
\underbrace{\sum a_{nb} \cdot u_{nb}'}_\mathbf{\mathcal{H}} -
\frac{1}{\rho}
\underbrace{\frac{p'_{i+1}-p'_{i-1}}{2\Delta x}}_{\nabla p'}

Because the matrix \mathbf{\mathcal{A}} only contains diagonal components, we can easily invert this matrix. That is, the diagonal matrix can be inverted as:

\mathbf{\mathcal{A}}^{-1}=
\begin{bmatrix}
\frac{1}{a_1} & 0 & \dots & 0\\[1em]
0 &\frac{1}{a_2} &  & 0\\[1em]
\vdots & & \ddots & \vdots \\[1em]
0 & & \dots & \frac{1}{a_n}
\end{bmatrix}

Therefore, we can multiply by this inverted matrix, which will provide us with the following correction equation:

\mathbf{U}'=\mathbf{\mathcal{A}}^{-1}\mathbf{\mathcal{H}}-\frac{1}{\rho}\mathbf{\mathcal{A}}^{-1}\nabla p'

For the SIMPLE algorithm, we made the simplification of dropping the neighbouring contributions, that is, we drop the first term on the right-hand side. This results in the following simplified equation:

\mathbf{U}'=-\frac{1}{\rho}\mathbf{\mathcal{A}}^{-1}\nabla p'

We also said that a correction \phi' can be expressed as \phi'=\phi^c - \phi^*, i.e. it is the result from subtracting the guessed value \phi^* from the correct value \phi^c. So, rewriting the above equation in terms of correct and guessed velocity results in:

\mathbf{U}^c=\mathbf{U}^*-\frac{1}{\rho}\mathbf{\mathcal{A}}^{-1}\nabla p'

The discrete equation we obtained was:

u_i^c=u^*_i-\frac{1}{\rho a_p}\frac{p'_{i+1}-p'_{i-1}}{2\Delta x}

In the literature, though, the term \mathbf{\mathcal{A}}^{-1}\mathbf{\mathcal{H}} isn’t dropped from the derivation. Instead, we keep it in the derivation and then if we want to have the pure SIMPLE algorithm, we simply set the coefficients in \mathbf{\mathcal{H}} to zero. If we would like to use the SIMPLEC alternative, for example, those coefficients are not zero, and so we have a unified way of expressing both the SIMPLE algorithm and its derivatives.

In the next step, we insert the equation for the corrected velocity into the continuity equation \nabla\cdot \mathbf{U}. This results in:

\nabla\cdot\mathbf{U}'=\nabla\cdot(\mathbf{\mathcal{A}}^{-1}\mathbf{\mathcal{H}})-\frac{1}{\rho}\nabla\cdot (\mathbf{\mathcal{A}}^{-1}\nabla p')=0

Placing the pressure on the left-hand side and the velocity on the right-hand side, we have the following pressure Poisson equation:

\frac{1}{\rho}\nabla\cdot (\mathbf{\mathcal{A}}^{-1}\nabla p')=\nabla\cdot(\mathbf{\mathcal{A}}^{-1}\mathbf{\mathcal{H}})

After we have solved this equation for the corrected pressure, we have a value for the correct pressure through p^c=p^*+p'. Then, we use this correct pressure field to get the correct velocity field. To do that, we rewrite our correction equation, i.e.

\mathbf{U}'=\mathbf{\mathcal{A}}^{-1}\mathbf{\mathcal{H}}-\frac{1}{\rho}\mathbf{\mathcal{A}}^{-1}\nabla p'

in terms of correct quantities:

\mathbf{U}^c=\mathbf{U}^{n+1}=\mathbf{\mathcal{A}}^{-1}\mathbf{\mathcal{H}}-\frac{1}{\rho}\mathbf{\mathcal{A}}^{-1}\nabla p^c

This is the SIMPLE algorithm in matrix form, and you see that both this form and the discrete equation form solve one and the same thing.

Regardless of which form you implement, you will see that the SIMPLE algorithm doesn’t work, at least not until we have made some modifications. You see, we made crude approximations to the momentum equation, such as linearisation and dropping some terms entirely. While we derived the equation for an unsteady problem, the resolution of the SIMPLE algorithm in time isn’t great due to all the simplifications made. Coupled with convergence issues, it is best used for steady state problems only.

We need to introduce under-relaxation at this point, which is the process of updating the flow field only partially with the newly obtained quantities. If we did not use under-relaxation, the inaccurate SIMPLE algorithm would overshoot the correct solution slightly in each iteration/time step, and due to the non-linear nature of our equations, these inaccuracies would grow exponentially in time, ultimately causing our simulation to diverge.

For the generic quantity \phi, we can apply under-relaxation as:

\phi^{n+1}=(1-\alpha)\phi^n+\alpha\phi^{n+1},\quad\quad 0<\alpha\le 1

This will ensure that we are only slightly updating velocity and pressure for each iteration/time step. The pressure projection method of Chorin does not require under-relaxation, though in some cases, it may help to stabilise the simulation.

How to go from the SIMPLE algorithm to the PISO algorithm

I have mentioned that modifications to the matrix \mathbf{\mathcal{H}} will determine which version of the SIMPLE algorithm we are using. But there is one other algorithm which is commonly used, and that is the PISO algorithm, which stands for Pressure Implicit with Splitting of Operators (or PIwSoO, though PISO is probably more catchy).

The PISO algorithm is essentially the SIMPLE algorithm, but we simply insert an additional loop. We solve the same equations and start with the momentum equation (I’ll be using the matrix notation here, which is just a bit more compact):

\mathbf{\mathcal{M}}\mathbf{U}^{n+1}=\mathbf{RHS}(\mathbf{b}(\mathbf{U}^n), \nabla p^n)

Now, we start a loop that is not present in the SIMPLE algorithm. We say that we solve the following equations within a loop:

\frac{1}{\rho}\nabla\cdot (\mathbf{\mathcal{A}}^{-1}\nabla p')=\nabla\cdot(\mathbf{\mathcal{A}}^{-1}\mathbf{\mathcal{H}})

and

\mathbf{U}'=\mathbf{\mathcal{A}}^{-1}\mathbf{\mathcal{H}}-\frac{1}{\rho}\mathbf{\mathcal{A}}^{-1}\nabla p'

Why? Because the matrix \mathbf{\mathcal{H}} depends on the velocity \mathbf{U}. We use the first equation to obtain a new pressure field, and then we use the second equation to correct the velocity field based on this new pressure. As soon as we update the velocity, \mathbf{\mathcal{H}} changes as well since it depends on the velocity, which just changed. Thus, we loop a few times over these two equations until we satisfy them.

We don’t do that in the SIMPLE algorithm, and this is the reason why we had to introduce under-relaxation here because of this inaccuracy here.

Thus, the PISO algorithm extends the SIMPLE algorithm by adding an additional loop. If you ask me, that is an extension to the SIMPLE algorithm rather than a new method in its own right, but I suppose Patankar and Spalding set the precedence by making the SIMPLE method its own algorithm instead of an extension to Chorin’s pressure projection method, so I suppose this modification counts as a new method.

Bonus: the PIMPLE algorithm

If you have never used OpenFOAM (which means you must still have your sanity; look after it before it is gone!), then you will likely never have heard about the PIMPLE algorithm. It is a combination of PISO and SIMPLE, hence the name PIMPLE (great choice, especially when you are googling it. It is only beaten by the carbuncle problem in CFD, though I would strongly suggest not to google carbuncle, or at least not while you are eating!)

It is an unpublished algorithm that extends the PISO algorithm by, well, who would have thought, yet another loop! As we have established, adding a new loop warrants the introduction of a new method, and so we have the PIMPLE algorithm.

It can be summarised as:

• Start loop
• Solve PISO algorithm
• Finish loop

Yep, that’s PIMPLE for you. If you run PIMPLE with just one iteration, it will even tell you that it is running in PISO mode. Why would you want to do that? We can belittle the PIMPLE algorithm all we want, but it does remove one essential shortcoming of both the SIMPLE and PISO algorithm, and that is the problem of the semi-implicit nature of the pressure.

We said that because the pressure is evaluated at the previous time level n and not at time level n+1, we are limited in the CFL numbers we can use before the simulation starts to diverge. The PIMPLE algorithm introduces this additional loop to stabilise the solution algorithm, and, as a result, we can run simulations at much higher CFL numbers.

It is an expensive trade-off, but if we can increase the CFL number by a factor of, say 5 times, but the additional computation only takes 3 times longer, then we are saving computational cost, as we can now integrate our solution much more efficiently in time.

But it is an unnecessary one. If we look back at Chorin’s pressure projection method, it is fully implicit (we have dropped the pressure gradient from the equation). As a result, we can use CFL numbers that are as high as we want. We are only limited by the temporal resolution we want to achieve. From a numerical point, there is nothing holiding us back to use a high CFL number here.

But, there are times when we may want to use an additional loop as we do in the PIMPLE algorithm. This idea was introduced in 1991 by Jameson and it is known these days as dual time stepping. All we do is introduce a second time derivative in the momentum equation. This can be written as:

\frac{\partial \mathbf{u}}{\partial t}+\frac{\partial \mathbf{u}}{\partial \tau} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}

The way that this works is to solve the inner part of the equation, i.e.:

\frac{\partial \mathbf{u}}{\partial \tau} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}

as a steady-state problem. We advance this solution in time until the time derivative becomes zero. If this is the case, then we have obtained the solution at the next time level n+1 for the velocity. Since we solve this equation iteratively, we are also solving the pressure correction equation within each iteration, and so it does not matter if the pressure is evaluated at time level n, as we will iteratively satisfy this momentum equation.

With sufficient iterations, the above equation can be solved, and then we can update the velocity through the real-time derivative. The procedure isn’t complicated, but it takes some time to wrap your head around. If you are not fully understanding the concept or have more questions, I get it. It took me some time to understand as well how to implement it.

But the takeaway message is that we can use dual time stepping to cheaply and efficiently solve the inner iterations until we get a steady state result, which we can then use to integrate in time with a large CFL number. This is another good idea for an article I have to write in the future, I have placed it on my ever growing to do list.

This is essentially what the PIMPLE algorithm does, i.e. it solves the momentum equation iteratively. We can see it as an extension of the PISO algorithm using a dual-time stepping procedure.

The verdict

OK, I want to summarise the stupidity that is the SIMPLE algorithm and all of its derivatives. I really want to get the point across that SIMPLE and co are all sub-optimal implementations of a perfectly fine method that has been around for longer than SIMPLE and any of its friends. We shouldn’t be afraid of using the pressure projection algorithm.

Let’s start with the pressure projection method of Chorin. It is fully implicit, rooted in an exact mathematical derivation and treatment, and as a result it does not require any under-relaxation to make the method work.

If we want to mildly annoy ourselves, then we add the pressure gradient back into the equation (which is wrong from a mathematical standpoint), and as a result, we lose the fully implicit nature. We fix this by adding a dual time stepping procedure, which allows us to have arbitrarily large CFL numbers again, but we have just increased the computational cost compared to Chorin’s original method.

Next, we drop the dual time stepping so we have a more efficient solution, but we are now limited in the maximum CFL number that we can use. This is the PISO algorithm. If we further butcher the pressure Poisson equation by removing non-linear contributions, we are losing accuracy and solving a related but not exact discretised form of the Navier-Stokes equation. As a result, we have to introduce under-relaxation and typically want to limit ourselves to using the SIMPLE algorithm for steady-state flows only.

I think I have made my point.

Summary

OK, let’s collect ourselves, there was a lot of information in one article. So let’s summarise the main points here:

• Incompressible and compressible flows require different solution algorithm due to the density changing its behaviour. We differentiate between density-based descriptions, where we solve for the density and then use an equation of state to link it to the pressure, and pressure-based descriptions, where the density is may be constant, and we link pressure and velocity through some iterative procedure.
• We saw that we have conservative and non-conservative (primitive) variable formulations, and we prefer conservative variables for compressible flows, especially when we have discontinuouties in our flow (e.g. shock waves). For smooth data, the variable description does not matter and we prefer primitive variables for incompressible flows. They make life easier.
• Compressible flows do not change the Navier-Stokes equations, and we solve them as they are. We augment them with an equation of state to link density and pressure, which closes the system.
• For incompressible flows, the density is considered constant, and so we can no longer solve the equations as they are. Instead, we have to construct new equations that derive from either the momentum or continuity equation, that will iteratively link the pressure field to the velocity.
• There are two main families of algorithms that we have looked at; preconditioning like the Artificial Compressibility method, and pressure projection methods, like Chorin’s exact projection, and approximate projections like SIMPLE, PISO, and PIMPLE.

It is common for people specialising in CFD to pick on flow regime (i.e. incompressible or compressible) and then sticking with it. To the best of my judgement, this is due to differences in variable formulation and solution algorithms, as well as differences in numerical schemes (e.g. shock-capturing schemes for compressible flows vs. upwind-based schemes for incompressible flows). It almost feels like both incompressible and compressible flows have almost nothing in common.

I really hope you will be able to see how to solve both sides, and that you will not be scared writing a solver for either incompressible or compressible flows. Now that we have a good idea of how to discretise the Navier-Stokes equations, which schemes to use, how to treat boundaries, and what algorithms to use to solve incompressible and compressible flows, we are ready to face the biggest problem in CFD; how to deal with turbulent flows. This will be discussed in the next article, I’ll see you over there!

---

Tom-Robin Teschner is a senior lecturer in computational fluid dynamics and course director for the MSc in computational fluid dynamics and the MSc in aerospace computational engineering at Cranfield University.