The SIMPLE Algorithm

for Pressure-Velocity Coupling

ME 448/548 Notes

Gerald Recktenwald
Portland State University
Department of Mechanical Engineering
gerry@pdx.edu

February 10, 2020

ME 448/548 SIMPLE

Overview

1. Basic ideas
2. Mesh ideas: structured versus unstructured, and staggered
3. SIMPLE: a classic segregated solution method for incompressible flows
4. Convergence checking
5. Summary and recommendations

ME 448/548 SIMPLE page 1

 Basic Ideas: High-level View of CFD

CFD modeling involves these steps

1. Define the geometry of the fluid volume

2. Locate boundary regions and assign values to boundary condition parameters
3. Specify thermophysical properties and physics, e.g., buoyancy, turbulence, radiation
4. Create the mesh
5. Define parameters that control the solution
6. Solve the flow model
7. Vizualize the results

Commercial CFD packages have preprocessors (graphical user interfaces) for setting up
the model and viewing the results. The order of preprocessing steps may be di↵erent for
di↵erent CFD packages.
Users control step 6 by specifying parameters the a↵ect how the solver updates the
velocity, pressure and other scalar fields. That control is indirect, at best.

ME 448/548 SIMPLE page 2

Meshes:
Structured,
Unstructured,
and Staggered

ME 448/548 SIMPLE page 3

 Structured and Unstructured Meshes

The SIMPLE algorithm can be developed for non-staggered, unstructured meshes. We

will use the case of a staggered, structured mesh because the nomenclature is simpler. It
is also the historical origin of the method.

Structured Meshes: Cells (nodes) are arranged in rows and columns

Lx
ny

∆y
Ly
...

j=1
y i=1 2 3 ... nx
∆x
x

Each cell has the same number of neighbors. The connectivity is regular and uniform.

ME 448/548 SIMPLE page 4

Structured and Unstructured Meshes

Unstructured Meshes: Cells and nodes are not arranged in rows and columns.

y y

x x

ME 448/548 SIMPLE page 5

 Structured and Unstructured Meshes

Trimmer and Polyhedral volume meshes in STAR-CCM+

ME 448/548 SIMPLE page 6

Structured and Unstructured Meshes

Unstructured Meshes: Cells and nodes can have arbitrary connectivity.

ME 448/548 SIMPLE page 7

 Structured and Unstructured Meshes

Unstructured Meshes: Connectivity can vary with position in the domain.

ME 448/548 SIMPLE page 8

Staggered Mesh

Main characteristics
Main control volume
• Velocities are computed at cell faces N N
υn
• Pressure is computed at cell centers uw P ue W P E
W υs E υs

S S

Control volume and neighbors of υs

N

uw P
W E

Control volume and neighbors of uw

ME 448/548 SIMPLE page 9

 General Polyhedral Mesh Elements in STAR-CCM+

See Theory ! Numerical Flow Solution ! Finite Volume Discretization in the on-line
STAR-CCM+ documentation (“Help” menu).

ME 448/548 SIMPLE page 10

Basic Discretization – 2D (1)

∆x

N
υn
For a two-dimensional, incompressible flow, applying the finite uw P ue
∆y
υs
volume method on a staggered mesh leads to a set of three equations W E

for each interior control volume. S

u u u u u pi 1,j pi,j (u)

aS ui,j 1 a W ui 1,j + aP ui,j aE ui+1,j aN ui,j+1 = + bo
xw
v v v v v pi,j 1 pi,j (v)
aS vi,j 1 a W vi 1,j + aP vi,j aE vi+1,j aN vi,j+1 = + bo
ys
ui+1,j ui,j vi,j+1 vi,j
+ =0
x y

where b(u)
o and b(v)
o are source terms due to non-uniform viscosity, gravity, and other
e↵ects.

ME 448/548 SIMPLE page 11

 Basic Discretization – 2D (2)

These equations are nonlinear. For example, if the central di↵erence scheme is used,
u
aE = F (ui,j , ui+1,j , ⌫, geometry)

The solution process – even when a coupled solver is used – requires these basic steps

1. Make a guess at the current flow field

2. Compute the coefficients of the (now) linearized momentum and continuity equations
3. Solve the equations

These steps are repeated until the guessed flow field does not result in a change in the
coefficients.

ME 448/548 SIMPLE page 12

Basic Discretization – 2D (3)

Recall the structure of the coefficient matrix for the finite-volume approximation to the
2D Poisson equation.

𝜑k–nx

𝜑k–1
𝜑k
–aS –aW aP –aE –aN

𝜑k+1
= bk

𝜑k–nx

