The finite volume method for

convection-diffusion problems
TMMV08 – Assignment 1
Division of Applied Thermodynamics and Fluid Mechanics
Department of Management and Engineering, Linköping University
SE-581 83, Linköping

January 8, 2024

Introduction
The dynamics of fluid flow are controlled by various transport processes that can be mathematically
described by the governing Navier-Stokes equations. A generalization of these processes can be
established for an arbitrary scalar quantity ϕ (e.g., u, v, w and T ) by introducing the generic
transport equation (here on conservative form):

∂(ρϕ)
+ div(ρϕu) = div(Γ grad ϕ) + Sϕ (1)
| ∂t
{z } | {z } |
convective term
{z } |{z}
diffusive term source term
transient term

where u is the velocity vector, ρ fluid density and Γ the diffusion coefficient. Γ can vary depending
on the variable of choice, e.g. u → Γ = dynamics viscosity or T → Γ = thermal conductivity.
This differential equation shows that the property ϕ can be altered by changes over time (transient
term), presence of a velocity field (convection term), due to spatial gradients (diffusive term) and
anything else that may impact the scalar distribution (e.g. gravity) but cannot accounted in the
former terms (source term). In order to solve this partial differential equation for a given problem
a numerical method can be adopted.
One of the most common numerical methods used in CFD is the so called finite volume method
(FVM), which integrates all terms in the equation(s) over a specific control volume (i.e. grid or
mesh) where also gradients are approximated by different discretization schemes (e.g. the up-
wind differencing, central differencing or QUICK scheme). These numerical schemes have different
inherent properties that influence the accuracy of the results depending on the flow problem.
Diffusive-like transport can also appear artificially in the solution, commonly known as false
diffusion (or numerical diffusion). This unwanted phenomenon is influenced by the choice of nu-
merical method strategy, e.g. mesh type/quality as well as choice of discretization schemes, and
can have detrimental effects on the CFD predictions. As a practitioner it is therefore essential
to understand the underlying theory associated to these numerical approaches and how they may
impact the current solution.

1
Assignment
This assignment consist of a range of tasks/questions that should be addressed in a written report
(see Examination). This will be done by performing different simulation exercises in MATLAB
and ANSYS CFX, derived from the course book Versteeg and Malalasekera (2007) and reading
suggested literature. In the first two parts you are to solve different steady-state one-dimensional
(1D) diffusion problems (Part I) and convection-diffusion problems (Part II) using MATLAB. In the
final part (Part III)you are going to used the commercial solver ANSYS CFX to study the impact
from artificial diffusion for different simulation strategies that can contaminate the accuracy of the
results.

Part I – Diffusion
In this problem the scalar transport process is only driven by pure diffusion under stationary
conditions, i.e. in Eq. 1 the transient, convective and source term are all vanished:

div(Γ grad ϕ) = 0 (2)

Here the transport of property ϕ in Eq. 2 will be solved numerically using the FVM for a given
set of 1D problems given in Versteeg and Malalasekera (2007, Ch. 4). Start by reading Ch. 2.5
(p. 24–26), Ch. 4 (p. 115–133) and relevant lecture notes thoroughly to make sure that you have
understood the origin of Eq. 1; the difference between differential and integral form, and the basic
concepts of the the FVM. Thereafter continue to study Example 4.1 (p. 118–121). A MATLAB
script that solves this problem for one specific setup will be provided. However, to be able to answer
the questions related to this part of the assignment you need to understand this implementation
and the related theory.

Part II – Convection & Diffusion

In this problem the scalar transport process is governed by both convection and diffusion (stationary
and source-free):

div(ρϕu) = div(Γ grad ϕ) (3)

