CYCLONE

Fluid Dynamics

Dolfyn Developers Guide

DRAFT

CFD-10xxxx
20th May 2010

Cyclone Fluid Dynamics B.V.

Sweelincklaan 4 Tel.: +31-40-22 30 491 Web: www.cyclone.nl
NL–5583 XM Waalre Fax.: +31-40-22 30 490 email: info@cyclone.nl
 CYCLONE
Fluid Dynamics

Dolfyn Developers Guide

DRAFT

Cyclone Fluid Dynamics BV

Author:
H.W. Krüs

CFD-10xxxx
Cyclone Fluid Dynamics B.V.
20th May 2010

Copyright
c Cyclone Fluid Dynamics B.V., 2010. All rights reserved.

Cyclone Fluid Dynamics B.V.

Sweelincklaan 4 Tel.: +31-40-22 30 491 Web: www.cyclone.nl
NL–5583 XM Waalre Fax.: +31-40-22 30 490 email: info@cyclone.nl
 CYCLONE
Fluid Dynamics

ii
 CYCLONE
Fluid Dynamics

Preface

This is just a collection of notes, thoughts and ideas about the dolfyn code. Do
not expect a decent structure nor an in depth methodology manual. Consult for
the latter the various appropriate text books.

Thanks to Yousef Saad of University of Minnesota, Department of Computer

Science and Engineering for Sparsekit2 (www-users.cs.umn.edu/∼saad/) and
the people behind Lapack (www.netlib.org/lapack/) and countless others like
those involved in the Hypre and VisIt packages at Lawrence Livermore
(www.llnl.gov/CASC/linear_solvers/ respectively www.llnl.gov/visit/)

A special thanks goes to Franco Magagnato of the Department of Fluid Machin-

ery, University of Karlsruhe, for sharing his Sparc code years ago (Structured
PArallel Research Code). Sparc showed and proved the concept of f95 and
MPICH on a heterogeneous network consisting of different machines and op-
erating systems. Within a day the code was up and running over a network with
Alpha/Unix and PC/Linux workstations. A quality statement for the three pieces
of software involved (Sparc, f95 and MPICH).

Thanks also to those involved in for example TEX, LATEX and pdfTEX, xfig,
NEdit, the Gimp, Perl, Gnu and Linux. Intel made it’s contribution by releasing
their Intel Fortran95 compiler (ifc/ifort) for free to the Linux community. And
of course, Andy Vaught for his amazing work on g95 (www.g95.org). Thanks
to IBM for Data Explorer or OpenDX as it is now called (www.opendx.org). An
impressive tool. Also CiteSeer at Penn State was, and is, very useful.

Simplicity and accuracy are the keywords. This means for example that there are
no 2D and axisymmetric modes of the code. Also readability is more important
than speed (rely instead on the quality of today’s compilers and on the array
syntax of f95).

Dolfyn is based on the ideas, algorithms, and procedures presented in the book
by Joel Ferziger & Milovan Períc ‘Computational Methods for Fluid Dynamics’
[13]. This book is also acts as Theoretical Manual.

Dolfyn contains contributions by Thomas Geenen, Johan Jacobs, Shibo Kuang,

Max Staufer, Bouke Tuinstra.

Enjoy!

CFD-10xxxx iii
 CYCLONE
Fluid Dynamics

iv
 CYCLONE
Fluid Dynamics

Contents

Preface iii

Contents v

List of Figures vii

1 Introduction 1
1.1 Fortran 9x 3
1.2 Further reading 4

2 A transported scalar 9
2.1 Conservation of momentum 12
2.2 Diffusion term and non-orthogonal meshes 14
2.3 Calculation of wall shear stresses 14
2.3.1 Aerodynamically smooth walls 15
2.3.2 Calculation of turbulent wall shear stresses 16
2.3.3 Rough walls 17
2.4 Residuals 18

3 Convective term and differencing schemes 23

3.1 Standard upwind differencing, UD 27
3.2 Second order central differencing, CD 27
3.2.1 Linear interpolation, CD1 27
3.2.2 Gradient based averaging, CD2 27
3.2.3 Simple averaging, CD3 28
3.3 Blended differencing 28
3.3.1 Constant blending 28
3.3.2 Variable blending, hybrid scheme 28
3.3.3 Variable blending, power law scheme 28
3.4 Second order linear upwind differencing, LUD 29
3.5 Third order linear upwind differencing, QUICK 29
3.6 Second order switching schemes 30
3.6.1 MinMod differencing, MinMod 30
3.6.2 Gamma differencing, GAM 30
3.6.3 Third order monotone upstream, MUSCL 30
3.7 VOF scheme(s) 30
3.8 Discusion 30

4 Calculation of gradients 35
4.1 Gauss 35
4.2 Least squares 35
4.3 Gradient limiters 37

CFD-10xxxx v
 CYCLONE
Fluid Dynamics

4.3.1 Barth and Jespersen 37

4.3.2 Van Albada 38
4.3.3 Venkatakrishnan 41
4.3.4 Polynomial function 42
4.3.5 Uniform regions 43
4.4 Discussion 44

5 Heat transfer 47
5.1 Laminar heat transfer 49
5.2 Turbulent heat transfer 49

6 Lagrangian particles 53
6.1 Basic equations 53
6.1.1 Drag force 53
6.1.2 Virtual mass force 54
6.1.3 Body forces 54
6.1.4 Adding all forces 54
6.2 Integrating particles in time 57
6.3 Particle within a cell 58
6.4 Intersecting a cell face 59
6.5 Interpolating velocities 60
6.6 Examples 62

Appendices 69

A Geometry input files 71

B Mathematical toolbox 73
B.1 Gauss’ theorem 73
B.2 Sum of two vectors 73
B.3 Multiplying a vector 74
B.4 Length of a vector and vector products 74

Bibliography 79

Index 83

vi
 CYCLONE
Fluid Dynamics

List of Figures

1.1 Dolfyn data and result files 2

1.2 Example Fortran95 code 3

2.1 An arbitrary control volume 10

2.2 Correction 11
2.3 Over relaxed non-orthogonal treatment 15
2.4 Implementation of N.Scholz’s universal velocity profile [39] 18

3.1 Definition of upwind, central and downwind nodes 24

3.2 Convection Boundness Criterion with UDS, CDS and LUDS 25
3.3 31
3.4 Smart, MinMod and Gamma-schemes 31

4.1 Average as a result of fraction x 39

4.2 Influence of the bias ε2 on the van Albada function 40
4.3 Venkatarishnan limiter vs Barth Jespersen and van Albada. 42
4.4 Transition function s(x) 44

5.1 Heat transfer near a wall 48

6.1 Characteristic response of a particle to forces 55

6.2 Example of particles passing through cells 58
6.3 Example of particle inside and outside of a cell (red) 59
6.4 Particle intersects face 59
6.5 Free falling particles, influence of particle Courant number 62
6.6 Free falling particles, details showing steps within cells 63
6.7 Particles in a stagnation flow 64
6.8 Particles in a stagnation flow, influence of particle Courant num-
ber 65
6.9 Small particles/streamlines in a stagnation flow 66
6.10 Small particles/streamlines in an 45 angle test 67

A.1 Cell definitions 71

B.1 Sum of two vectors 74

B.2 Inner product of two vectors (hashed area of outer product) 75
B.3 Vector components near a wall 77

CFD-10xxxx vii
 CYCLONE
Fluid Dynamics

viii
 CYCLONE
Fluid Dynamics

The purpose of computing is insight, not numbers.

R. Hamming, 1962

CFD-10xxxx ix
 CYCLONE
Fluid Dynamics

x
 CYCLONE
Fluid Dynamics

1 Introduction

Heat and mass transfer, transport phenomena, and fluid dynamics play an im-
portant role in various products, machines and in our daily life.

Examples are:

• Automotive industry

– External flow around vehicles

– Fouling
– Thermal passenger comfort

• Heating and cooling

– Heat exchangers
– Cooling of electronics

• Flows in pipes and tubes

– Industrial installations
– Air conditioning

• Wind in the built environment

– Wind forces
– Wind nuisance for pedestrians and cyclists

• Process technology

dolfyn is a face based implicit Finite Volume Method code, employing primitive
variables on 3D unstructured polyhedral meshes targeted towards these indus-
trial type of problems.

dolfyn is accompanied by a preprocessor. The preprocessor writes a geometry

file in a format suitable for dolfyn. The input is a set of three files which describe
the cells, vertices en the boundaries (see Appendix A). Dolfyn reads the geome-
try file and a separate input file (format is described in a separate document). In
this input file the user sets numerical and model data, boundary conditions etc.
The input file can be edited with any simple editor. The file structure is shown
in Figure 1.1.

The code uses a ‘segregated solver’ which means that the transport equations
are solved sequentially. Because the coupled non-linear equations have been

CFD-10xxxx 1
 CYCLONE 1 Introduction
Fluid Dynamics

Hypermesh
ICEM/CFD others ...

geometry data (*.cel *.vrt *.bnd)

standard editor
grid preprocessor

control data deck (*.din)

grid data (*.geo)
DOLFYN
result translator
restart data file (*.rst)

output files (*.txt)

restart file

output
opendx data file (*.odx)

OpenDX GMV others ...

Figure 1.1: Dolfyn data and result files

linearised implicitly, several iterations are needed in order to get a converged

solution. Using the currently known fluid properties and mass fluxes the three
momentum transport equations are solved. At this stage the continuity equation
may not be satisfied and a pressure correction is set up. This correction is solved
to obtain the correct pressure and velocity fields and the mass fluxes at the faces
satisfy the continuity equation. After this stage the transport equations for the
turbulence models, energy or species is solved. Finally the fluid properties are
adjusted and the process is repeated until convergence.

The grid or mesh used by dolfyn is based on ‘face based unstructured polyhedral
cells’. Unstructured does not mean tetrahedral cells only (in 3D, or triangles in
2D), but refers to the way the topology is implemented. In using polyhedra any
type of mesh can be used, even ‘structured hexahedra’s’. It only depends on the
kind of preprocessor, or ‘mesher’, or ‘grid generator’ one is using. The same
applies to the postprocessor to visualise the final results; a postprocessor which
only can handle hexahedral cells is of no use when one employs tetrahedral cells.
On the other hand results of hexahedral cells can be converted to a tetrahedral
structure and displayed (this technique is partly used for OpenDX for example).
However, ‘real polyhedra’ still pose a problem (except, maybe, for GMV, but
GMV is not Open Source). Thus dolfyn is never the limiting code in the process
and it all depends on the kind of, favourite, pre- and postprocessors at hand.

2
1.1 Fortran 9x CYCLONE
Fluid Dynamics

1.1 Fortran 9x

Dolfyn uses Fortran? But isn’t that hopelessly old fashioned?

Lots of people who have not encountered the latest Fortran compilers react in
this way. Fortran 77 is indeed a bit ‘old fashioned’ and had some serious limita-
tions (however not it’s speed).

The language was developed way back in 1954 at IBM and has been updated
since. Well known versions are Fortran IV (or Fortran 66), Fortran 77 (the well
known old standard), Fortran 90 and soon followed by Fortran 95. In 2003
further enhancements were made, notably enhancements of derived types, asyn-
chronous transfer, stream access and procedure pointers and, very important,
much improved interoperability with the C programming language. If number
crunching is what you are up to then Fortran (95+) still is the best horse for your
course. With g95 and f03gl you can even do OpenGL!

Modern Fortran is also highly portable. It is relatively easy to write an entirely

portable program in Fortran, even without recourse to a preprocessor. Nor does
one needs tools like the configure and build system (‘automake’, ‘autoconf’,
and ‘configure’ etcetera). Hunting down ‘include files’ all over the place and
file name case sensivity are absent (something one has to take care of when
developing a multiplatform solution). And finally, not to mention Fortran’s vec-
torisation and parallelisation capalilities are excellent.

If the before mentioned arguments are not enough. Well, the most important
argument for Fortran is it’s readability. Even people who hardly know the lan-
guage can read and understand the code (see the example in Figure 1.2). “An-
other important point is that the time required to learn Fortran 90 is much less
than the time to learn either C or C++. ... But the issue is what freshmen should
learn as their first language and for that I would recommend Fortran 90 hands
down.”1 Readability should not be underestimated since it is the core of the
peer-reviewed human eye parallel debugging device.

.
.
Xn = Face ( i )%n ! surface vector
c a l l n o r m a l i s e ( Xn ) ! normal v e c t o r w i t h l e n g t h 1

Up = Up − Uw ! v e l oc i t y difference vector

dp = d o t _ p r o d u c t ( Up , Xn ) ! d o t p r o d u c t o f t h e two v e c t o r s
Un = dp ∗ Xn ! normal v e l o c i t y component
Ut = Up − Un ! t a n g e n t i a l v e l o c i t y component
.
.

Figure 1.2: Example Fortran95 code

1
"Fortran 90 vs C++ - An educational perspective", by Dr. John K Prentice, Quetzal Com-
putational Associates, 1996. Note that Fortran 90 is addressed and advocated, imagine Dr.
Prentice’s enthusiasm about Fortran 95 or even Fortran 2003.

CFD-10xxxx 3
 CYCLONE 1 Introduction
Fluid Dynamics

