UFVM V1.0 U FVM User Guide
User Manual:
Open the PDF directly: View PDF 
.
Page Count: 56
| Download | |
| Open PDF In Browser | View PDF |
9/6/2017
uFVM v1.0
CFD Academic Tool developed in
MATLAB® – User Manual
Computational Mechanics Lab
AMERICAN UNIVERSITY OF BEIRUT
1 CONTENTS
2
Introduction ............................................................................................................................. 3
2.1
2.1.1
About The Manual ................................................................................................... 3
2.1.2
Terms of Use ............................................................................................................ 3
2.2
3
4
Using The Manual ............................................................................................................ 3
About The Code Development and Authors ................................................................... 3
2.2.1
What’s uFVM? ......................................................................................................... 3
2.2.2
The Place and Time of Development....................................................................... 3
2.2.3
Developers ............................................................................................................... 3
2.3
Range of Applicability ...................................................................................................... 4
2.4
A Glance into uFVM’s Discretization and Solution Methods .......................................... 4
Structure of The Code ............................................................................................................. 6
3.1
Source Files ...................................................................................................................... 6
3.2
Tutorials ........................................................................................................................... 6
Introductory Test Case ............................................................................................................ 7
4.1
Main File: Running a Case ............................................................................................... 7
4.2
Example ........................................................................................................................... 7
4.2.1
Header ..................................................................................................................... 7
4.2.2
Setup Solver Class .................................................................................................... 8
4.2.3
Read OpenFOAM Files ............................................................................................. 8
4.2.4
Setup Time ............................................................................................................... 9
4.2.5
Setup Equations ....................................................................................................... 9
4.2.6
Initialize Case ........................................................................................................... 9
4.2.7
Running the Case ..................................................................................................... 9
5
Reading and Plotting OpenFOAM Mesh................................................................................ 10
6
An Introductory Presentation of uFVM Functions ................................................................ 15
6.1
Application Class............................................................................................................ 15
6.2
Reading OpenFOAM Files and Translation to uFVM ..................................................... 15
6.3
Fixing Time Settings ....................................................................................................... 16
6.4
Interpreting the Model/Equations ................................................................................ 16
6.5
Running the Case ........................................................................................................... 20
6.5.1
1
Time Loop/Convergence Loop............................................................................... 20
6.5.2
6.5.2.1
Transient Term .................................................................................................. 23
6.5.2.2
Convection Term ............................................................................................... 24
6.5.2.3
Diffusion Term ................................................................................................... 24
6.5.2.4
Source Term....................................................................................................... 25
6.5.2.5
‘mdot_f’ Term .................................................................................................... 26
6.5.2.6
Pressure Correction Treatment for Free Surface Flow...................................... 28
6.5.2.7
False Transience ................................................................................................ 28
6.5.2.8
Gradient Computation....................................................................................... 29
6.5.2.9
Implicit Under-Relaxation.................................................................................. 31
6.5.2.10
Residual Form of the Equation ...................................................................... 32
6.5.2.11
Residual Computation ................................................................................... 32
6.5.2.12
Assembling Fluxes to Global Assembly Matrix .............................................. 34
6.5.3
2
Solving the Equation Algebraic System ................................................................. 39
6.5.3.1
SOR Solver ......................................................................................................... 40
6.5.3.2
ILU Solver ........................................................................................................... 40
6.5.3.3
AMG Linear Solver ............................................................................................. 40
6.5.4
7
Assembling Equation Terms .................................................................................. 21
Correcting Equation Solution ................................................................................ 41
6.6
Convergence .................................................................................................................. 43
6.7
Post-processing ............................................................................................................. 44
Tutorials ................................................................................................................................. 45
7.1
Basic ............................................................................................................................... 45
7.2
Incompressible .............................................................................................................. 48
7.3
Heat Transfer ................................................................................................................. 53
2 INTRODUCTION
2.1 USING THE MANUAL
2.1.1 About The Manual
This manual provides an overview about uFVM code. The manual describes with sufficient
details the structure of the code, the range of applicability and who may use it. The way through
which a CFD case is prepared is then described and plenty of tutorials are accordingly provided.
In other words, this manual presents a complete set of instructions for the user to follow in
order to setup a CFD problem.
2.1.2 Terms of Use
uFVM is an academic CFD tool made for learning purposes. It provides a package of libraries and
algorithms that the user can comfortably follow up. Handling, distributing or modifying is fully
permissible; the user has the full permission to add any piece of code or modify an existing one.
2.2 ABOUT THE CODE DEVELOPMENT AND AUTHORS
2.2.1 What’s uFVM?
The name of the code presents an abbreviation letters of the finite volume method (FVM) that
the code is based on. The “u” at the beginning of the name points for a fluid flow. The code is
developed in Matlab® environment because it is assumed that the majority of interested people
are familiar with this environment.
2.2.2 The Place and Time of Development
The code is developed in the computational mechanics lab at the American University of Beirut,
Beirut, Lebanon. The development has started in 2003 and was built and updated gradually
through years. Lots of versions were made each of them had a different structure but
necessarily the same theoretical background.
2.2.3 Developers
The code is a direct accomplishment of the CFD group at the American University of Beirut. The
CFD group is a team of professors, graduate and undergraduate students. Their main objective is
to build computational knowledge and work on plenty of related topics in both tracks
development and application. The team is very familiar with Ansys Fluent® and OpenFOAM®;
they utilize these packages for many purposes. uFVM for example was validated in reference to
these packages. The group’s accomplishments, research topics and published work are posted
on their website https://feaweb.aub.edu.lb/research/cfd. The major contributor to the code is
3
Professor Marwan Darwish1. A co-contributor to the code is Professor Fadl Moukalled2. The
other contributors to the code are Master and PhD students who accomplished their theses and
dissertations from the computational mechanics lab at the American University of Beirut.
2.3 RANGE OF APPLICABILITY
uFVM works for incompressible fluid flows with any type of mesh (structured and unstructured).
For consistency, the code doesn’t include any geometry modeling or meshing capabilities. It
accepts mesh files of OpenFOAM format.
It is worth mentioning that uFVM is a solver which solves the conservation equations (transport
equations) where the user is able to investigate any physical quantity of interest which is
transported by means of a physical phenomenon like convection and diffusion. It is possible also
to set any form of an explicit term into the conservation equation like external forces; these are
treated as source terms. A transient treatment is also included. The domain usually assumes a
fluid flow. However, the code may still apply to solid domains if the user seeks certain transport
quantities in a solid, obviously the temperature distribution.
The code may also handle multi-phase flows. The current distributed version of uFVM doesn’t
include compressible and multi-phase solvers as they are still under construction and revision.
For further information about the status of multi-phase solver, the user may contact any of the
contributors.
It is not the purpose of this code to provide a CFD tool for conducting fluid flow simulations for
heavy/complex applications. There are two issues to raise here. First, this code is made for those
who are mainly learning CFD and/or interested in CFD code development. This code provides a
very useful and helpful means for those people. Second, the user should necessarily realize that
Matlab is a highly user friendly language; this makes it very convenient for learning issues much
more than it is for conducting real life engineering applications. However, this friendly user
specialty had an expense at the computational time; Matlab is slower than other lower level
languages.
2.4 A GLANCE INTO UFVM’S DISCRETIZATION AND SOLUTION METHODS
Only pressure-based methods are available in uFVM with SIMPLE method as the default scheme.
The default convection scheme is the Upwind scheme, but however, different convection
schemes are also available like (SOU, QUICK, SMART, etc.). Gradient computation is based on
the first-order Gauss approximation whether it is cell-based or nodal-based. The “ILU” and
“SOR” solvers are available along with a multi-grid (AMG) solver which utilizes any of the fixed
cycles (V, F, or W Cycle). Under-relaxation factors for any given quantity are treated implicitly in
1
Professor Darwish has gained his PhD from Brunel University in 1985. His major research topics are
computational fluid dynamics, classical fluid mechanics and material sciences. He’s currently acting as a
full time professor at the American University of Beirut.
2
Professor Moukalled has gained his PhD from Louisiana University in 1987. His major research topics are
compressible flows, computational fluid dynamics, modeling energy systems. He’s currently acting as a
full time professor at the American University of Beirut.
4
the equations. All these solution methods and controls are presented in a more detailed
framework in later parts of the manual.
5
3 STRUCTURE OF THE CODE
The uFVM directories are distributed into sources, which are the routines that make up the
code, and tutorials that include the test cases.
3.1 SOURCE FILES
The source code is available in the ‘ufvm/src’ directory. It is distributed to different folders,
‘fvm’ and ‘utilities’. ‘fvm’ contains the finites volume methods, algorithms and solution. ‘utilities’
contains all auxiliary functions.
3.2 TUTORIALS
There are set of tutorials that allows the user to work with ufvm. These tutorials are classified as
basic, incompressible, compressible, multiphase and heatTransfer. Cases that are to be
simulated are to be made in OpenFOAM format. OpenFOAM cases include 3 main directories:
‘0’, ‘system’ and ‘constant’. The ‘0’ directory is where initial and boundary conditions are
specified. The ‘system’ directory includes the solution methods, the finite volume schemes and
the time and write controls of the simulation. The ‘constant’ directory includes the mesh files,
the fluid properties (transport and/or thermophysical), gravity properties and turbulence
properties. For further information, refer to the tutorials.
6
4 INTRODUCTORY TEST CASE
4.1 MAIN FILE: RUNNING A CASE
The code runs from a main script, usually called ‘run’. The ‘run’ file has to be located in an
OpenFOAM case as mentioned earlier. It represents the case study or the problem definition.
This is the only file of importance for the user. The main file contains a set of functions that build
up the model. However, the user has to add the path of uFVM source files by typing the
following command into the command window:
addpath(genpath('~location of uFVM~/uFVM/src'));
4.2 EXAMPLE
The incompressible ‘elbow’ example will be presented here.
%-------------------------------------------------------------------------%
% written by the CFD Group @ AUB, 2017
% contact us at: cfd@aub.edu.lb
%==========================================================================
% Case Description:
%
In this test case a water flow in an elbow is simulated
%-------------------------------------------------------------------------% Setup Case
cfdSetupSolverClass('incompressible');
% Read OpenFOAM Files
cfdReadOpenFoamFiles;
% Setup Time Settings
cfdSetupTime;
% Setup Equations
cfdDefineEquation('U', 'ddt(rho*U) + div(rho*U*U) = laplacian(mu*U) + div(mu*transp(grad(U)))
- grad(p) + rho*g'); % Momentum
cfdDefineEquation('p', 'div(U) = 0'); % Continuity
cfdDefineEquation('T', 'ddt(rho*Cp*T) + div(rho*Cp*U*T) = laplacian(k*T)'); % Energy
% Initialize case
cfdInitializeCase;
% Run case
cfdRunCase;
Listing 1 - Run file of the case named ‘elbow’ in the incompressible tutorials directory
4.2.1 Header
The script starts by few words expressing a header and providing a summary of the case:
%-------------------------------------------------------------------------%
% written by the CFD Group @ AUB, 2017
% contact us at: cfd@aub.edu.lb
%==========================================================================
% Case Description:
%
In this test case a water flow in an elbow is simulated
%--------------------------------------------------------------------------
7
4.2.2 Setup Solver Class
The user has to set the class of the application. The application can be one of the following:
basic-incompressible-compressible-multiphase-heatTransfer. These are the standard
OpenFOAM tutorial classes.
4.2.3 Read OpenFOAM Files
OpenFOAM files are imported at this level to Matlab and stored in a global variable called
‘Domain’. In fact, all OpenFOAM input information are stored in a structure called foam within
the global Domain variable. The table below presents the settings that the user must or may
include in the OpenFOAM files:
Option/Setting Name
residualControl
relaxationFactors
ddtSchemes
Description
The equations’ tolerances
(convergence criteria)
The under-relaxation factors
to be applied to the
equations
The transient scheme (steady
or transient). The options
are:
Location
system/fvSolution
system/fvSolution
system/fvSchemes
steadyState
Euler
divSchemes
The convection scheme. The
options are:
system/fvSchemes
Gauss upwind
Gauss linear
startFrom
This tells the solver from
where to start. The options
are:
system/controlDict
firstTime
startTime
latestTime
startTime
stopAt
The initial time of the
simulation
This tells the solver until
when to stop. The options
are:
system/controlDict
system/controlDict
endTime
endTime
deltaT
writeControl
The final time of the
simulation
The time step of the
simulation
The criterion that the solver
writes the results based on.
The options are:
system/controlDict
system/controlDict
system/controlDict
timeStep
writeInterval
8
Write results every specified
number of time steps (if the
system/controlDict
writeControl is set to
timeStep)
4.2.4 Setup Time
The function cfdSetupTime sets time controls and other conditions for the simulation. Other
conditions of the simulation are related to multiphase flow model which will not be considered
in the current uFVM version. The following table presents the valid arguments of the function:
4.2.5
Setup Equations
The user has to define the equations that are to be solved. The equations that are being solved
in the above example are described below.
Momentum:
𝜕(𝜌𝒗)
+ ∇. (𝜌𝒗𝒗) = 𝜇∇2 𝒗 + ∇. {𝜇(∇𝒗)𝑇 } − ∇𝑝 + 𝜌𝒈
𝜕𝑡
Continuity:
∇. 𝒗 = 0
Energy:
𝜕(𝜌𝑐𝑝 𝑇)
+ ∇. (𝜌𝑐𝑝 𝒗𝑇) = 𝑘∇2 𝑇
𝜕𝑡
4.2.6 Initialize Case
The case is initialized after that by evaluating all the fields with their corresponding values on
the mesh. For example, the field variables (variables to be solved for), have initial values that are
imported from the OpenFOAM ‘0’ directory. The properties are evaluated based on their values
dictated by the ‘transportProperties’ file or ‘thermophysicalProperties’ file in case of a
compressible flow; these are found in the ‘constant’ directory.
4.2.7 Running the Case
The case is made to run at this level, iterating over the equations until convergence. The
convergence criteria can be found in the ‘fvSolution’ file in the ‘system’ directory within
OpenFOAM folders.
9
5 READING AND PLOTTING OPENFOAM MESH
When in the case, a user may read the OpenFOAM mesh and visualize it before running the
case. The following commands are all what the user needs to do so:
To read the mesh, the function cfdReadPolyMesh has to be called. The uFVM will go to the
‘constant/polyMesh’ directory, read the files there, and store the mesh information in the
database. To plot the mesh, the user has to use the function cfdPlotMesh; the mesh will then
be plotted:
Figure 1 - Plotted mesh
To witness the details of the read mesh, the user has to call for the mesh from the database by
writing:
mesh = cfdGetMesh
10
The mesh will be printed out such as:
mesh =
struct with fields:
nodes: [1×4851 struct]
numberOfNodes: 4851
caseDirectory: 'constant/polyMesh'
numberOfFaces: 12800
numberOfElements: 4000
faces: [1×12800 struct]
numberOfInteriorFaces: 11200
boundaries: [1×3 struct]
numberOfBoundaries: 3
numberOfPatches: 3
elements: [1×4000 struct]
numberOfBElements: 1600
numberOfBFaces: 1600
cconn: {4000×1 cell}
csize: [4000×1 double]
The information printed out above are structures and values that store the quantitative and
qualitative details of the mesh. For instance, see the number of elements (numberOfElements:
4000) and the number of patches (numberOfPatches: 3). The information of any element or face
can be witnessed by accessing it from the corresponding structure. If I want to see the details of
element number 100 in the mesh, I can simply write mesh.elements(100)and the output is:
11
>> theMesh.elements(100)
ans =
struct with fields:
index: 100
iNeighbours: [80 99 120 300]
iFaces: [235 291 294 295 12205 12595]
iNodes: [104 335 336 105 125 356 357 126]
volume: 0.1250
faceSign: [-1 -1 1 1 1 1]
numberOfNeighbours: 4
OldVolume: 0.1250
centroid: [3×1 double]
The previous commands can also be applied to faces and nodes. The user may also visualize the
elements, faces and patches. In order to plot patch #1 for instance, the function
cfdPlotPatches(1) has to be called. The output figure is:
12
Figure 2 - Plotted mesh and patch number 1
Elements and faces may also be visualized using the functions cfdPlotElements and
cfdPlotFaces. If more than an entity are to be plotted, say 3 entities of indecis 1, 2 and 3, and
considering plotting the elements, the user has to call cfdPlotElements([1 2 3]). The
output figure is:
13
Figure 3 - Visualizing elements 1, 2 and 3
14
6 AN INTRODUCTORY PRESENTATION OF UFVM FUNCTIONS
6.1 APPLICATION CLASS
Different applications give rise to different ways of treating the fluid flow equations. Thus a user
has to set an application class by choosing any of the following applications (basic,
incompressible, compressible, heatTransfer and multiphase). It is significantly important to set
the right application class because it is the basis of many things inside the code.
6.2 READING OPENFOAM FILES AND TRANSLATION TO UFVM
As mentioned previously, uFVM cases are prepared in OpenFOAM format. This choice was made
because OpenFOAM is a widely used CFD library (C++ based). Preparing an OpenFOAM case is
quite straight forward for each application class. Refer to the tutorials and check the way in
which each application is being prepared.
uFVM’s function cfdReadOpenFoamFiles reads all the files that exist in the 3 standard
OpenFOAM directories ‘0’, ‘constant’ and ‘system’. It searches for specific entries and blocks
within the files and stores them in the uFVM’s data base. See the following table for details
about the read files:
OpenFOAM
Case
Directory
fvSolution
'system'
Directory
block: SIMPLE/PISO/ALGORITHM
block: relaxationFactors
fvSchemes
controlDict
polyMesh
'constant'
Directory
transportProperties/thermophysicalProperties
turbulenceProperties
g
U
'0'
Directory
block: solvers
p
T
block: ddtSchemes
block: divSchemes
block: laplcaianSchemes
Other Finite Volume Operators
key: startFrom
key: startTime
key: stopAt
key: endTime
other control options
boundary
points
faces
owner
neighbour
Others
Collect all properties and their attributes (values,
dimensions, models, expressions)
key: simulationType i.e. RAS
key: RASModel
key: turbulence
key: printCoeffs
Collect gravitational acceleration vale i.e. (0 0 -9.8)
class
internalField
boundaryField
class
internalField
boundaryField
class
internalField
boundaryField
Other fields
The above FOAM details are imported and stored in the data base ‘Domain’ under the structure
‘foam’ as such:
15
6.3 FIXING TIME SETTINGS
Of the imported control settings, the time controls which exist in the controlDict structure in the
global Domain variable, are stored again in a more useful way in another structure called time.
The structure ‘time’ contains the following info for the arbitrary case of elbow:
>> Domain.time
ans =
struct with fields:
startFrom: 'latestTime'
startTime: 0
stopAt: 'endTime'
endTime: 10000
deltaT: 10
fdt: 1.0000e+09
type: 'Transient'
The above settings will be used later on to set the start and end of the time loop or convergence
loop. If the case is transient, ‘deltaT’ in the above structure is the time step.
6.4 INTERPRETING THE MODEL/EQUATIONS
The equations in uFVM are interpreted in such a way that each term of the equation has its own
attributes. Generally, any conservation equation has the following form:
16
𝜕(𝜌𝜙)
+⏟
∇. (𝜌𝒗𝜙) = ⏟
∇. (𝛤 𝜙 ∇𝜙) + ⏟
𝑄𝜙
⏟𝜕𝑡
𝑆𝑜𝑢𝑟𝑐𝑒
𝑇𝑟𝑎𝑛𝑠𝑖𝑒𝑛𝑡
𝐶𝑜𝑛𝑣𝑒𝑐𝑡𝑖𝑜𝑛
𝐷𝑖𝑓𝑓𝑢𝑠𝑖𝑜𝑛
Note: The diffusion term ∇. (𝛤 𝜙 ∇𝜙) can only be written as 𝛤∇2 𝜙 when the diffusion coefficient
𝛤 is constant. However, when implementing the diffusion term in uFVM, it is always written as
laplacian(tau*phi) whether 𝛤 is constant or not.
The general conservation equation above is expressed in uFVM as follows:
cfdDefineEquation('phi', 'ddt(rho*phi) + div(rho*U*phi) = laplacian(tau*phi) + Q');
𝑄 in the above implementation has to be an explicit expression that includes fields and/or
properties. Inside uFVM, each of the above terms is interpreted in a different way. Each term
consists of an explicit part and an implicit one; the explicit part is the part which holds the
previous iteration value, and the implicit part is kept as an unknown to be assembled with other
implicit parts of other terms in a global matrix. The explicit part is roughly all the parameters
that are multiplied with 𝜙 and the implicit part is 𝜙. In the case of a general conservation
equation, the explicit part of the transient term
𝜕(𝜌𝜙)
𝜕𝑡
is 𝜌 and the implicit part is 𝜙. In addition,
the explicit part of the transient term, is called the rho field. For the convection and diffusion
terms, the explicit parts are 𝜌𝒗 and 𝛤 respectively, and the names of the latter explicit fields are
the psi field and the gamma field respectively. The source term is treated explicitly; it is
evaluated based on available field values.
The table below summarizes the treatment of the terms of a general conservation equation:
𝜕(𝜌𝜙)
𝜕𝑡
𝛻. (𝜌𝒗𝜙)
The term type: Transient Term
the rho field: 𝜌
the variable field: 𝜙
The term type: Convection Term
the psi field: 𝜌𝒗
the variable field: 𝜙
𝛻. (𝛤 𝜙 𝛻𝜙)
The term type: DiffusionTerm
the gamma field: 𝛤 𝜙
the variable name: 𝜙
𝑄𝜙
17
The term type: Source term
Example1:
Take the momentum equation in the elbow case as an example:
𝜕(𝜌𝒗)
+ ∇. (𝜌𝒗𝒗) = 𝜇∇2 𝒗 + ∇. {𝜇(∇𝒗)𝑇 } − ∇𝑝 + 𝜌𝒈
𝜕𝑡
It is written in uFVM in the following form:
'ddt(rho*U) + div(rho*U*U) = laplacian(mu*U) + div(mu*transp(grad(U))) - grad(p) + rho*g'
The momentum equation above includes 6 terms, each of them is interpreted in uFVM as shown
here:
ddt(rho*U)
The term type: Transient Term
Finite Volume Operation: Implicit
the rho field: rho
the variable field: U
term sign: 1
div(rho*U*U)
The term type: Convection Term
Finite Volume Operation: Implicit
the psi field: rho*U
the variable field: U
term sign: 1
laplacian(mu*U)
The term type: Stress Term
Finite Volume Operation: Implicit
the gamma field: mu
the variable name: U
term sign = -1
div(mu*transp(grad(U)))
The term type: Source term
Finite Volume Operation: Explicit
explicit operators available: grad
term sign: -1
grad(p)
The term type: Source term
Finite Volume Operation: Explicit
explicit operators available: grad
term sign: 1
rho*g
The term type: Source term
Finite Volume Operation: Explicit
explicit operators available: none
term sign: -1
Example2:
Take the Energy equation in the elbow case as another example:
𝜕(𝜌𝑐𝑝 𝑇)
+ ∇. (𝜌𝑐𝑝 𝒗𝑇) = 𝑘∇2 𝑇
𝜕𝑡
It is written in uFVM in the following form:
'ddt(rho*Cp*T) + div(rho*Cp*U*T) = laplacian(k*T)'
The Energy equation above includes 3 terms, each of them is interpreted in uFVM as shown
here:
18
ddt(rho*Cp*T)
The term type: Transient Term
Finite Volume Operation: Implicit
the rho field: rho*Cp
the variable field: T
term sign: 1
div(rho*Cp*U*T)
The term type: Convection Term
Finite Volume Operation: Implicit
the psi field: rho*Cp*U
the variable field: T
term sign: 1
laplacian(k*T)
The term type: Stress Term
Finite Volume Operation: Implicit
the gamma field: k
the variable name: T
term sign = -1
However, a special case arises for the continuity equation which is usually called also the
pressure equation ‘p’. The ‘p’ equation is treated in a different way. The pressure equation for
an incompressible application class is the same for any fluid flow problem within the
incompressible category. Thus, the user may ignore the pressure equation expression keeping
only the equation name (‘p’). The reason that the pressure equation expression is included in
the Matlab listing above and we are repeating it here is just for showing the exact models that
this code is solving:
cfdDefineEquation('p', 'div(U) = 0');
While it is sufficient to write:
cfdDefineEquation('p');
Automatically, all the terms of the continuity equation are treated in a single term assembly
process, this term is called ‘mdot’. For incompressible fluid flows, the continuity equation stated
above is theoretically converted into another equation called the pressure equation, and that’s
why the name ‘p’ corresponds to the continuity equation.
A summary of how the pressure equation is derived from the continuity equation is presented
here:
The continuity equation for an incompressible flow:
∇. 𝒗 = 𝟎
19
→
∑
𝑚̇𝑓 = 0
𝑓~𝑛𝑏(𝐶)
Dividing 𝑚̇𝑓 into a predicted and correction components:
𝑚̇𝑓 = 𝑚̇𝑓∗ + 𝑚̇𝑓′
So,
∑ 𝑚̇𝑓′ +
𝑓~𝑛𝑏(𝐶)
∑ 𝑚̇𝑓∗ = 0
𝑓~𝑛𝑏(𝐶)
The Rhie-Chow interpolation for 𝑚̇𝑓∗ is
̅̅̅̅̅̅̅
(𝑛)
(𝑛)
𝑚̇𝑓∗ = 𝜌𝒇 ̅̅̅
𝒗𝑓∗ . 𝑺𝑓 − 𝜌𝒇 ̅̅̅̅
𝑫𝑓𝒗 (∇𝑝𝑓 − ∇𝑝𝑓 ) . 𝑺𝑓
and for 𝑚̇𝑓′ is
𝑚̇𝑓′ = −𝜌𝒇 ̅̅̅̅
𝑫𝑓𝒗 ∇𝑝𝑓′ . 𝑺𝑓
Therefore, the equation of the incompressible pressure equation is:
∑ −𝜌𝑓 ̅̅̅̅
𝑫𝑓𝒗 (∇𝑝′ 𝑓 ). 𝑆𝑓 + ∑ 𝑚̇𝑓∗ = 0
𝑓~𝑛𝑏(𝐶)
𝑓~𝑛𝑏(𝐶)
⏟
𝑚𝑑𝑜𝑡_𝑓
For a compressible continuity equation, the resulting pressure equation is a follows (you may
find the proof in the book):
𝑚̇𝑓∗
𝑉𝐶 𝐶𝜌 ′
(𝜌𝐶∗ − 𝜌𝐶° )
𝑝𝐶 + ∑ 𝐶𝜌,𝑓 ( ∗ ) 𝑝𝑓′ + ∑ −𝜌𝑓 ̅̅̅̅
𝑫𝑓𝒗 (∇𝑝′ 𝑓 ). 𝑆𝑓 +
+ ∑ 𝑚̇𝑓∗ = 0
∆𝑡
𝜌𝑓
∆𝑡
𝑓~𝑛𝑏(𝐶)
𝑓~𝑛𝑏(𝐶)
𝑓~𝑛𝑏(𝐶)
⏟
𝑚𝑑𝑜𝑡_𝑓
6.5 RUNNING THE CASE
6.5.1 Time Loop/Convergence Loop
The time loop and the convergence loop are located in the cfdRunCase where the time
settings are utilized. The following listing shows the time and convergence loops for a transient
simulation:
20
time = cfdGetTime;
startTime = time.startTime;
endTime = time.endTime;
deltaT = time.deltaT;
% Time Loop: Loop until the final time
timeIter = 1;
cumulativeIter = 1;
for t=startTime:deltaT:endTime
% Time settings
currentTime = t + deltaT;
cfdSetCurrentTime(currentTime);
% Update previous time step fields
cfdTransientUpdate;
cfdPrintCurrentTime(currentTime);
% Convergence Loop: Loop until convergence for the current time step
for iter=1:20
cfdPrintIteration(cumulativeIter);
cfdPrintResidualsHeader;
%
cfdUpdateFields;
%
for iEquation=1:theNumberOfEquations
% Assemble the current equation and correct it
[rmsResidual, maxResidual, lsResBefore, lsResAfter] = cfdAssembleAndCorrectEquation(theEquationNames{iEquation});
% Print the equation residuals
cfdPrintResiduals(cfdGetBaseName(theEquationNames{iEquation}),rmsResidual,maxResidual,lsResBefore,lsResAfter);
% If multigrid solver is assigned, print the AMG solver
% settings
theEquation = cfdGetModel(theEquationNames{iEquation});
isMultigrid = theEquation.multigrid.isActive;
if isMultigrid
cfdPrintLinearSolver(theEquationNames{iEquation});
if iEquationnumberOfInteriorFaces)))
localElementIndices = theNode.iElements;
for iElement = localElementIndices
theElement = fvmElements(iElement);
C = theElement.centroid;
d = cfdMagnitude(N-C);
localPhi = phi(iElement);
localPhiNode = localPhiNode + localPhi/d;
localInverseDistanceSum = localInverseDistanceSum + 1/d;
end
else
localBFacesIndices = theNode.iFaces(theNode.iFaces>numberOfInteriorFaces);
for iBFace = localBFacesIndices
theFace = fvmFaces(iBFace);
C = theFace.centroid;
iBElement = numberOfElements+(iBFace-numberOfInteriorFaces);
d = cfdMagnitude(N-C);
localPhi = phi(iBElement);
localPhiNode = localPhiNode + localPhi/d;
localInverseDistanceSum = localInverseDistanceSum + 1/d;
end
end
localPhiNode = localPhiNode/localInverseDistanceSum;
%
%
phi_n(iNode) = localPhiNode;
%
end
Listing 13 - Interpolating the cell values to the nodes
Then the following code calculates the node values to the faces, and can be found in the
function cfdInterpolateFromNodesToFaces:
30
for iFace=1:numberOfFaces
theFace = fvmFaces(iFace);
iNodes = theFace.iNodes;
C = theFace.centroid;
localSumOfInverseDistance=0;
localPhi=0;
for iNode=iNodes
theNode = fvmNodes(iNode);
N=theNode.centroid;
d=cfdMagnitude(C-N);
localPhi = localPhi + phiNodes(iNode)/d;
localSumOfInverseDistance = localSumOfInverseDistance+1/d;
end
localPhi = localPhi/localSumOfInverseDistance;
phi_f(iFace) = localPhi;
end
Listing 14 - Interpolating the node values to the faces
Finally, we calculate the gradient according to the Green-Gauss method which can be found in
the function cfdComputeGradientNodal:
for iFace=1:theNumberOfInteriorFaces
%
theFace = fvmFaces(iFace);
%
iElement1 = theFace.iOwner;
iElement2 = theFace.iNeighbour;
%
Sf = theFace.Sf;
%
%
phiGrad(:,iElement1) = phiGrad(:,iElement1) + phi_f(iFace)*Sf;
phiGrad(:,iElement2) = phiGrad(:,iElement2) - phi_f(iFace)*Sf;
end
%=====================================================
% BOUNDARY FACES contribution to gradient
%=====================================================
for iBPatch=1:theNumberOfBElements
%
iBFace = theNumberOfInteriorFaces+iBPatch;
iBElement = theNumberOfElements+iBPatch;
theFace = fvmFaces(iBFace);
%
iElement1 = theFace.iOwner;
%
Sb = theFace.Sf;
phi_b = phi(iBElement);
%
phiGrad(:,iElement1) = phiGrad(:,iElement1) + phi_b*Sb;
end
%----------------------------------------------------% Get Average Gradient by dividing with element volume
%----------------------------------------------------for iElement =1:theNumberOfElements
theElement = fvmElements(iElement);
phiGrad(:,iElement) = phiGrad(:,iElement)/theElement.volume;
end
Listing 15 - Calculation of the Green-Gauss gradient based on the face values calculated from node values
6.5.2.9 Implicit Under-Relaxation
As suggested by Patankar’s approach, a relaxation factor 𝜆𝜙 is introduced into the algebraic
system.
31
𝑎𝐶
𝜙 +
𝜆𝜙 𝐶
∑
𝑎𝐹 𝜙𝐹 = 𝑏𝐶 +
𝐹~𝑁𝐵(𝐶)
(1 − 𝜆𝜙 )
𝑎𝐶 𝜙𝐶∗
𝜆𝜙
However, for an equation in the correction form, an implicit under-relaxation is made as follows:
𝑎𝐶 ′
𝜙 +
𝜆𝜙 𝐶
∑
𝑎𝐹 𝜙𝐹′ = 𝑏𝐶 − (𝑎𝐶 𝜙𝐶∗ +
𝐹~𝑁𝐵(𝐶)
∑
𝑎𝐹 𝜙𝐹∗ )
𝐹~𝑁𝐵(𝐶)
In uFVM, it is done in the function cfdApplyURF:
theEquation = cfdGetModel(theEquationName);
urf = theEquation.urf;
theCoefficients = cfdGetCoefficients;
theCoefficients.ac = theCoefficients.ac/urf;
cfdSetCoefficients(theCoefficients);
Listing 16 - Introducing under-relaxation to the algebraic equation
6.5.2.10 Residual Form of the Equation
In fact, uFVM assembles a residual or correction form of the equation instead of the direct
equation. This correction makes use of the previous iteration values of 𝜙𝐶 such that:
𝜙𝐶 = 𝜙𝐶∗ + 𝜙𝐶′
Satisfying in the standard equation form:
𝑎𝐶 (𝜙𝐶∗ + 𝜙𝐶′ ) +
∑
𝑎𝐹 (𝜙𝐹∗ + 𝜙𝐹′ ) = 𝑏𝐶
𝐹~𝑁𝐵(𝐶)
→ 𝑎𝐶 𝜙𝐶′ +
∑
𝑎𝐹 𝜙𝐹′ = 𝑏𝐶 − (𝑎𝐶 𝜙𝐶∗ +
𝐹~𝑁𝐵(𝐶)
∑
𝑎𝐹 𝜙𝐹∗ )
𝐹~𝑁𝐵(𝐶)
Once the algebraic system is calculated and 𝜙𝐶′ is determined, the exact value 𝜙𝐶 is updated as:
𝜙𝐶 = 𝜙𝐶∗ + 𝜙𝐶′
6.5.2.11 Residual Computation
The residuals are criteria upon which the user decides to consider the results correct enough.
𝜙
The residual 𝑅𝐶 of the algebraic equation at element 𝐶 is calculated as:
𝜙
𝑅𝐶 = |𝑏𝐶 − (𝑎𝐶 𝜙𝐶 +
∑
𝐹~𝑁𝐵(𝐶)
The maximum residual over the cells is calculated as:
32
𝑎𝐹 𝜙𝐹 )|
𝜙
𝜙
𝑅𝐶,𝑚𝑎𝑥 = max (𝑅𝐶 )
And the root-mean-squared of the residuals over the cells is:
2
𝑅𝐶,𝑟𝑚𝑠 = √
∑𝐶~𝑎𝑙𝑙 𝑐𝑒𝑙𝑙𝑠 (𝑅𝐶𝜙 )
𝜙
𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑒𝑙𝑒𝑚𝑒𝑛𝑡𝑠
Normalized residuals provide better insight of the convergence. They are calculated as follows:
𝜙
𝜙
𝑅𝐶,𝑠𝑐𝑎𝑙𝑒𝑑 =
𝑅𝐶
max(𝑎𝐶 𝜙𝐶 )
After that, you we calculate the maximum and root-mean square scaled residual.
However, since in uFVM the equations are assembled in the residual (correction) form as
mentioned earlier, the residual at element 𝐶 is simply:
𝑎𝐶 𝜙𝐶′ +
∑
𝑎𝐹 𝜙𝐹′ = 𝑏𝐶 − (𝑎𝐶 𝜙𝐶∗ +
𝐹~𝑁𝐵(𝐶)
∑
𝑎𝐹 𝜙𝐹∗ ) =
𝐹~𝑁𝐵(𝐶)
So,
𝜙
𝑅𝐶 = |𝑏𝐶,𝑟𝑒𝑠 |
The following listing presents the residual calculation; it is available in the function
cfdComputeNormalizedResidual:
% Loop over elements and calculate residual at each element
Rc = abs(bc);
% Residuals. Calculate for convenience. Otherwise, they are not used
Rc_max = max(Rc);
Rc_rms = sqrt(sum(Rc.^2)/theNumberOfElements);
% Get phi scale from data base.
% phi_scale = max(abs(phi)). And if phi is zero, phi_scale is set to 1
phi_scale = cfdGetScale(theEquationUserName);
% Normalized Residuals
Rc_scaled = Rc / (max(abs(ac))*phi_scale);
Rc_max_scaled = max(Rc_scaled);
Rc_rms_scaled = sqrt(sum(Rc_scaled.^2)/theNumberOfElements);
MAXResidual = Rc_max_scaled;
RMSResidual = Rc_rms_scaled;
Listing 17 - Residual calculation
A special case arises with the pressure correction equation as it is by default in the correction
form as shown here:
33
𝑎𝐶 𝑝𝐶′ +
𝑝
𝑎𝐹 𝑝𝐹′ = 𝑏𝐶
∑
𝐹~𝑁𝐵(𝐶)
However, the residual of this equation is determined as a continuity criterion
∇. 𝒗 = 0
→
∑
𝑚̇𝑓 = 0
𝑓~𝑛𝑏(𝐶)
Thus, a quantity called 𝐷𝑖𝑣𝑒𝑓𝑓𝑒𝑐𝑡𝑖𝑣𝑒 = ∑𝑓~𝑛𝑏(𝐶) 𝑚̇𝑓 is calculated at each iteration to judge the
convergence of the continuity equation, so, we have for the continuity equation:
𝜙
𝑅𝐶 = 𝐷𝑖𝑣𝑒𝑓𝑓𝑒𝑐𝑡𝑖𝑣𝑒
The corresponding implementation can be found in the function:
cfdComputeEffectiveDivergence:
% Interior Faces Contribution
theNumberOfInteriorFaces = cfdGetNumberOfInteriorFaces;
iFaces = 1:theNumberOfInteriorFaces;
owners = [theMesh.faces(iFaces).iOwner]';
neighbours = [theMesh.faces(iFaces).iNeighbour]';
for iFace=1:theNumberOfInteriorFaces
iOwner = owners(iFace);
iNeighbour = neighbours(iFace);
%
effDiv(iOwner)
= effDiv(iOwner)
+ mdot_f(iFace);
effDiv(iNeighbour) = effDiv(iNeighbour) - mdot_f(iFace);
end
% Boundary Faces Contribution
theNumberOfPatches = cfdGetNumberOfPatches;
for iPatch=1:theNumberOfPatches
theBoundary = theMesh.boundaries(iPatch);
numberOfBFaces = theBoundary.numberOfBFaces;
% cfdGetBoundaryIndex
iFaceStart = theBoundary.startFace;
iFaceEnd = iFaceStart+numberOfBFaces-1;
iBFaces = iFaceStart:iFaceEnd;
owners = [theMesh.faces(iBFaces).iOwner]';
mdot_b = theMdotField.phi(iBFaces);
for iBFace=1:numberOfBFaces
iOwner = owners(iBFace);
effDiv(iOwner) = effDiv(iOwner) + mdot_b(iBFace);
end
end
Listing 18 - Calculating the effective divergence as a residual criterion for the continuity equation
6.5.2.12 Assembling Fluxes to Global Assembly Matrix
After calculating the face or/and element fluxes from each term, these fluxes are to be summed
to construct the algebraic system coefficients.
6.5.2.12.1 Algebraic Systems Representation
The algebraic system 𝐴𝜙 ′ = 𝑏 is constructed:
𝑎11
[ ⋮
𝑎𝑁1
34
… 𝑎1𝑁 𝜙1′
𝑏1
⋱
⋮ ][ ⋮ ] = [ ⋮ ]
′
… 𝑎𝑁𝑁 𝜙𝑁
𝑏𝑁
The coefficient matrix 𝐴 is a highly sparse matrix, so it is stored as an array for the diagonal
coefficients in addition of a data structure containing the non-zeros at each row (off-diagonal
coefficients). The diagonal array is named in uFVM as ac and the off-diagonal coefficient data
structure as anb.
6.5.2.12.2 Assembling Algebraic System coefficients
If the discretized term gives rise to element fluxes like the transient and source terms, the fluxes
𝑓𝑙𝑢𝑥𝐶, 𝑓𝑙𝑢𝑥𝐶 ° , 𝑓𝑙𝑢𝑥𝑉 and 𝑓𝑙𝑢𝑥𝑇 are assembled as follows given that the algebraic equation is
in the residual (correction) form:
𝑎𝐶 = 𝑎𝐶 + 𝑓𝑙𝑢𝑥𝐶
𝑎𝐶° = 𝑎𝐶° + 𝑓𝑙𝑢𝑥𝐶 °
𝑏𝐶 = 𝑏𝐶 − 𝑓𝑙𝑢𝑥𝑇
If the term includes faces fluxes (𝑓𝑙𝑢𝑥𝐶𝑓, 𝑓𝑙𝑢𝑥𝐹𝑓, 𝑓𝑙𝑢𝑥𝑉𝑓 and 𝑓𝑙𝑢𝑥𝑇𝑓) like the convection and
diffusion terms, the assembling is done such as:
𝑎𝐶 = 𝑎𝐶 +
∑ 𝑓𝑙𝑢𝑥𝐶𝑓
𝑓~𝑛𝑏(𝐶)
𝑎𝐹 = 𝑎𝐹 + 𝑓𝑙𝑢𝑥𝐹𝑓
𝑏𝐶 = 𝑏𝐶 −
∑ 𝑓𝑙𝑢𝑥𝑇𝑓
𝑓~𝑛𝑏(𝐶)
Recall that for element fluxes:
𝑓𝑙𝑢𝑥𝑇 = 𝑓𝑙𝑢𝑥𝐶𝜙𝐶∗ + 𝑓𝑙𝑢𝑥𝐶 ° 𝜙𝐶° + 𝑓𝑙𝑢𝑥𝑉
And for face fluxes,
𝑓𝑙𝑢𝑥𝑇𝑓 = 𝑓𝑙𝑢𝑥𝐶𝑓𝜙𝐶∗ + 𝑓𝑙𝑢𝑥𝐹𝑓𝜙𝐹∗ + 𝑓𝑙𝑢𝑥𝑉𝑓
Therefore, the algebraic equation has the general form:
(𝑓𝑙𝑢𝑥𝐶 + ∑ 𝐹𝑙𝑢𝑥𝐶𝑓 ) 𝜙𝐶′ + ∑ 𝐹𝑙𝑢𝑥𝐶𝑓
𝜙𝐹′ = − (𝑓𝑙𝑢𝑥𝑇 + ∑ 𝐹𝑙𝑢𝑥𝑇𝑓)
⏟
𝑎𝐹
𝐹~𝑁𝐵(𝐶)
𝑓~𝑛𝑏(𝐶)
𝑓~𝑛𝑏(𝐶)
⏟
⏟
𝑎𝐶
𝑓~𝑛𝑏(𝐶)
𝑏𝐶
So,
𝑎𝐶 = 𝑓𝑙𝑢𝑥𝐶 +
∑ 𝐹𝑙𝑢𝑥𝐶𝑓
𝑓~𝑛𝑏(𝐶)
𝑎𝐹 = 𝐹𝑙𝑢𝑥𝐶𝑓
35
𝑏𝐶 = − (𝑓𝑙𝑢𝑥𝑇 +
∑ 𝐹𝑙𝑢𝑥𝑇𝑓)
𝑓~𝑛𝑏(𝐶)
The pressure equation again has a special case since it is already in the residual form. The
discretized incompressible pressure correction equation:
( ∑ 𝜌𝑓 𝒟𝑓 ) 𝑝𝐶′ + ∑ ⏟
−𝜌𝑓 𝒟𝑓 𝑝𝐹′ = − ∑ 𝑚̇𝑓∗
𝐹~𝑁𝐵(𝐶)
⏟ 𝑓~𝑛𝑏(𝐶)
𝑎𝐹
⏟𝑓~𝑛𝑏(𝐶)
𝑓~𝑛𝑏(𝐶)
𝑎𝐶
𝑏𝐶
So, the fluxes here are to be calculated as:
1) Diagonal coefficient 𝑎𝐶 :
𝑎𝐶 = 𝑓𝑙𝑢𝑥𝐶 +
∑ 𝐹𝑙𝑢𝑥𝐶𝑓 =
𝑓~𝑛𝑏(𝐶)
∑ 𝜌𝑓 𝒟𝑓
𝑓~𝑛𝑏(𝐶)
→ 𝑓𝑙𝑢𝑥𝐶 = 0 and 𝐹𝑙𝑢𝑥𝐶𝑓 = 𝜌𝑓 𝒟𝑓
2) Off-diagonal coefficients 𝑎𝐹 :
𝑎𝐹 = 𝐹𝑙𝑢𝑥𝐹𝑓 = −𝜌𝑓 𝒟𝑓
3) Right-hand-side:
𝑏𝐶 = − (𝑓𝑙𝑢𝑥𝑇 +
∑
𝐹𝑙𝑢𝑥𝑇𝑓) = − ∑
𝑓~𝑛𝑏(𝐶)
𝑚̇𝑓∗
𝑓~𝑛𝑏(𝐶)
→ 𝑓𝑙𝑢𝑥𝑇 = 0 and 𝐹𝑙𝑢𝑥𝑇𝑓 = 𝑚̇𝑓∗
For the compressible pressure correction equation, the discretized equation is:
(𝜌𝐶∗ − 𝜌°𝐶 )
𝐶𝜌,𝑓
𝐶𝜌,𝑓
𝑉𝐶 𝐶𝜌
+ ∑ ( ∗ ‖𝑚̇𝑓∗ , 0‖) + ∑ 𝜌𝑓∗ 𝒟𝑓 ) 𝑝𝐶′ + ∑ (− ∗ ‖−𝑚̇𝑓∗ , 0‖ − 𝜌𝑓∗ 𝒟𝑓 ) 𝑝𝐹′ = − ∑ 𝑚̇𝑓∗ −
(
∆𝑡
𝜌𝑓
𝜌𝑓
∆𝑡
𝐹~𝑁𝐵(𝐶) ⏟
𝑓~𝑛𝑏(𝐶)
𝑓~𝑛𝑏(𝐶)
⏟ 𝑓~𝑛𝑏(𝐶)
⏟
𝑎𝐶
𝑓~𝑛𝑏(𝐶)
𝑎𝐹
So, the fluxes of the compressible pressure correction equation are calculated as:
1) Diagonal coefficient 𝑎𝐶 :
36
𝑏𝐶
𝑎𝐶 = 𝑓𝑙𝑢𝑥𝐶 +
∑ 𝐹𝑙𝑢𝑥𝐶𝑓 =
𝑉𝐶 𝐶𝜌
𝑓~𝑛𝑏(𝐶)
→ 𝑓𝑙𝑢𝑥𝐶 =
𝑉𝐶 𝐶𝜌
∆𝑡
∆𝑡
+
∑ (
𝑓~𝑛𝑏(𝐶)
and 𝐹𝑙𝑢𝑥𝐶𝑓 =
𝐶𝜌,𝑓
𝜌∗𝑓
𝐶𝜌,𝑓
𝜌∗𝑓
‖𝑚̇ ∗𝑓 , 0‖) + ∑
∗
𝜌𝑓 𝒟𝑓
𝑓~𝑛𝑏(𝐶)
‖𝑚̇ ∗𝑓 , 0‖ + 𝜌∗𝑓 𝒟𝑓
2) Off-diagonal coefficients 𝑎𝐹 :
𝑎𝐹 = 𝐹𝑙𝑢𝑥𝐹𝑓 = − 𝐶𝜌𝜌,𝑓∗ ‖−𝑚̇𝑓∗ , 0‖ − 𝜌𝑓∗ 𝒟𝑓
𝑓
3) Right-hand-side:
𝑏𝐶 = − (𝑓𝑙𝑢𝑥𝑇 +
∑ 𝐹𝑙𝑢𝑥𝑇𝑓) = − (
𝑓~𝑛𝑏(𝐶)
→ 𝑓𝑙𝑢𝑥𝑇 =
(𝜌∗𝐶 −𝜌°𝐶 )
∆𝑡
(𝜌∗𝐶 − 𝜌°𝐶 )
∆𝑡
+
∑ 𝑚̇ ∗𝑓 )
𝑓~𝑛𝑏(𝐶)
and 𝐹𝑙𝑢𝑥𝑇𝑓 = 𝑚̇𝑓∗
The listing below shows the calculation of the fluxes of the pressure correction equation:
37
%
% assemble term I
%
rho_f [v]_f.Sf
%
U_bar_f = dot(vel_bar_f(:,:)',Sf(:,:)')';
FLUXVf = FLUXVf + rho_f.*U_bar_f;
%
% Assemble term II and linearize it
%
- rho_f ([DPVOL]_f.P_grad_f).Sf
%
DUSf = [DU1_f.*Sf(:,1),DU2_f.*Sf(:,2),DU3_f.*Sf(:,3)]; % S'f
eDUSf = [DUSf(:,1)./cfdMagnitude(DUSf),DUSf(:,2)./cfdMagnitude(DUSf),DUSf(:,3)./cfdMagnitude(DUSf)];
DUEf =
[cfdMagnitude(DUSf).*eCN(:,1)./dot(eCN(:,:)',eDUSf(:,:)')',cfdMagnitude(DUSf).*eCN(:,2)./dot(eCN(:,:)',eDUSf(:,:)'
)',cfdMagnitude(DUSf).*eCN(:,3)./dot(eCN(:,:)',eDUSf(:,:)')'];
geoDiff = cfdMagnitude(DUEf)./cfdMagnitude(CN);
DUTf = DUSf - DUEf;
FLUXCf = FLUXCf + rho_f.*geoDiff;
FLUXFf = FLUXFf - rho_f.*geoDiff;
FLUXVf = FLUXVf - rho_f.*dot(p_grad_f(iFaces,:)',DUTf(:,:)')';
%
% assemble term III
%
rho_f ([P_grad]_f.([DPVOL]_f.Sf))
%
FLUXVf = FLUXVf + rho_f.*dot(p_grad_bar_f(iFaces,:)',DUSf(:,:)')';
%
% assemble terms VIII and IX
%
(1-URF)(U_f -[v]_f.S_f)
%
FLUXVf = FLUXVf + (1.0 - mdot_f_URF)*(mdot_f_previous - rho_f.*U_bar_f);
%
% compute Rhie-Chow interpolation of mdot_f and updated it in the data base
%
mdot_f = FLUXCf .* pressureC + FLUXFf .* pressureN + FLUXVf;
theMdotField.phi(iFaces) = mdot_f;
cfdSetMeshField(theMdotField);
%
% assemble total flux
%
FLUXTf = mdot_f;
%
% assemble terms X (for compressible flow)
%
applicationClass = cfdGetApplicationClass;
if strcmp(applicationClass, 'compressible')
theDrhodpField = cfdGetMeshField('C_rho');
C_rho = theDrhodpField.phi;
C_rho_f = cfdInterpolateFromElementsToFaces('Average', C_rho);
C_rho_f = C_rho_f(iFaces);
FLUXCf = FLUXCf + (C_rho_f ./ rho_f) .* max(mdot_f, 0);
FLUXFf = FLUXFf - (C_rho_f ./ rho_f) .* max(-mdot_f, 0);
% Add transient contribution
if isTransient
deltaT = cfdGetDt;
volume = [theMesh.elements.volume]';
FLUXCE = volume(iElements) .* C_rho(iElements) / deltaT;
FLUXTE = (rho(iElements) - density_old(iElements)) .* volume(iElements) / deltaT ;
end
end
Listing 19 - Calculating fluxes of the pressure correction equation
In uFVM, the assembling of element fluxes is shown here in the listing below retrieved from the
function cfdAssembleIntoGlobalMatrixElementFluxes:
38
% Call coefficients
ac = theCoefficients.ac;
ac_old = theCoefficients.ac_old;
bc = theCoefficients.bc;
% Assemble element fluxes
for iElement = 1:numberOfElements
ac(iElement)
= ac(iElement)
+ theFluxes.FLUXCE(iElement);
ac_old(iElement) = ac_old(iElement) + theFluxes.FLUXCEOLD(iElement);
bc(iElement)
= bc(iElement)
- theFluxes.FLUXTE(iElement);
end
% Store updated coefficients
theCoefficients.ac = ac;
theCoefficients.ac_old = ac_old;
theCoefficients.bc = bc;
Listing 20 - Assembling element fluxes
The assembly of face fluxes is shown here, and it is retrieved from the function
cfdAssembleIntoGlobalMatrixFaceFluxes:
% Call coefficients
ac = theCoefficients.ac;
anb = theCoefficients.anb;
bc = theCoefficients.bc;
%
% Assemble fluxes of interior faces
%
for iFace = 1:numberOfInteriorFaces
theFace
= theMesh.faces(iFace);
iOwner
= theFace.iOwner;
iOwnerNeighbourCoef
= theFace.iOwnerNeighbourCoef;
iNeighbour
= theFace.iNeighbour;
iNeighbourOwnerCoef
= theFace.iNeighbourOwnerCoef;
%
% assemble fluxes for owner cell
%
ac(iOwner)
= ac(iOwner)
+ theFluxes.FLUXC1f(iFace);
anb{iOwner}(iOwnerNeighbourCoef) = anb{iOwner}(iOwnerNeighbourCoef) + theFluxes.FLUXC2f(iFace);
bc(iOwner)
= bc(iOwner)
- theFluxes.FLUXTf(iFace);
%
% assemble fluxes for neighbour cell
%
ac(iNeighbour)
= ac(iNeighbour)
- theFluxes.FLUXC2f(iFace);
anb{iNeighbour}(iNeighbourOwnerCoef) = anb{iNeighbour}(iNeighbourOwnerCoef) - theFluxes.FLUXC1f(iFace);
bc(iNeighbour)
= bc(iNeighbour)
+ theFluxes.FLUXTf(iFace);
end
%
% assemble fluxes of boundary faces
%
for iBFace=numberOfInteriorFaces+1:numberOfFaces
theBFace = theMesh.faces(iBFace);
iOwner
= theBFace.iOwner;
%
% assemble fluxes for owner cell
%
ac(iOwner) = ac(iOwner) + theFluxes.FLUXC1f(iBFace);
bc(iOwner) = bc(iOwner) - theFluxes.FLUXTf(iBFace);
end
% Store updated coefficients
theCoefficients.ac = ac;
theCoefficients.anb = anb;
theCoefficients.bc = bc;
Listing 21 - Assembling face fluxes
6.5.3 Solving the Equation Algebraic System
The algebraic system is solved 𝐴𝜙 ′ = 𝑏 iteratively. Direct solving of the system using Gaussian
elimination is very expensive because the matrix 𝐴 is usually large and highly sparse. There are
plenty of solvers which iteratively try to approximate the solution of the system. In uFVM, two
iterative solvers are implemented, Successive Over-relaxation (SOR) and Incomplete Lower
Upper (ILU).
39
6.5.3.1 SOR Solver
The SOR is a Gauss-Seidal solver except that it includes a factor 𝜔 which enhances the progress
of the solver. Check the function cfdSORSolver:
for iElement=1:numberOfElements
cconn = theCoefficients.cconn{iElement};
local_dphi = bc(iElement);
for iLocalNeighbour = 1:length(cconn)
iNeighbour = cconn(iLocalNeighbour);
local_dphi = local_dphi - anb{iElement}(iLocalNeighbour)*dphi(iNeighbour);
end
dphi(iElement) = local_dphi/ac(iElement);
end
for iElement=numberOfElements:-1:1
cconn = theCoefficients.cconn{iElement};
local_dphi = bc(iElement);
for iLocalNeighbour = 1:length(cconn)
iNeighbour = cconn(iLocalNeighbour);
local_dphi = local_dphi - anb{iElement}(iLocalNeighbour)*dphi(iNeighbour);
end
dphi(iElement) = local_dphi/ac(iElement);
end
Listing 22 - SOR solver
6.5.3.2 ILU Solver
Incomplete factorization of the matrix 𝐴 is an efficient preconditioner which allows for an
accelerated convergence rate of the solver. The corresponding function is: cfdILUSolver
for i1=1:numberOfElements
dc(i1) = ac(i1);
end
for i1=1:numberOfElements
dc(i1) = 1.0/dc(i1);
rc(i1) = bc(i1);
i1NbList = theCoefficients.cconn{i1};
i1NNb = length(i1NbList);
if(i1~=numberOfElements-1)
% loop over neighbours of iElement
j1_ = 1;
while(j1_<=i1NNb)
jj1 = i1NbList(j1_);
% for all neighbour j > i do
if((jj1>i1) && (jj1<=numberOfElements))
j1NbList = theCoefficients.cconn{jj1};
j1NNb = length(j1NbList);
i1_= 0;
k1 = -1;
% find _i index to get A[j][_i]
while((i1_<=j1NNb) && (k1 ~= i1))
i1_ = i1_ + 1;
k1 = j1NbList(i1_);
end
% Compute A[j][i]*D[i]*A[i][j]
if(k1 == i1)
dc(jj1) = dc(jj1) - anb{jj1}(i1_)*dc(i1)*anb{i1}(j1_);
else
disp('the index for i in j is not found');
end
end
j1_ = j1_ + 1;
end
end
end
Listing 23 - ILU solver
6.5.3.3 AMG Linear Solver
An algebraic multigrid solver remove low-frequency error components. In the context of
multigrid solvers, the direct iterative solvers like ILU solver are regarded as smoothers. Below is
40
a glance to the AMG code, it is advisable to access the code to see much more details about it.
Refer to cfdApplyAMG.
gridLevel = 1;
nCycle = 1;
if(strcmp(cycleType,'V-Cycle'))
while ((nCycle<=maxCycles)&&(finalResidual>rrf*initialResidual))
%
% Apply V-Cycle
%
finalResidual = cfdApplyVCycle(gridLevel,smootherType,maxLevels,preSweep,postSweep,rrf);
%
nCycle = nCycle + 1;
end
elseif(strcmp(cycleType,'F-Cycle'))
while ((nCycle<=maxCycles)&&(finalResidual>rrf*initialResidual))
%
% Apply F-Cycle
%
finalResidual = cfdApplyFCycle(gridLevel,smootherType,maxLevels,preSweep,postSweep,rrf);
%
nCycle = nCycle + 1;
end
elseif(strcmp(cycleType,'W-Cycle'))
while ((nCycle<=maxCycles)&&(finalResidual>rrf*initialResidual))
%
% Apply W-Cycle
%
finalResidual = cfdApplyWCycle(gridLevel,smootherType,maxLevels,preSweep,postSweep,rrf);
%
nCycle = nCycle + 1;
end
end
Listing 24 - Multigrid solver
Note: uFVM utilizes the AMG solver by default for the pressure equation only, while applies
direct iterative solvers for all other equations.
6.5.4 Correcting Equation Solution
All the equations are corrected just after solving the algebraic system 𝐴𝜙 ′ = 𝑏 such that:
𝜙𝐶 = 𝜙𝐶∗ + 𝜙𝐶′
However, after solving the pressure correction equation, the other fields have also to be
corrected:
For incompressible flow:
𝒗
∗
′
𝒗∗∗
𝐶 = 𝒗𝐶 − 𝑫𝐶 (∇𝑝 )𝐶
(𝑛)
𝑝𝐶∗ = 𝑝𝐶 + 𝜆𝑝 𝑝𝐶′
𝑚̇𝑓∗∗ = 𝑚̇𝑓∗ − 𝜌𝑓∗ ̅̅̅̅
𝑫𝑓𝒗 ∇𝑝𝑓′ . 𝑺𝑓
For compressible flow:
𝒗
∗
′
𝒗∗∗
𝐶 = 𝒗𝐶 − 𝑫𝐶 (∇𝑝 )𝐶
(𝑛)
𝑝𝐶∗ = 𝑝𝐶 + 𝜆𝑝 𝑝𝐶′
𝜌𝐶∗∗ = 𝜌𝐶∗ + 𝜆𝜌 𝐶𝜌 𝑝𝐶′
41
𝑚̇𝑓∗∗ = 𝑚̇𝑓∗ − 𝜌𝑓∗ ̅̅̅̅
𝑫𝑓𝒗 ∇𝑝𝑓′ . 𝑺𝑓 + (
𝑚̇𝑓∗
𝜌𝑓∗
) 𝐶𝜌,𝑓 𝑝𝑓′
Listing of correcting 𝑚̇𝑓∗ :
applicationClass = cfdGetApplicationClass;
if strcmp(applicationClass, 'compressible')
% Get density field
theDensityField = cfdGetMeshField('rho');
rho = theDensityField.phi;
% Get the convected density at the faces
pos = zeros(size(mdot_f));
pos(mdot_f>0) = 1;
rho_f = rho(iOwners).*pos(iFaces) + rho(iNeighbours).*(1 - pos(iFaces));
% Get drhodp field
theDrhodpField = cfdGetMeshField('C_rho');
C_rho = theDrhodpField.phi;
C_rho_f = cfdInterpolateFromElementsToFaces('Average', C_rho);
C_rho_f = C_rho_f(iFaces);
% Correct bt adding compressible contribution
mdot_f(iFaces) = mdot_f(iFaces) + (mdot_f(iFaces) ./ rho_f) .* C_rho_f .* pp(iOwners);
end
% Correct
mdot_f(iFaces) = mdot_f(iFaces) + FLUXC1f(iFaces).*pp(iOwners) + FLUXC2f(iFaces).*pp(iNeighbours);
Listing 25 - Correcting 𝑚̇𝑓∗
Listing of correcting 𝒗∗𝐶 :
thePPField = cfdGetMeshField('PP');
ppGrad = thePPField.phiGradient;
%
DUPPGRAD = [DU1.*ppGrad(iElements,1),DU2.*ppGrad(iElements,2),DU3.*ppGrad(iElements,3)];
%
theVelocityField = cfdGetMeshField('U');
vel = theVelocityField.phi;
%
vel(iElements,:) = vel(iElements,:) - DUPPGRAD(iElements,:);
Listing 26 – Correcting 𝒗∗𝐶
Listing of correcting 𝜌𝐶∗ :
thePressureCorrectionField = cfdGetMeshField('PP');
pp = thePressureCorrectionField.phi;
theDensityField = cfdGetMeshField('rho');
rho = theDensityField.phi;
theDrhodpField = cfdGetMeshField('C_rho');
C_rho = theDrhodpField.phi;
theEquation = cfdGetModel('rho');
URFRho = theEquation.urf;
rho = rho + 0.7 .* C_rho .* pp;
Listing 27 – Correcting 𝜌𝐶∗
Further details of the implementation of field corrections can be found at
cfdCorrectEquation.
42
6.6 CONVERGENCE
Once the case is run, the solution will be updated at each iteration while at the mean time a
real-time plot will be displayed showing the residuals of the equations. A sample is shown in the
figure below:
Figure 4 - Real-time residuals monitor
Once the case is run until convergence, or the maximum number of iterations are reached in a
steady state simulation, a phrase ‘Solution is converged!’ will show up on the screen notifying
the user.
In case of divergence, a pop-up message box will be displayed notifying the user that the
program has detected divergence. The box is presented below:
Figure 5 - Divergence notification
43
6.7 POST-PROCESSING
The user may plot the results on figures. The corresponding functions to plot the resulting fields
is cfdPlotField. The function takes an argument the name of the field to be plotted.
Considering that the user wants to plot the velocity field, they have to call the following in the
command window:
cfdPlotField('U')
To plot the velocity vectors, the function to be used is cfdPlotVelocity. The function takes
as arguments the vector scale, the transparency of the faces and the vector skipping criterion. If
the user wants to plot the vector field with vector scale of 1, full transparency and 10 vector
skips (skip a vector every 10 vectors), they have to call the following in the command window:
cfdPlotVelocity('vectorScale',1,'faceAlpha',0,'vectorSkip',0);
44
7 TUTORIALS
In this chapter, some tutorials will be provided from each application class. 5 application classes
can be simulated so far within uFVM: basic, incompressible, compressible, Heat Transfer and
multiphase.
7.1 BASIC
The ‘stepProfile’ case is considered here. The mesh is a square domain with 20 by 20 structured
cells. A constant velocity profile (2𝑖̂ + 1𝑗̂) is assumed over the domain as shown in the figure
below. We attempt to solve a pure convection problem in order to evaluate the quality of
convection schemes. The convected quantity is 𝜙 and has dirichlet boundary conditions at two
patches shown below in the figure, while on the other part of the boundary, the quantity 𝜙 is
set to zero gradient.
𝜕𝜙
=0
𝜕𝑛
𝜙=0
𝜙=1
The equation that is to be solved is
∇. (𝜌𝒗𝜙) = 0
where 𝒗 is the velocity field set as constant field (Refer to the 0 directory in the tutorials and
look at the file named ‘U’). The corresponding run file is as shown here:
45
%-------------------------------------------------------------------------%
% written by the CFD Group @ AUB, 2017
% contact us at: cfd@aub.edu.lb
%==========================================================================
% Case Description:
%
In this test case the square cavity problem is considered with a
%
uniform velocity profile throughout the domain. The objective is to
%
investigate the convection schemes (the default now is set to first
%
order upwind).
%-------------------------------------------------------------------------% Setup Case
cfdSetupSolverClass('basic');
% Read OpenFOAM Files
cfdReadOpenFoamFiles;
% Setup Time Settings
cfdSetupTime;
% Setup Equations
cfdDefineEquation('phi', 'div(rho*U*phi) = 0'); % Convection
% Initialize case
cfdInitializeCase;
% Run case
cfdRunCase;
In addition, we set the default option of the divSchemes in the ‘fvShemes’ directory to ‘Gauss
upwind’ which corresponds to first order upwind. Running the case generates the following
residuals plot:
Figure 6 - Residuals monitor for the upwind convection of the quantity 𝜙
46
The contour of 𝜙 is as follows:
Figure 7 - Contour plot of phi subject to first order upwind convection
47
7.2 INCOMPRESSIBLE
The ‘elbow’ case is considered here. The following is the mesh with inlet velocities shown:
16
64
1 m/s
4
3 m/s
This case simulates a water flow in the elbow at steady state conditions with no gravitational
acceleration. The governing equations are:
Momentum:
∇. (𝜌𝒗𝒗) = 𝜇∇2 𝒗 + ∇. {𝜇(∇𝒗)𝑇 } − ∇𝑝
Continuity:
∇. 𝒗 = 0
Energy:
∇. (𝜌𝑐𝑝 𝒗𝑇) = 𝑘∇2 𝑇
The corresponding run case is:
48
%-------------------------------------------------------------------------%
% written by the CFD Group @ AUB, 2017
% contact us at: cfd@aub.edu.lb
%==========================================================================
% Case Description:
%
In this test case a water flow in an elbow is simulated at steady state
%-------------------------------------------------------------------------% Setup Case
cfdSetupSolverClass('incompressible');
% Read OpenFOAM Files
cfdReadOpenFoamFiles;
% Setup Time Settings
cfdSetupTime;
% Setup Equations
cfdDefineEquation('U', 'div(rho*U*U) = laplacian(mu*U) + div(mu*transp(grad(U))) - grad(p)'); % Momentum
cfdDefineEquation('p', 'div(U) = 0'); % Continuity
cfdDefineEquation('T', 'div(rho*Cp*U*T) = laplacian(k*T)'); % Energy
% Initialize case
cfdInitializeCase;
% Run case
cfdRunCase;
Listing 28 - Run file of the 'elbow' case
Refer to the ‘elbow’ case in the tutorials directory. Running the case will solve the problem
where it converges after 105 iterations, and the figure below shows the residuals history of the
equations (U, p, and T):
Figure 8 - Residuals history of the elbow case
49
We call the following functions:
cfdPlotField('U')
cfdPlotField('p')
cfdPlotField('T')
cfdPlotVelocity('vectorScale',1,'faceAlpha',0,'vectorSkip',0);
The results for velocity ‘U’, pressure ‘p’ and temperature ‘T’ in addition to the velocity vector
field will be plotted. Presented below are the results:
Figure 9 - Velocity magnitude contour (m/s)
50
Figure 10 - Pressure contour (pa)
51
Figure 11 - Temperature contour (k)
52
Figure 12 - Velocity vectors
7.3 HEAT TRANSFER
A hot room case will be presented here. A buoyancy driven air flow is simulated in a room which
has different temperature values on its floor and ceiling while the walls are adiabatic. The mesh
is shown here:
Ceiling 𝑇 = 290 𝐾
Floor 𝑇 = 320 𝐾
53
In this problem, Boussinesq approximation is implemented to account for buoyancy effects due
to temperature gradients.
The governing equations are:
Momentum:
∇. (𝒗𝒗) = 𝜈∇2 𝒗 + ∇. {𝜈(∇𝒗)𝑇 } − ∇𝑝 − 𝑔𝛽(𝑇 − 𝑇𝑟𝑒𝑓 )
Continuity:
∇. 𝒗 = 0
Energy:
∇. (𝒗𝑇) = 𝛼∇2 𝑇
The corresponding run case is:
%-------------------------------------------------------------------------%
% written by the CFD Group @ AUB, 2017
% contact us at: cfd@aub.edu.lb
%==========================================================================
% Case Description:
%
In this test case a hot room is simulated with boussinesq
%
approximation
%-------------------------------------------------------------------------% Setup Case
cfdSetupSolverClass('heatTransfer');
% Read OpenFOAM Files
cfdReadOpenFoamFiles;
% Setup Time Settings
cfdSetupTime;
% Define new properties
cfdSetupProperty('alpha', 'model', 'nu/Pr');
% Setup Equations
cfdDefineEquation('U', 'div(U*U) = laplacian(nu*U) + div(nu*transp(grad(U))) - grad(p) - g*beta*(T - TRef)'); %
Momentum
cfdDefineEquation('p', 'div(U) = 0'); % Continuity
cfdDefineEquation('T', 'div(U*T) = laplacian(alpha*T)'); % Energy
% Initialize case
cfdInitializeCase;
% Run case
cfdRunCase;
Running the case, gives the following velocity field at an arbitrary section plane:
54
Figure 13 - Velocity field at a section plane
55
Source Exif Data:
File Type : PDF File Type Extension : pdf MIME Type : application/pdf PDF Version : 1.5 Linearized : No Page Count : 56 Language : en-US Tagged PDF : Yes Title : uFVM v1.0 Author : Computational Mechanics Lab Subject : CFD Academic Tool developed in MATLAB® – User Manual Creator : Microsoft® Word 2013 Create Date : 2017:10:09 17:22:33+03:00 Modify Date : 2017:10:09 17:22:33+03:00 Producer : Microsoft® Word 2013EXIF Metadata provided by EXIF.tools