PhD course: Implementation of a finite volume

solver for incompressible, two-dimensional

laminar and turbulent flow
Lars Davidson
Division of Fluid Dynamics
Dept. of Mechanics and Maritime Sciences
Chalmers University of Technology
SE-412 96 Göteborg
http://www.tfd.chalmers.se/~lada, lada@chalmers.se
June 30, 2019

1 CFD course
1.1 Background
The traditional method for CFD in industry and universities is Reynolds-Averaged
Navier-Stokes (RANS). It is a fast method and mostly rather accurate. In every
industry which is working on fluid mechanics, CFD is used.

1.2 Aim
Unfortunately, most engineers and many researchers have limited knowledge of
what a CFD CFD code is doing. The object of this course is to close that
knowledge gap.

1.3 When?
This PhD course is given every second year (2015, 2017, 2019 . . . ) in Study
period 3. If you want to take part, register on the course by sending me an
Email. The course page can be found at PhD course in CFD

1.4 Lectures
There will be no lectures. There will be one meeting (1-2 hours) per week. In
these meetings the students have the opportunity for questions and discussion.
The first meeting will take place in week 1 (Study period 3).

1
2 Part I
Write a 2D Navier-Stokes solver for incompressible, laminar flow. Compute the
flow in a square lid-driven cavity at Re = Uwall L/ν = 1000; set Uwall = L = 1
and ν = 1/Re. Compare with numerical simulations in [1]. A mesh with 12 × 12
nodes (constant size) should be used. Compare with data on the course www
page. Write a small report (max 10 pages). This gives 4 ETC credits.

• Collocated grid arrangement, i.e non-staggered

• The convection-diffusion can be treated as in assignment K2 in the MTF072
CFD course

The convection scheme for first-order upwind is derived in the Appendix A

A flow chart could look like this. You should first look at the book [5] as
well as the reports for CALC-BFC [3] and TEACH3D [4]. The reports can be
downloaded here. You should also look at the lecture notes (Chapter 2-9).

1. Compute all geometrical quantities such as areas, volume, interpolation

factors. I recommend that you put the boundary nodes at the boundaries,
see Chapter 4, p. 10 in Lecture notes1
2. Set all variables (u, v, p) as well as the convections, ṁe and ṁn , to zero.
3. Compute the coefficients, aW , aE , aS and aN . Use ṁe and ṁn in Eq. 2
(from the previous iteration) when computing the convective fluxes; the
first iteration they are zero. The coefficients are the same for u and v.
Initially you should use 1st-order upwind. Later on, you could try central
differencing or second-order upwinding.
4. Compute the source terms of u and v, i.e. −∂p/∂x and −∂p/∂y and put
them in SU ; for u is reads SU = −∂p/∂x∆V . The boundary condition for
the pressure is ∂p/∂n = 0 at all boundaries.
5. Introduce implicit under-relaxation (see pp. 52-53, Chapter 6 in Lecture
notes1 ), of say 0.5, and compute modified SUmod and amod
P .
6. Solve u and v (recall that SU is different for u and v) . Use either TDMA
or Gauss-Seidel (G-S) If you use TDMA (see Chapter 7), use it along both
x (I lines) and y (J lines) directions. Note that you must probably make
a couple of G-S sweeps for p′ (10 or more) in order to make the equations
converge, i.e.

do nn=1,nsweep(nphi)
do i= 2,nim1
do j= 2,njm1
1 http://www.tfd.chalmers.se/~lada/comp fluid dynamics/lecture notes.html

2
 rhs =ae(i,j)*phi(i+1,j,nphi)+aw(i,j)*phi(i-1,j,nphi)
. +an(i,j)*phi(i,j+1,nphi)+as(i,j)*phi(i,j-1,nphi)
. +su(i,j)
phi(i,j,nphi)= rhs/ap(i,j)
end do
end do
end do

where nsweep(nphi) > 3 where nphi=u, v or p’.

7. Compute ṁ∗e = 0.5(uE + uP )Ae ρ and ṁ∗n with central differencing
8. Add Rhie-Chow as
dp = pEE − 3pE + 3pP − pW
dv = dp Ae /4/auP e (1)
ṁe = ṁ∗e + ρAe dv

see p. 55 in Chapter 6 (collocated grid) in Lecture notes1 . auP e denotes

aP for the velocity interpolated to face e; Ae denotes the area of the east
face.
9. Compute the coefficients for the p′ equation. The boundary condition at
all boundaries is ∂p′ /∂n = 0. This is achieved by setting the corresponding
coefficient to zero (e.g. at the east boundary aE = 0).
10. Compute the source term in the p′ equation (i.e. net mass flux) as SU =
−(ṁe − ṁw + ṁn − ṁs )
11. The initial conditions of p′ is zero, i.e. set p′ = 0 in all points (at every
iteration).
12. Solve p′ using TDMA
13. Since we have Neumann boundary conditions on p′ at all boundaries, the
level of p′ it not defined. Choose a cell (e.g. cell (2, 2)), and make p′ = 0
(i.e. p′new = p′ − p′I=J=2 ).
14. Correct ṁe and ṁn . For ṁe it reads

ṁnew
e = ṁe + aE (p′P − p′E ) (2)

where aE is the coefficient in the p′ equation which is computed as

auP e = 0.5(auP,I+1 + auP,I ), aE = ρA2e /auP e

15. Correct the velocities, u and v. For u, it reads

p′w − p′e
unew
P = uP + ∆V
∆xauP

3
 16. Correct p as (αP ′ ≃ 0.5)

pnew
P = pp + αP ′ p′P

where αP is the under-relaxation coefficient (αP ≃ 0.5).

17. Compute residuals for u and v as RHS minus LHS for each cell, i.e.
I=N
X I−1 J=N
X J−1
R= |aE uE + aW uW + aN uN + aS uS + SU − aP uP )|
I=2 J=2