The current July 2008 version of dolfyn only consists of about 25,000 lines of
Fortran code. About 60% of it, in the order of 15,000 lines of code, is actually
related to dolfyn itself. 16% of the lines only deal with reading the din control
file and setting the control parameters (easy of use always tends to increase the
‘interactive part’ of a code). Another 10% are subroutines for the various out-
put interfaces. Finally 18% accounts for the Fortran77 SparseKit2 and Lapack
subroutines.

1.2 Further reading

Of course the basics of dolfyn can be found in Ferziger and Períc’ ‘Computa-
tional Methods for Fluid Dynamics‘ [12][13]. Older classic textbooks are by
Roache [37] and Patankar [32], the influence of the latter is not to be under-
estimated. Today ‘Patankar’ is useful as an introduction but ‘proves’ that a
co-located code, like dolfyn, is not feasible. The breakthrough came with the
Rhie and Chow paper [35] in 1983. But ‘staggered grids’ today still have their
supporters like Wesseling [52][50].

Dolfyn relies heavily on the work and theses by people like Spalding, Gosman,
Issa, Jasak, Weller and many others at Imperial College and in Britain (examples
are [17][33][25]). Imperial College is at the roots of various leading commer-
cial CFD codes (from the seventies, early eighties onwards). Spalding started
in mid 70’s with Concentration Heat and Momentum Ltd (CHAM) in Wimble-
don and in 1981 (single block hexahedral structured) Parabolic Hyperbolic Or
Elliptic Numerical Integration Code Series, or PHOENICS, appeared. Fluent has
it’s (ancient) roots at Sheffield and saw it’s first launch in 1983 by Creare Inc2 .
Computational Dynamics was founded in 1987 in London by Gosman and Issa
and they launched Star-CD, from version two onwards hexahedral unstructured.
Down the road in Harwell/Oxford the UK Atomic Energy Authority (AEA) re-
leased the multi-block structured Flow3D, later known as CFX3 .

Another stronghold is the T-3 group at Los Alamos [27] were k-ε was born in
1967 [19][20][18] and the Volume-Of-Fluid, VOF, method in the eighties. Los
Alamos has a long tradition in open source codes (the Marker-and-Cell, MAC,
method by Harlow and Welch, the first successful code with primitive variables
thanks to a novel staggered scheme [14][21][49][5]. Followed by ‘SOLA’ [24],
‘SOLA-VOF’ by Hirt and Nichols [31], Arbitrary Lagrangian-Eulerian or ALE
[23] and other codes, and the impressive ‘Kiva’ series of combustion codes by
Hirt, Amsden and others [36][8][7][15][6][3][4].

Yet another world of CFD applications originates naturally from the aerospace
world. Here too the vast majority of the codes were multi-block structured
codes see for example the interesting introductions as Anderson’s ‘Computa-
2
Fluent acquired Fluid Dynamics International (FDI) and it’s Finite Element Method based
package Fidap. Fluent was bought by Ansys in 2006.
3
Later the privatised part of AEA, called AEA Technology, aquired Canadian Advanced Sci-
entific Computing (ASC) and it’s coupled solver TASCflow in 1997. CFX was bought by
Ansys in 2003.

4
1.2 Further reading CYCLONE
Fluid Dynamics

tional Fluid Dynamics: The Basics with Applications’ [9] or ‘Computational

Fluid Mechanics and Heat Transfer’ by Tannehill, Anderson and Pletcher [42].
Tetrahedral unstructured codes were rare and some of the techniques used for
these kind of codes now start to appear in the general purpose codes.

More on boundary layers can be found in Schlichting [39], and on turbulence

start with Tennekes and Lumley [43] and then Hinze [22]. Turbulence modelling
for CFD is dealt with by Wilcox [51]. Another book with much attention to
modelling is Pope’s ‘Turbulent flows’ [34].

On general working with CFD one can use Versteeg and Malalasekera’s ‘An In-
troduction to Computational Fluid Dynamics’ (use the latest 2007 edition [48]
instead of the 1995 edition [47]). Shaw’s ‘Using Computational Fluid Dynam-
ics’ [40] is even older (1992) and is today not so useful anymore.

Do not start with the numerics without proper background in fluid dynamics.
Consult your local popular textbooks. And keep Van Dyke’s ‘Album of fluid
motion’ [45] and other books like ‘Visualized Flow’ by Nakayama [30] at hand.

CFD-10xxxx 5
 CYCLONE 1 Introduction
Fluid Dynamics

6
1.2 Further reading CYCLONE
Fluid Dynamics

Just as with buses, so with Open Source compilers

for Fortran95: you wait ages for one, then two arrive at once.

Clive Page

CFD-10xxxx 7
 CYCLONE 1 Introduction
Fluid Dynamics

8
 CYCLONE
Fluid Dynamics

2 A transported scalar

The basic ideas of how a transport equation can be solved will be illustrated in
this chapter.

The basis is formed by the assumption that a variable varies within a control
volume according to
φ(~x) = φP + (~x − ~xP ) · (∇φ)P . (2.1)
With the centroid ~xP defined as:
Z
(~x − ~xP ) dV = 0. (2.2)
VP
Now when equation (2.1) is integrated over volume V results in
Z Z
φ(~x) dV = φP + (~x − ~xP ) · (∇φ)P dV
VP VP
Z Z
= φP dV + [ (~x − ~xP ) dV ] · (∇φ)P
VP V
| P {z }
=0
= φP V P (2.3)
This means that equation (2.2) is an essential geometrical proposition. Because
~xP is in the centre and the variable
R varies linearly; the value of φ at P is equal
to the mean value, and therefore V φ(~x)dV is equal to φP VP with second order
accuracy.

Imagine a known flow field (which might be at rest as well) in which the scalar
quantity φ like a temperature or a concentration is transported. For some control
volume one can write down the following conservation equation:
∂
(ρ φ) + div( ρ ~vr φ ) = div( Γφ grad φ ) + qφ . (2.4)
|∂t {z } | {z } | {z } |{z}
unsteady term convective term diffusive term source term

Where ~vr is the relative velocity (~v − ~vc ) between the fluid moving with ~v and
the sweeping coordinate velocity ~vc .

It simply states that the change in time of φ is the result of the balance between
the transport (convection) and diffusion in and out of the the control volume,
and, when present, the result of a source (or sink).

Equation (2.4) for an arbitrary moving volume V with surface S transforms into
Z Z Z Z
d ~= ~+
ρ φ dV + (ρ ~vr φ) · dS ( Γφ grad φ) · dS qφ dV (2.5)
dt V S S V

CFD-10xxxx 9
 CYCLONE 2 A transported scalar
Fluid Dynamics

y
N faces = 7










f=2
  
 

f=3
f=1
P
  

f=4
  

  

  

  

  

  

  

  

  

  

  

f=7
  

  

f=5
f=6
  

  

  

  

  
  

  

  

  

Figure 2.1: An arbitrary control volume

~ is the surface vector and the relative velocity is now related to the
where S
moving surface.

For steady state cases without sources (or sinks) this equation reduces to a simple
convection/diffusion equation:

div( ρ ~v φ ) − div( Γφ grad φ ) = 0 (2.6)

or using Gauss’ theorem

Z Z
ρφ~v · ~n dS − Γ grad φ · ~n dS = 0 (2.7)
S S

and can be casted into:

Nfaces
X Nfaces
X
Ffc − Ffd = 0
f =1 f =1

or
Nfaces
X
(Ffc − Ffd ) = 0
i=f

with Z
Ffc = ρφ~v · ~n dS ≈ ṁf φf (2.8)
S
and ṁf as the mass flux through face ‘f’:
Z
ṁf = ρ~v · ~n dS ≈ (ρ~v · ~n)f Sf (2.9)
Sf

As the face can be arbitrarily oriented all the velocity components can contribute
to the mass flux. Thus each velocity component is multiplied by the correspond-
ing surface vector component (see 8.6.1).

ṁf = ρf (S x ux + S y uy + S z uz )f (2.10)

Now assume that the fluid properties ρ and Γ are known (eg. from a previous
iteration). Then the value of φ and its gradient normal to the cell face is needed.

10
 CYCLONE
Fluid Dynamics












xn
x p x f’
  

  

  

  











N

∆2
  

  

P

∆1

xf
∆1
λ=
∆1 + ∆ 2
x

Figure 2.2: Correction

Recall that the normal of a cell face is positive directed outwards (in Dolfyn
defined by the Face(i)%cell1 and Face(i)%cell2, the first one is by definition
node ‘p’ and the second, if present, the ‘next’ node ‘n’).

The values at the faces are needed at the face centre. The second order linear
interpolation of two φ’s is on a equidistant grid:
φP + φN
φf =
2
and for non-equidistant grids:

φf = λf φN + (1 − λf )φP (2.11)

with λ defined as:

xf − xP
λf = (2.12)
xN − xP
The factor λ can be calculated and stored in advance. This is the ‘central differ-
encing scheme’ (CDS).

For non-smooth grids however the before mentioned accuracy will probably be
lost. The value is not evaluated at ~xf but at ~xf 0 which is defined by the ratio of
the original lengths (λ) and lies on the straight line connecting the two centres.
Part of the accuracy can be restored if one evaluates the gradient of φ at ~xf 0
(= (∇φ)f 0 ) and travels from ~xf 0 to ~xf . This is shown in Figure 2.2.

For equidistant grids, and some smooth expanding grids, the correction is not
necessary (~xf 0 − ~xf = 0). For highly distorted grids the procedure might be re-
peated several times (hence the variable Ngrad in subroutine GradientPhiGauss).

There is also a simpler first order method called ‘upwind differencing scheme’
(UDS), in which the upstream nodal value is used:

φP if (~v · ~n)f > 0
ṁf = (2.13)
φN if (~v · ~n)f < 0

The upwind differencing scheme is robust and never results in under- or over-
shoots, nor is it oscillatory. However, it is inaccurate especially when the flow is

CFD-10xxxx 11
 CYCLONE 2 A transported scalar
Fluid Dynamics

oblique trough the cell face. Thus it is possible to combine both scheme’s into a
blended or hybrid scheme:

Ff = (1 − β) Ffuds + β Ffcds (2.14)

or the using a deffered correction approach

 old
Ff = Ffuds + γ Ffcds + Ffuds (2.15)
|{z} | {z }
implicit explicit

a hybrid scheme in which the higher order scheme is evaluated explicitely as a

corrector step to the stable lower order implicit scheme.

Thus the flux can be approximated by:


c max(ṁf , 0) φP + min(ṁf , 0) φN for UDS
Ff = (2.16)
ṁf (1 − λf ) φP + ṁf λf φN for CDS

For the diffusive fluxes one needs the gradient vector at the cell face centre. The
CDS scheme is normally used:
Z
d
Ff = Γ grad φ · ~n dS ≈ (Γgrad φ · ~n)f Sf (2.17)
Sf

For simple grids one can simply use:

 
∂φ φN − φP
= (∇φ)f ace · ~xf ace ≈
∂n f ace |~xN − ~xP |

As both values at P and N show up; this relationship can be used simply in an
implicit form:

However this form is too simple for complex meshes. A correction term is
discussed in paragraph XXX.

So for one control volume there are a number of fluxes, each of them related
between the central node ‘P’ and the neighbour defined by the face.

surface area components sx , sy and sz ==> normal

2.1 Conservation of momentum

Consider the conservation equations for mass and momentum in the form of
equation (2.4) for an arbitrary moving volume V with surface S:
Z Z
d ~=0
ρ dV + (ρ ~vr ) · dS
dt V S

and Z Z Z Z
d ~=
ρui dV + (ρui~vr ) · dS T~ · dS
~+ ρbi dV
dt V S S V

12
2.1 Conservation of momentum CYCLONE
Fluid Dynamics

With bi the ‘body force’ in the direction of the Cartesian coordinate xi (for
example gravity, centrifugal and Coriolis forces, electromagnetic forces etc.).

The stress tensor T~ is the molecular rate of transport of momentum and incor-
porates the surface forces (pressure, normal and shear stresses, surface tension
etc.).

~ = ~ndS results
Introducing the outwards pointing unit normal ~n to surface dS
in: Z Z
d
ρ dV + (ρ ~vr ) · ~n dS = 0 (2.18)
dt V S

and
Z Z Z Z
d
ρui dV + (ρui~vr ) · ~n dS = t~i · ~n dS + ρbi dV (2.19)
dt V S S V

For isotropic Newtonian fluids the components τij , of the viscous stress tensor
are:  
∂ui ∂uj 2
τij = −p δij + µ + − δij µ div~v (2.20)
∂xj ∂xi 3
with the dynamic fluid viscosity µ and the Kronecker delta δij .

The stress tensor  

τxx τxy τxz
Π =  τyx τyy τyz 
τzx τzy τzz
with τxx = σx etc, and τxy = τyx , τxz = τzx and τyz = τzy only contains six
elements and is symmetric
 
σx τxy τxz
Π =  τxy σy τyz  (2.21)
τxz τyz σz

and the elements are

σx = −p − 32 µ div ~v + µ( ∂u ∂u
∂x + ∂x )
2 ∂v ∂v
σy = −p − 3 µ div ~v + µ( ∂y + ∂y ) (2.22)
2 ∂w ∂w
σz = −p − 3 µ div ~v + µ( ∂z + ∂z )

∂v
τxy = τyx = µ( ∂x + ∂u
∂y )
∂w
τyz = τzy = µ( ∂y + ∂v
∂z ) (2.23)
∂u ∂w
τzx = τxz = µ( ∂z + ∂x )
With p as the local thermodynamic pressure. Pleas note that for incompressible
flows div ~v = 0.

Because the momentum equation R (2.19) is a vector equation, the corresponding

∂ ∂ui
diffusive term is ∂xj (µ ∂xj ) or S µ grad ui · ~n dS (see FP 7.1.1).

One can extract the pressure p from equation (2.22)

CFD-10xxxx 13
 CYCLONE 2 A transported scalar
Fluid Dynamics

2.2 Diffusion term and non-orthogonal meshes

The discretisation of the diffusion term is:

Z Nfaces
X
∇ · (ρΓφ ∇φ) dV = ~ · (ρΓφ ∇φ)f
S
VP f =1
Nfaces
X
= ~ · (∇φ)f .
(ρΓφ )f S (2.24)
f =1

If the mesh is orthogonal then the diffusion term is simply

S ~ φN − φP .
~ · (∇φ)f = |S| (2.25)
~
|d|

Another method would be calculating the cell centered gradients

Nfaces
1 X ~
(∇φ)P = Sφf (2.26)
VP
f =1

and interpolating to the cell face using (2.11) again:

(∇φ)f = λf (∇φ)N + (1 − λf )(∇φ) P. (2.27)

However this approach cannot be used for non-orthogonal meshes [25].

~ · (∇φ)f is split into two parts

The product S

~ · (∇φ)f = ~s · (∇φ)f ) +
S ~t · (∇φ)f ) . (2.28)
| {z } | {z }
orthogonalpart non−orthogonal corr.

