9.

PRE- AND POST-PROCESSING SPRING 2021

9.1 Stages of a CFD analysis

9.2 The computational mesh
9.3 Boundary conditions
9.4 Flow visualisation
Appendix: Summary of vector calculus
Examples

9.1 Stages of a CFD Analysis

Pre-Processing

This stage includes:

• creating geometry;
• generating a grid;
• specifying the equations to be solved;
• specifying boundary conditions.

It may depend upon:

• the desired outputs (e.g. force coefficients, heat transfer, …);
• the capabilities of the solver.

Solving

With commercial CFD codes the solver is often operated as a “black box”. Nevertheless, user
intervention is needed to specify solution methods and choose discretisation schemes.
Convergence must be checked and under-relaxation factors changed if necessary.

Post-Processing

The raw output of the solver is the set of values of each field variable (𝑢, 𝑣, 𝑤, 𝑝, …) at each
node. Relevant data must be extracted and manipulated. For example, a subset of surface
pressures, shear stresses and cell-face areas is required in order to compute a drag coefficient.

Most commercial CFD vendors supplement their flow solvers with grid-generation and flow-
visualisation tools, all operated from a graphical user interface (GUI), which simplifies the
setting-up of a CFD simulation and the subsequent extraction and manipulation of data.

Popular commercial CFD codes include:

• STAR-CCM+ (Siemens)
• Fluent (ANSYS)
• FLOW-3D (Flow Science)
• PowerFLOW (Dassault Systèmes)
• COMSOL (part of the COMSOL multiphysics suite)

CFD 9–1 David Apsley

