Difference between revisions of "CFD"

From CUSF Wiki
Jump to navigation Jump to search
(discussion of continuity eqn)
(nearly finished crash course)
Line 16: Line 16:
[todo] Topics list:
[todo] Topics list:


* Meshing
* Meshing /
* Continuity equation 1D
* Continuity equation 1D /
* Momentum equation 1D
* Momentum equation 1D
* Energy equation 1D
* Energy equation 1D
Line 27: Line 27:
Every CFD solver works by splitting up the region of interest (a process called "meshing" the "domain") into a bunch of small boxes (called the "cells" in the "mesh"), and solving a set of equations in each of those cells, taking into account how the flow variables (velocity, temperature, pressure) are related between neighbouring cells. Most commonly:
Every CFD solver works by splitting up the region of interest (a process called "meshing" the "domain") into a bunch of small boxes (called the "cells" in the "mesh"), and solving a set of equations in each of those cells, taking into account how the flow variables (velocity, temperature, pressure) are related between neighbouring cells. Most commonly:


* '''The continuity equation''': mass is not created or destroyed.
* '''The continuity equation''': mass is not created or destroyed (conservation of mass).
* '''The momentum equation''': "F=ma" but for fluids.
* '''The momentum equation''': "F=ma" but for fluids (conservation of momentum).
* '''The energy equation''': energy is not created or destroyed.
* '''The energy equation''': energy is not created or destroyed (conservation of energy).
 
=== The continuity equation ===
 


To clarify this, consider the example below of a simple two-dimensional pipe, image left. Now, let us "mesh" the "domain" into four "cells" (split the pipe into 4 boxes), image right. <gallery widths="687" heights="336" perrow="2">
To clarify this, consider the example below of a simple two-dimensional pipe, image left. Now, let us "mesh" the "domain" into four "cells" (split the pipe into 4 boxes), image right. <gallery widths="687" heights="336" perrow="2">
Line 36: Line 39:
</gallery>
</gallery>
Consider an arbitrary cell within this mesh labelled "C", with its western and eastern neighbours labelled "W" and "E" respectively. The faces between W & C and C & E are labelled WC and CE. Position, velocity, and pressure are abbreviated with X, U, and p respectively. Any variable with a subscript corresponds to the value of that variable at that location. E.g. U<sub>W</sub> is the velocity of the fluid at W.
Consider an arbitrary cell within this mesh labelled "C", with its western and eastern neighbours labelled "W" and "E" respectively. The faces between W & C and C & E are labelled WC and CE. Position, velocity, and pressure are abbreviated with X, U, and p respectively. Any variable with a subscript corresponds to the value of that variable at that location. E.g. U<sub>W</sub> is the velocity of the fluid at W.
[[File:Cell-neighbours.png|none|frame|An arbitrary cell C neighboured by cells W and E on the west and east sides respectively]]


The continuity equation tells us that mass is not created or destroyed - said in another way, rate of change of mass at any point = rate of mass flow in - rate of mass flow out. For an arbitrary cell C, we can write [todo eqn]
The continuity equation tells us that mass is not created or destroyed - said in another way, rate of change of mass at any point = rate of mass flow in - rate of mass flow out. For an arbitrary cell C, we can write [todo eqn]


where [eqn rho] is the density at a given location. If we assume a constant density in the fluid (e.g. water), we can eliminate the density term. If we assume a "steady" flow, then nothing changes with time, i.e. all time derivatives are zero. So the equation becomes [todo eqn]
where [eqn rho] is the density at a given location. If we assume an "incompressible" (constant density) fluid, e.g. water, we can eliminate the density term. If we assume a "steady" flow, then nothing changes with time, i.e. all time derivatives are zero. The equation becomes [todo eqn]
 
Take a moment to predict what our solution to this velocity equation will look like, given these assumptions.