The two vectors ~s and ~k should give the surface again

~ = ~s + ~t.
S (2.29)

Where ~s simply is parallel to d~ (= ~xP − ~xN ) and equation (2.25) can still be
used.

In [25] the over-relaxed approach is chosen. The contribution of φP and φN is

increased when the non-orthogonality is increased.

2.3 Calculation of wall shear stresses

The viscous stresses near a wall are ([13] eq. (8.75)):

 
∂vn
τnn = 2µ =0, (2.30)
∂n wall

14
2.3 Calculation of wall shear stresses CYCLONE
Fluid Dynamics

d
N
P s

k
S

Figure 2.3: Over relaxed non-orthogonal treatment

and  
∂vt
τnt = µ . (2.31)
∂n wall
In which the coordinate t is in the direction of the velocity vector near the wall.
The third coordinate s is orthogonal to t and the normal n.

The surface integral of τnt results in a force:

Z
~
fwall = ~t τnt dS ≈ (~t τnt S)S (2.32)
S

in which f~wall points in the direction of ~t (with unit length) and magnitude τnt
multiplied by the surface area S.

2.3.1 Aerodynamically smooth walls

High Reynolds number flows normally require ‘wall functions’ to represent the
fast varying distributions of velocity, turbulence et cetera. Wall functions are
based upon the following assumptions:

• Changes are predominantly normal to the wall, leading to one-dimensional

behaviour (L. Prandtl’s great contribution).
• Pressure gradients are neglected and uniform shear stress in the boundary
layer.
• Shear stress and velocity vectors are aligned and do not twist in the layer.
• Turbulence energy prdoction is balanced by dissipation.
• Linear variation of the length scale (connected to Von Karman’s constant).

The standard formulation of ‘the law of the wall’ for the turbulent boundary
layer along a smooth wall reads:
1
u+ = ln(y + ) + Bsmooth (2.33)
κ

CFD-10xxxx 15
 CYCLONE 2 A transported scalar
Fluid Dynamics

with:

u+ = (u − uw )/uτ , the dimensionless velocity.

u, the velocity near the wall.
yw , wall velocity (in the plane of the wall).
y + , dimensionless distance from the wall defined by y + = yuτ /ν.
κ, Von Karman’s constant (κ = 0.42).
uτ , the shear velocity (defined by τw = ρu2τ ).

In the numerical world the following formulation for the complete boundary
layer is often used:
y+ for y + ≤ ym
+ ,

+
u = 1 + + + (2.34)
κ ln(Ey ) for y > ym .

with E as the wall roughness coefficient.

+ is defined by
The matching point ym
+ 1
ym − ln(Ey + ) = 0 , (2.35)
κ
which is a trancedental equation and cannot be solved analytically.

The value of E for an aerodynamically smooth wall can be derived by equating

(2.33) and (2.34).
1 1
ln y + + Bsmooth = ln E y +



κ κ
E = eκBglad , (2.36)
with Bsmooth = 5.5 and κ = 0.4 follows E = 9.0.

For the turbulence variables in the boudary layer it is assumed that production is
balanced by dissipation.
kρ
k+ = ( = Cµ−1/2 balanced) (2.37)
τw
and
y 1
+ = 3/2
= Cµ3/4 (2.38)
k κ
with:

k + , dimensionless turbulent kinetic energy k.

ρ, density.
1/2
τw , the shear stress (= ρ Cµ k).
Cµ , an emprical constant (Cµ ' 0.09).
+ , dimensionless turbulent dissipatiion .

2.3.2 Calculation of turbulent wall shear stresses

In subroutine FluxUVW the effective wall shear stress is used. Yet in the law
of the wall the shear stress is calculated differently. Start with the laminar

16
2.3 Calculation of wall shear stresses CYCLONE
Fluid Dynamics

formulation:  
∂U
τw = µef f (2.39)
∂n wall
or
Ut
τw = µef f (2.40)
∆n
Rewriting
τw
µef f =
Ut /∆n
ρ u2τ
=
Ut /∆n
ρ u2τ ∆n
=
Ut
ρ uτ ∆n µ uτ
=
µ Ut
ρ uτ ∆n µ
=
µ Ut /uτ
+
y µ
= (2.41)
u+

Finally the minimum viscosity possible is the laminar viscosity. Thus a simple
test on y + /u+ ≥ 1 is needed.

2.3.3 Rough walls

The law of the wall for a rough wall reads [39]:

 
+ 1 y
u (y) = ln + Brough , (2.42)
κ ks

with ks as Nikuradse’s sand roughness height (which is an absolute height), and

Brough an empirical constant equal to 8.5 [39].

Another description (commonly used for atmospheric boundary layers) is:

 
1 y − yd
u+ (y) = ln , (2.43)
κ y0

with yd as the displacement thickness (simple translation) and y0 , or also known

as z0 , an aerodynamic roughness parameter (which is not a physical height).

Set yd = 0 then the relationship between y0 and ks is:

z0 = ks /eκBrough . (2.44)

With Brough = 8.5 and κ = 0.4 follows y0 = 1/30 ks , or ks = 30y0 .

CFD-10xxxx 17
 CYCLONE 2 A transported scalar
Fluid Dynamics

Universal velocity profile, smooth and rough walls

40

Smooth (E = 9.0) k+
k+ = 7.5
k+ = 12.5
70
k+ = 20.0
k+ = 35.0
30 gh
rou 300
to
th
oo
sm 1 000
ally
y+ = 70 u lic
dra 3 000
hy
u+ [-]

m
20 n fro 10 000
tio fully rough
ansi
f tr 30 000
no
gio
Re 100 000
y+ = 5

10

0
0 1 2 3 4 5 6
log10(y+) [-]

Figure 2.4: Implementation of N.Scholz’s universal velocity profile [39]

Now extent equation (2.34) to cater for extra roughness1 .

1
u+ (y, ks ) = ln(Ey + ) − ∆B(ks ) , (2.45)
κ

which means that roughness retards the velocity at y + . Recast equation (2.45)
in the form:
1 1
u+ (y, ks ) = ln(Ey + ) − ln(∆B 0 (ks )) . (2.46)
κ κ
With ∆B 0 written as
∆B 0 (ks ) = 1 + f (ks ) (2.47)

automatically lets ∆B vanish when ks diminishes.

uτ ks
aerodynamically smooth: ν < 5
uτ ks
transition: 5< ν < 70 (2.48)
uτ ks
completely rough: ν > 70

2.4 Residuals

Consider the equation:

A ~x = B .
1
Adjusting E alone is sufficient, as higher values of E mean rougher walls, and even lower
values than 9.0 are numerically possible but physically useless. The method presented here
allows for a more direct interpretation.

18
2.4 Residuals CYCLONE
Fluid Dynamics

After some iterations the approximate solution is ~x n which is inexact and leaves
a non-zero residual:
A ~x n = B − ρ~n .
or
~ n = B − A ~x n .
ρ (2.49)
This is an absolute residue and does not give any clue of the degree of conver-
gence.

In the subroutine ’Set_Normalisation_Factors’ the residuals are normalised us-

ing the mass flux entering the domain:
X
Σṁin = ρi (u~i · Ai ) (≡ MP ) (2.50)
inlets

For the momentum equations the residual is normalised with the total momen-
tum m~v entering the domain:
X p
Mv = ṁi u~i · u~i (2.51)
inlets

The turbulence energy:

X
Mk = ṁi (u~i · u~i ) (2.52)
inlets

The turbulent dissipation:

Mk
Mε = (2.53)
L/V̄inlet
Note that the latter contains an arbitrary length scale L (and is thus not a very
reliable measure).

All the above mentioned normalisation factors depend on the mass flux flowing
into the domain. In case there is no inflow then ....

CFD-10xxxx 19
 CYCLONE 2 A transported scalar
Fluid Dynamics

20
2.4 Residuals CYCLONE
Fluid Dynamics

Big whirls have little whirls

That feed on their velocity;
And little whirls have
smaller whirls
And so on to viscosity.

L.F. Richardson, 1926

CFD-10xxxx 21
 CYCLONE 2 A transported scalar
Fluid Dynamics

22
 CYCLONE
Fluid Dynamics

3 Convective term and differencing

schemes

The discretisation of the convection term is:

Z Nfaces
X
∇ · (ρ~vr φ) dV = ~ · (ρ~vr φ)f
S
VP f =1
Nfaces
X
= ~ · (ρ~vr )f φf
S
f =1
Nfaces
X
= F φf (3.1)
f =1

were one recognises F as the mass flux through the face. Assume that the mass
flux is known and the process reduces to interpolating φ on the face.

One possibility was already introduced as equation (2.11) and it is called the
‘Central Differencing Scheme’:

φf = λf φN + (1 − λf )φP (3.2)

with λ defined as:

xf − xP
λf = . (3.3)
xN − xP
For non-smooth grids a correction can be applied. This scheme is called ‘CD’
in dolfyn.

Another central differencing scheme can be defined using the central nodes, the
gradients and the distances towards the face center. The last step is averaging
the two found values:
1
φf = (φP + ∇φP (~xf − ~xP ) + φN + ∇φN (~xf − ~xN )). (3.4)
2
This scheme is available as ‘CD2’.

Finally a very crude, but simple, method is avaible (just for testing purposes)
and that is simply averaging the values of φP and φN :
1
φf = (φP + φN ). (3.5)
2
Of course this scheme is second order for equidistant grids. This scheme is
named ‘CD3’. This scheme is also used in [28] and is the basis of the gradients
a major CFD code.

CFD-10xxxx 23
 CYCLONE 3 Convective term and differencing schemes
Fluid Dynamics

φD

flow direction

φf
φC
φU

U C D
face f

Figure 3.1: Definition of upwind, central and downwind nodes

All the central differencing schemes can cause unphysical, i.e. unbounded, re-
sults.

The classic scheme which is guaranteed to be bounded is the ‘Standard First

Order Upwind Differencing Scheme’:

φP if (~v · ~n)f > 0
φf = . (3.6)
φN if (~v · ~n)f < 0

See the classic texts and critics on the Upwind Scheme, Upstream Scheme,
Donor Cell method, or ‘Pigpen method’1 of Roache in 1972 [37], or Patankar in
1980 [32].

Next a combination of UDS and CDS is possible by blending the two schemes,
or blending any second order scheme with UDS:

φf = (1 − γ)φlow
f
order
+ γφhigh
f
order
. (3.7)

were γ is a fixed constant for the domain.

Although in [13] it is claimed that a blend between UDS and CDS should be
suffient, sometimes other bounded high resolution schemes are desired. Fol-
lowing [25][26] and [11] the Normalised Variable Approach (NVA) is intro-
duced. Under certain circumstances boundness can be achieved by introducing
a small amount of first order numerical ‘diffusion’. First start with the Convec-
tion Boundness Criterion for unstructured meshes.
1
“D.B. Spalding is purported to have referred to upwind differencing as the “pigpen method,”
the idea being that if ζ is taken as a concentration fraction, then we should smell the pigpen
when we’re on the downwind side, and not on the upwind side (except for diffusion effects).”
[37]

24
 CYCLONE
Fluid Dynamics

UDS
LUDS
~
φf
Quick CDS

3/4

1/2
3/8

~
1/2 1 φC

Figure 3.2: Convection Boundness Criterion with UDS, CDS and LUDS

The nomalised variable is defined as (see Figure 3.1):

φ − φU
φ̃ = . (3.8)
φD − φU

The value at face f is then a function of φ at U (upwind node), C (the central

node), and D (the downwind node): φ̃f = f (φ̃C ) with

φC − φU
φ̃C = (3.9)
φD − φU

and
φf − φU
φ̃f = . (3.10)
φD − φU
Now it is required that φC and φf are locally bounded between φU and φD , or

φU ≤ φC ≤ φD (3.11)

or
φU ≥ φC ≥ φD . (3.12)
The Convection Boundness Criterion is then simply:

0 ≤ φ̃C ≤ 1. (3.13)

In the Normalised Variable Diagram (NVD) the boundness criteria can be pre-
sented, see Figure 3.2, with φ̃f as a function of φ̃C . Within the grey region and
on the line φ̃f = φ̃C the scheme is bounded.

In the CBC the upwind value of φC is needed for the definition of φ̃C . Following
[25] again one can argue that φC is bounded by φU and φD and therefore also

CFD-10xxxx 25
 CYCLONE 3 Convective term and differencing schemes
Fluid Dynamics

by the interpolated values at φ− +

f and φf (see Figure 3.3). The new definition of
φ̃C is
φC − φ−f φ+f − φC
φ̃C = + − = 1 − −. (3.14)
φf − φf φ+f − φf

Recall Equation (2.11):

φf = λf φN + (1 − λf )φP (3.15)
with λ defined as:
xf − xP
λf = . (3.16)
xN − xP
Thus φ+
f − φC is (with φN = φD and φP = φC )

φ+
f − φC = λf φD + (1 − λf )φC − φC
= λf φD − λf φC
= λf (φD − φC )
φD − φC
= λf (xN − xC )
xN − xC
= λf (∇φ)f · d̂ (xN − xC ). (3.17)
With the unit vector d̂ parallel with CD
~d
d̂ = . (3.18)
|d|
−
For the two faces φ+
f and the virtual upstream face φf the process is analogous:
−
φ+
f − φf
φ+
f − φ−
f = − (x+ −
f − xf )
x+
f − xf
−
= (∇φ)C · d̂ (x+
f − xf ). (3.19)
Now Equation 3.14 turns into
φ+
f − φC
φ̃C = 1− −
φ+
f − φf

λf (∇φ)f · d̂ (xD − xC )
= 1− −
(3.20)
(∇φ)C · d̂ (x+
f − xf )

Using the definition of λ which is defined as the distance from point C to φ+ f

divided by CD) and the fact that C is located in the middle of the cell gives
−
x+
f − xf x+
f − xC
=2 = 2λf (3.21)
xD − xC xD − xC
and Eqn. (3.20) transforms into
(∇φ)f · d̂
φ̃C = 1−
2 (∇φ)C · d̂
(∇φ)f · d
= 1−
2 (∇φ)C · d
φD − φC
= 1− (3.22)
2 (∇φ)C · d