Popular non-commercial CFD solvers include
• OpenFOAM ([Link]
• CodeSaturne ([Link]
These are open-source codes and you would be expected to have a reasonably good knowledge
of C++ or Fortran to use them.

An excellent web portal for all things CFD is [Link]

9.2 The Computational Mesh

9.2.1 Mesh Structure

The purpose of the grid generator is to decompose the flow domain into control volumes.

The primary outputs are:

• cell vertices;
• connectivity information.

Various arrangements of nodes and vertices may be employed.

u p

cell-centred storage cell-vertex storage staggered velocity mesh

The shapes of control volumes may depend on the capabilities of the solver. Structured-grid
codes use quadrilaterals in 2D and hexahedra in 3D. Unstructured-grid solvers often use
triangles (2D) or tetrahedra (3D), but newer codes can use arbitrary polyhedra.

hexahedron tetrahedron

In all cases it is necessary to specify connectivity; that is, which cells are adjacent to each other,
and which faces and vertices they share. For structured grids with (𝑖, 𝑗, 𝑘) numbering this is
straightforward, but for unstructured grids advanced data structures are required.

CFD 9–2 David Apsley

9.2.2 Areas and Volumes

To calculate a flux requires the vector area A of a cell face. This has magnitude equal to the
area and direction normal to it. Its orientation may also distinguish the outward direction.
𝑚𝑎𝑠𝑠 𝑓𝑙𝑢𝑥 = ρu • A (1)

𝑑𝑖𝑓𝑓𝑢𝑠𝑖𝑣𝑒 𝑓𝑙𝑢𝑥 = −Γ∇ϕ • A (2)

To find the total amount of some property in a cell requires its volume 𝑉:
𝑎𝑚𝑜𝑢𝑛𝑡 = ρ𝑉ϕ (3)

Areas

The fundamental unit of surface is a triangle.

Triangles.
The vector area of a triangle with side vectors s1 and s2 is the cross product s2
1
A
A = s1 ∧ s2 (4)
2 s1
The orientation depends on the order of vectors in the cross product. A

Quadrilaterals
4 (or more) points do not, in general, lie in a plane. However, since the sum of the outward
vector areas from any closed surface is zero, i.e.
r4
∑ A𝑓 = 0 r3
∮ dA = 0 or (5)
∂𝑉 faces

any surface spanning the same set of points has the same vector area. Hence,
it can be divided into planar triangles and their vector areas summed. A

By adding vector areas of, e.g., triangles 123 and 134, the vector area of any r2
r1
surface spanned by these points and bounded by the side vectors is found
(see the Examples) to be half the cross product of the diagonals:
1 1
A= d13 ∧ d24 = (r3 − r1 ) ∧ (r4 − r2 ) (6)
2 2

General polygons
The vector area of an arbitrary polygonal face may be found by breaking
it up into triangles and summing the individual vector areas. Because it
spans the same set of vertices, the overall vector area is independent of
how it is decomposed.

CFD 9–3 David Apsley

Volumes

If r ≡ (𝑥, 𝑦, 𝑧) is the position vector, then

∂𝑥 ∂𝑦 ∂𝑧
∇•r= + + =3
∂𝑥 ∂𝑦 ∂𝑧
Integrating over an arbitrary control volume:

∫ ∇ • r d𝑉 = 3𝑉
𝑉

Dividing by 3 and using the divergence theorem (see Appendix) gives the volume of a cell:
1
𝑉= ∮ r • dA (7)
3 ∂𝑉

If a polyhedral cell has plane faces this can be evaluated as

1
𝑉= ∑ r𝑓 • A𝑓 (8)
3
faces

where A𝑓 is the vector area of face f and r𝑓 the position vector of any point on that face. It
doesn’t matter which point is chosen since, for any other two vectors r1 and r2 on that face,
r2 • A𝑓 − r1 • A𝑓 = (r2 − r1 ) • A𝑓 =0
This is zero because r2 – r1 lies in the plane and hence is perpendicular to A𝑓 .

If the cell faces are not planar, then the volume depends on how each face is
broken down into triangles. A consistent approach is to connect the vertices
surrounding a particular face with a common central point: for example, the
rf
average of the vertices:
1
r𝑓 = ∑ r𝑖 (9)
𝑁
vertices

The (possibly non-planar) face is then the assemblage of triangles connecting to r𝑓 .

s3
Using the general formula, it is readily shown that the volume of a s1
tetrahedron formed from side vectors s1 , s2 , s3 (taken in a right-handed
sense) is s2
1
𝑉= s • s ∧ s3 (10)
6 1 2

CFD 9–4 David Apsley

2-d Cases

In 2-d cases consider 3-d cells to be of unit depth. The “volume” of the cell is then
V
numerically equal to its planar area.

The side “area” vectors are most easily obtained from projected areas:
Δ𝑦
A=( ) (11) CELL s y
−Δ𝑥 y
0
where edge Δs = (Δ𝑥, Δ𝑦) is taken anticlockwise around the cell. A

x
0
9.2.3 Cell-Averaged Derivatives -x

The average value of any function 𝑓 over a cell is

𝑡𝑜𝑡𝑎𝑙 1
𝑓𝑎𝑣 ≡ = ∫ 𝑓 d𝑉 (12)
𝑣𝑜𝑙𝑢𝑚𝑒 𝑉 𝑉

By applying the divergence theorem (see the Appendix) to (ϕ, 0,0) = ϕe𝑥 ,
∂ϕ 1 ∂ϕ
( ) ≡ ∫ d𝑉
∂𝑥 𝑎𝑣 𝑉 𝑉 ∂𝑥
1
= ∫ ∇ • (ϕ, 0,0) d𝑉
𝑉 𝑉
1
= ∫ ∇ • (ϕe𝑥 ) d𝑉
𝑉 𝑉
1
= ∮ ϕe𝑥 • dA
𝑉 ∂𝑉
1
= ∮ ϕ d𝐴𝑥
𝑉 ∂𝑉

Considering all three derivatives, and assuming polyhedral cells this becomes
∂ϕ/ ∂𝑥
1
(∂ϕ/ ∂𝑦) = ∑ ϕ𝑓 A𝑓 (13)
𝑉
∂ϕ/ ∂𝑧 𝑓𝑎𝑐𝑒𝑠
𝑎𝑣

If the cell is cartesian this reduces to the expected form since

area A

∂ϕ ϕ𝑒 − ϕ𝑤 ϕ𝑒 𝛢 − ϕ𝑤 𝛢 ϕ𝑒 𝛢𝑒𝑥 + ϕ𝑤 𝛢𝑤𝑥
= = = V
∂𝑥 Δ𝑥 𝐴Δ𝑥 𝑉
Here, the only outward area vectors with non-zero 𝑥 components are
A𝑒 = (𝐴, 0,0), A𝑤 = (−𝐴, 0,0). x

CFD 9–5 David Apsley

Classroom Example 1

A tetrahedral cell has vertices at A(2, –1, 0), B(0, 1, 0), C(2, 1, 1) and D(0, –1, 1).

(a) Find the outward vector areas of all faces. Check that they sum to zero.

(b) Find the volume of the cell.

(c) If the values of ϕ at the centroids of the faces (indicated by their vertices) are
ϕ𝐵𝐶𝐷 = 5, ϕ𝐴𝐶𝐷 = 3, ϕ𝐴𝐵𝐷 = 4, ϕ𝐴𝐵𝐶 = 2
find the volume-averaged derivatives (∂ϕ/ ∂𝑥)𝑎𝑣 , (∂ϕ/ ∂𝑦)𝑎𝑣 , (∂ϕ/ ∂𝑧)𝑎𝑣 .

Classroom Example 2

In a 2-dimensional unstructured mesh, one cell is a pentagon. The coordinates of the vertices
are as shown in the figure, whilst the average values of a scalar ϕ on edges a to e are:
ϕ𝑎 = −7, ϕ𝑏 = 8, ϕ𝑐 = −2, ϕ𝑑 = 5, ϕ𝑒 = 0
(1,3)
b
c
(5,1)
(-3,0)

d a

(-2,-3)
e
(2,-4)
Find:
(a) the area of the pentagon;

(b) the cell-averaged derivatives (∂ϕ/ ∂𝑥)𝑎𝑣 , and (∂ϕ/ ∂𝑦)𝑎𝑣 .

CFD 9–6 David Apsley

Classroom Example 3 (Exam 2016)

The figure shows the vertices of a single triangular cell in a 2-d unstructured finite-volume
mesh. The accompanying table shows the pressure 𝑝 and velocity (𝑢, 𝑣) on the faces marked
a, b, c at the end of a steady, incompressible flow simulation.

(3,5)

edge 𝑝 𝑢 𝑣
c a 3 5 –1
b b 5 7 –5
c 2 uc –6
y, v
(1,1)
a
x, u (4,0)

Find:
(a) the area of the triangle;

(b) the net pressure force (per unit depth) on the cell;

(c) the outward volume flow rate (per unit depth) for all faces;

(d) the missing velocity component 𝑢𝑐 ;

(e) the cell-averaged velocity gradients 𝜕𝑢/𝜕𝑥, 𝜕𝑢/𝜕𝑦, 𝜕𝑣/𝜕𝑥, 𝜕𝑣/𝜕𝑦.

(f) Define, mathematically, the acceleration (material derivative) Du /D𝑡. If the velocity
at the centre of the cell is u = (16/3, −4), use this and the gradients from part (e) to
calculate the acceleration.

Data. The cell-averaged derivative of ϕ in volume 𝑉 is given by

∂ϕ 1
( ) = ∑ ϕ𝑓 𝐴𝑓,𝑖
∂𝑥𝑖 𝑎𝑣 𝑉
𝑓𝑎𝑐𝑒𝑠 𝑓

CFD 9–7 David Apsley

9.2.4 Classification of Meshes

Structured meshes are those whose control volumes can be indexed by (𝑖, 𝑗, 𝑘) for 𝑖 = 1, . . . , 𝑛𝑖 ,
𝑗 = 1, . . , 𝑛𝑗 , 𝑘 = 1, . . . , 𝑛𝑘 , or by a group of such blocks (multi-block structured meshes).

Structured grids can be cartesian (lines parallel to 𝑥𝑦𝑧 axes) or curvilinear (usually curved to
fit boundaries). The grid is orthogonal if all grid lines cross at 90; examples include cartesian,
cylindrical-polar or spherical-polar grids. Some flows can be treated as axisymmetric, and in
these cases, the flow equations can be expressed in terms of polar coordinates (𝑟, θ, 𝑧), rather
than cartesian coordinates, with minor modifications.

Cartesian mesh Curvilinear mesh

Unstructured meshes are useful for complex geometries. However, data structures, solution
algorithms, grid generators and plotting routines for such meshes are very complex.

Unstructured cartesian mesh Unstructured triangle mesh

9.2.5 Fitting Complex Boundaries With Structured Meshes

Blocking Out Cells

Some bluff-body flows can be computed on single-block

cartesian meshes by blocking out cells. Solid-surface boundary
conditions are applied to cell faces abutting the blocked-out
region, whilst values of flow variables are forced to zero within
the blocked-out region by a modification of the source term for
those cells. If the scalar-transport equations for a single cell are discretised as

ϕ𝑃 − ∑ 𝑎𝐹 ϕ𝐹 = 𝑏𝑃 + 𝑠𝑃 ϕ𝑃

then the source-term coefficients are modified:

𝑏𝑃 → 0, 𝑠𝑃 → −(𝑙𝑎𝑟𝑔𝑒 𝑛𝑢𝑚𝑏𝑒𝑟) (14)
Rearranging for ϕ𝑃 , this gives
∑ 𝑎𝐹 ϕ𝐹
ϕ𝑃 = (15)
𝑎𝑃 + 𝑙𝑎𝑟𝑔𝑒 𝑛𝑢𝑚𝑏𝑒𝑟
ensuring that the computational variable ϕ𝑃 is effectively forced to zero in these cells.
However, the computer still stores values and carries out operations for these points, so that it
is essentially performing a lot of redundant work. A better approach is to fit several structured-
mesh blocks around the body. Multi-block grids are discussed below.

CFD 9–8 David Apsley

Volume-of-Fluid (VOF) Approach

In the volume-of-fluid approach the fraction 𝑓 of f=0

each cell filled with fluid is stored: 𝑓 = 0 for cells
with no fluid, 𝑓 = 1 for cells completely within the f=1
interior of the fluid and 0 < 𝑓 < 1 for cells cutting 0<f<1
a boundary. At solid boundaries 𝑓 is determined by
the surface contour. At moving free surfaces a f=0
time-dependent equation is solved for 𝑓, with the
surface being reconstructed from its values.

A related technique is the level-set method where each phase-change surface (e.g. free surface
or solid bed) is a contour of an indicator variable 𝑓, which changes smoothly, rather than
discontinuously as in the VOF approach. For example, 𝑓 might vary between 0 at the bed and
1 at the free surface.

Body-Fitted Meshes

The majority of general-purpose codes allow body-fitted grids.

Mesh lines are distorted so as to fit curved boundaries. Accuracy
in turbulent-flow calculations demands a high density of grid
cells close to solid surfaces, but often the grid need only be
refined in the direction normal to the surface.

The use of non-cartesian grids has important consequences.

• It is necessary to store a lot of geometric data for each

control volume; for example, in our multi-block-structured code STREAM we need to
store (𝑥, 𝑦, 𝑧) components of the cell-face-area vectors for “east”, “north” and “top”
faces of each cell, plus the volume of the cell itself – a total of 10 arrays.
=
• Unless the mesh is orthogonal, the diffusive flux is no longer const.
dependent only on the nodes immediately on either side of a face;
e.g. for the east face: st.
con
=
∂ϕ ϕ𝐸 − ϕ𝑃 E
Γ 𝐴 ≠ −Γ ( )𝐴 P
∂𝑛 Δ𝑃𝐸
The derivative of ϕ normal to the face involves cross-derivative
terms parallel to the cell face and nodes other than P and E. These
extra terms are usually transferred to the source term and have to
be treated explicitly, which tends to reduce stability.

• Interpolated values of all three velocity components are needed v

to evaluate the mass flux through a single face. This requires
approximations in the pressure-correction equation that can u
slow down convergence.

CFD 9–9 David Apsley

Multi-Block Structured Meshes

In multi-block structured meshes the domain is decomposed into a small number of regions, in
each of which the mesh is structured.

Grid lines may match at the interfaces between blocks, so that the
edge vertices are common to two blocks. Alternatively, arbitrary
interfaces permit block boundaries where vertices do not match
2 3 4
and interpolation is required. An important example is a sliding
interface, used in rotating machinery such as pumps or turbines. 1 5
On each iteration of a scalar-transport equation the discretised equations are solved implicitly
within each block, with values from the adjacent blocks providing internal boundary conditions
which are updated explicitly at the end of the iteration. (In an arrangement where cell vertices
match at block boundaries and blocks also meet “whole-face to whole-face”, this can be done
more accurately by extending and overlapping adjacent blocks.)

It is generally desirable to avoid sharp changes in grid

2 3 4 direction and/or non-orthogonality (which lead to lower
1 6 7 8 5 accuracy and stability problems) in rapidly-changing
regions of the flow.

Many distributed-memory parallel-processing CFD implementations (e.g. STREAM, with

MPI parallelisation) use blockwise domain decomposition, with one or more blocks assigned
to each computer processor. The distribution of blocks amongst processors tries to:
• share load (roughly, cell count) as equally as possible; otherwise some processors will
be idly waiting for others to catch up at the end of an iteration;
• minimise the transfer of information between cores (which is inherently much slower
than in-core transfers).

Chimera (or Overset) Meshes

Some solvers, including STAR-CCM+, allow overlapping

blocks (chimera, or overset, meshes). Here, there is a
simple background mesh with one or more overset meshes
that exchange information with it. These are particularly
useful for moving objects, since the overset mesh can
move with the object to which it is attached, whilst the
background mesh is geometrically simple and stationary.

9.2.6 Disposition of Grid Cells

To resolve detail, more points are needed in rapidly-changing regions of the flow such as:
• solid boundaries;
• shear layers;
• separation, reattachment and impingement points;
• flow discontinuities (shocks, hydraulic jumps).

CFD 9 – 10 David Apsley

Simulations must demonstrate grid-independence; i.e. that a finer-resolution grid would not
significantly change the numerical solution. This usually requires a sequence of calculations
on successively finer grids.

Boundary conditions for turbulence models impose limitations on cell sizes near walls. Low-
Reynolds-number models (resolving the viscous sublayer) typically require that 𝑦 + < 1 for the
near-wall node, whilst high-Reynolds-number models (relying on wall functions) require the
near-wall node to lie in the log-law region, say 30 < 𝑦 + < 150. This is difficult to prescribe a
priori; many codes now permit 𝑦 + -dependent blending between the two approaches.

9.2.7 Multi-Grid Methods

Multiple grids – combining cells so that there are 1, 1/2, 1/4, 1/8, ... times the number of cells
in each direction compared with a reference fine grid – are used in multigrid methods. A
moderately accurate solution is obtained quickly on the coarsest grid, where the number of cells
is small and changes propagate rapidly across the domain. This is then interpolated and iterated
to solution on successively finer grids.

If two levels of grid are used then Richardson extrapolation may be used both to estimate the
error and refine the solution. If the basic discretisation is known to be of order 𝑛 and the exact
(but unknown) solution of some property ϕ is denoted ϕ∗ , then the error ϕ − ϕ∗ may be taken
as proportional to Δ𝑛 , where Δ is the mesh spacing. For solutions ϕΔ and ϕ2Δ , respectively, on
two grids with mesh spacing Δ and 2Δ,
ϕΔ − ϕ ∗ = 𝐶Δ𝑛
ϕ2Δ − ϕ ∗ = 𝐶(2Δ)𝑛
These two equations for two unknowns, ϕ∗ and 𝐶, can be solved to get a better estimate of the
exact solution:
2𝑛 ϕΔ − ϕ2Δ
ϕ ∗= (16)
2𝑛 − 1
and the error in the fine-grid solution:
ϕ2Δ − ϕΔ
𝐶Δ𝑛 = (17)
2𝑛 − 1

Classroom Example 4

A numerical scheme known to be second-order accurate is used to calculate a steady-state

solution on two regular cartesian meshes A and B, where the finer mesh A has half the grid
spacing of mesh B. The values of the solution ϕ at a particular point are found to be 0.74 using
mesh A and 0.78 using mesh B. Use Richardson extrapolation to:

(a) estimate an improved value of the solution at this point;

(b) estimate the error at this point using the mesh-A solution.

CFD 9 – 11 David Apsley

9.3 Boundary Conditions

INLET
• Velocity inlet
The values of transported variables are specified on the boundary, either by a
predefined profile or by doing an initial fully-developed-flow calculation.
• Stagnation (or reservoir) boundary
Total pressure and total temperature (in compressible flow) or total head (in
incompressible flow) fixed; this is a common inlet condition in compressible flow.

OUTLET
• Standard outlet
Zero normal gradient (𝜕ϕ/𝜕𝑛 = 0) for all variables.
• Pressure boundary
As for standard outlet, except fixed value of pressure; this is the usual outlet
condition in compressible flow if the exit flow is subsonic (Ma < 1) and in free-
surface flow if the exit flow is subcritical (Fr < 1).
• Radiation
Prevent reflection of wave-like motions at outflow boundaries by solving a
simplified first-order wave equation ∂ϕ/ ∂𝑡 + 𝑐 ∂ϕ/ ∂𝑛 = 0 with wave velocity 𝑐.

WALL
• Non-slip wall
The default case for solid boundaries (zero velocity relative to the wall; wall stress
computed by viscous law or wall functions).
• Slip wall
Only the velocity component normal to the wall vanishes; used if it is not necessary
to resolve a thin boundary layer on an unimportant wall boundary.

SYMMETRY PLANE
• 𝜕ϕ/𝜕𝑛 = 0, except for the velocity component normal to the boundary, which is zero.
• Used where there is a geometric plane of symmetry.
• Used also as a far-field boundary condition, because it ensures that there is no flow
through, nor viscous stresses on, the boundary.

PERIODIC
• Used in repeating flow; e.g. rotating machinery, regular arrays.

FREE SURFACE
• Pressure fixed (dynamic boundary condition);
• No net mass flux through the surface (kinematic boundary condition).

CFD 9 – 12 David Apsley

9.4 Flow Visualisation

CFD has a reputation for producing colourful output and, whilst some of it is promotional, the
ability to display results effectively may be an invaluable design tool.

9.4.1 Available Packages

Visualisation tools are often packaged with commercial CFD products. However, many
excellent stand-alone applications or libraries are available, including the following.
• TECPLOT ([Link]
• AVS ([Link]
• EnSight ([Link]
• ParaView ([Link]
• VisIt ([Link]
• Dislin ([Link]
• Gnuplot ([Link]

9.4.2 Types of Plot

x-y plots

Simple two-dimensional graphs can be drawn by

hand or by many plotting packages. They are the
most precise and quantitative way to present
numerical data and, since laboratory data is often
gathered by straight-line traverses, they are a
common way of comparing experimental and
computational data. Logarithmic scales allow the
identification of important effects occurring at very
small scales, particularly near solid boundaries.
Line graphs are widely used for profiles of flow
variables and for plots of surface quantities such as
pressure or skin-friction.

One way of visualising flow development

is to use several successive profiles.

CFD 9 – 13 David Apsley

Line or Shaded Contour Plots

On a surface, a contour line (isoline) is a

curve along which some property is constant.
(The corresponding entity in 3-d space is an
isosurface.)

Any field variable may be contoured. In

contrast to line graphs, contour plots give a
global view of the flow field, but are less
useful for reading precise numerical values.
Detail occurring in small regions is often
obscured.

If contour intervals are equal then clustering

of lines indicates rapid changes in flow
quantities (e.g. shocks and discontinuities).

Colour is an excellent medium for conveying

information and good for on-screen and
presentational analysis of data. Simple
packages flood the region between isolines
with a fixed colour for that interval. More
advanced packages allow a pixel-by-pixel
gradation of colour between isolines.

Vector Plots

Vector plots display vector quantities (usually velocity, occasionally stress) with arrows whose
orientation indicates direction and whose size (and/or colour) indicates magnitude. They are a
popular and informative means of illustrating the flow field in two dimensions, although if grid
densities are high then either
interpolation to a uniform grid or a
reduced set of output positions is
necessary to prevent large localised
numbers of arrows obscuring plots.
It can be difficult to select a scale for
arrow lengths when large velocity
differences are present, especially in
recirculation zones where the mean
flow speed is low. In three
dimensions, vector plots can be
deceptive because of the angle from
which they are viewed.

CFD 9 – 14 David Apsley

Streamline Plots

Streamlines are everywhere parallel to the

local mean velocity vector. In 3D they must be
obtained by integration:
dx
=u
d𝑡
In advanced software they can be coloured –
usually by velocity magnitude, but sometimes
by other scalars such as temperature.
2
In 2-d incompressible flow a more accurate method is to contour the
streamfunction ψ. This function is defined by fixing ψ arbitrarily at
volume flux
one point and then setting the change in ψ between two points as the
volume flow rate (per unit depth) across any curve joining them:
2
ψ2 − ψ1 = 𝑣𝑜𝑙𝑢𝑚𝑒 𝑓𝑙𝑢𝑥 = ∫ u • n d𝑠
1 1

(The sense used here is clockwise about the start point, although the opposite sign convention
is equally valid). ψ is well-defined in incompressible flow because, by continuity, the flow rate
across any curve connecting two points must be the same. Contours of ψ are streamlines
because a curve of constant ψ is, by definition, one across which there is no flow.

For a path consisting of small increments d𝑥 and d𝑦 parallel to the axes, 2

dψ = 𝑢 d𝑦 − 𝑣 d𝑥 (18)
dy u
Conversely, the velocity components are related to ψ by
∂ψ dx
∂ψ 1
𝑢= , 𝑣=− (19) -v
∂𝑦 ∂𝑥

Computationally, ψ is stored at cell vertices. This is convenient 3

because the flow rates across cell edges are already stored as part of
the calculation. In the cartesian cell shown: ue
ψ2 = ψ1 + 𝑞𝑠 = ψ1 − 𝑣𝑠 Δ𝑥
1 2
ψ3 = ψ2 + 𝑞𝑒 = ψ2 + 𝑢𝑒 Δ𝑦 -vs

If isolines are at equal ψ increments, then

clustering of lines corresponds to high speed,
whilst regions where streamlines are further apart
signify low velocities. As with vector plots, this
makes it difficult to visualise the flow pattern in
low-velocity regions such as recirculation zones
and a smaller increment in ψ is needed here,

CFD 9 – 15 David Apsley

Classroom Example 5

(a) Two adjacent cells in a 2-dimensional cartesian mesh are shown below, along with the
cell dimensions and some of the velocity vectors (in m s–1) normal to cell faces. The
value of the stream function ψ at the bottom left corner is ψ𝐴 = 0. Find the value of the
stream function at the other vertices B – F. (You may use either sign convention).
D E 2 F

2 3 1 0.1 m
5
A B C

0.3 m 0.2 m

(b) Sketch the pattern of streamlines across the two cells in part (a).

CFD 9 – 16 David Apsley

Mesh Plots

The computational mesh is usually visualised by plotting the edges of control volumes. It can
be difficult to visualise fully-unstructured 3-d meshes, and usually only surface meshes or mesh
projections onto a plane section are portrayed.

Composite Plots

To maximise information it is common to combine plots of the above types, emphasising the
behaviour of several important quantities in a single picture. (In STAR-CCM+ parlance:
multiple “displayers” in a single “scene”.)

Complicated 3-d plots are often enhanced by the use of perspective, together with lighting and
other special effects such as translucency.

CFD 9 – 17 David Apsley

Appendix: Summary of Vector Calculus (Non-Examinable)

The operator
∂ ∂ ∂
∇≡ ( , , )
∂𝑥 ∂𝑦 ∂𝑧
(called del or nabla) is both a vector and a differential operator.

It is used to define:
∂ϕ ∂ϕ ∂ϕ
gradient: 𝑔𝑟𝑎𝑑 ϕ ≡ ∇ϕ ≡( , , ) acting on scalar field 
∂𝑥 ∂𝑦 ∂𝑧

∂𝑓𝑥 ∂𝑓𝑦 ∂𝑓𝑧

divergence: 𝑑𝑖𝑣 f ≡ ∇ • f ≡ + + acting on vector field f = (𝑓𝑥 , 𝑓𝑦 , 𝑓𝑧 )
∂𝑥 ∂𝑦 ∂𝑧

𝐢 𝐣 𝐤
∂ ∂ ∂
curl: 𝑐𝑢𝑟𝑙 f ≡ ∇ ∧ f ≡ |∂𝑥 ∂𝑦 ∂𝑧
| acting on vector field f = (𝑓𝑥 , 𝑓𝑦 , 𝑓𝑧 )
𝑓𝑥 𝑓𝑦 𝑓𝑧

Two important (and related) things about 𝑔𝑟𝑎𝑑 𝜙 are worth noting:
• ∇ϕ has direction that of the fastest rate of increase of increase of ϕ;
• ∇ϕ is perpendicular to the surfaces ϕ = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡

2
∂2 ϕ ∂2 ϕ ∂2 ϕ
div(grad ϕ) ≡ ∇ • ∇ϕ ≡∇ ϕ ≡ 2+ 2+ 2 is called the Laplacian.
∂𝑥 ∂𝑦 ∂𝑧

Integral Theorems dA

Gauss’s Divergence Theorem

For arbitrary closed volume 𝑉, with bounding surface 𝜕𝑉:
V
∫ ∇ • f d𝑉 = ∮ f • dA V
𝑉 ∂𝑉

This converts between a volume integral and a surface integral.

Stokes’s Theorem
For arbitrary open surface 𝐴, with bounding curve 𝜕𝐴:
A ds
∫ ∇ ∧ f • dA = ∮ f • ds
𝐴 ∂𝐴
A
This converts between a surface integral and a line integral.

CFD 9 – 18 David Apsley

Examples

Q1.
(a) Explain what is meant, in the context of computational meshes, by:
(i) structured;
(ii) multi-block structured;
(iii) unstructured;
(iv) chimera meshes.

(b) Define the following terms when applied to structured meshes:

(i) cartesian;
(ii) curvilinear;
(iii) orthogonal.

Q2.
(a) For wind-loading calculations a CFD calculation of airflow
is to be performed about the building complex shown.
Sketch a suitable computational domain, indicating the
specific types of boundary condition that are applied at each
boundary of the fluid domain. For each boundary type
summarise the mathematical conditions that are applied to
each flow variable (velocity, pressure, turbulent scalar).

(b) Define the drag coefficient for an object and explain how it would be calculated from
the raw data obtained in a CFD simulation.

r4
Q3. r3
Show that the vector area of a (possibly non-planar) quadrilateral is half
the cross product of its diagonals; i.e.
1 1 A
A= d13 ∧ d24 = (r3 − r1 ) ∧ (r4 − r2 )
2 2
r2
r1

Q4. (*** Advanced / optional ***)

Use Gauss’s Divergence Theorem to derive the following formulae for the volume of a cell
and the cell-averaged derivative of a scalar field:
1
(a) 𝑉= ∮ r • dA
3 ∂𝑉

∂ϕ 1
(b) ( ) = ∮ ϕ d𝐴𝑥
∂𝑥 𝑎𝑣 𝑉 ∂𝑉
where r is a position vector and dA a small element of (outward-directed) area on the closed
surface 𝜕𝑉.

CFD 9 – 19 David Apsley

Q5.
The figure shows part of a non-cartesian 2-d mesh. A single N

quadrilateral cell is highlighted and the coordinates of the corners W (0,4)

marked. The values of a variable ϕ at the cell-centre nodes (3,3)
(labelled geographically) are: (-2,2) P

ϕ𝑃 = 2, ϕ𝐸 = 5, ϕ𝑊 = 0, ϕ𝑁 = 3, ϕ𝑆 = 1. E
(2,0)

S
Find:
(a) the area of the cell;
(b) the cell-averaged derivatives (𝜕ϕ/𝜕𝑥)𝑎𝑣 and (𝜕ϕ/ ∂𝑦)𝑎𝑣 , assuming that cell-face
values are the average of those at the nodes either side of that cell face.

Q6.
One quadrilateral face of a hexahedral cell in a finite-volume mesh has vertices at
(2, 0, 1), (2, 2, – 1), (0, 3, 1), (– 1, 0, 2)
(a) Find the vector area of this face (in either direction).
(b) Determine whether or not the vertices are coplanar.

Q7.
(a) The vertices of a particular tetrahedral cell in a finite-volume mesh are
(0,0,0), (4,0,0), (1,4,0), (1,2,4)
Find:
(i) the outward vector area of each face;
(ii) the volume of the cell.

(b) For this cell a scalar ϕ has value 6 on the largest face, 2 on the smallest face and 3 on
the other two faces. Find the cell-averaged derivatives (𝜕ϕ/𝜕𝑥, 𝜕ϕ/𝜕𝑦, 𝜕ϕ/𝜕𝑧)𝑎𝑣 .

(c) Determine whether the point (1,2,1) lies inside, outside or on the boundary of the cell.

CFD 9 – 20 David Apsley

Q8. (Exam 2015)
A tetrahedral cell has vertices A(2, 2, 0), B(1, 5, 3), C(4, –1, 3), D(–1, –1, 4), as shown.
z
D(-1,-1,4)

C(4,-1,3)
B(1,5,3)

A(2,2,0)
x y

(a) Find the four outward face-area vectors and the volume of the cell.

(b) The average pressures on the four faces of the cell (indicated by vertices) are:
𝑃𝐵𝐶𝐷 = 1, 𝑃𝐴𝐶𝐷 = 3, 𝑃𝐴𝐵𝐷 = 3, 𝑃𝐴𝐵𝐶 = 4
Find the net pressure force vector on the cell.

(c) The average velocity vectors on three of the faces are

4 6 4
u𝐴𝐶𝐷 = (2), u𝐴𝐵𝐷 = ( 0), u𝐴𝐵𝐶 = (−1)
0 −1 1
Assuming incompressible flow, find the outward volume flux on all cell faces.

(d) The average velocity vector on the last face is

5
u𝐵𝐶𝐷 = (0 )
𝑤
Find the unknown velocity component 𝑤.

Data.
The volume 𝑉 of a polyhedral cell with plane faces is given by
1
𝑉= ∑ r𝑓 • A𝑓
3
𝑓

where r𝑓 is any position vector on cell face 𝑓 and A𝑓 is the vector area of cell face 𝑓.

CFD 9 – 21 David Apsley

Q9.
The figure below shows a quadrilateral cell, together with the coordinates of its vertices and
the velocity components on each face. If the value of the stream function ψ at the bottom left
corner is 0, find:
(a) the values of ψ at the other vertices;
(b) the unknown velocity component 𝑣𝑛 .
(3,4)
n face velocity
(0,3)
u v
w 3 –3
cell e
w s 0 2
e 3 1
y,v n 1 vn
(0,0) s (2,0)

x,u

Q10.
The figure right shows a region of a 2-d structured mesh. The
coordinates of the vertices of one cell are marked, together
12
with the pressure at nearby nodes.
(5,6)
(a) Assuming that pressure on an edge is the average of (1,5)
that at nodes on either side of the edge, find the (𝑥, 𝑦) 4 10 16
components of the pressure force (per unit depth) on
each edge of the shaded cell and hence the net pressure
(1,2) (6,1)
force vector on the cell.
14
(b) Find the area of the cell and the cell-averaged pressure
gradients 𝜕𝑝/𝜕𝑥 and 𝜕𝑝/𝜕𝑦.

(c) The figure right shows the velocity components on each (8,4) A
edge of the same cell. Find the outward volume flow rate
B
(per unit depth) from each edge and the missing velocity
component 𝑢. (u,2)
(7,2)
(d) If the value of the streamfunction ψ at vertex A is 3,
calculate the streamfunction at the other vertices. (You C
may use either sign convention.) (9,0) D

CFD 9 – 22 David Apsley

Q11. (Exam 2017)
The figure shows a single quadrilateral cell in a 2-d finite-volume mesh. Coordinates of the
vertices are given in the figure, whilst edge-averaged values of velocity (𝑢, 𝑣) and pressure 𝑝
(in consistent units) are given in the adjacent table. The flow is incompressible.
(-1,4)
n
(5,3) edge 𝑝 𝑢 𝑣
e 1 9 2
w n 5 7 –5
e w 19 1 –4
s 21 2 ?
(0,0)
s
(3,-1)
Find (in consistent units):
(a) the area of the cell;
(b) the net pressure force (per unit depth) on the cell;
(c) the outward volume flow rate (per unit depth) for all faces;
(d) the missing velocity component 𝑣𝑠 ;
(e) the value of the stream function ψ at all vertices, assuming it to be zero at the vertex
(0,0). (You may use either sign convention.)

Q12. (Exam 2018)

The figure shows a single tetrahedral cell in a finite-volume mesh. Coordinates of vertices B
and C are given in the figure, and vertex O is at the coordinate origin. Outward cell-face area
vectors are given in the accompanying table (as column vectors), along with face-averaged
velocity and pressure in consistent units.
Face Outward Velocity Pressure
z vector area
C(0,0,2) −6 2
OBC ( 2) (−12) 5
0 2

O −1 2
A OAC (−2) ( 2) 140
0 3
x y 3 4
B(2,6,0) OAB (−1) (−10) 40
−7 1
Find: 4
(a) the outward cell-face area vector for face ABC; ABC ? ( ?) 32
(b) the volume of the cell; 3
(c) the outward volume flow rate through each face;
(d) the missing 𝑦-velocity component on face ABC;
(e) the average pressure gradient;
(f) the coordinates of point A.

CFD 9 – 23 David Apsley

Q13. (Exam 2019)
The figure shows a single quadrilateral cell in a finite-volume mesh. Coordinates of its vertices
are given in the figure. Cell-face-averaged velocity (𝑢, 𝑣) and pressure, 𝑝, are given in the
table, in consistent units. The flow is incompressible.
(2,4)
Face 𝑢 𝑣 𝑝
n
e 5 –5 9
e
(-1,2) n 1 –6 6
w –1 –2 10
w (3,1)
s s 3 ? –9
(0,0)
(a) Find:
(i) the area of the cell;
(ii) the net pressure force on the cell (per unit depth);
(iii) the missing v-velocity component on face 𝑠.

(b) Explain:
(i) the main principles of the finite-volume method for CFD;
(ii) the SIMPLE pressure-correction method;
(iii) the main features of flux-limited advection schemes;
(iv) what is meant by a RANS turbulence closures and, specifically, by an eddy-
viscosity turbulence model.

Q14. (Exam 2020 – part)

The figure shows a tetrahedral cell in an unstructured mesh. The coordinates of the vertices are
marked in situ and the cell-face values of a scalar ϕ are given below the figure. Find:
(i) the volume of the cell;
(ii) the cell-averaged derivatives 𝜕ϕ/𝜕𝑥, 𝜕ϕ/𝜕𝑦, 𝜕ϕ/𝜕𝑧.
z
A(1,1,2)

y
C(-1,1,0)

D(1,1,1)

B(2,-1,0) x

ϕ𝐵𝐶𝐷 = 2, ϕ𝐴𝐶𝐷 = 9, ϕ𝐴𝐵𝐷 = 12, ϕ𝐴𝐵𝐶 = 4

CFD 9 – 24 David Apsley

Q15. (*** Advanced / optional ***)
(a) Give a physical definition of the stream function ψ in a 2-d incompressible flow, and
show how it is related to the velocity components.

(b) Show that if the flow is irrotational then ψ satisfies Laplace’s equation. If the flow is
not irrotational how is 2 ψ related to the vorticity?

The figure below and accompanying table show a 2-d quadrilateral cell, with coordinates at the
vertices and velocity components on the cell edges. The flow is incompressible.
(4,4) face velocity
n
(0,3) 𝑢 𝑣
e w 3 2
w s 0 2
y
(1,1) s (3,1) e 3 –1
x n 1 vn

(c) If the value of the stream function at the bottom left (sw) corner is 0, calculate the stream
function at the other vertices.

(d) Find the 𝑦-velocity component, 𝑣𝑛 , whose value is not given in the table.

(e) By calculating ∮ u • ds for the quadrilateral cell estimate the cell vorticity.

Q16. (*** Advanced / optional ***)

A quadrilateral cell in a 2-d finite-volume mesh is shown in the figure below. The coordinates
of the vertices are marked in the figure; the velocities on the edges are given in the table.

(a) Find the area 𝐴 of the quadrilateral.

(b) Find the cell-averaged derivatives (∂𝑢/ ∂𝑥)𝑎𝑣 , (∂𝑢/ ∂𝑦)𝑎𝑣 , (∂𝑣/ ∂𝑥)𝑎𝑣 , (∂𝑣/ ∂𝑦)𝑎𝑣 .

(c) By calculating the line integral from the velocity and geometric data, confirm that both
the cell-averaged velocity derivatives and the discrete form of Stokes’s Theorem,
1
(ω𝑧 )𝑎𝑣 = ∮ u • ds
𝐴 ∂𝐴
give the same cell-averaged vorticity, (ω𝑧 )𝑎𝑣 . (This is, in fact, true in general.)
(-1,3)
edge 𝑢 𝑣
n w 9 –6
w
(2.5,2) s 1 0
e 10 12
e
n 13 6
y (-2,0)
s
x (0,-1)

CFD 9 – 25 David Apsley