Here the transport of property ϕ in Eq. 3 will also be affected by a velocity field (which in this
case already is known, see Ch. 6.1, p. 179–180 for further explanation). Read Ch. 5 (p. 134–178)
in Versteeg and Malalasekera (2007), which is about the procedure of solving convection-diffusion
problems. After you have obtained a good understanding about the different discretization schemes
and their properties, focus on solving the examples given in the Ch. 5.1–5.6 and Ch. 5.9 (i.e.
Example 5.1, 5.2 and 5.4) using the corresponding schemes (central -, upwind -, and QUICK
differencing scheme) you have read about. A MATLAB script that solves Example 5.2 (Case 1)
using the upwind differencing (UD) scheme will be provided. Here you need to modify the current
script in order to implement the other schemes. Note that you may not be able to use the script
unless you have good understanding about the different methods and involving parameters.

2
Comments

• For large system of algebraic equations (involving large matrices) an iterative solver is typically
used to obtain the final solution, read Versteeg and Malalasekera (2007, Ch. 7.6, p. 223–229).
In iterative methods, the solution accuracy will depend on the level of convergence of the
involved equations. For this assignment the Gauss-Seidel method have been provided with as
a function to your current scripts. Open and study the script in relation to the theory.

• The term accuracy can refer to i) the order of accuracy of the discretization scheme, or
ii) the accuracy of the final solution. The order of accuracy of the numerical scheme will
naturally contribute to the final solution accuracy, however, you need to be careful within
your discussion and distinguish between ”order of accuracy” and ”solution accuracy”.

Part III – False Diffusion

False diffusion will to some level always be present in a discretized solution of a partial differential
equation due to numerical modeling inaccuracies (such as numerical scheme and type of mesh).
This will result in erroneous smearing of the solution variables (analog to adding extra viscosity)
with general loss of accuracy. In this part you are to investigate the level of false diffusion in some
commonly used numerical schemes. This phenomena is described in Versteeg and Malalasekera
(2007, p. 150–151). Consider a two-dimensional (2D) domain with the a mesh-misaligned flow
(Fig. 1a and Fig. 5.14 Versteeg and Malalasekera (2007, p. 150)) using pure quadrilateral elements.

Figure 1: Fluid domain with quadrilateral mesh and boundary conditions for studying the effect of false
diffusion. The solid line indicates the flow direction. The exact solution along the X–X diagonal (dashed
line) exhibit a step change from ϕ = 100 to ϕ = 0 when crossing the solid diagonal.

For you convenience a complete setup of these cases have been prepared in ANSYS Workbench
with CFX, and are at your disposal. A description on how to setup these cases are given in the
Appendix, which we encourage you to follow. Identify the settings of the current setup and make
the necessary adjustments to perform the tasks and answer the provided questions.

3
Numerical model details:

• The domain dimensions are 1 x 1 x 0.01 m

• As CFX only can handle three-dimensional (3D) domains, one cell was used in the thickness
direction
• An extra transport equation for the convective term was added (Additional Variable 1) in
order to solve the scalar field along with the standard fluid flow solver
• As we intend to only investigate the numerical diffusion, i.e. not accounting for any physical
diffusion, the Kinematic Diffusivity was set to zero (Flow Analysis 1>Fluid Models>Additional
Variable 1 )
• Boundary condition are:
– Inlet: Velocity components with addition of the scalar value condition
– Outlet: Static pressure (relative pressure of zero, i.e. 1 atm)
– Symmetry condition on the renaming faces
– Free-slip wall on the top and lower wall in Case 2
• The scheme for solving the convective term (Advection Scheme) of the Additional Variable
is set to Upwind (first-order). Here you also need to change to the other requested schemes
• A root-mean-square (RMS) convergence criteria of 10−6 is used for all equation residuals
• No turbulence modeling is used

4
Tasks and Questions
Perform and answer the following tasks and questions by reading the relevant literature and do the
simulation exercises addressed for each part of this assignment. For TQ2–TQ4 present your results
in graphs against the provided analytical solutions.