Applying the finite-volume method to the 2D incompressible flow equations give three
matrices with this structure. However, the unknowns for each of these matrices are
coupled. (Actually, there’s a bit of a problem with the pressure equation.)

ME 448/548 SIMPLE page 13

 Basic Discretization – 2D (4)

The coupled system of equations looks like this for a staggered, structured, 2D mesh.

u bou

v = bou

p 0

Diagonal lines in the coefficient matrix represent the position of non-zero entries.

ME 448/548 SIMPLE page 14

Coupled or Segregated Solution

We have a two basic choices

u bou
1. Solve the large sparse system as a
single, linearized system
2. Use sequential, a.k.a. segregated v = bou
method to separately solve each
linearized system for the dependent
vectors u, v , and p.
p 0

SIMPLE uses the segregated approach

ME 448/548 SIMPLE page 15

 Refer to the STAR-CCM+ documentation
Theory ! Numerical Flow Solution
! Finite Volume Discretization
! Segregated Flow Solver

ME 448/548 SIMPLE page 16

Segregated Solution Procedure (1)

SIMPLE: Semi-Implicit Method for Pressure-Linked Equations

Preview: For laminar flow in a 2D domain

Unknowns: u, v , p: x-direction velocity, y -direction velocity and pressure

1. Solve u equation assuming v and p are frozen
2. Solve v equation assuming u and p are frozen
3. Solve a pressure correction equation for p that enforces continuity
4. Update u, v and p from the pressure correction

References: Ferziger and Perić, §,7.5.1, pp. 188–196

Patankar
Tu, Yeoh, and Liu, §4.3.3

ME 448/548 SIMPLE page 17

 Segregated Solution Procedure (1)

Apply the finite-volume discretization separately to each momentum equation

• For the u equation, v and p are frozen
(u) (u)
A u=b

• For the v equation, u and p are frozen

(v) (v)
A v=b

Each of these matrix equations is a linearized version of the finite-volume model.

The momentum source terms, b(u) and b(v), consist of pressure terms and other terms
due to velocity gradients and non-uniform viscosity.
(u) (u) (u) (v) (v) (v)
b = bo bp b = bo bp
⇤
b(u) (v)
p and bp are contributions due to the frozen pressure field (p ).

ME 448/548 SIMPLE page 18

Segregated Solution Procedure (2)

Discretization of momentum equations results in a non-linear, coupled system of

equations for u and v (in 2D).
(u) (u) (v) (v)
A u=b A v=b

• A(u) is the coefficient matrix, and b(u) is the source term (vector) for the u equation
when v and p are frozen.
• A(v) is the coefficient matrix, and b(v) is the source term (vector) for the v equation
when u and p are frozen.
• Momentum equations are nonlinear, so A(u) and A(v) depend on u and v .
• Nonlinearity requires iterative solution: A(u) and A(v) are updated with new values of
u and v , and the systems are solved again.

What about the continuity equation?

Where is the equation for pressure?

ME 448/548 SIMPLE page 19

 Pressure in Momentum Source Terms

Pressure is the dominant source term in the momentum equations

(u) (u) (u) (u) (u) (u)
A u=b where b = bP + bo ⇠ du(pW p P ) + . . . + bo

(v) (v) (v) (v) (v) (v)

A v=b where b = bP + bo ⇠ dv (pS p P ) + . . . + bo

The du and dv are coefficients that depend on the geometry of of the control volume.

Note: The ⇠ symbol is meant to convey the basic features of the terms, and allows a
degree of notational abuse. To make these expressions more precise, read the
⇠ symbol as “contains terms like”.

Pressure does not have its own natural

equation for incompressible flow.

ME 448/548 SIMPLE page 20

Pressure Correction and Continuity (1)

Solution to the linearized momentum equations ∆x

(u) (u) (v) (v) N

A u=b A v=b
υn
uw P ue
gives updated u and v fields. ∆y
W υs E

In general these new u and v fields will not satisfy the continuity equation for each cell

uw y ue y + vs x vn x 6= 0

ME 448/548 SIMPLE page 21

 Pressure Correction and Continuity (2)
Assume that we can find a velocity correction at each node that adjusts the velocity field
so that it satisfies the continuity equation.
0
uw = ũw + uw (1)
0
vs = ṽs + vs (2)

• uw and vs are the velocity components that satisfy the continuity equation.
• ũw and ṽs are the velocity components obtained by satisfying the momentum
equations.
• u0w and vs0 are the corrections to the velocity components needed to satisfy the
continuity equation.