26
3.1 Standard upwind differencing, UD CYCLONE
Fluid Dynamics

because the product (∇φ)f · d is equal to φD − φC .

Using this framework all kinds of differencing schemes can be defined based on
the value of φ̃C .

3.1 Standard upwind differencing, UD

In this case simply use the upwind value.

φf = φC . (3.23)

It is called ‘standard (first order) upwind differencing’ as it is the oldest upwind

scheme and (still) often used as the default differencing scheme. It is very dissi-
pative, especially when the (hexahedral) grids do not follow, or are aligned with,
the main flow direction. However it is very robust which makes it an attractive
scheme for the startup phase of a simulation.

3.2 Second order central differencing, CD

As has been noted before the central differencing scheme can cause unphysical
oscillations (which can have desastrous effects for example when the turbulent
kinetic energy becomes negative). However, central differencing is important to
LES simulations as it preserves the turbulent kinetic energy better where upwind
schemes tend to decay the energy faster.

All three methods presented below more or less act in the same way (for example
with respect to oscillating results).

3.2.1 Linear interpolation, CD1

The standard method proposed by [13] is the linearly weighted central differ-
encing scheme. Including the non-orthogonal correction it reads
φf = λf φD + (1 − λf )φC + (∇φ)f (~xf 0 − ~xf ). (3.24)

3.2.2 Gradient based averaging, CD2

Another CDS method is based on averaging the results including both the up-
and downstream gradients as in
1
φf = (φC + ∇φC (~xf − ~xC ) + φD + ∇φD (~xf − ~xD )) . (3.25)
2
The attractivness of this formulation is that it automatically deals with non-
orthogonal, or skewed, cells.

CFD-10xxxx 27
 CYCLONE 3 Convective term and differencing schemes
Fluid Dynamics

3.2.3 Simple averaging, CD3

The most simple method, discarding all the geometrical information.

1
φf = (φC + φD ). (3.26)
2

3.3 Blended differencing

3.3.1 Constant blending

Using a fixed blending, by an amount of γ, between UD and CDS is the proposed

method in [13]:
φf = γφCD + (1 − γ)φUD . (3.27)

3.3.2 Variable blending, hybrid scheme

Instead of using a fixed blending factor in the domain one could differentiate on
the basis of the local Peclet number P e. The Peclet number is defined as
(ρ u)f ∆xf
Pe = . (3.28)
Γf

If the viscosity is used for the diffusion coefficient Γ it is immedeately obvious

that the Peclet number remsembles the Reynolds number Re. Hence the term
‘cell Reynolds number’ is often used as well.

In [13] and [32] it is shown that the oscillations depend on the value of P e.
When P e ≤ 2 at every grid node, no oscillations occur. ‘This is a suffient, but
not necessary condition for boundness of CDS solution’ [13]. Spaldings hybrid
scheme [41] was designed to switch from CDS to UD at any node at which
P e ≥ 2. ‘This is too restrictive and reduces the accuracy’ [13] as it is only
needed in regions where the gradient is high. Another way of looking at it is to
use the viscosity and regard it as the cell Reynolds number; even for very mild
air flow the Peclet number would jump up and UD would be used everywhere.
For a further discusion see [13].

The hybrid scheme is therefore not implemented.

3.3.3 Variable blending, power law scheme

Patankar introduced a scheme based on the exact one-dimensional solution which

should provide better results than the hybrid scheme [32]. It switches to UD
when P e ≥ 10.

28
3.4 Second order linear upwind differencing, LUD CYCLONE
Fluid Dynamics

The original formulation is of the form [32]

φ(x) − φ0 exp(P e Lx ) − 1
= (3.29)
φL − φ0 exp(P e) − 1

hence the name ‘power law’.

In [48] an equivalent expression is given:

φf = φC − βf (φD − φC ) (3.30)

for 0 < P e < 10 with

(1 − 0.1P ef )5
βf = . (3.31)
P ef

Especially the original formulation should not be used as the numerous calls
to the exponentional function would be too costly. The power law scheme, al-
though it performs better than the hybrid scheme, inherits the same defficiencies.
The power law scheme is not implemented.

3.4 Second order linear upwind differencing, LUD

The original formulation includes a second upwind node in order to find a gra-
dient and extrapolate to the face. As gradients are available the same idea can
be accomplished using:

φf = φC + ∇φC (~xf − ~xC ). (3.32)

Two forms are implemented: LUD using the CBC framework (and reverting to
UD outside the range) and the unbounded version LUX.

When the gradient ∇φC is limited –using a separate slope limiter– then the CBC
framework is superflous. In such (limited) cases the second order linear upwind
differencing scheme combines relative simplicity, robustness and second order
accuracy. Albeit that slope limiting might deteriorate convergence behavior; the
only downside of slope limited LUD.

3.5 Third order linear upwind differencing, QUICK

Historically the third order linear upwind differencing ‘Quadratic Upwind Inter-
polation for Convective Kinematics’, or QUICK is important. In principle this
scheme involves an additional upstream node (or φU U ). The search for such a
‘virtual’ node is not useful. Moreover the results in the past of QUICK were
rarely superior.

CFD-10xxxx 29
 CYCLONE 3 Convective term and differencing schemes
Fluid Dynamics

In the NVD diagram, see Figure 3.2, the QUICK scheme is found to be equal to
[2]
3 3
φf = + φC . (3.33)
8 4

The QUICK scheme is not implemented.

3.6 Second order switching schemes

3.6.1 MinMod differencing, MinMod

Following Alves et.al. [1]

3.6.2 Gamma differencing, GAM

3.6.3 Third order monotone upstream, MUSCL

The third order ‘Monotone Upstream-centered Scheme for Conservation Laws‘,

or MUSCL, is a blending of a CDS scheme and the LUD scheme. In short

φf = γφCDS + (1 − γ)φLUD . (3.34)

It has not (yet) been implemented. It can be easily done, either with or without
the CBC framework.

3.7 VOF scheme(s)

Ubbink en High Resolution Interface Capturing HRIC

3.8 Discusion

The schemes are not complete yet. On very skewed, really bad or dirty, meshes
instabilities still occur. This point needs further attention as a remedy has not
been found yet.

30
3.8 Discusion CYCLONE
Fluid Dynamics

φD
flow direction

φf+
φC
φU
φf−

U C D
−
face f face +
f

Figure 3.3:

minmod
~ LUDS

φf gamma
CDS

smart
1

Quick

~
β
1/3 1/2 1 φC
m

Figure 3.4: Smart, MinMod and Gamma-schemes

CFD-10xxxx 31
 CYCLONE 3 Convective term and differencing schemes
Fluid Dynamics

32
3.8 Discusion CYCLONE
Fluid Dynamics

The last decent thing written in C

was Schuberts 5th Sympony.

Ron Bell

CFD-10xxxx 33
 CYCLONE 3 Convective term and differencing schemes
Fluid Dynamics

34
 CYCLONE
Fluid Dynamics

4 Calculation of gradients

Currently only two methods have been used: Gauss and Least Squares. Both are
called from subroutine GradientPhi.

Both depend on the number of neighbours; the more the merrier (and the accu-
racy improves). And here tetrahedra, especially near boundaries, are suffering
from the reduced number of neighbours.

4.1 Gauss

Consider Gauss’ divergence theorem applied to a gradient:

Z Z
grad φ dV = φ d~s
V S
nf
1 X
(grad φ)P0 ≈ φj ~sj (4.1)
VP0
j=1

with φj the value at the centre of face j.

When the mesh is orthogonal only one pass is needed as the data is available to
calculate (reconstruct) the gradient. In case of non-orthogonal meshes at least
one extra pass is necessary for the correction.

This is done in subroutine GradientPhiGauss.

4.2 Least squares

Another possibility is based on the following thought: one knows the values at
the cell centre surrounding point P0 . Thus for one pair

d~j · (grad φ)P0 = φPj − φP0

or the set
d~j · (grad φ)P0 = φPj − φP0 (j = 1, . . . , n) (4.2)

CFD-10xxxx 35
 CYCLONE 4 Calculation of gradients
Fluid Dynamics

with d~j = ~xPj − ~xP0 the vector from P0 to Pj . This results in general to more
sets or pairs than there are unknowns. For example

δφ δφ δφ
dj+11 + dj+12 + dj+13 = φPj+1 − φP0
δx1 δx2 δx3
δφ δφ δφ
dj+21 + dj+22 + dj+23 = φPj+2 − φP0
δx1 δx2 δx3
δφ δφ δφ
dj+31 + dj+32 + dj+33 = φPj+3 − φP0
δx1 δx2 δx3
δφ δφ δφ
dj+41 + dj+42 + dj+43 = φPj+4 − φP0
δx1 δx2 δx3
..
.
δφ δφ δφ
dj+n1 + dj+n2 + dj+n3 = φPj+n − φP0
δx1 δx2 δx3

δφ δφ δφ
With the unknowns δx1 , δx2 , and δx3 . Now multiply equation (4.2) with the
transpose of d~j :

δφ
 
δx1
d~jT d~j ·  δφ
= d~jT ∆φPj,0 (j = 1, . . . , n) (4.3)
 
δx2 
δφ
δx3 P0

And the sum of (4.3) is:

δφ
 
nf δx1 nf
X X
d~jT d~j ·  δφ
= d~jT ∆φPj,0 (j = 1, . . . , n) (4.4)
 
δx2 
j=1 δφ j=1
δx3 P0

or
  δφ

 2  
Xnf δx 1 δx2 δx1 δx3 δx1 δx1 nf
X δx1 ∆φPj,0
 δx1 δx2 δx2 2 δx3 δx2  ·  δφ  δx2 ∆φPj,0 
=

 δx2 
j=1 δx1 δx3 δx2 δx3 δx3 2 j
δφ j=1 δx3 ∆φPj,0 j
δx3 P0

This looks like the classic:

A ~x = B

With
nf
X
A= d~jT d~j (4.5)
j=1

and
nf
X
B= d~jT (φPj − φP0 ) (4.6)
j=1

and the solution is ~x = (grad φ)P0 .

This is done in subroutine GradientPhiLeastSquares.

36
4.3 Gradient limiters CYCLONE
Fluid Dynamics

Note that A is a symmetric 3 x 3 matrix with only geometrical coefficients (the

distances from P0 to Pj ). So inverting A gives:
nf
X
−1
(grad φ)P0 = A d~jT (φPj − φP0 ) (4.7)
j=1

A−1 can be stored once and the gradient is simply calculated every time from
equation (4.7). Of course storing will need some extra memory (6 elements per
cell).