At the moment, we are using values of velocity at the faces between the cells. But in a CFD solver, we only store the values at the cell centroid. The process of calculating face values using cell values on either side of the face is called interpolation. The most obvious interpolation "scheme" (method) is linear interpolation - that is, assume that all variables vary linearly between two points. So for example, halfway between W and C the value of U would be 1/2 Uw + 1/2 Uc [todo eqn]. The continuity equation becomes [todo eqn]
At the moment, we are using values of velocity at the faces between the cells. But in a CFD solver, we only store the values at the cell centroid. The process of calculating face values using cell values on either side of the face is called interpolation. The most obvious interpolation "scheme" (method) is linear interpolation - that is, assume that all variables vary linearly between two points. So for example, halfway between W and C the value of U would be 1/2 Uw + 1/2 Uc [todo eqn]. The continuity equation becomes [todo eqn]


We can now apply this equation for each cell in the mesh. [todo 4 eqns]
We can now apply this equation for cells 2 and 3. However, we can't apply this equation to cells 1 and 4, since they are both missing a neighbour (as they are on the boundary of the domain). Without knowing what the value of U is at the boundaries, we cannot solve (or even fully write down) the continuity equation. Therefore we must apply some special treatment at the boundaries (called boundary conditions) in order to solve this issue. For the moment, let the flow velocity on the left boundary be fixed at 5 m/s (called a Dirichlet boundary - mnemonic: kind of sounds like "directly") and let the flow velocity on the right boundary be the same as the flow velocity at cell 4 (also known as a zero gradient condition, or von Neumann boundary). Applying these conditions, we can now write down equations for cells 1 and 4, as well as the equations for cells 2 and 3 as found before. [todo 4 eqns]
 
We now have a set of four simultaneous equations to solve for four unknowns (the velocity at each cell centroid), and we can solve this now, and find that the velocity everywhere is 5 m/s. Hooray!
 
Let us simplify the notation a bit using matrices and vectors. [todo matrix eqn] MU = B
 
We could now invert the matrix M, and calculate the solution vector U using U=M^-1 * B [todo eqn]. The matrix M always has size n x n, where n is the number of cells in the mesh (in this case 4). But typical CFD cases will have n in the thousands or millions - inverting matrices of this size becomes very time consuming. So instead, CFD solvers use an iterative approach on the equation MU = B, where they start with a guess for U, and improve their guess with successive iterations of an algorithm until the solution vector U is "close enough" (this threshold is determined by the CFD user). The error in a solution for a variable is characterised by the "residual". A small residual means the solution has only a small error, and vice versa. The CFD user tells the software what residual value qualifies as "close enough".
 
Important note, the solution vector U does not store the components of velocity, i.e. Ux, Uy, Uz, but instead stores the velocities at each cell centroid, i.e. U1, U2, U3, U4. In 1D this might be more clear, but in higher dimensions it is important not to get confused.
 
=== The momentum equation ===
In a similar process to the continuity equation, we can convert "F=ma but for fluids" into an actual equation, with the intermediate equation: rate of change of momentum at any point = rate of momentum flow in - rate of momentum flow out. For an arbitrary cell C experiencing pressure forces and an external force F, we can write [todo eqn]
 
Again we apply the assumptions of steady incompressible flow, and arrive at [todo eqn]
 
Note the term p/rho [todo eqn] is sometimes called the "kinematic pressure" with units [todo dims], and is sometimes what the CFD software (e.g. in OpenFOAM) will confusingly call p. So when finding the pressure drop when using an incompressible fluid, you may have to multiply by density. Also note that the absolute pressure doesn't matter here - it is only the differences in pressure that matter (only true for an incompressible flow - for compressible flow, density is a function of absolute pressure).
 
For a more complete equation, you would also include the effects of viscosity and turbulence in this equation (added to the force F) - these would be called the Navier-Stokes equations. For simplicity, we will neglect viscosity and turbulence.
 
By applying the boundary conditions of zero pressure gradient on the left boundary, and a fixed pressure of 0 on the right boundary, we arrive at these equations [todo 4 eqns]


We now have a set of four simultaneous equations to solve for four unknowns (the pressure at each cell centroid), and we can solve this now, and find that the pressure everywhere is 0 Pa. Hooray!


The continuity equation tells us that mass is not created or destroyed. And if we assume the flow through the pipe is "steady" (i.e. it does not change with time, so fluid does not accumulate within any cell), we should get the result that the flow velocity in every cell is the same. But for the purpose of understanding what a CFD solver does, let us continue with the example.
Again, CFD solvers actually work with these equations in matrix form, and use an iterative approach to solve them, until the residual is below a threshold value set by the user.