Part I & II
TQ1. Briefly explain (in word only) the main steps for solving the transport equation problem
numerically using the finite volume method in the domain of interest. Clarify the key steps moving
from the governing equation to its integral and discretized form. What is the Gauss divergence
theorem used for? How is a large system of nodal points typically solved? What is the role of the
boundary condition and where do they enter the problem?

TQ2. Why do we need to care about convergence in numerical analysis? Explain how the con-
vergence level is defined in the provided scripts. What does the ”residual ” term in the script
quantify/measure in this case, and how is this related to the convergence criteria? Evaluate how
the solution accuracy is influenced by comparing three different convergence levels (10−1 , 10−2 and
10−6 ) using Example 5.2 Case 2 in the book (p. 149).

TQ3. How is the solution accuracy influenced by the mesh (node) density? Here you need to
perform a so called mesh sensitivity study. Compare three different mesh resolutions of the domain
(i.e. N, 2N vs. 4N, where N is the total number of nodal points). When does the solution start to
get mesh insensitive? Evaluate Example 5.2 Case 2 in the book (p. 149). Also, explain why the
solution for the pure diffusion case (Example 4.1, p. 118) is mesh insensitive.

TQ4. Investigate the solution accuracy and boundedness features of the central -, upwind -, and
quadratic-upwind (QUICK) differencing scheme by adjusting the required parameter (Peclet num-
ber and mesh density) in the concerned examples. Briefly discuss the properties: order of accuracy,
conservativeness, boundedness, transportiveness for each scheme in relation to your findings. High-
light different results that show the strength and weakness of each scheme.

Part III
TQ5. Investigate the impact of false diffusion in miss-aligned flow using a quadrilateral mesh in
ANSYS CFX using the following CFX schemes: Upwind, QUICK, and High Resolution, for three
(3) different mesh resolutions: 10 x 10, 50 x 50 vs. 100 x 100. Identify the level of false diffusion
qualitatively using surface plots (post-processing hints in the Appendix ) and/or across the X–X
diagonal (see Fig. 5.15, p. 150 in the book) and discuss the schemes properties in relation to these
findings.

5
Examination
Intended Learning Outcomes

Under successful completion of this assignment you should be able:

• to demonstrate the general procedures involved in a finite volume method for solving simpler
transport equations governing fluid flow problems
• to understand the basic properties of common numerical methods (i.e. numerical schemes,
grid, convergence, etc.) and how these choices may influence the accuracy of the solution
• to apply appropriate analyzing techniques in order to facilitate interpretation and discussion
of results in relation to a set of aims (questions) of an assignment

Written report. The answer to each task/question (TQ1–5) in this assignment should be pre-
sented in a written report by the group. You are encouraged to be inspired by relevant literature
and have discussions with you fellow students, but the results and written text should be your own
work.
You need to pay special attention to present your results in a clear, concise and informative
way, that best conveys the message you intend to send to the reader. The report is limited to 5
pages, not including cover page and reference list, using the provided formatting guidelines (see
provided task/question report template in Word or LaTeX format). If this limit is violated or if the
figures/tables are unclear (e.g. not interpretable due to poor quality), the report will be sent back
for immediate revision. Thus, try to present the content as compact and clear as possible. The
report should be uploaded on the designated assignment folder on the LISAM course homepage.

MATLAB script. In addition, you should also provide a working MATLAB script that solves
Example 5.4 in Versteeg and Malalasekera (2007, p. 159–162) and reproduce Fig. 5.19 (p. 162).
The code should consist of one single file (besides the provided Gauss-Seidel function) and work
by simply running the script. The working script should be uploaded to the designated assignment
folder on the LISAM course homepage.

Grading
Pass (2P). To be able to pass this assignment the written report needs to be performed appro-
priately with clear and concise response to the related tasks and questions. Each part should be
addressed with enough details (e.g. in form of relevant theory, results, and discussion) such that
we are convinced that you have satisfied the intended learning outcomes of this assignment. If the
report lack substantial parts or entail major errors and/or incorrect statements, the report will be
sent back for revision. In addition, the MATLAB script should produce the correct output.