and scale it with an appropriate value (see Chapter 4, p. 15 in Lecture

notes2 )
The residual for the continuity equation is simply the magnitude of the
net mass flux in each cell (which is the source term in the p′ equation).
18. If the residuals are larger than, say, 10−3 , go to Item 3 and repeat (next
iteration)

Compare your results with data on the course www page.

3 Part II
Here you should implement the two-equation k − ω model of Wilcox [6]. It reads
  
∂v̄i k ∂ µt ∂k
ρ = Pk + µ+ − ρβ ∗ kω
∂xi ∂xj σk ∂xj
  
∂v̄i ω ω ∂ µt ∂ω
ρ = Cω1 Pk + µ+ − Cω2 ρω 2 (3)
∂xi k ∂xj σω ∂xj
 
ρk ∂v̄i ∂v̄i ∂v̄j
µt = , Pk = µt +
ω ∂xj ∂xj ∂xi

where β ∗ = 0.09, cω1 = 5/9, cω2 = 3/40 and σk = σω = 2. The momentum

equation reads
 
∂ 1 ∂ p̄B ∂ ∂v̄i
(v̄i v̄j ) = − + (ν + νt ) (4)
∂xj ρ ∂xi ∂xj ∂xj

where the cross-diffusion term has been neglected. The turbulent kinetic energy
in the Boussinesq assumption is also neglected (or it can be assumed to be
incorporated in the pressure, i.e. p̄B = p̄ + 2k/3, see [2]).
The boundary conditions for the turbulent quantities are k = 0 at the wall
and ωw = 6ν/(cω2 y 2 ) for the cells adjacent to the walls. The ω boundary
condition is set by using SP and SU as SP = −1E10, SU = 1E10 · ωw . You
need to use some reasonable initial conditions on k and ω, otherwise the solution
2 http://www.tfd.chalmers.se/~lada/comp fluid dynamics/lecture notes.html

4
may easily diverge. I recommend to use kinit = 10−2 and ωinit = 100. Under-
relaxate turbulent viscosity as νtnew = ανt +(1−α)νtold . Note that you may need
to use rather strong under-relaxation on all variables (even down to α = 0.1).
The Reynolds number is 100 000. A mesh with 66 × 66 nodes are used with
15% stretching from the walls to the center.
Compare your results with data on the course www page.
Write a short report (max 10 pages). This part gives an additional 3.5 ETC
credits.

A First-order upwind
The convective term reads
∂vΦ
=0 (5)
∂y
We integrate and get
Qe − Qw = 0 (6)
where the fluxes are Qn (north face) and Qs (south face). First-order upwind
gives

Qs = max(Cs , 0)ΦS − max(−Cs , 0)ΦP

(7)
Qn = max(Cn , 0)ΦP − max(−Cn , 0)ΦN

where Cs = (vA)s and Qs = (vAΦ)s . We get

Qs − Qn
= max(Cs , 0)ΦS − max(−Cs , 0)ΦP − max(Cn , 0)ΦP + max(−Cn , 0)ΦN (8)
= max(Cs , 0)ΦS + max(−Cn , 0)ΦN − (max(−Cs , 0) + max(Cn , 0))ΦP

which can be written

Qs − Qn = Cs ΦS − Cn ΦP , Cs , Cn > 0
(9)
Qs − Qn = Cn ΦN − Cs ΦP , Cs , Cn < 0

(note that − max(−Cn , 0) = Cn if Cn < 0)

Add the continuity equation times ΦP , i.e (Cs − Cn )ΦP , to the last line of
Eq. 8 so that

Qs − Qn
= max(Cs , 0)ΦS + max(−Cn , 0)ΦN − (max(Cs , 0) + max(−Cn , 0))ΦP (10)
= max(Cs , 0)(ΦS − ΦP ) + max(−Cn , 0)(ΦN − ΦP )

The first line is obtained because

max(−x, 0) + x = max(x), max(x, 0) − x = max(−x) (11)

5
Equation 10 (using Eq. 6) can now be written

Qs − Qn = aS (ΦS − ΦP ) + aN (ΦN − ΦP ) = 0 (12)

where
aS = max(Cs , 0)
(13)
aN = max(−Cn , 0)

Equation 12 can be written on usual SIMPLE finite volume form as

a P ΦP = a N ΦN + a S ΦS
(14)
aP = aS + aN

References
[1] O. Botella and R. Peyret. Benchmark spectral results on the lid-driven
cavity flow. Computers & Fluids, 27(4):421–433, 1998.
[2] L. Davidson. Fluid mechanics, turbulent flow and turbulence modeling. Di-
vision of Fluid Dynamics, Dept. of Applied Mechanics, Chalmers University
of Technology, Gothenburg, 2014.
[3] L. Davidson and B. Farhanieh. CALC-BFC: A finite-volume code employ-
ing collocated variable arrangement and cartesian velocity components for
computation of fluid flow and heat transfer in complex three-dimensional
geometries. Rept. 95/11, Dept. of Thermo and Fluid Dynamics, Chalmers
University of Technology, Gothenburg, 1995.
[4] L. Davidson and P. Hedberg. A general computer program for transient,
three-dimensional, turbulent recirculating flow. Technical report, Dept. of
Applied Thermo and Fluid Dynamics, Chalmers University of Technology,
Gothenburg, 1986.
[5] H. K. Versteegh and W. Malalasekera. An Introduction to Computational
Fluid Dynamics - The Finite Volume Method. Longman Scientific & Tech-
nical, Harlow, England, 1995.
[6] D. C. Wilcox. Reassessment of the scale-determining equation. AIAA Jour-
nal, 26(11):1299–1310, 1988.