4.3 Gradient limiters

4.3.1 Barth and Jespersen

The definition of the cell centroid in eq. (2.2) is also crucial to the multi-
dimensional limiter of Barth and Jespersen (see [10]). In a cell the variation
of a variable can be described according to:

φ(~x) = φcentroid + grad φ · ∆~r (4.8)

where grad φ represents an estimate of the solution gradient in the cell from
surrounding centroid data (using Gauss or least squares) and ~r the position inside
the cell (~r = ~x − ~xcentroid ). Consider a ‘limited’ form of the reconstructed field
of φ in cell i:
φi (~x) = φi + Φ grad φi · (~x − ~xi ) (4.9)
with 0 ≤ Φ ≤ 1. Then find the largest admissible Φ not exceeding the maximum
and minimum of the neighbouring centroid values. The procedure is as follows.
First compute
φmin = min(φi , φneighbors ) (4.10)
and
φmax = max(φi , φneighbors ) (4.11)
then require
φmin ≤ φi (~x) ≤ φmax . (4.12)

For linear reconstructions, extrema in φi (~x) occur at the nodes, or vertices, of the
faces. The original procedure then visits each node k to determine the limited
value Φ which satisfies (4.12):
φmax −φk

 min(1, φk −φi ) if φk − φi > 0

min −φ
Φk = min(1, φφk −φ i
k
) if φk − φi < 0 (4.13)

1 if φk − φi = 0


with Φ = min(Φ1 , Φ2 , . . . , Φk ).

This procedure guarantees that the variables satisfy the monotonicity principle
when evaluated anywhere within a face (when using the nodes for the test). A

CFD-10xxxx 37
 CYCLONE 4 Calculation of gradients
Fluid Dynamics

less restrictive procedure would be satisfying at the face centres. In [10] it is said
that “numerical evidence shows that both methods work well”, however this is
not true for 2D models consisting of one layer of cells in the third dimension.

Finally using a different notation which follows the implementation. First find
the neighboring maximum and minimum values, then define the maximum change
of a neighbor of cell i as

∆max = φmax − φi
∆min = φmin − φi (4.14)

and the change at each node, centroid or face, k as

∆k = φk − φi . (4.15)

next define the ratio:

∆max


 ∆k if ∆k > 0
rk = ∆min (4.16)
∆k if ∆k < 0

1 if ∆k = 0


and the final limiter is Φ = min(r1 , r2 , . . . , rk ).

The Barth and Jespersen limiters come in three flavours: centroid, face, and
node based, with the face based version as the default. See subroutine Gradient-
PhiLimiterBarthJespersen.

4.3.2 Van Albada

The Barth and Jespersen method relies on the non-differentiable min-function.

Such functions can kick in or not, also called ’switching behaviour’, and as a
result of this non-smooth behaviour convergence (to steady state) is severely
hampered.

Using the (simple) arithmetic average which is defined as

n
1 1X
x̄ = (x1 + x2 + . . . + xn ) = xi . (4.17)
n n
1

Thus the average of two values a and b is

1
ave arithmetic (a, b) = (a + b). (4.18)
2
Which is equivalent to the central difference of a value at a face centre (using an
equidistant mesh).

Now suppose that b is related to a according to

b = xa, (4.19)

38
4.3 Gradient limiters CYCLONE
Fluid Dynamics

arithmetic(x)
1.4 harmonic(x)
van Albada
Barth Jespersen

1.2

0.8

y
0.6

0.4

0.2

0
0 0.5 1 1.5 2 2.5 3 3.5 4
x

Figure 4.1: Average as a result of fraction x

with x as the fraction of a, and thus eq. (4.18) becomes

1 1
ave arithmetic (a, xa) = (a + xa) = (1 + x) a. (4.20)
2 2

Another average is the harmonic average, typically appropriate for situations

when the average of rates is desired. Defined as
n n
H= 1 1 1 1 = Pn 1 , xi > 0 for all i. (4.21)
x1 + x2 + x3 ... + xn 1 xi

The harmonic average of two values a and b is thus

2
ave harmonic (a, b) = 1 1 (4.22)
a + b

which can be recasted into

|b|a + |a|b
ave harmonic (a, b) = (4.23)
|a| + |b|

a form of an earlier version of the formula used by van Leer [44]. Note that this
form has a risk when a and b approach zero.

Assume (4.19) again and insert it in equation (4.22) and the result is:

2x
ave harmonic (a, xa) = a, (4.24)
x+1
with limx→∞ = 2a.

CFD-10xxxx 39
 CYCLONE 4 Calculation of gradients
Fluid Dynamics

van Albada
1.4 50%
25%
10%

1.2

0.8
y

0.6

0.4

0.2

0
0 0.5 1 1.5 2 2.5 3 3.5 4
x

Figure 4.2: Influence of the bias ε2 on the van Albada function

Van Albada et al. [44] choose the following averaging variant:

(b2 + ε2 )a + (a2 + ε2 )b
ave van Albada (a, b) = , (4.25)
a2 + b2 + 2ε2
with ε2 as a very small non-vanishing bias.

Setting ε2 to zero yields

b2 a + a2 b
ave van Albada (a, b) = , (4.26)
a2 + b2
and using b = xa gives

(xa)2 a + a2 (xa)
ave van Albada (a, xa) =
a2 + (xa)2
x2 a3 + xa3
=
a2 + x2 a2
x2 + x
= a. (4.27)
1 + x2
In a similar way one can study the influence of ε. Assume ε = a, or the bias as
a fixed fraction of a:
((xa)2 + (a)2 )a + (a2 + (a)2 )(xa)
ave van Albada (a, xa) =
a2 + (xa)2 + 2(a)2
x + x + 2 + 2 x
2
= a. (4.28)
1 + x2 + 22
The result is shown in Figure 4.2. As long as ε is small, say less than 1% of a,
its influence is minimal.

40
4.3 Gradient limiters CYCLONE
Fluid Dynamics

The result of (4.25) tends to 12 (a + b) if a and b are finite differences of a smooth

solution, but tend to the smallest value where the solution is not smooth (if
b = xa then limx→∞ = a). Overall the form resembles the harmonic variant.
See also Figure 4.1.

Because the shape, or behavior, of (4.25), resembles the characteristics of the

Barth and Jespersen function (min(1, x)) one can use it as a limiter. Start with
finding the neighboring maximum and minimum changes ∆max and ∆min of
cell i as in (4.14), and the change ∆k at each node, centroid or face, k as in
(4.15), with the ratios rk of (4.16) and the function
rk2 + rk
φk = . (4.29)
1 + rk2
The final limiter is Φ = min(φ1 , φ2 , . . . , φk ). See also subroutine Gradient-
PhiLimiterAlbada.

Note that in this form of the limiter, i.e. (4.29), there is no need for a small bias
ε as a division by zero is not possible unlike the original form as in eq. (4.26).

4.3.3 Venkatakrishnan

Venkatakrishnan took van Albadas concept further with a function of the ratio
of two second degree polynomials (with boundary conditions as φ(∞) = 1 and
φ(0) = 0):
x2 + 2x
φ(x) = 2 . (4.30)
x +x+2

Originally his final form looks like [46]:

(∆2+ + 2 )∆− + 2∆2− ∆+
 
1
Φ= (4.31)
∆− ∆2+ + 2∆2− + ∆− ∆+ + 2

with ∆+ = ∆max and ∆− = ∆k with 2 again as a bias.

Setting 2 = 0 and recasting (4.31) into

∆2+ + 2∆− ∆+
Φ =
∆2+ + 2∆2− + ∆− ∆+
∆2+ (1 + 2 ∆ −
∆+ )
=
∆2+ (1 + 2( ∆ − 2
∆+ ) +
∆−
∆+ )

∆−
and with r0 = ∆+
2r0 + 1
Φ= . (4.32)
2(r0 )2 + r0 + 1
∆+ ∆max
Using the previous definitions r = 1/r0 = ∆− = ∆k

2( r1k ) + 1 rk2 + 2rk

φk = = (4.33)
2( r1k )2 + ( r1k ) + 1 rk2 + rk + 2

CFD-10xxxx 41
 CYCLONE 4 Calculation of gradients
Fluid Dynamics

van Albada
1.4 Venkatarishnan
polynomial (x<1.5)
Barth Jespersen

1.2

0.8
y

0.6

0.4

0.2

0
0 0.5 1 1.5 2 2.5 3 3.5 4
x

Figure 4.3: Venkatarishnan limiter vs Barth Jespersen and van Albada.

and the limiter is again Φ = min(φ1 , φ2 , . . . , φk ). The result compared to Barth

Jespersen and van Albada is shown in Figure 4.3. The code can be found in
GradientPhiLimiterVenkatarishnan.

4.3.4 Polynomial function

Clearly Figure 4.3 suggests that a polynomial function of the form

φ(x) = a3 x3 + a2 x2 + a1 x + a0 (4.34)
can be considered with either two or three sections.

For the case with two sections the boundary conditions are:
dφ(0)
φ(0) = 0 with =1
dx
and at x = xc
dφ(xc )
φ(xc ) = 1 with = 0.
dx

For φ(0) = 0 it simply follows that a0 = 0. Differentiating (4.34) results in

dφ(x) = ( 3a3 x2 + 2a2 x + a1 ) dx, (4.35)
dφ(0)
thus dx = 1 simply means that a1 = 1.

dφ(xc )
At xc is by definition φ(xc ) = 1 and dx = 0. With
dφ(xc )
= 3a3 x2c + 2a2 xc + 1 = 0
dx

42
4.3 Gradient limiters CYCLONE
Fluid Dynamics

follows
1 + 3a3 x2c
a2 = − . (4.36)
2xc
Finally
1 + 3a3 x2c 2
a3 x3c − xc + xc = 1
2xc
or
1 3
a3 x3c − xc − a3 x3c + xc = 1
2 2

1 1
− a3 x3c + xc = 1
2 2
so
x−2
a3 = . (4.37)
x3c

The case of xc = 1.5 is plotted in Figure 4.3.

4.3.5 Uniform regions

Assume a uniform flow, field distribution, where ∆min ≈ ∆max , then it is de-
sirable to circumvent the limiter. Note that in such cases the ratio’s rk are ap-
proximately unity and in case the Venkatarishnan limiter is used the gradient
is reduced by 25% (in equation (4.33) is φ(1) = 3/4); by contrast the Barth
Jespersen and van Albada limiters do not alter the gradient in that point.

Michalak and Ollivier-Gooch suggest to ‘switch the limiter off’ when [29]

3
δφ ≡ (∆max − ∆min ) < (k∆x) 2 . (4.38)

Switching the limiter off should not introduce new switching behavior. This can
be accomplished by:
Φ = σ + (1 − σ)Φ̃ (4.39)

where Φ̃ is the active limiter and σ is

δφ2 ≤ (k∆x)3


  1 
δφ2 −(k∆x) 3
σ= s (k∆x)3
(k∆x)3 < δφ2 < 2(k∆x)3 (4.40)

0 δφ2 ≥ 2(k∆x)3


with the transition function s(x) is defined by, see Figure 4.4

s(x) = 2x3 − 3x2 + 1. (4.41)

CFD-10xxxx 43
 CYCLONE 4 Calculation of gradients
Fluid Dynamics

1
transition function s(x)

0.8

0.6
y

0.4

0.2

0
0 0.2 0.4 0.6 0.8 1
x

Figure 4.4: Transition function s(x)

4.4 Discussion

Gradient limiters can be used for any CBC scheme and even for the definition
of φ̃C . However some care is needed as the risk of extra dispersion is there.

The gradient limiters seem useful for the gradient based second order central
differencing scheme (CD2) and the unbounded form of the second order linear
upwind differencing scheme (the unbounded form as the gradient limiter already
ensures boundness).

44
4.4 Discussion CYCLONE
Fluid Dynamics

Users of CFD codes should especially be careful,

as the optimism of salesman is legendary.

J.H. Ferziger, M. Períc.

CFD-10xxxx 45
 CYCLONE 4 Calculation of gradients
Fluid Dynamics

46
 CYCLONE
Fluid Dynamics

5 Heat transfer

For non-isothermal flows the energy equation reads:

∂
(ρ cp T ) + div( ρ ~vr cp T ) = div( λ grad T ) + qT (5.1)
|∂t {z } | {z } | {z } |{z}
unsteady term convective term diffusive term source terms

The primitive variable is the temperature T thus in order to adhere to the standard
scalar transport equation the heat capacity is rearranged to:
∂ λ qT
(ρ T ) + div( ρ ~vr T ) = div( grad T ) + (5.2)
∂t cp cp

Where λ/cp appears as the thermal diffusivity ΓT .

Near solid walls either the wall temperature Tw or the heat flux qw through the
wall can be specified. In both cases the heat transfer coefficient is needed.

The heat transfer coefficient h is a simple way to lump the heat transfer in one
empirical formula. Starting from

Q̇ = h A (Tw − T )

and rearranging
Q̇/A qw
h= = (5.3)
Tw − T Tw − T
Thus once the heat flux qw and the temperatures Tw and T are known the heat
transfer coefficient is defined by equation (5.3). Likewise if instead of the tem-
perature in the cell adjacent to the wall T some other reference or bulk tem-
perature Tbulk is used, another value of h is found. This latter value should
correspond to the values found in empirical handbooks of heat transfer.