Well done (3P). In additions to above, the majority of the written report should be well executed,
with clear and appropriate presentation to each task and question. If the report lack substantial
parts or manifests major errors and/or incorrect statements, this grade will not be rewarded upon
revision.

6
References

References

Versteeg, H.K., Malalasekera, W., 2007. An introduction to computational fluid dynamics: the
finite volume method. Pearson Education.

7
Appendix
Part III
Case Setup

• Setup the geometry in DesignModeler. As CFX only can handle 3D domains, one can mimic
this by only considering the geometry as a thin plate. For simplicity, set the dimensions to
1 x 1 x 0.01 m
• Use symmetry boundary conditions on the 1 x 1 m (large) faces
• In CFX, a transport equation for an arbitrary scalar can be added by:
– Expressions, Functions and Variables>Additional Variables>right-click>Insert>Additional
Variable>Apply with the default settings (specific, unit-less, scalar)
– Default Domain>Fluid Models>Additional Variable 1>Option>Transport Equation
– To not account for any physical diffusion, set the Kinematic Diffusivity to zero
• More information for about additional variables and transport equations can be found in the
CFX-Solver Modeling Guide, which can be accessed via the ANSYS Help manual using the
path: // Modeling Guide // 1. Basic Capabilities Modeling // 1.2. Physical Models //
1.2.12. Additional Variables
• For the outlet boundaries, simply use Outlet with relative static pressure equal to zero (i.e.
1 atm). Opening type boundary condition may also work in these cases

Change numerical schemes

• For the flow equations you don’t need to change the default numerical scheme (i.e. keep the
High Resolution scheme). The objective is to evaluate the schemes impact on the transport
of the additional variable and not on the flow itself (which here is known).
• Make sure set the appropriate advection (convection) scheme for the transport equation.
This can manually be done by: Solver Control>Equation Class Settings>Equation Class
Settings>Additional Variable>Advection Scheme

Note! The QUICK scheme is not accessible in the manual list in CFX. This scheme can be
activated by:
– Right-click>Solver Control>Edit in Command Editor>Equation class: av>Advection
Scheme:>Option = Upwind (in this example)
– Change the current scheme name ”Upwind” to ”QUICK” and confirm by pressing Pro-
cess
– Open the Command Editor again and make sure that the correct settings have been
applied

• In order to minimize the convergence impact on the solution, set a strict convergence criteria
for the equation residuals. The criteria for the Additional Variable can be found at:

8
 – Solver Control>Equation Class Settings>Additional Variable>Advection scheme>Residual
Target
• The turbulence model can be deactivated by:
– Default Domain>Fluid Models>Turbulence>Option>None (Laminar)
• Information about the specific schemes in CFX can be found in the ANSYS Help manual using
the path: // Theory Guide // 11. Discretization and Solution Theory // 11.1. Numerical
Discretization // 11.1.1. Discretization of the Governing Equations // Modeling Guide //
15. Advice on Flow Modeling // 15.5. Advection Scheme Selection (Step forward to view
15.5.1-15.5.4)

Post-processing

• To only get the nodal values across a line in CFD-post, choose Line Type>Cut (These points
can be highlighted by Location>Point Cloud> at the specific line). Note, be aware of using
interpolated values as the result could be misleading for coarse number of cells
• Surface plots could also be a good way to highlight the overall scalar distribution. In CFD-
Post you can Insert>Location>Plane (or simply use one of the symmetry planes), and change
the Color to variable
• Insert>Chart to plot variables (xy-plots) against each other, e.g. the Additional variable
versus the line diagonal length (L). This is a good way to preview results prior to MATLAB
• A new variable can be computed for the length of the X-X diagonal line by Variable>Solution>Right-
click>New...>Name = L>Expression = sqrt(xˆ2+(y-1 [m])ˆ2)