Equations (1) and (2) define u0w , and vs0 values for each cell in the domain. The ũw and
ṽs fields are obtained by solving the momentum equations. We need to develop another
procedure for computing the u0w and vs0 fields. Ultimately the u0w and vs0 fields are
obtained by requiring each cell to satisfy the continuity equation.

ME 448/548 SIMPLE page 22

Pressure Correction and Continuity (3)

Assume that the velocity corrections can be derived from suitable pressure corrections
0
pP = p̃P + pP (3)

• pP is the desired pressure at the cell center, i.e., the one that brings the velocity field
into balance so that continuity equation is satisfied for each cell.
• p̃P is the current guess at the pressure at the cell center.
• p0P is the change in p̃P needed to satisfy the continuity equation.

ME 448/548 SIMPLE page 23

 Pressure Correction and Continuity (4)

Unifying assumption:

Assume that the velocity corrections are uniquely determined by the pressure corrections.
0 0 0
uw ⇡ cw (pW pP ) (4)
0 0 0
vs ⇡ cs(pS pP ) (5)

where cw and cs are coefficients that depend on terms from the momentum equation.

(We’ve skipped a couple of steps here.)

ME 448/548 SIMPLE page 24

Pressure Correction and Continuity (5)

Tying the velocity corrections to the pressure corrections allows us to close the system of
equations. An equation for the p0 is obtained by substituting the Equations (1), (2), (4),
and (5) into the continuity equation. The result is

(p0 ) 0 (p0 )
A p =b

where
(p0 )
b ⇠ ũw y ũe y + ṽs x ṽn x

In words: The source term for the p0 field is the local error in the continuity equation
that results from updating the u and v fields from the previous iteration to
get to ũ and ṽ .

ME 448/548 SIMPLE page 25

 The Pieces (in 2D)

The momentum equations

(u) (u) (v) (v)
A u=b A v=b ( ?)

The pressure correction equation

(p0 ) 0 (p0 )
A p =b (??)

Momentum and pressure corrections

0 0 0
u = ũ + u v = ṽ + v p = p̃ + p

ũ and ṽ are the solutions to Equations (?) obtained when the previous field variables (u,
v , and p) are used to compute A(u), b(u), A(v) and b(v). p̃ is the pressure field from the
previous global iteration.
0
p0 is the solution to Equation (??) when b(p ) is computed from the newly obtained ũ
and ṽ .

ME 448/548 SIMPLE page 26

SIMPLE Algorithm

Repeat the following steps until convergence:

1. Compute A(u) and b(u), and solve A(u)u = b(u) to obtain ũ, a tentative guess at the
x-direction velocity field.
2. Compute A(v) and b(v), and solve A(v)v = b(v) to obtain ṽ , a tentative guess at the
y -direction velocity field. Note that A(v) and b(v) use the most recently available ũ
field, and the v field from the previous iteration
0 0 0 0
3. Compute A(p ) and b(p ), and solve A(p )p0 = b(p ) to obtain the pressure correction
field p0.
4. Correct the velocities and pressures.
0 0 0
pP p̃P + pP uw ũw + uw vs ṽs + vs

5. Test for convergence, and return to Step 1 if necessary.

The SIMPLE algorithm is described by Patankar [2], and Ferzier and Perić [1, § 7.5.1],
and Tu et al.[3, §4.3.3].

ME 448/548 SIMPLE page 27

 Inner and Outer Iterations

The iteration number in the STAR-CCM+ residual plot is the number of outer iterations.

1. Linearize u equation with v and p frozen

Iteratively solve u equation Inner
iteration

2. Linearize v equation with u and p frozen

Iteratively solve v equation
Outer Inner
iteration
iteration
3. Assemble pressure correction equation
Iteratively solve the pressure
Inner
correction equation iteration

4. Correct u, v and p

Solution of the linearized equations in steps 1, 2 and 3 also use iterative methods (as
opposed to direct matrix solution). Those iterations are called “inner iterations”.

ME 448/548 SIMPLE page 28

Convergence Checking (1)

Each outer iteration of the SIMPLE loop, A(u) and b(u) are updated with the newest
information about the u, v and p fields.

Before solving
(u) (u)
A u=b (6)
compute the residual
(u) (u) (u)
r =b A u (7)
Note that the u in these equations is from the previous iteration.

If r (u) is small, then the old u is a good approximation to the system of equations defined
by the newly updated coefficients in A(u) and b(u).

ME 448/548 SIMPLE page 29

 Convergence Checking (2)

Similarly,
(v) (v) (v)
r =b A v
is a measure of how well the old v satisfies the linearized equations based on the newly
updated coefficients in A(v) and b(v)

Each of the unknown fields, u, v , p (in 2D) will have a residual.