Recall that λ/cp is used as the thermal diffusivity in the fluid domain. Introduc-
ing the (laminar) Prandtl number:
cp µ
Pr = (5.4)
λ
hence
λ µ
= (5.5)
cp Pr
which links and scales the laminar viscosity directly to the thermal diffusivity.

For an given wall temperature Tw and heat transfer coefficient h the heat source
Q̇ = qw A = h A (Tw − T ) near the boundary cell is known and splits in two

CFD-10xxxx 47
 CYCLONE 5 Heat transfer
Fluid Dynamics

Ts Tb Tw
T
α δn
λ q

h x

Figure 5.1: Heat transfer near a wall

parts: h A Tw and h A T , the former enters the source term ST for the boundary
cell and the latter enters the diagonal term Ap as h A.

When the heat flux qw is used as a boundary condition then the wall temperature
is first calculated. This value is then used in the source term.

The temperature boundary condition Tw does not need to be specified at the

solid fluid interface. Assuming one dimensional heat conduction through a wall
of s thick and a heat conduction of λs gives (see Figure 5.1):

λs
q= (Tb − Tw ) (5.6)
s
In the surroundings with Ts and an assumed heat transfer coefficient of αs results
in:
q = αs (Ts − Tb ) (5.7)
The heat flux through all three components must be equal thus

λs
q = αs (Ts − Tb ) = (Tb − Tw ) = h (Tw − T ). (5.8)
s
The overall temperature difference is then:

Ts − T = (Ts − Tb ) + (Tb − Tw ) + (Tw − T )

q q q
= + +
αs λ/s h
1 1 1
= q( + + ) (5.9)
αs λ/s h
P
and can be regarded as a chain of resistances as in q = ∆T / Ri or

∆T
q= 1 P . (5.10)
h + Ri

s
R=
λ

48
5.1 Laminar heat transfer CYCLONE
Fluid Dynamics

5.1 Laminar heat transfer

In the case of laminar flow the heat transfer is dictated by heat conduction from
the centre of the wall cell to the wall over a distance of δn. The heat flux simply
reads:
(Tw − T )
q = λ
δn
µ cp (Tw − T )
=
Pr δn
1
= P r δn
(Tw − T )
µ cp
1
= 1 P (Tw − T ). (5.11)
h + Ri
With the local heat transfer coefficient h
µ cp
h= . (5.12)
P r δn
Note that λ/cp is the thermal diffusivity ΓT in the transport equation. Also the
heat source (or sink) for the boundary cell is needed. In the code itself

Q̇ = H (Tw − T ) (5.13)

is used, with
1
H= 1 P A. (5.14)
µ/(P rδn) + cp Ri
However for output purposes H has to be transferred into h again
cp
h=H (5.15)
A
with a heat flux in W/m2
qw = h (Tw − T ) (5.16)
and a total heat balance of
X
Q= qwi Ai . (5.17)

5.2 Turbulent heat transfer

The turbulent heat transfer is basically modelled analog to the viscous boundary
layer. Where the law of the wall profile u+ = f (y + ) is scaled by the turbulent
Prandtl number P rt (or σT ). In short dimensionless form:

T + = P rt u+ (5.18)

However, in general a sublayer resistance factor (commonly denoted by P ) by

C. Jayatilleke is accounted for as well:
h  P r 3/4 ih  
−0.007 P r i
P = 9.24 − 1 1 + 0.28 exp (5.19)
P rt P rt

CFD-10xxxx 49
 CYCLONE 5 Heat transfer
Fluid Dynamics

which brings eqn. (5.18) to

T + = P rt (u+ + P ) (5.20)

The dimensionless temperature T + is defined by:

ρ µt (Tw − T )
T+ = (5.21)
qw

Using the definition of the heat transfer coefficient (5.3), and eqns. (5.20) and
(5.21) gives
qw ρ Cp µt
h= = (5.22)
Tw − T P rt (u+ + P )
One of the great consequences of eqn. (5.22) is that one does not need to perform
a thermal simulation in order to find the local heat transfer coefficients. On the
other hand, a local heat transfer might become high, but when the temperature
difference Tw − T is negligible, there will hardly be any heat flux and thus no
cooling (or warming up).

50
5.2 Turbulent heat transfer CYCLONE
Fluid Dynamics

Programming is all about readability.

M. Sahami, Stanford, 2009

CFD-10xxxx 51
 CYCLONE 5 Heat transfer
Fluid Dynamics

52
 CYCLONE
Fluid Dynamics

6 Lagrangian particles

6.1 Basic equations

Multi-phase flows often contain particles or droplets. Let ~u be the velocity of

the continuous phase, ~vp the current droplet velocity and ~xp the droplet position.

Then the Newtonian equation of the motion of a droplet with mass mp is

d~vp X ~
mp = Fi . (6.1)
dt
i

After summing all the forces the equation of motion (6.1) can be solved using

d~xp
= ~vp . (6.2)
dt

6.1.1 Drag force

The most important force is the drag force

1
F~d = Cd (Re) ρf |~u − ~vp |(~u − ~vp ) Sp (6.3)
2

with the drag coefficient Cd as a function of the droplet Reynolds number Re.
Sp is the frontal area of the particle.

The droplet Reynolds number is defined by

ρf |~u − ~vp | Dp
Rep = . (6.4)
µ

The drag coefficient is



 24, 000. , Re ≤ 0.001


 24/Re , 0.001 ≤ Re ≤ 0.2
Rep 24 (1 + 0.15 Rep0.687 )/Rep , 0.2 < Re ≤ 1000 (6.5)
0.44 , 1000 < Re ≤ 200 000.




0.10 , Re > 200 000.


CFD-10xxxx 53
 CYCLONE 6 Lagrangian particles
Fluid Dynamics

6.1.2 Virtual mass force

For accelerating or decelerating particles the surrounding carrier fluid needs to

be accelerated or decelerated too. Potential flow theory, also know as irrotational
flow (without viscosity), shows that an accelerating particle give rise to a force
of:
d(~vp − ~u)
F~v = −Cv ρf Ωp . (6.6)
dt
Using potential flow theory the ‘virtual mass coefficient’ Cv is found to be equal
to 0.5.

For a steady state flow equation (6.6) reduces to

d~vp
F~v = −Cv ρf Ωp . (6.7)
dt
The mass of a particle mp then simply increases by Cv ρf Ωp to (ρp + Cv ρf )Ωp .
Thus for a particle with density ρp of 1000. kg/m3 in air (ρf ≈ 1.2 kg/m3 ) the
increase is neglible (just 0.06%).

6.1.3 Body forces

Gravity, or any other external force, can also act on a particle. It is incorpo-
rated into the ‘body force’. When gravity acts on a particle it will also let a
particle float or the buoyant force is equal to the weight of the displaced fluid
(Archimedes’s principle). The net result is
F~b = mp ~g = (ρp − ρf )Ωp ~g . (6.8)

6.1.4 Adding all forces

Starting with equation (6.1) and using Sp = 3

2 Ωp /Dp , ρ0p = ρp + Cv ρf and
~vs = ~u − ~vp :
d~vp X
mp = F~i
dt
i
d~vp
ρp Ωp = F~d + F~v + F~b
dt
d~vp 1 3 Ωp
(ρp + Cv ρf )Ωp = Cd ρf k~vs k (~u − ~vp ) + (ρp − ρf ) Ωp ~g
dt 2 2 Dp
d~vp 1 3 1
ρ0p = Cd ρf k~vs k (~u − ~vp ) + (ρp − ρf ) ~g
dt 2 2 Dp
d~vp 3 ρf ρp − ρf
= Cd k~vs k (~u − ~vp ) + ~g (6.9)
dt 4 ρ0p Dp ρ0p

Equation (6.9) can be solved analytically provided some terms are assumed to
be constant:
d~vp
= C1 (~u − ~vp ) + C2 ~g (6.10)
dt

54
6.1 Basic equations CYCLONE
Fluid Dynamics

τ
1 V1

v/V1

0/0 1 2 3 4 5 6
t /τ

Figure 6.1: Characteristic response of a particle to forces

with
3 ρf
C1 = Cd k~vs k (6.11)
4 ρ0p Dp
and
ρp − ρf
C2 = . (6.12)
ρ0p
Rearanging 6.10 to
d~vp
= dt
C1 (~u − ~vp ) + C2 ~g
and integrating
−1
ln( C1 (~u − ~vp ) + C2 ~g ) = t + c
C1
C1 (~u − ~vp ) + C2 ~g = e−C1 t+c
C1 (~u − ~vp ) + C2 ~g = C e−C1 t
or
C2 C −C1 t
~vp = (~u + ~g ) − e (6.13)
C1 C1

At t = 0 the particle has an initial velocity of ~vp0 thus C is:

C
~vp0 = C1 (~u − ~vp ) + C2 ~g −
C1
or
C = ( C1 ~u + C2~g ) − C1~vp0 .
The particle velocity is then
C2 C2
~vp = ( ~u + ~g ) − ( ~u + ~g − ~vp0 ) e−C1 t (6.14)
C1 C1
| {z } | {z }
~v1 ~v2

and is illustrated in Figure 6.1. In a quiescent fluid (~u = 0) and after some time
a particle will drop subjected to the action of gravitational force ~g
C2
~vp∞ = ~g
C1

CFD-10xxxx 55
 CYCLONE 6 Lagrangian particles
Fluid Dynamics

or
(ρp − ρf ) Dp ~g
~vp∞ = 3 .
4 Cd ρf k~vs k
For small particles the drag coefficient equals to 24/Re. Thus with k~vs k =
k~u − v~p k = kv~p k
(ρp − ρf )
~vp∞ = 24 3 ρf ~g
Re 4 Dp k~vp k
(ρp − ρf )
= 24µ 3 ρf
~g
ρf k~vp kDp 4 Dp k~
vp k
(ρp − ρf )
= 18µ 1 ~g
Dp Dp
(ρp − ρf ) 2
= Dp ~g . (6.15)
18µ
Which is the well known particle settling or sedimentation velocity for small
particles or ‘Stokes’ law’.

The characteristic time in the response of equation (6.14) can be found using

dv~p 
~vp∞ = ~vpo + ∆t
dt  t=0
or
~vp∞ = ~vpo + (C1 (~u − ~vp ) + C2 ~g ) ∆t

and consider again a quiescent fluid (~u = 0) and a starting velocity ~vp0 = 0.

The characteristic time τmomentum or τM is then

C2
~g = ~vpo + (C1 (~u − ~vp ) + C2 ~g ) ∆t
C1
C2
~g = C2 ~g ∆t
C1
1
∆t =
C1
4 ρ0p Dp
τM = (6.16)
3 Cd ρf k~vs k
which meaning becomes clear in Figure 6.1. Note that the gravitational effects
are cancelled out and the characteristic time is only based on the drag force.

For light or small particles in the Stokes’ range equation (6.16) reduces to
4 ρ0p Dp
τM = 24
3 Re ρf k~vs k
4 ρ0p Dp
= 24µ
3 ρf k~vs k Dp ρf k~vs k
ρ0p Dp2
= (6.17)
18 µ
note that the division of |~u − ~vp |/|~u − ~vp | which goes to 0/0 when the slip
velocity drops to zero has been eliminated.

56
6.2 Integrating particles in time CYCLONE
Fluid Dynamics

6.2 Integrating particles in time

The most common method to integrate equation (6.2) is to use a simple explicit
method known as ‘explicit Euler’ or ‘forward Euler’. It is so popular because it
relies only on the currently known values at time level tn :
φn+1 = φn + f (φn , tn ) ∆t (6.18)
It is only first order accurate but more important it has restrictions in the size of
∆t for stability (set by τM ).

Equation (6.17) also shows what happens when very light and small particles are
going to be used. Consider a tiny (Dp = 1 10−6 m) droplet of water (ρp = 1000.
kg/m3 ) in air (µ = 18 10−6 Pa s) when τM will drop to ≈ 3 10−6 s. Therefore
explicit Euler has it’s limitations when any kind of droplet has to be catered for.

An alternative is to use ‘implicit Euler’ or ‘backward Euler’ method

φn+1 = φn + f (φn+1 , tn+1 ) ∆t (6.19)
in which there is, in principle, no restriction to the time step size ∆t. The dis-
advantage is that the values of f (φn+1 , tn+1 ) at the new time tn+1 are unknown
and extra work is necessary.

When equation (6.9) is casted in the form of (6.10)

d~vp
= C1 (~u − ~vp ) + C2 ~g
dt
with only ~vp and t as variables then it can be solved analytically. The solution
is:
~vp = ~v1 + ~v2 e−C1 t (6.20)
with
C2
~v1 = ~u + ~g
C1
and
C2 ~ − ~vp0
~v2 = ~u + ~g − ~vp0 = v1
C1
and the constants C1 and C2 defined in 6.11 and 6.12 and the characteristic time
τM = 1/C1 .

The solution of ~vp in equation (6.20) approaches ~v1 rapidly (see Figure 6.1).
How close ~vp is to the final solution ~v1 depends on it’s initial velocity ~vp0 , or
the ratio of ~vp0 /~v1 . The closer the ratio is to 1 the larger the time step can be
allowed to be.