=== The energy equation ===
In case we also wanted to find the temperature of the fluid throughout the domain, we could solve the energy equation: rate of change of energy at any point = rate of energy flow in - rate of energy flow out. I will not go through the derivation here (although it's fairly simple) - it is enough to know that the CFD solver performs an iterative approach to solve this equation for temperature until the residual is low enough.


Then, let us set the velocity at the start of the pipe to be 5 m/s to the right (called a "boundary condition").  
=== Higher dimensions ===
In 1D, we have been solving for conservation of momentum in just one direction - the x axis. In 2D, we must also solve for momentum in the y axis to obtain the component of velocity in the y direction. Similarly in 3D, we need to solve for momentum in all three axes: x, y, and z to obtain the components of velocity in each direction. The momentum equation can be generalised to work in arbitrary dimensions


Why do we mesh? The reason we split a domain into cells is because a flow variable (e.g. velocity) could vary continuously along the length of the pipe, meaning you'd need an infinite amount of memory to store it. But we don't have infinite memory, so instead we store the flow variables at each cell centroid (here marked with circles), and approximate the variables as varying linearly between each cell centroid.
Why do we mesh? The reason we split a domain into cells is because a flow variable (e.g. velocity) could vary continuously along the length of the pipe, meaning you'd need an infinite amount of memory to store it. But we don't have infinite memory, so instead we store the flow variables at each cell centroid (here marked with circles), and approximate the variables as varying linearly between each cell centroid.

Revision as of 00:04, 4 January 2022

This page is currently a work in progress by Nik Lebedenko - if you want to ask questions / see a mistake / request more details to be added to this page, feel free to email me or message me on slack.

Computational Fluid Dynamics (CFD) is the process of solving fluid mechanics equations with the use of computers. Think F=ma, but like... more. (I will elaborate further)

Crash course on the Finite Volume Method (FVM)

This section is incomplete

There are many resources for learning how CFD works properly (hopefully linked below, TODO), but you're not here for that. This section is designed for people who are doing a CFD simulation for just a one-time thing, and they can treat it like a black box that gives them answers (that is how I started). So, here are the bare minimum requirements for understanding roughly what is going on in any CFD solver. Nonetheless, I will assume you know:

  • Some basic high-school physics about gases (e.g. what is pressure, how does it change with temperature)
  • What a partial derivative is,
  • What simultaneous equations are,
  • How to write simultaneous equations as matrices, and
  • What a vector is.

[todo] Topics list:

  • Meshing /
  • Continuity equation 1D /
  • Momentum equation 1D
  • Energy equation 1D
  • Turbulence
  • Higher dimensions
  • Numerical schemes and Stability


Every CFD solver works by splitting up the region of interest (a process called "meshing" the "domain") into a bunch of small boxes (called the "cells" in the "mesh"), and solving a set of equations in each of those cells, taking into account how the flow variables (velocity, temperature, pressure) are related between neighbouring cells. Most commonly:

  • The continuity equation: mass is not created or destroyed (conservation of mass).
  • The momentum equation: "F=ma" but for fluids (conservation of momentum).
  • The energy equation: energy is not created or destroyed (conservation of energy).

The continuity equation

To clarify this, consider the example below of a simple two-dimensional pipe, image left. Now, let us "mesh" the "domain" into four "cells" (split the pipe into 4 boxes), image right.

  • 2D pipe with a flow from left to right

  • 2D pipe split into cells, with cell centroids marked with circles

Consider an arbitrary cell within this mesh labelled "C", with its western and eastern neighbours labelled "W" and "E" respectively. The faces between W & C and C & E are labelled WC and CE. Position, velocity, and pressure are abbreviated with X, U, and p respectively. Any variable with a subscript corresponds to the value of that variable at that location. E.g. UW is the velocity of the fluid at W.

An arbitrary cell C neighboured by cells W and E on the west and east sides respectively

The continuity equation tells us that mass is not created or destroyed - said in another way, rate of change of mass at any point = rate of mass flow in - rate of mass flow out. For an arbitrary cell C, we can write [todo eqn]

where [eqn rho] is the density at a given location. If we assume an "incompressible" (constant density) fluid, e.g. water, we can eliminate the density term. If we assume a "steady" flow, then nothing changes with time, i.e. all time derivatives are zero. The equation becomes [todo eqn]

Take a moment to predict what our solution to this velocity equation will look like, given these assumptions.

At the moment, we are using values of velocity at the faces between the cells. But in a CFD solver, we only store the values at the cell centroid. The process of calculating face values using cell values on either side of the face is called interpolation. The most obvious interpolation "scheme" (method) is linear interpolation - that is, assume that all variables vary linearly between two points. So for example, halfway between W and C the value of U would be 1/2 Uw + 1/2 Uc [todo eqn]. The continuity equation becomes [todo eqn]

We can now apply this equation for cells 2 and 3. However, we can't apply this equation to cells 1 and 4, since they are both missing a neighbour (as they are on the boundary of the domain). Without knowing what the value of U is at the boundaries, we cannot solve (or even fully write down) the continuity equation. Therefore we must apply some special treatment at the boundaries (called boundary conditions) in order to solve this issue. For the moment, let the flow velocity on the left boundary be fixed at 5 m/s (called a Dirichlet boundary - mnemonic: kind of sounds like "directly") and let the flow velocity on the right boundary be the same as the flow velocity at cell 4 (also known as a zero gradient condition, or von Neumann boundary). Applying these conditions, we can now write down equations for cells 1 and 4, as well as the equations for cells 2 and 3 as found before. [todo 4 eqns]

We now have a set of four simultaneous equations to solve for four unknowns (the velocity at each cell centroid), and we can solve this now, and find that the velocity everywhere is 5 m/s. Hooray!

Let us simplify the notation a bit using matrices and vectors. [todo matrix eqn] MU = B

We could now invert the matrix M, and calculate the solution vector U using U=M^-1 * B [todo eqn]. The matrix M always has size n x n, where n is the number of cells in the mesh (in this case 4). But typical CFD cases will have n in the thousands or millions - inverting matrices of this size becomes very time consuming. So instead, CFD solvers use an iterative approach on the equation MU = B, where they start with a guess for U, and improve their guess with successive iterations of an algorithm until the solution vector U is "close enough" (this threshold is determined by the CFD user). The error in a solution for a variable is characterised by the "residual". A small residual means the solution has only a small error, and vice versa. The CFD user tells the software what residual value qualifies as "close enough".

Important note, the solution vector U does not store the components of velocity, i.e. Ux, Uy, Uz, but instead stores the velocities at each cell centroid, i.e. U1, U2, U3, U4. In 1D this might be more clear, but in higher dimensions it is important not to get confused.

The momentum equation

In a similar process to the continuity equation, we can convert "F=ma but for fluids" into an actual equation, with the intermediate equation: rate of change of momentum at any point = rate of momentum flow in - rate of momentum flow out. For an arbitrary cell C experiencing pressure forces and an external force F, we can write [todo eqn]

Again we apply the assumptions of steady incompressible flow, and arrive at [todo eqn]

Note the term p/rho [todo eqn] is sometimes called the "kinematic pressure" with units [todo dims], and is sometimes what the CFD software (e.g. in OpenFOAM) will confusingly call p. So when finding the pressure drop when using an incompressible fluid, you may have to multiply by density. Also note that the absolute pressure doesn't matter here - it is only the differences in pressure that matter (only true for an incompressible flow - for compressible flow, density is a function of absolute pressure).

For a more complete equation, you would also include the effects of viscosity and turbulence in this equation (added to the force F) - these would be called the Navier-Stokes equations. For simplicity, we will neglect viscosity and turbulence.

By applying the boundary conditions of zero pressure gradient on the left boundary, and a fixed pressure of 0 on the right boundary, we arrive at these equations [todo 4 eqns]

We now have a set of four simultaneous equations to solve for four unknowns (the pressure at each cell centroid), and we can solve this now, and find that the pressure everywhere is 0 Pa. Hooray!

Again, CFD solvers actually work with these equations in matrix form, and use an iterative approach to solve them, until the residual is below a threshold value set by the user.

The energy equation

In case we also wanted to find the temperature of the fluid throughout the domain, we could solve the energy equation: rate of change of energy at any point = rate of energy flow in - rate of energy flow out. I will not go through the derivation here (although it's fairly simple) - it is enough to know that the CFD solver performs an iterative approach to solve this equation for temperature until the residual is low enough.

Higher dimensions

In 1D, we have been solving for conservation of momentum in just one direction - the x axis. In 2D, we must also solve for momentum in the y axis to obtain the component of velocity in the y direction. Similarly in 3D, we need to solve for momentum in all three axes: x, y, and z to obtain the components of velocity in each direction. The momentum equation can be generalised to work in arbitrary dimensions

Why do we mesh? The reason we split a domain into cells is because a flow variable (e.g. velocity) could vary continuously along the length of the pipe, meaning you'd need an infinite amount of memory to store it. But we don't have infinite memory, so instead we store the flow variables at each cell centroid (here marked with circles), and approximate the variables as varying linearly between each cell centroid.

Fluid-structure interaction (FSI)

IMPORTANT: Although I hope that these instructions will be easily understood by anyone, this is not a beginner's project! If you want a proper introduction to CFD, this is not it. (I am in the process of creating a beginner's guide which will be available soon™)

This technique allows complex interactions between flexible structures and fluid flows. This has been used to predict fin flutter at hypersonic speeds for Griffin 1.

Simulations of this kind are achievable with the combination of the following free open-source software packages:

  • OpenFOAM (website, install on Ubuntu) - a C/C++ library for the solution of fluid mechanics problems formulated with the Finite Volume Method (FVM),
    • Important: there are two main versions of OpenFOAM, provided at OpenFOAM.com and OpenFOAM.org. I have only used the .com version, but everything here should also be possible with the .org version.
  • FEniCS (website, install on Ubuntu) - a Python interface to the DOLFIN C++ library for the solution of solid mechanics problems formulated with the Finite Element Method (FEM), and
  • preCICE (website, install on Ubuntu) - a very powerful library which allows the coupling of arbitrary solvers using human-readable syntax.
    • Note: in order to use OpenFOAM and FEniCS with preCICE, it is also necessary to download their respective adapters. See this page for OpenFOAM and this page for FEniCS from the preCICE website for further information.

Step-by-step instructions

IMPORTANT: Although I include the links to the official websites for Ubuntu installation instructions, I recommend you follow my instructions for installation (unless you are a masochist), since the official instructions contain some mistakes and omissions.

Install Ubuntu

If you are running a Windows machine (Windows 10 or above), I recommend installing Windows Subsystem for Linux (WSL - use WSL1, discussion of WSL2 below) instead of setting up a dual-boot partition.

  • Note: there is a choice between WSL1 and WSL2, with the main difference being that WSL2 has significantly improved read-write speeds compared with WSL1 while staying within the Linux file system, but WSL2 has significantly worse read-write speeds when moving between the Linux and Windows file systems. In practical terms, this means that the Linux-based software (OpenFOAM, FEniCS) will read and write faster, at the expense of Windows-based ParaView installations loading results much slower when post-processing. Since OpenFOAM spends significantly more time solving equations than reading and writing files, and gigabytes of data need to be loaded across file systems for post-processing in Windows, I recommend sticking with WSL1. If you prefer WSL2, you can install ParaView on Ubuntu, installing XMing to view the GUI (see this page and this page for further information)

Install OpenFOAM on Ubuntu

The Ubuntu installation instructions as of 31/12/2021 are summarised below (adapted from this website):

curl -s https://dl.openfoam.com/add-debian-repo.sh | sudo bash
sudo apt-get install openfoam2012-default

Ensure that your ~/.bashrc file has the following line appended (with XXXX replaced with your version number):

source /usr/lib/openfoam/openfoamXXXX/etc/bashrc

Restart your Ubuntu. Now check the installation completed correctly by running a tutorial case:

cd ~
mkdir -p OpenFOAM-sims/tutorials
cp -r $FOAM_TUTORIALS/incompressible OpenFOAM-sims/tutorials
cd OpenFOAM-sims/tutorials/incompressible/icoFoam/cavity/cavity/
blockMesh > log.blockMesh
icoFoam > log.icoFoam

Note: I use OpenFOAMv2012, but later versions should also work.

Install FEniCS on Ubuntu

The Ubuntu installation instructions as of 31/12/2021 are summarised below (adapted from this website):

sudo apt-get install software-properties-common
sudo add-apt-repository ppa:fenics-packages/fenics
sudo apt-get update
sudo apt-get install fenics

Install preCICE on Ubuntu

The Ubuntu installation instructions as of 31/12/2021 are summarised below (adapted from this website): 

sudo apt update
sudo apt install build-essential cmake libeigen3-dev libxml2-dev libboost-all-dev petsc-dev python3-dev python3-numpy
cd ~
mkdir FSI
cd FSI

Now either install the Ubuntu package directly with the instructions below (which did not work for me)...

wget https://github.com/precice/precice/releases/download/v2.3.0/libprecice2_2.3.0_focal.deb
sudo apt install ./libprecice2_2.3.0_focal.deb

... or install from source:

tar -xzf v2.3.0.tar.gz
cd precice-2.3.0
mkdir build
cd build
cmake -DBUILD_SHARED_LIBS=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=~/FSI/precice-2.3.0 -DPRECICE_MPICommunication=ON .. ### IMPORTANT! INCLUDE THE TWO DOTS!
make -j 4

Test that the cmake command worked as expected:

cd ~/FSI/precice-2.3.0/build
ctest --output-on-failure

If no errors were shown, install the software and test it worked:

make install
make test_install

Note: I used preCICE v2.3.0, but later versions should also work.

Install the OpenFOAM-preCICE adapter

The Ubuntu installation instructions as of 31/12/2021 are summarised below (adapted from this website):

cd ~/FSI
wget https://github.com/precice/openfoam-adapter/releases/download/v1.0.0/openfoam-adapter_v1.0.0_OpenFOAMv1812-v2012.tar.gz
tar -xzf openfoam-adapter_v1.0.0_OpenFOAMv1812-v2012.tar.gz
cd openfoam-adapter_v1.0.0_OpenFOAMv1812-v2012/
export LD_LIBRARY_PATH=/usr/lib/x86_64-linux-gnu:$LD_LIBRARY_PATH ### This line may-or-may-not be needed. YMMV.
./Allwmake

Note: this is for OpenFOAMv2012. Please see the preCICE website for information on adapters for later versions.

Install the FEniCS-preCICE adapter

The Ubuntu installation instructions as of 31/12/2021 are summarised below (adapted from this website). This requires Python3, which should already be installed if you followed my previous instructions: 

(optional, should not need this line) python3 -m pip install scipy
python3 -m pip install --user fenicsprecice

Test the installation ran correctly

... by running this tutorial case: 

cd ~/FSI
git clone --branch=master --depth 1 https://github.com/precice/tutorials.git
cd ~/FSI/tutorials/perpendicular-flap/fluid-openfoam
./run.sh

The command line output should have paused after this output. If it did not pause, something went wrong. 

---[precice]  I am participant "Fluid"
---[precice]  Setting up master communication to coupling partner/s

Now open another Ubuntu terminal, and run these commands:

cd ~/FSI/tutorials/perpendicular-flap/solid-fenics
./run.sh

Both terminals should now have log outputs flying. After around 1 minute of computation, both terminals should finish at (roughly) the same time. You can view the results with the following commands, assuming you have installed ParaView on Windows and added the folder containing its .exe file to the Windows PATH environment variable:

cd ~/FSI/tutorials/perpendicular-flap/fluid-openfoam
paraview.exe fluid-openfoam.foam

Congratulations! You should now have a functional FSI installation.

Note that as of 31/12/2021, the unmodified FEniCS-preCICE adapter only works in 2D, so you will have to edit the adapter (or download my edited version, TODO add link here) in order to add 3D FSI functionality. Or just wait a few months until I get around to submitting a pull request for my version.

TODO:

  • ParaView install instructions
  • Resources for learning CFD and FEM
  • A complete line-by-line breakdown of simpleFoam
  • rhoCentralFoam resources
  • rhoCentralFoam rotating reference frame
  • Tips and Tricks from experience
  • How to do 3D hypersonic fin flutter
Retrieved from ‘https://wiki.cusf.co.uk/index.php?title=CFD&oldid=706’