The residual plot in StarCCM+ shows the residuals of each of the dependent vector
variables.

ME 448/548 SIMPLE page 30

Convergence Checking (3)

r (u) and r (v) are residual vectors: there is a residual for each interior cell in the domain.

To check convergence, we need to look at a scalar value, not an entire vector.

Let be a normalized scalar value of the residual.

kruk krv k
u = , v =
kru,ref k krv,ref k

The values of kru,ref k and krv,ref k are chosen so that u ⇠ 1 and v ⇠ 1 on the first
iteration, or close to the first iteration.

ME 448/548 SIMPLE page 31

 Convergence Checking (4)

For the pressure correction equation, the residuals are not as informative as the source
term
(p0 )
b ⇠ ũw y ũe y + ṽs x ṽn x
0
If b(p ) is small, then the current guess at the velocity field comes close to satisfying the
continuity equation before the pressure and velocity corrections are applied.

To check convergence of the continuity equation, monitor c

0
kb(p )k1
c =
kb(p0)k1,ref

ME 448/548 SIMPLE page 32

Monitoring Residuals

Convergence history for a well-behaved, laminar flow simulation with STAR-CCM+:

Residuals

0.1

0.01
Continuity
Residual

X−momentum
Y−momentum
Z−momentum
0.001

1E−4

1E−5
0 100 200 300 400 500
Iteration

Remember: The iteration number in this plot refers to the outer iterations.

ME 448/548 SIMPLE page 33

 Monitoring Residuals

Convergence history for a less well-behaved, turbulent flow simulation with STAR-CCM+:
Residuals

1000

100

10

Continuity
1
X−momentum
Residual

Y−momentum
0.1
Z−momentum
Tke
0.01
Tdr

0.001

1E−4

1E−5
0 200 400 600 800 1000 1200
Iteration

ME 448/548 SIMPLE page 34

Monitoring Residuals

Non-convergence for a laminar flow that is unsteady, but the user is trying to use a steady
solution method. A flow with lower inlet velocity was used for an initial guess.

ME 448/548 SIMPLE page 35

 Summary

• SIMPLE is a segregated solution method: The u, v , w, and p fields are solved

separately. Coupling between these field variables is achieved via velocity and pressure
corrections.
• Convergence occurs when the residuals for each of the equations is reduced to a value
below a tolerance.
• The “residuals” are normalized scalars. The normalization is chosen so that ⇠ 1 on
the first iteration.
• c is the largest local (cell-wise) error in the continuity equation before the
pressure-based corrections to the velocity fields are applied. After the velocity field is
corrected, local mass conservation is obtained for all cells. Thus, at the end of each
iteration of SIMPLE, both global and local mass conservation is guaranteed.
• By default, StarCCM+ does not use a stopping criterion other than limiting the
simulation 10 1000 outer iterations. Turning on convergence checking requires user
action.

ME 448/548 SIMPLE page 36

STAR-CCM+ Gives You a Choice

When specifying the Physics Continuum:

Coupled flow solves the large system of equations for u, v , w and p simultaneously (with
an iterative method).

Segregated flow uses SIMPLE for steady problems or PISO for transient problems.

ME 448/548 SIMPLE page 37

 Adding a Convergence Criterion in Star-CCM+

To add a convergence criterion in Star-CCM+

1. Right click on Stopping Criteria
2. Select Create New Criterion!From
Monitor. . . and select one of the field
variables.
I recommend using Continuity as a baseline
convergence criterion.

For hard problems, one of the turbulence field

variables or temperature might also be used as a
convergence criterion.

The same convergence criteria can be used with

coupled and segregated solvers.

ME 448/548 SIMPLE page 38

Recommendations

• Always inspect the residual plot after each run.

• If the residual tolerances have not been met, continue the solution by increasing the
maximum allowable number of steps. This approach continues the solution from the
last iteration.
• For some hard problems, other adjustments to the solution parameters might be
necessary in order obtain convergence.
• STAR-CCM+ also has a coupled flow solver, which may or may not converge better
(and faster) than the segregated solver. The coupled flow solver uses more memory
than the segregated solver.

ME 448/548 SIMPLE page 39

References

[1] Joel H. Ferziger and Milovan Perić. Computational Methods for Fluid Dynamics.
Springer-Verlag, Berlin, third edition, 2001.
[2] S.V. Patankar. Numerical Heat Transfer and Fluid Flow. Hemisphere, Washington
D.C., 1980.
[3] Jiyuan Tu, Guan Heng Yeoh, and Chaoqun Liu. Computational Fluid Dynamics: A
Practical Approach. Butterworth-Heinemann, 2008.

ME 448/548 SIMPLE page 40