The velocity 6.20 can be integrated once again to obtain the new position:
~v2
~xp = ~v1 t + (1 − e−C1 t ) (6.21)
C1

xxxxx

The basic scheme is:

CFD-10xxxx 57
 CYCLONE 6 Lagrangian particles
Fluid Dynamics

Figure 6.2: Example of particles passing through cells

1. choose a suitable time step ∆t

2. use the current particle velocity Vp , or fluid velocity U (whichever is

higher), to get an estimate of the time needed to pass through a cell

3. calculate the new particle velocity ~vp

4. use equation (6.21) to find the new particle position using the new particle
velocity ~vp

5. check if the particle has left the cell

The starting time step ∆t is estimated using the Courant number defined by

U δt
Co = . (6.22)
δx
As there are three dimensions the choose the lowest of
Co ∆x Co ∆y Co ∆z
∆tc = min( , , ). (6.23)
ux uy uz

And a value of Co ≈ 0.34 in order to step in two or three steps through a typical
cell (see an example in Figure 6.2).

6.3 Particle within a cell

One of the major tasks is to find out whether a particle or droplet is within or
outside of a cell. The easiest way to find out for an arbitrary cell is to use the
cell face normals (positive pointing outwards, see Figure 6.3) and a normalised
position vector from the particle to the cell face centre

x̂pf = ~xcell face − ~xp . (6.24)

Next simply calculate the dot product using the cell face normal ~xn

x̂pf · ~xn

58
6.4 Intersecting a cell face CYCLONE
Fluid Dynamics

particle

particle outside cell

Figure 6.3: Example of particle inside and outside of a cell (red)

particle velocity

λ>0
λ>0
particle
λ<0

λ<0

λ<0 



 
 

 
 



Figure 6.4: Particle intersects face

The result is

 x̂pf · ~xn > 0 particle inside
x̂ · ~xn = 0 particle on the face
 pf
x̂pf · ~xn < 0 particle outside

Simply loop over all faces and check if the dot product does not become nega-
tive.

6.4 Intersecting a cell face

Once a particle is within a cell it will travel in some direction. It is then nec-
essary to find out on which cell face the particle will leave the cell. When this
cell face is known then the neighbouring cell is by virtue of dolfyns unstruc-
tured addressing method also known. Hence there is no need for special particle
searching algorithms. Consider a particle in a cell at point ~xp with a velocity ~vp .
Any point on the line, including an intersection point ~xi made with this velocity
is defined as
~xi = ~xp + λ ~vp . (6.25)

CFD-10xxxx 59
 CYCLONE 6 Lagrangian particles
Fluid Dynamics

A cell face is defined by its cell face centre and normal. For any point on the
cell face, including an intersection point ~xi , the normal is perpendicular to the
vector from the cell face centre to the intersection point. In such a case the dot
product is equal to zero. Hence

(~xi − ~xcell face ) · ~xn = 0 . (6.26)

Insert 6.25 into 6.26 and solve for λ

(~xp + λ ~vp − ~xcell face ) · ~xn = 0

(~xp − ~xcell face ) · ~xn + λ ~vp · ~xn = 0

or
(~xp − ~xcell face ) · ~xn
λ=− (6.27)
~vp · ~xn
The lowest positive λ found is the intersection point. The point itself is then
found using equation (6.25).

When the dot product ~vp · ~xn = 0 then ~vp is parallel to the face defined by ~xn .
A special check is needed in case two, or more, faces share the same normal
vector.

6.5 Interpolating velocities

The Lagrangian procedure described here requires that the local fluid velocity
at some point within a cell is available. One possible method is to start directly
from the cell centred velocity components and interpolate from all neighbouring
cells. This path will lead to large interpolating molecules surrounding the cell
the particle is currently in.

Another method is based on interpolating the velocities on to the corner nodes of

the cells first (using the cell centred values and the velocity gradients) and then
use these velocities as a starting point. Next the velocity inside the cell can be
interpolated using for example shape functions. Such functions, however, need
to take into account the cell topology in order to be able to map parameters into
a computational domain. This method requires the cell topology to be used in
dolfyn which is not desirable.

The currently used interpolation method is based on weighting the (interpolated)

velocity components on the nodes as in

NX
nodes

d= Wi di . (6.28)
i=1

Furthermore it is a requirement that when the interpolation point coincides with

a node its value should be equal to the nodal value.

60
6.5 Interpolating velocities CYCLONE
Fluid Dynamics

Finally the interpolation method should yield values no smaller than the mini-
mum di and no larger than the maximum di . Thus
NX
nodes

Wi = 1 , with 0 ≤ Wi ≤ 1 . (6.29)
i=1

For any point within a cell the next weighting function is used

( r1 )2
Wi = P i 1 2 (6.30)
( ri )

with
ri = k~xnode − x~p k .

As the weight function is known the fluid velocity ~up at a particle position ~xp is
then simply
NX
nodes

~up = Wi ~unodei . (6.31)

i=1

CFD-10xxxx 61
 CYCLONE 6 Lagrangian particles
Fluid Dynamics

6.6 Examples

Figure 6.5: Free falling particles, influence of particle Courant number

62
6.6 Examples CYCLONE
Fluid Dynamics

Figure 6.6: Free falling particles, details showing steps within cells

CFD-10xxxx 63
 CYCLONE 6 Lagrangian particles
Fluid Dynamics

Figure 6.7: Particles in a stagnation flow

64
6.6 Examples CYCLONE
Fluid Dynamics

Figure 6.8: Particles in a stagnation flow, influence of particle Courant number

CFD-10xxxx 65
 CYCLONE 6 Lagrangian particles
Fluid Dynamics

Figure 6.9: Small particles/streamlines in a stagnation flow

66
6.6 Examples CYCLONE
Fluid Dynamics

Figure 6.10: Small particles/streamlines in an 45 angle test

CFD-10xxxx 67
 CYCLONE 6 Lagrangian particles
Fluid Dynamics

68
 CYCLONE
Fluid Dynamics

Appendices

CFD-10xxxx 69
 CYCLONE 6 Lagrangian particles
Fluid Dynamics

70
 CYCLONE
Fluid Dynamics

A Geometry input files

You will need three geometry files casename.cel, casename.vrt, and casename.bnd,
where casename as the name of your case. The preprocessor only needs these
simple files:

• A file with cell definitions (*.cel)

• A file with the vertices (*.vrt)

• And the boundaries (*.bnd)

The results are stored in casename.geo.

The definitions of the cells in *.cel are:

nr , i 1 , i 2 , i 3 , i 4 , i 5 , i 6 , i 7 , i 8 , i d

With the cell number ‘nr’, and the 8 corner nodes of a hexaeder (i1 t/m i8). Each
cell has a pointer to an ‘id’ in a table. This ‘id’ is reserved for the future (for the
time being set it to ‘1’).

The accompanying ‘format statement’ is a ‘free format’ (*).

Cell shapes are defined as (see Figure A.1):

Hexaeder i1, i2, i3, i4, i5, i6, i7, i8

Prism or wedge i1, i2, i3, i3, i5, i6, i7, i7

Pyramid i1, i2, i3, i4, i5, i5, i5, i5

4 3 6 5 6 5
4 3
7 8 7 8

8 8
7 7
1 2 1 2 1 1
2

5 6 5 6 3 4 2 34
hexahedral prism or wedge pyramid tetrahedral

Figure A.1: Cell definitions

CFD-10xxxx 71
 CYCLONE A Geometry input files
Fluid Dynamics

Tetraeder i1, i2, i3, i3, i5, i5, i5, i5

The vertices in *.vrt are:

nr , x , y , z

With the number ‘nr’ and it’s position in 3D (x, y, z). The accompanying ‘format
statement’ is: i9,6x,3g16.9.

Finally the boundaries in *.bnd are:

nr , i 1 , i 2 , i 3 , i 4 , i d

With number ‘nr’ and the 4 corner nodes of a face of a hexaeder (i1 t/m i4). The
boundaries point to an ‘id’ which values will be set in the control file (*.din).
All boundary faces without an ‘id’ will be given the default id of ‘0’.

The accompanying ’format statement’ is in ’free format’ (*). Two variants are
defined:

Quad i1, i2, i3, i4

Tri i1, i2, i3, i3

72
 CYCLONE
Fluid Dynamics

B Mathematical toolbox

B.1 Gauss’ theorem

Let V be a volume in R3 with a smooth and closed surface S with an outwards

facing normal ~n = (n1 , n2 , n3 ). Further let w
~ = (w1 , w2 , w3 ) be a differentiable
vector field then
ZZZ ZZ
div(w)~ dx dy dz = ~ · ~n dS
w (B.1)
V S

RR RR H
~ · ~n dS (or S w
The integral S w ~ · ~n dΩ, or S w
~ · ~n dΩ) is called the flux of
w
~ through S.

B.2 Sum of two vectors

Graphically, a vector is represented by an arrow, defining the direction, and the

length of the arrow defines the vector’s magnitude. If one end of the arrow is
the origin O and the tip of the arrow is A. Then the vector may be represented
algebraically by OA (also as ~a, a, or a (bold)) with magnitude ||~a||.

Two vectors, ~a and ~b are equal if they have the same magnitude and direction,
regardless of whether they have the same initial points. A vector having the
same magnitude as ~a but in the opposite direction to ~a is denoted by −~a.

The sum of two vectors is defined by:

~c = ~a + ~b (B.2)

or
ci = ai + bi
when ~a + ~b = ~0 then ~b = −~a or bi = −ai . Following properties hold:

1. ~a + ~b = ~b + ~a (commutative property)

2. ~a + (~b + ~c) = (~a + ~b) + ~c (associative property)

3. ~a + ~0 = ~a

4. ~a + (−~a) = ~0

CFD-10xxxx 73
 CYCLONE B Mathematical toolbox
Fluid Dynamics

y
C
A

a c

b
B
x

Figure B.1: Sum of two vectors

Vector subtraction is defined as the difference of two vectors, ~a − ~b, is a vector

~c that is, ~c = ~a − ~b, or ~c = ~a + (−~b).Thus vector subtraction can be represented
as a vector addition.

B.3 Multiplying a vector

Is λ  R, (λ 6= 0), and ~a a vector, ~a 6= ~0, then the length of ||λ~a|| = |λ| ||~a||.
For λ = 0, or ~a = 0, is λ~a = ~0.

Following multiplication properties hold:

1. λ(µ~a) = (λµ)~a (associative property)

2. 1~a = ~a

3. λ(~a + ~b) = λ~a + λ~b (distributive property)

4. (λ + µ)~a = λ~a + µ~a (distributive property)

B.4 Length of a vector and vector products

Is λ  R, (λ 6= 0), and ~a a vector, ~a 6= ~0, then the length of ||λ~a|| = |λ| ||~a||.
For λ = 0, or ~a = 0, is λ~a = ~0.

A unit vector is one which has a magnitude of 1. Is ~a a vector (~a 6= ~0), then the
unit vector in the direction of ~a is â. Suppose â = λ~a (λ > 0), then

1 = ||â|| = λ~a (B.3)

and so
1
λ=
||~a||

74
 CYCLONE
Fluid Dynamics

B
c b
φ C
x

Figure B.2: Inner product of two vectors (hashed area of outer product)

The unit vector in the direction of ~a is

~a
â = (B.4)
||~a||

Replacing ~a by â is called normalising of ~a (not defined for ~a = ~0).

The breaking up of a vector into it’s component parts is known as resolving a

vector. Notice that the representation of ~a by it’s components, eg. ax and ay is
not unique. Depending on the orientation of the coordinate system with respect
to the vector in question, it is possible to have more than one set of components.

The inner product, scalar, or dot product is defined by:

X
~a · ~b = ai bi (B.5)
i

which results in a scalar.

The scalar product of two vectors, ~a and ~b denoted by ~a · ~b, is defined as the
product of the magnitudes of the vectors times the cosine of the angle between
them, as illustrated in Figure B.2.

~a · ~b = ||~a|| ||~b|| cos φ (B.6)

or
~a · ~b
cos φ = (B.7)
||~a|| ||~b||

And in particular when â · â = b̂ · b̂ = 1, since the angle between a vector and

itself is 0 and the cosine of 0 is 1.

Alternatively, when â · b̂ = 0, since the angle between ~a and ~b is 90◦ and the
cosine of 90◦ is 0.

In general then, if â · b̂ = 0 and neither the magnitude of ~a nor ~b is 0, then ~a and

~b must be perpendicular.

CFD-10xxxx 75
 CYCLONE B Mathematical toolbox
Fluid Dynamics

Following properties hold:

1. ~a · ~b = ~b · ~a (commutative property)

2. ~a · (~b + ~c) = ~a · ~b + ~a · ~c (distributive property)

3. (λ~a) · ~b = ~a · (λ~b) = λ(~a · ~b)

4. When ~a · ~a ≥ 0 and ~a · ~a = 0 then ~a = 0.

Let a vector normal (a vector perpendicular to a plane) be denoted as ~n, a posi-

tion vector of some point in a plane denoted as ~a and a typical point on the plane
denoted as ~r. Then the vector ~b = ~r −~a lies within the plane and is perpendicular
to ~n. Thus the vector equation for a plane is:

(~r − ~a) · ~n = 0 (B.8)

The cross or outer product ~a × ~b of two vectors ~a and ~b in R3 is defined as:

 
a2 b3 − a3 b2
~a × ~b =  a3 b1 − a1 b3  (B.9)
a1 b2 − a2 b1

which is a vector with length

||~a × ~b|| = ||~a|| ||~b|| sin φ (B.10)

where φ is the angle between ~a and ~b (0 ≤ φ ≤ π).

Following properties hold:

1. ~a × ~b = −~b × ~a

2. (~a + a~0 ) × ~b = ~a × ~b + a~0 × ~b

3. ~a × (~b + b~0 ) = ~a × ~b + ~a × b~0

4. (λ~a) × ~b = ~a × (λ~b) = λ(~a × ~b)

Using the elementary trigonometric definition of the cosine in Figure B.3:

||~un ||
cos α =
||~uP ||
From the definition of the dot product:

~a · ~b
cos φ =
||~a|| ||~b||
or here with the surface normal ~xn0 facing inwards:
~xn0 · ~uP
cos φ =
||~xn0 || ||~uP ||

76
 CYCLONE
Fluid Dynamics

y
uP

un
α ut
P

x n’

face
x

xn

Figure B.3: Vector components near a wall

The length of the normal ~xn0 is by definition 1, thus

||~un || ~xn0 · ~uP

=
||~uP || ||~uP ||
or
||~un || = ~xn0 · ~uP (B.11)
Now only the length of ~un is known (||~un ||) without a direction. Vector ~un
points in the direction of ~xn0 thus multiplying the normal with the length results
in
||~un || ~xn0 = (~xn0 · ~uP ) ~xn0
or simply
~un = (~xn0 · ~uP ) ~xn0
Next recall that the normal should face outwards, and with ~xn0 = (−~xn ) follows

~un = ((−~xn ) · ~uP ) (−~xn )

Pulling −1 out
~un = (~xn · ~uP )(−1)(−1)~xn
results in
~un = (~xn · ~uP ) ~xn (B.12)
With other words the resulting vector ~un is independent of the choice of the
normal (pointing inwards or outwards).

The tangential component then follows from an subtraction of the now two
known vectors (~uP and ~un ):

~ut = ~uP − ~un (B.13)

CFD-10xxxx 77
 CYCLONE B Mathematical toolbox
Fluid Dynamics

78
Bibliography CYCLONE
Fluid Dynamics

Bibliography

[1] M.A. Alves, P. Cruz, A. Mendes, , F.D. Magalh aes, and P.J. Pinho,
F.T.Oliveira. Adaptive multiresolution approach for solution of hyperbolic
PDEs. Comput. Methods Appl. Mech. Engrg., 191:3909–3928, 2002.
[2] M.A. Alves, P.J. Oliveira, and F.T. Pinho. A convergent and universally
bounded interpolation scheme for the treatment of advection. Int. J. Numer.
Methods Fluids, 41:47–75, 2003.
[3] A.A. Amsden. KIVA-3: A KIVA program with block-structured mesh for
complex geometries. Technical Report LA-12503-MS, Los Alamos Nat.
Lab., March 1993.
[4] A.A. Amsden. KIVA-3V: A block-structured KIVA program for engines
with vertical or canted valves. Technical Report LA-13313-MS, Los
Alamos Nat. Lab., July 1997.
[5] A.A. Amsden and F.H. Harlow. The SMAC method: A numerical tech-
nique for calculating incompressible fluid flows. Technical Report LA-
4370, Los Alamos Nat. Lab., February 1970.
[6] A.A. Amsden, P.J. O’Rourke, and T.D. Butler. KIVA-II: A computer pro-
gram for chemically reactive flows with sprays. Technical Report LA-
11560-MS, Los Alamos Nat. Lab., February 1989.
[7] A.A. Amsden, J.D. Ramshaw, L.D. Cloutman, and P.J. O’Rourke. Im-
provements and extensions to the KIVA computer program. Technical Re-
port LA-10534-MS, Los Alamos Nat. Lab., October 1985.
[8] A.A. Amsden, J.D. Ramshaw, P.J. O’Rourke, and J.K. Dukowicz. KIVA:
A computer program for two and three-dimensional fluid flows with chem-
ically reactions and fuel sprays. Technical Report LA-10245-MS, Los
Alamos Nat. Lab., February 1985.
[9] J.D. Anderson. Computational Fluid Dynamics, the basics with applica-
tions. McGraw-Hill, 1995.
[10] T.J. Barth and D.C. Jespersen. The design and application of upwind
schemes on unstructured meshes. In Proceedings of the AIAA 27th
Aerospace Sciences Meeting, Reno, Nevada, January 1989.
[11] M.S. Darwish and F. Moukalled. Normalized variable and space formula-
tion methodology for high-resolution schemes. Num. Heat Transfer, Part
B, 30:217–237, 1994.
[12] J.H. Ferziger and M. Períc. Computational Methods for Fluid Dynamics.
Springer, Berlin, 1996.

CFD-10xxxx 79
 CYCLONE
Fluid Dynamics

[13] J.H. Ferziger and M. Períc. Computational methods for fluid dynamics. 3rd
rev. ed. Springer, Berlin, 2002.
[14] J.E. Fromm and F.H. Harlow. Numerical solution of the problem of vortex
street development. Physics of fluids, 6, July 1963.
[15] R.A. Gentry, B.J. Daly, and A.A. Amsden. KIVA-COAL: A modified ver-
sion of the KIVA program for calculating the combustion dynamics of a
coal-water slurry in a Diesel engine cylinder. Technical Report LA-11045-
MS, Los Alamos Nat. Lab., October 1987.
[16] U. Ghia, K. N. Ghia, and C. T. Shin. High-re solutions for incompressible
flow using the navier-stokes equations and a multigrid method. Journal of
Computational Physics, 48:387–411, 1982.
[17] A.D Gosman, W.M. Pun, A.K. Runchal, B.D. Spalding, and M. Wolfh-
stein. Heat and mass transfer in recirculating flows. Academic Press,
1969.
[18] F.H. Harlow and C.W. Hirt. Generalized transport theory of anisotropic
turbulence. Technical Report LA-4086, Los Alamos Nat. Lab., May 1969.
[19] F.H. Harlow and P.I. Nakayama. Turbulence transport equations. Physics
of fluids, 10, November 1967.
[20] F.H. Harlow and P.I. Nakayama. Transport of turbulence energy decay rate.
Technical Report LA-3854, Los Alamos Nat. Lab., February 1968.
[21] F.H. Harlow and J.E. Welch. Numerical calculation of time-dependent
viscous incompressible flow of fluid with free surface. Phys. Fluids, 8,
1965.
[22] J.O. Hinze. Turbulence. McGraw-Hill, 1975.
[23] C.W. Hirt, A.A. Amsden, and J.L. Cook. An arbitrary lagrangian-eulerian
method for all flow speeds. J. Comput. Phys., 14:227–253, 1974.
[24] C.W. Hirt, B.D. Nichols, and N.C. Romero. SOLA: A numerical solution
algorithm for transient fluid flows. Technical Report LA-5852, Los Alamos
Nat. Lab., April 1975.
[25] H. Jasak. Error analysis and estimation in the Finite Volume method with
applications to fluid flows. PhD thesis, Imperial College, University of
London, 1996.
[26] H. Jasak, H.G. Weller, and A.D. Gosman. High resolution NVD differ-
encing scheme for arbitrarily unstructured meshes. Int. J. Numer. Methods
Fluids, 31(2):431–449, 1999.
[27] N.L. Johnson. The legacy and future of CFD at Los Alamos.
[28] S.R. Mathur and J.Y. Murthy. A pressure-based method for unstructured
meshes. Num. Heat Transfer, Part B, fundamentals, 31:2:195–215, 1997.
[29] K. Michalak and C. Ollivier-Gooch. Limiters for unstructured higher-order
accurate solutions of the euler equations. In Forty-sixth Aerospace Sciences
Meeting, 2008.

80
Bibliography CYCLONE
Fluid Dynamics

[30] Y. Nakayama. Visualized flow; fluid mootion in basic and engineering

situations revealed by flow visualization. Pergamon Press, 1988.

[31] B.D. Nichols, C.W. Hirt, and R.S. Hotchkiss. SOLA-VOF: A solution
algorithm for transient fluid flow with multiple free boundaries. Technical
Report LA-8355, Los Alamos Nat. Lab., August 1980.

[32] S.V. Patankar. Numerical heat transfer and fluid flow. Hemisphere Publ.
Corp.; McGraw-Hill, 1980.

[33] M. Períc. A finite volume method for the prediction of three-dimensional

fluid flow in complex ducts. PhD thesis, Imperial College, University of
London, August 1985.

[34] S.B. Pope. Turbulent flows. Cambridge University Press, 2000.

[35] C.M. Rhie and W.L. Chow. Numerical study of the turbulent flow past an
airfoil with trailing edge separation. AIAA Journal 21 (11), 1983.

[36] W.C. Rivard, O.A. Framer, and T.D. Butler. RICE: A computer program for
multicomponent chemically reactive flows at all speeds. Technical Report
LA-5812, Los Alamos Nat. Lab., March 1975.

[37] P.J. Roache. Computational fluid dynamics. Hermosa Publishers, Albu-

querque NM, 1972.

[38] P.J. Roache. Verification and Validation in Computational Science and

Engineering. Hermosa Publishers, Albuquerque NM, 1998.

[39] H. Schlichting. Boundary-layer theory. Transl. by J. Kestin. 7th ed.

McGraw-Hill, 1979.

[40] C.T. Shaw. Using computational fluid dynamics. Prentice Hall, 1992.

[41] B.D. Spalding. A novel finite difference formulation for differential ex-
pressions involving both first and second derivates. Int.J. for numerical
methods in engineering, 4:551–559, 1972.

[42] J. C. Tannehill, D. A. Anderson, and R. H. Pletcher. Computational Fluid

Mechanics and Heat Transfer. Taylor and Francis, Washington, D.C.,
1997.

[43] H. Tennekes and J. L. Lumley. A First Course in Turbulence. MIT Press,

Cambridge, MA, 1983.

[44] G.D. van Albada, B. van Leer, and W.W. Roberts. A comparative study
of computational methods in cosmic gas dynamics. Astronomy and Astro-
physics, 108:76–84, 1982.

[45] M. Van Dyke. An Album of Fluid Motion. The Parabolic Press, 1982.

[46] V. Venkatakrishnan. On the accuracy of limiters and convergence to steady

state solutions. In Proceedings of the AIAA 31st Aerospace Sciences Meet-
ing, Reno, Nevada, January 1993.

[47] H.K. Versteeg and W. Malalasekera. An Introduction to Computational

Fluid Dynamics. Longman, 1995.

CFD-10xxxx 81
 CYCLONE
Fluid Dynamics

[48] H.K. Versteeg and W. Malalasekera. An Introduction to Computational

Fluid Dynamics, The Finite Volume Method. Pearson, 2nd edition, 2007.

[49] J.E. Welch, F.H. Harlow, J.P. Shannon, and B.J. Daly. The MAC method: A
cumputing technique for solving viscous, incompressible, transient fluid-
flow problems involving free surfaces. Technical Report LA-3425, Los
Alamos Nat. Lab., March 1966.

[50] P. Wesseling. Principles of computational fluid dynamics. Springer, Berlin,

2000.

[51] D.C. Wilcox. Turbulence Modeling for CFD. DCW Industries, 1993.

[52] M. Zijlema, A. Segal, and P. Wesseling. Finite volume computation of 2D

incompressible turbulent flows in general coordinates on staggered grids.
Technical Report 94-24, TU Delft, 1994.

82
 CYCLONE
Fluid Dynamics

Index

Least Squares, 33
blending, 22 lid driven cavity, 69, 85
body force, 13 Los Alamos, 4
boundary LUD, 27
definition, 68 LUX, 27
shapes, 68
boundness, 22 non-orthogonal diffusion, 14
Normalised Variable Approach, 22
CBC, 23 Normalised Variable Diagram, 23
CD, 21 NVA, 22
CD2, 21 NVD, 23
CD3, 21
CDS, 11, 21 OpenDX, iii
cell over-relaxed approach, 14
definitions, 67
preprocessor, 67
shapes, 68
prism, 68
centroid, 9
pyramid, 68
combustion, 4
Convection Boundness Criterion, 23 Richardson Extrapolation, 69
deffered correction, 12 Sparc, iii
differencing scheme Sparsekit2, iii
central, 11, 21 stress tensor, 13
central average, 21 subroutine
central w. gradients, 21 FluxUVW, 16
fixed blending, 22, 26 GradientPhi, 33
gamma scheme, 28 GradientPhiGauss, 11, 33
hybrid scheme, 26 GradientPhiLeastSquares, 34
linear upwind, 27
minmod scheme, 28 tests, 85
power law, 26 tetraeder, 68
quick scheme, 27
upwind, 11, 22, 25 UDS, 11, 22

FluxUVW, 16 vector
surface, 10
g95, iii verification, 69
Gamma scheme, 28 vertices, 68
Gauss, 33
GradientPhi, 33 wedge, 68
GradientPhiGauss, 11, 33
GradientPhiLeastSquares, 34

hexaeder, 68

83