Difference between revisions of "CFD"

From CUSF Wiki
Jump to navigation Jump to search
m
m (corrected Moukalled name)
 
(34 intermediate revisions by the same user not shown)
Line 1: Line 1:
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.
This page is maintained 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)
Computational Fluid Dynamics (CFD) is the process of solving fluid mechanics equations with the use of computers, usually with the Finite Volume Method.  
 
Essentially F=ma, but with many simultaneous equations.
 
[[OpenFOAM installation instructions|OpenFOAM installation instructions here]]
 
= Getting Started with CFD =
[[OpenFOAM installation instructions|OpenFOAM installation instructions here]]
 
Skip this section if you dislike broad-strokes philosophical discussion.
 
The path you should take for CFD differs depending on what your end goal is, and what budget you have. The main choice is between OpenFOAM and a commercial package (ANSYS Fluent / CFX, STAR-CCM, etc).
 
OpenFOAM Pros vs commercial packages:
 
* Free
* Open Source
* You can edit any part of the software, any setting, any way you want to.
* No limits, no "paid version only" features. (ANSYS limits you to 2 CPU threads and 100k cells for simulation)
 
OpenFOAM Cons vs commercial packages:
 
* Difficult to install
* Difficult to use
* Difficult to learn - less of a learning curve, more like a cliff.
* No graphical user interface
* Mesh generation software is terrible. IMO, the best way to generate a mesh in OpenFOAM, is not to generate a mesh in OpenFOAM
** Pro tip: use ANSYS meshing software (either the mechanical mesher, or the Fluent mesher) to generate CFD meshes. It will be 10x quicker (including setup time), take 5x less RAM, and be a higher quality than what you could achieve with OpenFOAM. Then import it into OpenFOAM with <code>fluentMeshToFoam <options> <.msh file></code>. ANSYS does not put a cap on the mesh size you can generate, only the mesh size you can simulate is limited.
** There are other free & open source meshing programs, e.g. Gmsh, but I have not tried them.
 
OpenFOAM throws you into the deep end, and forces you to learn nearly everything at once. If you stick through to the end of the learning cliff, you will understand much more about CFD than you would have otherwise learned with a commercial package.
 
If you just need to do one or two small sims, a commercial package e.g. ANSYS is your best choice. (Just not SOLIDWORKS CFD simulations).
 
The remainder of this wiki is dedicated to free and open source software, including OpenFOAM.


= Crash course on the Finite Volume Method (FVM) =
= 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:
There are many resources for learning how CFD works properly (linked below), but this is a crash course. This section is designed for people who would like to find out a bit more about how CFD works, beyond a "black-box" understanding, but without going too far in depth with all the details. So, here are the bare minimum requirements for understanding roughly what is going on under the hood in any CFD solver. I will assume you know:


* Some basic high-school physics about gases (e.g. what is pressure, how does it change with temperature)  
* Some basic high-school physics about gases (e.g. what is pressure, how does it change with temperature)  
Line 14: Line 47:
* What a vector is.
* What a vector is.


[todo] Topics list:
Sections marked "Optional" are not required for a basic level of understanding of how the software works, but give a bit more detail for those who want it.


* Meshing /
=== Overview ===
* Continuity equation 1D /
* Momentum equation 1D
* Energy equation 1D
* Turbulence
* Higher dimensions
* Numerical schemes and Stability


equation test?: <nowiki><html><math>E=mc^2</math></nowiki><nowiki></html></nowiki>
The most common method of doing CFD is with the Finite Volume Method (FVM). An informal derivation of some flow equations along with application of the FVM to them is in this section.


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:
''It is also possible to do CFD with the Finite Element Method (FEM), which is more typically associated with Structural Analysis / Elastodynamics problems, but we won't discuss that here.''
 
CFD solvers work 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 continuity equation''': mass is not created or destroyed (conservation of mass).
Line 33: Line 62:


=== The continuity equation ===
=== 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 uniform pipe, with a cross-section shown in the left image. Now, let us ''mesh'' the ''domain'' into four ''cells'' (split the pipe into 4 boxes), image right. <gallery widths="687" heights="336" perrow="2">
File:Pipe.png|2D pipe with a flow from left to right
File:Pipe.png|2D pipe cross-section with a flow from left to right
File:Pipe-meshed3.png|2D pipe split into cells, with cell centroids marked with circles
File:Pipe-meshed3.png|Pipe split into cells, with cell centroids marked with circles
</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. The following abbreviations for variables will be used:
Note: the cells here are 3D - we are looking at a 2D cross section. We will assume that the velocity is purely in the x direction.
Consider an arbitrary cell within this mesh labelled <math>P</math>, with its western and eastern neighbours labelled <math>W</math> and <math>E</math> respectively. The faces between <math>W</math> & <math>P</math> and <math>P</math> & <math>E</math> are labelled <math>w</math> and <math>c</math>. The following abbreviations for flow variables will be used:


* X: Position,
*<math>x</math>: Position,
*rho [eqn]: Density,
*<math>\rho</math> (rho): Density,
*p: Pressure,
*<math>p</math>: Pressure,
*U: Velocity (in the x-axis),
*<math>U</math>: Velocity (in the x-axis),
*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.
*<math>V</math>: Volume of the cell,
*<math>A</math>: Area of a face,
*Any variable with a subscript corresponds to the value of that variable at that location. E.g. <math>U_{\small{W}}</math> is the velocity of the fluid at <math>W</math>. The area of face <math>w</math> is <math>A_w</math>.
<gallery widths="687" heights="336">
<gallery widths="687" heights="336">
File:Cell-neighbours.png|An arbitrary cell C neighboured by cells W and E on the west and east sides respectively
File:Cell-neighbours.png|An arbitrary cell P neighboured by cells W and E on the west and east sides respectively
File:Linear-interp.png|A graph showing cell values (grey circles) and the linear interpolated face values (dotted blue circles)
File:Linear-interp.png|A graph showing cell values (grey circles) and the linear interpolated face values (blue circles)
</gallery>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]
</gallery>The continuity equation tells us that mass is not created or destroyed. For an arbitrary cell P, knowing that the mass flow rate through a given face is <math>{\rho}UA</math>, we can write
 
Rate of change of mass in <math>P</math> = Mass flow rate in - Mass flow rate out
 
<math>\frac{\partial (\rho V)_{\small{P}}}{\partial t} = ({\rho}UA)_{w} - ({\rho}UA)_{e}</math>


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]
If we assume an ''incompressible'' (constant density) fluid, e.g. water, we can divide by the density term <math>\rho</math> on both sides. If we assume a ''steady'' flow, then nothing changes with time, i.e. all time derivatives are zero. The equation becomes


Take a moment to predict what our solution to this velocity equation will look like, given these assumptions.
<math>0 = ({\rho}UA)_w - ({\rho}UA)_e</math>


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, as in the image above. The continuity equation becomes [todo eqn]
<math>0 = (UA)_w - (UA)_e</math> - '''(1)'''


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. Boundary conditions are specified by the CFD user. For the moment, we will use these boundary conditions:
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, as in the image above. Note: interpolation isn't unique to the variable <math>U</math> - it can be applied to all flow variables.


* Flow velocity on the left boundary is fixed at 5 m/s (called a Dirichlet boundary - mnemonic: kind of sounds like "directly").
In our case, the cells and faces are all perfectly uniformly spaced, so if we use linear interpolation to work out the face values,
* Flow velocity on the right boundary is 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]
<math>U_w = {1 \over 2} (U_{\small{W}} + U_{\small{P}}), \quad U_e = {1 \over 2} (U_{\small{P}} + U_{\small{E}})</math>


We now have a set of four simultaneous equations to solve for four unknowns (the velocity at each cell centroid). We can solve this now, and find that the velocity everywhere is 5 m/s. Hooray!
This is the same as saying that the face values are the average of the neighbouring cell values. Let's use a more general equation to represent linear interpolation.


Let us simplify the notation a bit using matrices and vectors. [todo matrix eqn] MU = B
<math>U_w = i_{\small{W}}U_{\small{W}} + (1-i_{\small{W}})U_{\small{P}}, \quad U_e = i_{\small{E}}U_{\small{P}} + (1-i_{\small{E}})U_{\small{E}}</math>


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 specifies what residual value qualifies as "close enough".
where <math>i_{\small{W}} = \frac{x_w - x_{\small{W}}}{x_{\small{P}} - x_{\small{W}}}</math> and <math> i_{\small{E}} = \frac{x_e - x_{\small{P}}}{x_{\small{E}} - x_{\small{P}}} </math>


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, since we are only concerned with x-velocity, but in higher dimensions it is important not to get confused.
''Don't be intimidated by the complicated-looking fractions!'' <math>i_{\small{W}}</math> simply represents how far the face <math>w</math> is from cell <math>W</math>, divided by the distance between <math>W</math> and <math>P</math>. This way, if the face is exactly half-way between <math>W</math> and <math>P</math>, then <math>i_{\small{W}} = {1 \over 2}</math>, giving the same equation as before.  
If instead, the face was closer to <math>W</math> than to <math>P</math>, then <math>i_{\small{W}} < {1 \over 2}</math>, and vice versa.


=== The momentum equation ===
After this point, the continuity equation will only work for cells that aren't on the boundary. We also factor out the areas of each face, since they are just constants. The continuity equation becomes
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]
<math>0 = (U_{\small{W}}i_{\small{W}} + U_{\small{P}}(1-i_{\small{W}})) \cdot A_w - (U_{\small{P}}i_{\small{E}} + U_{\small{E}}(1-i_{\small{E}})) \cdot A_e</math>


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).
We now expand and simplify the equation, combining the constants into single variables.


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.
<math>0 = U_{\small{W}}i_{\small{W}}A_w + U_{\small{P}}(1-i_{\small{W}})A_w - U_{\small{P}}i_{\small{E}}A_e - U_{\small{E}}(1-i_{\small{E}})A_e</math>


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]
<math>0 = U_{\small{W}}k_{\small{P,W}} + U_{\small{P}}k_{\small{P,P}} + U_{\small{E}}k_{\small{P,E}}</math>


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!
where
* <math>k_{\small{P,W}} = i_{\small{W}}A_w </math>
* <math>k_{\small{P,P}} = (1-i_{\small{W}})A_w - i_{\small{E}}A_e</math>
* <math>k_{\small{P,E}} = -(1-i_{\small{E}})A_e</math>.


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.
We can now apply this equation for cells 2 and 3. However, we can't apply this equation to cells 1 and 4, as they are on the boundary of the domain - the linear interpolation step does not work on the boundary. Without knowing what the value of U is at the boundaries, we cannot solve the continuity equation. Therefore we must apply some special treatment at the boundaries (called boundary conditions) in order to solve this issue. Boundary conditions are specified by the CFD user. We will use these boundary conditions:


=== The energy equation ===
*Flow velocity on the left boundary is fixed at 5 m/s (called a Dirichlet boundary - mnemonic: kind of sounds like "directly").
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 a given 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 for temperature using the conservation of energy until the residual is low enough.
**Plugging this condition into equation (1) gives the equation <math>0 = 5A_w - U_{\small{P}}i_{\small{E}}A_e - U_{\small{E}}(1-i_{\small{E}})A_e</math>
* The gradient of velocity <math>\frac {\partial U}{\partial x}</math> on the right boundary is zero - so the velocity on the right boundary is the same as the velocity at cell centroid 4.
**Plugging this condition into equation (1) gives the equation <math>0 = U_{\small{W}}i_{\small{W}}A_w + U_{\small{P}}(1-i_{\small{W}})A_w - U_{\small{P}}A_e</math>


In the case of an incompressible flow, this isn't necessary, since density is constant. But for a compressible flow, the density at a given location is a function of temperature, so it is necessary to find the temperature everywhere in the fluid, so that the density can be found.
* The constants are again combined into variables of the form <math>k_{a,b}</math>. Here is a summary of all the equations found so far:
*Cell 1: <math>0 = 5A_w + U_1k_{1,1} + U_2k_{1,2}</math>
*Cell 2: <math>0 = 0 + U_1k_{2,1} + U_2k_{2,2} + U_3k_{2,3}</math>
*Cell 3: <math>0 = 0 + U_2k_{3,2} + U_3k_{3,3} + U_4k_{3,4}</math>
*Cell 4: <math>0 = 0 + U_3k_{4,3} + U_4k_{4,4}</math>
 
We now have a set of four simultaneous equations to solve for four unknowns (the velocity at each cell centroid).
 
Let us simplify the notation a bit using matrices and vectors.
 
<math>
\begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} 5A_w \\ 0 \\ 0 \\ 0 \end{bmatrix} + \begin{bmatrix} k_{1,1} & k_{1,2} & 0 & 0 \\
                                                                                                                  k_{2,1} & k_{2,2} & k_{2,3} & 0 \\
                                                                                                                  0 & k_{3,2} & k_{3,3} & k_{3,4} \\
                                                                                                                  0 & 0 & k_{4,3} & k_{4,4} \end{bmatrix} \begin{bmatrix} U_1 \\ U_2 \\ U_3 \\ U_4 \end{bmatrix}
</math>
 
<math>
\begin{bmatrix} k_{1,1} & k_{1,2} & 0 & 0 \\
                k_{2,1} & k_{2,2} & k_{2,3} & 0 \\
                0 & k_{3,2} & k_{3,3} & k_{3,4} \\
                0 & 0 & k_{4,3} & k_{4,4} \end{bmatrix} \begin{bmatrix} U_1 \\ U_2 \\ U_3 \\ U_4 \end{bmatrix} = \begin{bmatrix} -5A_w \\ 0 \\ 0 \\ 0 \end{bmatrix}
</math>.
 
Abbreviated with a matrix of coefficients <math>\boldsymbol{M}</math>, a solution vector <math>\boldsymbol{U}</math>, and a right-hand-side vector <math>\boldsymbol{B}</math>:
 
<math> \boldsymbol{MU=B} </math>, where
 
<math> \boldsymbol{M}= \begin{bmatrix}
k_{1,1} & k_{1,2} & 0 & 0 \\
k_{2,1} & k_{2,2} & k_{2,3} & 0 \\
0 & k_{3,2} & k_{3,3} & k_{3,4} \\
0 & 0 & k_{4,3} & k_{4,4} \end{bmatrix} </math>,
<math>\boldsymbol{U} = \begin{bmatrix} U_1 \\ U_2 \\ U_3 \\ U_4 \end{bmatrix}</math>, and
<math>\boldsymbol{B} = \begin{bmatrix} -5A_w \\ 0 \\ 0 \\ 0 \end{bmatrix}</math>
 
At this point a CFD solver would calculate each of the coefficients of the matrix and RHS vector. Here we have a uniform cell spacing and a uniform face area of 1.
 
<math> \boldsymbol{M}= \begin{bmatrix}
-0.5 & -0.5 & 0 & 0 \\
0.5 & 0 & -0.5 & 0 \\
0 & 0.5 & 0 & -0.5 \\
0 & 0 & 0.5 & -0.5 \end{bmatrix} </math>
 
<math>\boldsymbol{B} = \begin{bmatrix} -5 \\ 0 \\ 0 \\ 0 \end{bmatrix}</math>
 
We could now invert the matrix <math>\boldsymbol{M}</math>, and calculate the solution vector <math>\boldsymbol{U}</math> using <math>\boldsymbol{U=M^{-1}B}</math>.
 
<math> \boldsymbol{M^{-1}}= \begin{bmatrix}
-1 & 1 & -1 & 1 \\
-1 & -1 & 1 & -1 \\
-1 & -1 & -1 & 1 \\
-1 & -1 & -1 & -1 \end{bmatrix} </math>
This yields a solution of U = 5 m/s everywhere. Hooray!
 
''Think: does this result make sense, given our assumptions?''
 
===Why matrix inversion rarely happens===
The matrix <math>\boldsymbol{M}</math> 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 <math> \boldsymbol{MU=B} </math>, where they start with a guess for <math>\boldsymbol{U}</math>, and improve their guess with successive iterations of an algorithm until the solution vector <math>\boldsymbol{U}</math> 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 specifies what residual value qualifies as "close enough".
 
===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, noting that
 
<math>UA = </math> the rate of flow of volume through a face, and
 
<math>\frac{mU}{V} = \frac{\rho U V}{V} = \rho U = </math> the momentum per unit volume of fluid, we can write
 
<math>\frac{\partial (\rho V U)_{\small{P}}}{\partial t} = ({\rho}UUA)_w - ({\rho}UUA)_e + (pA)_w - (pA)_e + F</math>


=== Higher dimensions ===
Next we apply the assumptions of steady incompressible flow and zero external force.


* In 1D, we have been solving for velocity in the x direction by applying conservation of momentum in just one direction - the x axis.
<math>0 = ({\rho}UUA)_w - ({\rho}UUA)_e + (pA)_w - (pA)_e + F</math>
* In 2D, we must also solve for velocity in the y direction. For this, we will apply conservation of momentum in the y axis (note that we will also solve the x-momentum equation for the x-velocity as well).
* Similarly in 3D, we need to solve the momentum equation 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 - the process isn't important for now, although you should know that the spatial derivatives (e.g. partial d/dx eqn todo) are written with different notation, resulting in some terms called the "gradient", the "divergence", and the "laplacian" of a variable. A nice summary of these operations can be found at this link [todo].


=== Viscosity and turbulence ===
<math>0 = (UUA)_w - (UUA)_e + \frac{(pA)_w}{\rho} - \frac{(pA)_e}{\rho} + \frac{F}{\rho}</math>
We ignored the effects of viscosity, but this is straightforward to include in the momentum equation, simply by adding an extra force to the equation to account for viscosity.


We ignored turbulence in the momentum equation for this example, but real flows often exhibit turbulence ("random" fluctuations), especially at high speeds and low viscosities. Since the effect of turbulence is to increase mixing between layers, it can be modelled using a modified viscosity. This modified viscosity is the sum of the ordinary viscosity and the "turbulent viscosity" - this turbulent viscosity is calculated at every location in the mesh by a turbulence model. The turbulence models calculate some variables which quantify the turbulence at every location, and use those variables to calculate the turbulent viscosity. Some common turbulence models include (along with some broad pros & cons):
<math>0 = (UUA)_w - (UUA)_e + \frac{(pA)_w}{\rho} - \frac{(pA)_e}{\rho}</math>


* K-epsilon model: best for use when the cells are quite large ("coarse") near the walls, and for preliminary simulations.
At this point, depending on the location of cell P, we would continue by applying either
* K-omega model: requires smaller (finer) cells near the walls in order to work best. Generally more accurate than K-epsilon, but more computationally expensive, and requires a finer mesh. Is also more sensitive to boundary conditions than K-epsilon.
* Boundary conditions for faces on the boundary, or
* K-omega SST (Shear Stress Transport) model: a combination of the K-epsilon and K-omega models, with an added shear stress transport model. Generally more accurate than K-omega and K-epsilon, but requires a finer mesh and is more computationally expensive than K-epsilon.
* An interpolation scheme (often linear) for internal faces.
Then we would calculate the coefficients of the matrix, and either invert the matrix or use an iterative method to solve it. For brevity, I will not go through this process, as the general process is similar with the continuity equation.


=== Stability tips ===
Optional exercise: try applying these steps with
* Zero pressure gradient on the left boundary,
* Fixed pressure of 0 on the right boundary,
* The solution found in the previous section for U.
You should find that the pressure everywhere is 0 Pa.


* CFD simulations can sometimes have variables blow up to infinity or show unphysical oscillations in the flow, if care is not taken with the setup process. The stability of a simulation (i.e. how likely it is to avoid blowing up) can be influenced by choice of numerical scheme, the quality of the mesh, and the selection of boundary conditions.
===Optional: The energy equation ===
* In OpenFOAM, each operation has its own numerical scheme selection - for example:
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 a given point = rate of energy flow in - rate of energy flow out. I will not go through the derivation here - it is enough to know that the CFD solver performs an iterative approach to solve for temperature using the conservation of energy until the residual is low enough.
** Interpolation (we used the linear scheme for this operation in the example above),
** Gradient,
** Divergence - this one is particularly influential to stability,
** Laplacian.
* Details on the selection of which scheme to use is beyond the scope of this section, but a good rule-of-thumb is that less accurate schemes (e.g. upwind interpolation) tend to be more stable, and vice versa. Gradient limiters can also be used to add stability at the expense of accuracy.
* Meshing is not as straightforward as shown here, especially with more complicated geometries. Sometimes, poor quality meshes can cause problems with the simulation. There are numbers that meshing software can calculate to represent the quality of a mesh (called a mesh quality metric).
** For example, meshes with orthogonal cells perform better than those with non-orthogonal cells. In our example, all the cells were perfectly orthogonal - that is, the faces were perpendicular to the lines between the adjacent cell centroids.
** How close the cells are to having perpendicular faces is called the "orthogonal quality" in ANSYS.
** In OpenFOAM, the opposite of this is given as the "non-orthogonality angle".
** Other examples include cell skewness (lower is better) and cell aspect ratio (closer to 1 is better). It is important to ensure these quality metrics are not too bad, otherwise the stability and accuracy of the simulation can be affected.


=== Further reading ===
In the case of an incompressible flow, this isn't necessary, since density is constant. But for a compressible flow, the density at a given location is a function of temperature, so it is necessary to find the temperature everywhere in the fluid, so that the density can be found.
Hopefully this section gives an understandable explanation of roughly what most CFD solvers do. If you want to explore the topics mentioned here in more detail, here are some resources that might help: [todo]


= Fluid-structure interaction (FSI) =
===Optional: Higher dimensions ===
'''''IMPORTANT:''''' Although I hope that these instructions will be easily understood by anyone, <u>this is not a beginner's project!</u> 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]].
* In this example, we have been solving for velocity in the x direction by applying conservation of momentum in just one direction - the x axis.
* In 2D, we must also solve for velocity in the y direction. For this, we will apply conservation of momentum in the y axis.
* Similarly in 3D, we need to solve the momentum equation 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 - the process isn't important for now, although you should know that the spatial derivatives (e.g. <math>\frac{\partial}{\partial x}</math>) are written with different notation, resulting in some terms called the ''gradient'', the ''divergence'', and the ''laplacian'' of a variable. This is called vector calculus - I recommend [https://www.khanacademy.org/math/multivariable-calculus this Khan Academy course], if this is unfamiliar to you.


Simulations of this kind are achievable with the combination of the following free open-source software packages:
===Optional: Viscosity and turbulence===
We ignored the effects of viscosity, but this is straightforward to include in the momentum equation, simply by adding an extra force to the equation to account for viscosity.


* OpenFOAM ([https://www.openfoam.com/ website], [https://develop.openfoam.com/Development/openfoam/-/wikis/precompiled/debian install on Ubuntu]) - a C/C++ library for the solution of fluid mechanics problems formulated with the Finite Volume Method (FVM),  
We ignored turbulence in the momentum equation for this example, but real flows often exhibit turbulence ("random" fluctuations), especially at high speeds and low viscosities. Since the effect of turbulence is to increase mixing between layers, it can be modelled using a modified viscosity. This modified viscosity is the sum of the ordinary viscosity and the ''turbulent viscosity'' - this turbulent viscosity is calculated at every location in the mesh by a turbulence model. The turbulence models calculate some variables which quantify the turbulence at every location, and use those variables to calculate the turbulent viscosity. Some common turbulence models include (along with some broad pros & cons):
** 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 ([https://fenicsproject.org/ website], [https://fenicsproject.org/download/ 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 ([https://precice.org/index.html website], [https://precice.org/quickstart.html 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 [https://precice.org/adapter-openfoam-overview.html this page for OpenFOAM] and [https://precice.org/adapter-fenics.html this page for FEniCS] from the preCICE website for further information.


== Step-by-step instructions ==
*K-epsilon model: best for use when the cells are quite large ("coarse") near the walls, and for preliminary simulations.
'''''IMPORTANT:''''' Although I include the links to the official websites for Ubuntu installation instructions, I recommend you <u>follow my instructions for installation</u> (unless you are a masochist), since the official instructions contain some mistakes and omissions.
*K-omega model: requires smaller (finer) cells near the walls in order to work best. Generally more accurate than K-epsilon, but more computationally expensive, and requires a finer mesh. Is also more sensitive to boundary conditions than K-epsilon.
*K-omega SST (Shear Stress Transport) model: a combination of the K-epsilon and K-omega models, with an added shear stress transport model. Generally more accurate than K-omega and K-epsilon, but requires a finer mesh and is more computationally expensive than K-epsilon.


=== Install Ubuntu ===
===Optional: Stability tips===
If you are running a Windows machine (Windows 10 or above), I recommend installing [https://docs.microsoft.com/en-us/windows/wsl/install-manual 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 [https://stackoverflow.com/questions/61110603/how-to-set-up-working-x11-forwarding-on-wsl2 this page] and [https://newbedev.com/how-to-set-up-working-x11-forwarding-on-wsl2 this page] for further information)


=== Install OpenFOAM on Ubuntu ===
*CFD simulations can sometimes have variables blow up to infinity or show unphysical oscillations in the flow, if care is not taken with the setup process. The stability of a simulation (i.e. how likely it is to avoid blowing up) can be influenced by choice of numerical scheme, the quality of the mesh, and the selection of boundary conditions.
The Ubuntu installation instructions as of 31/12/2021 are summarised below (adapted from [https://develop.openfoam.com/Development/openfoam/-/wikis/precompiled/debian this website]):<syntaxhighlight lang="bash">
*In OpenFOAM, each operation has its own numerical scheme selection - for example:
curl -s https://dl.openfoam.com/add-debian-repo.sh | sudo bash
**Interpolation (we used the linear scheme for this operation in the example above),
sudo apt-get install openfoam2012-default
**Gradient,
</syntaxhighlight>
**Divergence - this one is particularly influential to stability,
Ensure that your ~/.bashrc file has the following line appended (with XXXX replaced with your version number):<syntaxhighlight lang="bash">
**Laplacian.
source /usr/lib/openfoam/openfoamXXXX/etc/bashrc
*Details on the selection of which scheme to use is beyond the scope of this section, but a good rule-of-thumb is that less accurate schemes (e.g. upwind interpolation) tend to be more stable, and vice versa. Gradient limiters can also be used to add stability at the expense of accuracy.
</syntaxhighlight>
*Meshing is not as straightforward as shown here, especially with more complicated geometries. Sometimes, poor quality meshes can cause problems with the simulation. There are numbers that meshing software can calculate to represent the quality of a mesh (called a mesh quality metric).
Restart your Ubuntu. Now check the installation completed correctly by running a tutorial case:<syntaxhighlight lang="bash">
**For example, meshes with orthogonal cells perform better than those with non-orthogonal cells. In our example, all the cells were perfectly orthogonal - that is, the faces were perpendicular to the lines between the adjacent cell centroids.
cd ~
**How close the cells are to having perpendicular faces is called the ''orthogonal quality'' in ANSYS.
mkdir -p OpenFOAM-sims/tutorials
**In OpenFOAM, the opposite of this is given as the ''non-orthogonality angle''.
cp -r $FOAM_TUTORIALS/incompressible OpenFOAM-sims/tutorials
**Other examples include cell skewness (lower is better) and cell aspect ratio (closer to 1 is better). It is important to ensure these quality metrics are not too bad, otherwise the stability and accuracy of the simulation can be affected.
cd OpenFOAM-sims/tutorials/incompressible/icoFoam/cavity/cavity/
blockMesh > log.blockMesh
icoFoam > log.icoFoam
</syntaxhighlight>
Note: I use OpenFOAMv2012, but later versions should also work.


=== Install FEniCS on Ubuntu ===
===Optional: Common tripping points===
The Ubuntu installation instructions as of 31/12/2021 are summarised below (adapted from [https://fenicsproject.org/download/ this website]):<syntaxhighlight lang="bash">
The solution vector <math>\boldsymbol{U}</math> does not store the components of velocity, i.e. <math>U_x, U_y, U_z</math>, but instead stores the velocities at each cell centroid, i.e. <math>U_1, U_2, U_3, U_4</math>. Indeed, the solution vector <math>\boldsymbol{U}</math> summarises the x-velocity at every point in our domain, while hypothetically the scalars <math>U_x, U_y, U_z</math> would only specify the 3D velocity at one point. Since we are currently only concerned with x-velocity, this might be more clear, but in higher dimensions it is important not to get confused.
sudo apt-get install software-properties-common
sudo add-apt-repository ppa:fenics-packages/fenics
sudo apt-get update
sudo apt-get install fenics
</syntaxhighlight>


=== Install preCICE on Ubuntu ===
The term <math>p/\rho</math> is sometimes called the ''kinematic pressure''. When an incompressible fluid is used, it is sometimes what the CFD software (e.g. OpenFOAM) will confusingly call p. So when finding the pressure drop, you may have to multiply by density. The best way to check whether you need to multiply by density, is to look at the units of p that the CFD software gives you. Ordinary pressure has units <math>kgm^{-1}s^{-2}</math>, while kinematic pressure has units <math>m^{2}s^{-2}</math>.
The Ubuntu installation instructions as of 31/12/2021 are summarised below (adapted from [https://precice.org/quickstart.html this website]): <syntaxhighlight lang="bash">
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


===Optional: Further reading===
Hopefully this has been an understandable / intuitive explanation of roughly what most CFD solvers do. If you want to explore the topics mentioned here in more detail, here are some resources that might help:


</syntaxhighlight>
* '''RECOMMENDED: An Introduction to Computational Fluid Dynamics: The Finite Volume Method''' - commonly known as just "Versteeg" - https://books.google.co.uk/books/about/An_Introduction_to_Computational_Fluid_D.html?id=RvBZ-UMpGzIC
Now either install the Ubuntu package directly with the instructions below (which did not work for me)...<syntaxhighlight lang="bash">
** Absolutely brilliant. Please read either this or Moukalled before anything else.
wget https://github.com/precice/precice/releases/download/v2.3.0/libprecice2_2.3.0_focal.deb
**It is usually behind a paywall, but you can sometimes find an old freely available copy by simply [[google:versteeg+cfd+filetype:pdf|googling]] <code>versteeg cfd filetype:pdf</code>
sudo apt install ./libprecice2_2.3.0_focal.deb
* '''RECOMMENDED: The Finite Volume Method in Computational Fluid Dynamics''' - commonly known as just "Moukalled" - https://link.springer.com/book/10.1007/978-3-319-16874-6
</syntaxhighlight>
** Slower paced than Versteeg - recommended for those who are a little shaky on their fluid dynamics knowledge. Also contains some more OpenFOAM-specific knowledge.
... or install from source:<syntaxhighlight lang="bash">
**It is usually behind a paywall, but you can sometimes find an old freely available copy by simply [[google:moukalled+cfd+filetype:pdf|googling]] <code>moukalled cfd filetype:pdf</code>
tar -xzf v2.3.0.tar.gz
* Fluid Mechanics 101 - YouTube channel: https://www.youtube.com/channel/UCcqQi9LT0ETkRoUu8eYaEkg/
cd precice-2.3.0
** Has some very well explained and illustrated videos on a random selection of CFD topics.
mkdir build
* ANSYS Fluent Theory Guide - find it with [[google:ansys+fluent+theory+guide+filetype:pdf|google]].
cd build
** Usually a bit easier to find a copy of, but isn't ideal for someone who isn't already familiar with CFD and the finite volume method. See Versteeg for that.
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!
* OpenFOAM: A little user manual - https://github.com/ParticulateFlow/OSCCAR-doc/blob/master/openFoamUserManual_PFM.pdf
make -j 4
** More specific to OpenFOAM and its code syntax, but contains some discussion of the fundamental algorithms which are used everywhere, e.g. SIMPLE, PISO, PIMPLE, as well as some turbulence models.
</syntaxhighlight>
* OpenFOAM wiki - http://openfoamwiki.net/index.php/Main_Page
Test that the cmake command worked as expected:<syntaxhighlight lang="bash">
** Contains some useful derivations of equations as well as tips & tricks, but is a bit hit-and-miss in quality and ease of navigation.
cd ~/FSI/precice-2.3.0/build
ctest --output-on-failure
</syntaxhighlight>
If no errors were shown, install the software and test it worked:<syntaxhighlight lang="bash">
make install
make test_install
</syntaxhighlight>
Note: I used preCICE v2.3.0, but later versions should also work.


=== Install the OpenFOAM-preCICE adapter ===
= Fluid-structure interaction (FSI)=
The Ubuntu installation instructions as of 31/12/2021 are summarised below (adapted from [https://precice.org/adapter-openfoam-get.html this website]):<syntaxhighlight lang="bash">
[[FSI installation instructions|FSI installation instructions here]]
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
</syntaxhighlight>
Note: this is for OpenFOAMv2012. Please see the preCICE website for information on adapters for later versions.


=== Install the FEniCS-preCICE adapter ===
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]].
The Ubuntu installation instructions as of 31/12/2021 are summarised below (adapted from [https://precice.org/adapter-fenics.html this website]). This requires Python3, which should already be installed if you followed my previous instructions: <syntaxhighlight lang="bash">
(optional, should not need this line) python3 -m pip install scipy
python3 -m pip install --user fenicsprecice
</syntaxhighlight>


=== Test the installation ran correctly ===
Simulations of this kind are achievable with the combination of the following free open-source software packages:
... by running [https://precice.org/tutorials-perpendicular-flap.html this tutorial case]: <syntaxhighlight lang="bash">
cd ~/FSI
git clone --branch=master --depth 1 https://github.com/precice/tutorials.git
cd ~/FSI/tutorials/perpendicular-flap/fluid-openfoam
./run.sh
</syntaxhighlight>
The command line output should have paused after this output. If it did not pause, something went wrong. <syntaxhighlight>
---[precice]  I am participant "Fluid"
---[precice]  Setting up master communication to coupling partner/s
</syntaxhighlight>
Now open another Ubuntu terminal, and run these commands:<syntaxhighlight lang="bash">
cd ~/FSI/tutorials/perpendicular-flap/solid-fenics
./run.sh
</syntaxhighlight>
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:<syntaxhighlight lang="bash">
cd ~/FSI/tutorials/perpendicular-flap/fluid-openfoam
paraview.exe fluid-openfoam.foam
</syntaxhighlight>
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.
*OpenFOAM ([https://www.openfoam.com/ website], [[OpenFOAM installation instructions|OpenFOAM installation instructions here]]) - 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 ([https://fenicsproject.org/ website]) - a Python interface to the DOLFIN C++ library for the solution of solid mechanics problems formulated with the Finite Element Method (FEM), and
*preCICE ([https://precice.org/index.html website]) - 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 [https://precice.org/adapter-openfoam-overview.html this page for OpenFOAM] and [https://precice.org/adapter-fenics.html this page for FEniCS] from the preCICE website for further information.


=== TODO: ===
== Step-by-step instructions==
* ParaView install instructions
'''''IMPORTANT:''''' Although I hope that these instructions will be easily understood by anyone, <u>this is not a beginner's project!</u> If you want a better introduction to CFD, please read the crash-course section above, or any of the linked resources.
* 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


*
This section has moved to [[FSI installation instructions|a separate page]].

Latest revision as of 13:16, 12 February 2022

This page is maintained 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, usually with the Finite Volume Method.

Essentially F=ma, but with many simultaneous equations.

OpenFOAM installation instructions here

Getting Started with CFD

OpenFOAM installation instructions here

Skip this section if you dislike broad-strokes philosophical discussion.

The path you should take for CFD differs depending on what your end goal is, and what budget you have. The main choice is between OpenFOAM and a commercial package (ANSYS Fluent / CFX, STAR-CCM, etc).

OpenFOAM Pros vs commercial packages:

  • Free
  • Open Source
  • You can edit any part of the software, any setting, any way you want to.
  • No limits, no "paid version only" features. (ANSYS limits you to 2 CPU threads and 100k cells for simulation)

OpenFOAM Cons vs commercial packages:

  • Difficult to install
  • Difficult to use
  • Difficult to learn - less of a learning curve, more like a cliff.
  • No graphical user interface
  • Mesh generation software is terrible. IMO, the best way to generate a mesh in OpenFOAM, is not to generate a mesh in OpenFOAM
    • Pro tip: use ANSYS meshing software (either the mechanical mesher, or the Fluent mesher) to generate CFD meshes. It will be 10x quicker (including setup time), take 5x less RAM, and be a higher quality than what you could achieve with OpenFOAM. Then import it into OpenFOAM with fluentMeshToFoam <options> <.msh file>. ANSYS does not put a cap on the mesh size you can generate, only the mesh size you can simulate is limited.
    • There are other free & open source meshing programs, e.g. Gmsh, but I have not tried them.

OpenFOAM throws you into the deep end, and forces you to learn nearly everything at once. If you stick through to the end of the learning cliff, you will understand much more about CFD than you would have otherwise learned with a commercial package.

If you just need to do one or two small sims, a commercial package e.g. ANSYS is your best choice. (Just not SOLIDWORKS CFD simulations).

The remainder of this wiki is dedicated to free and open source software, including OpenFOAM.

Crash course on the Finite Volume Method (FVM)

There are many resources for learning how CFD works properly (linked below), but this is a crash course. This section is designed for people who would like to find out a bit more about how CFD works, beyond a "black-box" understanding, but without going too far in depth with all the details. So, here are the bare minimum requirements for understanding roughly what is going on under the hood in any CFD solver. 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.

Sections marked "Optional" are not required for a basic level of understanding of how the software works, but give a bit more detail for those who want it.

Overview

The most common method of doing CFD is with the Finite Volume Method (FVM). An informal derivation of some flow equations along with application of the FVM to them is in this section.

It is also possible to do CFD with the Finite Element Method (FEM), which is more typically associated with Structural Analysis / Elastodynamics problems, but we won't discuss that here.

CFD solvers work 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 uniform pipe, with a cross-section shown in the left image. Now, let us mesh the domain into four cells (split the pipe into 4 boxes), image right.

Note: the cells here are 3D - we are looking at a 2D cross section. We will assume that the velocity is purely in the x direction. Consider an arbitrary cell within this mesh labelled [math]\displaystyle{ P }[/math], with its western and eastern neighbours labelled [math]\displaystyle{ W }[/math] and [math]\displaystyle{ E }[/math] respectively. The faces between [math]\displaystyle{ W }[/math] & [math]\displaystyle{ P }[/math] and [math]\displaystyle{ P }[/math] & [math]\displaystyle{ E }[/math] are labelled [math]\displaystyle{ w }[/math] and [math]\displaystyle{ c }[/math]. The following abbreviations for flow variables will be used:

  • [math]\displaystyle{ x }[/math]: Position,
  • [math]\displaystyle{ \rho }[/math] (rho): Density,
  • [math]\displaystyle{ p }[/math]: Pressure,
  • [math]\displaystyle{ U }[/math]: Velocity (in the x-axis),
  • [math]\displaystyle{ V }[/math]: Volume of the cell,
  • [math]\displaystyle{ A }[/math]: Area of a face,
  • Any variable with a subscript corresponds to the value of that variable at that location. E.g. [math]\displaystyle{ U_{\small{W}} }[/math] is the velocity of the fluid at [math]\displaystyle{ W }[/math]. The area of face [math]\displaystyle{ w }[/math] is [math]\displaystyle{ A_w }[/math].

The continuity equation tells us that mass is not created or destroyed. For an arbitrary cell P, knowing that the mass flow rate through a given face is [math]\displaystyle{ {\rho}UA }[/math], we can write

Rate of change of mass in [math]\displaystyle{ P }[/math] = Mass flow rate in - Mass flow rate out

[math]\displaystyle{ \frac{\partial (\rho V)_{\small{P}}}{\partial t} = ({\rho}UA)_{w} - ({\rho}UA)_{e} }[/math]

If we assume an incompressible (constant density) fluid, e.g. water, we can divide by the density term [math]\displaystyle{ \rho }[/math] on both sides. If we assume a steady flow, then nothing changes with time, i.e. all time derivatives are zero. The equation becomes

[math]\displaystyle{ 0 = ({\rho}UA)_w - ({\rho}UA)_e }[/math]

[math]\displaystyle{ 0 = (UA)_w - (UA)_e }[/math] - (1)

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, as in the image above. Note: interpolation isn't unique to the variable [math]\displaystyle{ U }[/math] - it can be applied to all flow variables.

In our case, the cells and faces are all perfectly uniformly spaced, so if we use linear interpolation to work out the face values,

[math]\displaystyle{ U_w = {1 \over 2} (U_{\small{W}} + U_{\small{P}}), \quad U_e = {1 \over 2} (U_{\small{P}} + U_{\small{E}}) }[/math]

This is the same as saying that the face values are the average of the neighbouring cell values. Let's use a more general equation to represent linear interpolation.

[math]\displaystyle{ U_w = i_{\small{W}}U_{\small{W}} + (1-i_{\small{W}})U_{\small{P}}, \quad U_e = i_{\small{E}}U_{\small{P}} + (1-i_{\small{E}})U_{\small{E}} }[/math]

where [math]\displaystyle{ i_{\small{W}} = \frac{x_w - x_{\small{W}}}{x_{\small{P}} - x_{\small{W}}} }[/math] and [math]\displaystyle{ i_{\small{E}} = \frac{x_e - x_{\small{P}}}{x_{\small{E}} - x_{\small{P}}} }[/math]

Don't be intimidated by the complicated-looking fractions! [math]\displaystyle{ i_{\small{W}} }[/math] simply represents how far the face [math]\displaystyle{ w }[/math] is from cell [math]\displaystyle{ W }[/math], divided by the distance between [math]\displaystyle{ W }[/math] and [math]\displaystyle{ P }[/math]. This way, if the face is exactly half-way between [math]\displaystyle{ W }[/math] and [math]\displaystyle{ P }[/math], then [math]\displaystyle{ i_{\small{W}} = {1 \over 2} }[/math], giving the same equation as before. If instead, the face was closer to [math]\displaystyle{ W }[/math] than to [math]\displaystyle{ P }[/math], then [math]\displaystyle{ i_{\small{W}} \lt {1 \over 2} }[/math], and vice versa.

After this point, the continuity equation will only work for cells that aren't on the boundary. We also factor out the areas of each face, since they are just constants. The continuity equation becomes

[math]\displaystyle{ 0 = (U_{\small{W}}i_{\small{W}} + U_{\small{P}}(1-i_{\small{W}})) \cdot A_w - (U_{\small{P}}i_{\small{E}} + U_{\small{E}}(1-i_{\small{E}})) \cdot A_e }[/math]

We now expand and simplify the equation, combining the constants into single variables.

[math]\displaystyle{ 0 = U_{\small{W}}i_{\small{W}}A_w + U_{\small{P}}(1-i_{\small{W}})A_w - U_{\small{P}}i_{\small{E}}A_e - U_{\small{E}}(1-i_{\small{E}})A_e }[/math]

[math]\displaystyle{ 0 = U_{\small{W}}k_{\small{P,W}} + U_{\small{P}}k_{\small{P,P}} + U_{\small{E}}k_{\small{P,E}} }[/math]

where

  • [math]\displaystyle{ k_{\small{P,W}} = i_{\small{W}}A_w }[/math]
  • [math]\displaystyle{ k_{\small{P,P}} = (1-i_{\small{W}})A_w - i_{\small{E}}A_e }[/math]
  • [math]\displaystyle{ k_{\small{P,E}} = -(1-i_{\small{E}})A_e }[/math].

We can now apply this equation for cells 2 and 3. However, we can't apply this equation to cells 1 and 4, as they are on the boundary of the domain - the linear interpolation step does not work on the boundary. Without knowing what the value of U is at the boundaries, we cannot solve the continuity equation. Therefore we must apply some special treatment at the boundaries (called boundary conditions) in order to solve this issue. Boundary conditions are specified by the CFD user. We will use these boundary conditions:

  • Flow velocity on the left boundary is fixed at 5 m/s (called a Dirichlet boundary - mnemonic: kind of sounds like "directly").
    • Plugging this condition into equation (1) gives the equation [math]\displaystyle{ 0 = 5A_w - U_{\small{P}}i_{\small{E}}A_e - U_{\small{E}}(1-i_{\small{E}})A_e }[/math]
  • The gradient of velocity [math]\displaystyle{ \frac {\partial U}{\partial x} }[/math] on the right boundary is zero - so the velocity on the right boundary is the same as the velocity at cell centroid 4.
    • Plugging this condition into equation (1) gives the equation [math]\displaystyle{ 0 = U_{\small{W}}i_{\small{W}}A_w + U_{\small{P}}(1-i_{\small{W}})A_w - U_{\small{P}}A_e }[/math]
  • The constants are again combined into variables of the form [math]\displaystyle{ k_{a,b} }[/math]. Here is a summary of all the equations found so far:
  • Cell 1: [math]\displaystyle{ 0 = 5A_w + U_1k_{1,1} + U_2k_{1,2} }[/math]
  • Cell 2: [math]\displaystyle{ 0 = 0 + U_1k_{2,1} + U_2k_{2,2} + U_3k_{2,3} }[/math]
  • Cell 3: [math]\displaystyle{ 0 = 0 + U_2k_{3,2} + U_3k_{3,3} + U_4k_{3,4} }[/math]
  • Cell 4: [math]\displaystyle{ 0 = 0 + U_3k_{4,3} + U_4k_{4,4} }[/math]

We now have a set of four simultaneous equations to solve for four unknowns (the velocity at each cell centroid).

Let us simplify the notation a bit using matrices and vectors.

[math]\displaystyle{ \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} 5A_w \\ 0 \\ 0 \\ 0 \end{bmatrix} + \begin{bmatrix} k_{1,1} & k_{1,2} & 0 & 0 \\ k_{2,1} & k_{2,2} & k_{2,3} & 0 \\ 0 & k_{3,2} & k_{3,3} & k_{3,4} \\ 0 & 0 & k_{4,3} & k_{4,4} \end{bmatrix} \begin{bmatrix} U_1 \\ U_2 \\ U_3 \\ U_4 \end{bmatrix} }[/math]

[math]\displaystyle{ \begin{bmatrix} k_{1,1} & k_{1,2} & 0 & 0 \\ k_{2,1} & k_{2,2} & k_{2,3} & 0 \\ 0 & k_{3,2} & k_{3,3} & k_{3,4} \\ 0 & 0 & k_{4,3} & k_{4,4} \end{bmatrix} \begin{bmatrix} U_1 \\ U_2 \\ U_3 \\ U_4 \end{bmatrix} = \begin{bmatrix} -5A_w \\ 0 \\ 0 \\ 0 \end{bmatrix} }[/math].

Abbreviated with a matrix of coefficients [math]\displaystyle{ \boldsymbol{M} }[/math], a solution vector [math]\displaystyle{ \boldsymbol{U} }[/math], and a right-hand-side vector [math]\displaystyle{ \boldsymbol{B} }[/math]:

[math]\displaystyle{ \boldsymbol{MU=B} }[/math], where

[math]\displaystyle{ \boldsymbol{M}= \begin{bmatrix} k_{1,1} & k_{1,2} & 0 & 0 \\ k_{2,1} & k_{2,2} & k_{2,3} & 0 \\ 0 & k_{3,2} & k_{3,3} & k_{3,4} \\ 0 & 0 & k_{4,3} & k_{4,4} \end{bmatrix} }[/math], [math]\displaystyle{ \boldsymbol{U} = \begin{bmatrix} U_1 \\ U_2 \\ U_3 \\ U_4 \end{bmatrix} }[/math], and [math]\displaystyle{ \boldsymbol{B} = \begin{bmatrix} -5A_w \\ 0 \\ 0 \\ 0 \end{bmatrix} }[/math]

At this point a CFD solver would calculate each of the coefficients of the matrix and RHS vector. Here we have a uniform cell spacing and a uniform face area of 1.

[math]\displaystyle{ \boldsymbol{M}= \begin{bmatrix} -0.5 & -0.5 & 0 & 0 \\ 0.5 & 0 & -0.5 & 0 \\ 0 & 0.5 & 0 & -0.5 \\ 0 & 0 & 0.5 & -0.5 \end{bmatrix} }[/math]

[math]\displaystyle{ \boldsymbol{B} = \begin{bmatrix} -5 \\ 0 \\ 0 \\ 0 \end{bmatrix} }[/math]

We could now invert the matrix [math]\displaystyle{ \boldsymbol{M} }[/math], and calculate the solution vector [math]\displaystyle{ \boldsymbol{U} }[/math] using [math]\displaystyle{ \boldsymbol{U=M^{-1}B} }[/math].

[math]\displaystyle{ \boldsymbol{M^{-1}}= \begin{bmatrix} -1 & 1 & -1 & 1 \\ -1 & -1 & 1 & -1 \\ -1 & -1 & -1 & 1 \\ -1 & -1 & -1 & -1 \end{bmatrix} }[/math] This yields a solution of U = 5 m/s everywhere. Hooray!

Think: does this result make sense, given our assumptions?

Why matrix inversion rarely happens

The matrix [math]\displaystyle{ \boldsymbol{M} }[/math] 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 [math]\displaystyle{ \boldsymbol{MU=B} }[/math], where they start with a guess for [math]\displaystyle{ \boldsymbol{U} }[/math], and improve their guess with successive iterations of an algorithm until the solution vector [math]\displaystyle{ \boldsymbol{U} }[/math] 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 specifies what residual value qualifies as "close enough".

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, noting that

[math]\displaystyle{ UA = }[/math] the rate of flow of volume through a face, and

[math]\displaystyle{ \frac{mU}{V} = \frac{\rho U V}{V} = \rho U = }[/math] the momentum per unit volume of fluid, we can write

[math]\displaystyle{ \frac{\partial (\rho V U)_{\small{P}}}{\partial t} = ({\rho}UUA)_w - ({\rho}UUA)_e + (pA)_w - (pA)_e + F }[/math]

Next we apply the assumptions of steady incompressible flow and zero external force.

[math]\displaystyle{ 0 = ({\rho}UUA)_w - ({\rho}UUA)_e + (pA)_w - (pA)_e + F }[/math]

[math]\displaystyle{ 0 = (UUA)_w - (UUA)_e + \frac{(pA)_w}{\rho} - \frac{(pA)_e}{\rho} + \frac{F}{\rho} }[/math]

[math]\displaystyle{ 0 = (UUA)_w - (UUA)_e + \frac{(pA)_w}{\rho} - \frac{(pA)_e}{\rho} }[/math]

At this point, depending on the location of cell P, we would continue by applying either

  • Boundary conditions for faces on the boundary, or
  • An interpolation scheme (often linear) for internal faces.

Then we would calculate the coefficients of the matrix, and either invert the matrix or use an iterative method to solve it. For brevity, I will not go through this process, as the general process is similar with the continuity equation.

Optional exercise: try applying these steps with

  • Zero pressure gradient on the left boundary,
  • Fixed pressure of 0 on the right boundary,
  • The solution found in the previous section for U.

You should find that the pressure everywhere is 0 Pa.

Optional: 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 a given point = rate of energy flow in - rate of energy flow out. I will not go through the derivation here - it is enough to know that the CFD solver performs an iterative approach to solve for temperature using the conservation of energy until the residual is low enough.

In the case of an incompressible flow, this isn't necessary, since density is constant. But for a compressible flow, the density at a given location is a function of temperature, so it is necessary to find the temperature everywhere in the fluid, so that the density can be found.

Optional: Higher dimensions

  • In this example, we have been solving for velocity in the x direction by applying conservation of momentum in just one direction - the x axis.
  • In 2D, we must also solve for velocity in the y direction. For this, we will apply conservation of momentum in the y axis.
  • Similarly in 3D, we need to solve the momentum equation 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 - the process isn't important for now, although you should know that the spatial derivatives (e.g. [math]\displaystyle{ \frac{\partial}{\partial x} }[/math]) are written with different notation, resulting in some terms called the gradient, the divergence, and the laplacian of a variable. This is called vector calculus - I recommend this Khan Academy course, if this is unfamiliar to you.

Optional: Viscosity and turbulence

We ignored the effects of viscosity, but this is straightforward to include in the momentum equation, simply by adding an extra force to the equation to account for viscosity.

We ignored turbulence in the momentum equation for this example, but real flows often exhibit turbulence ("random" fluctuations), especially at high speeds and low viscosities. Since the effect of turbulence is to increase mixing between layers, it can be modelled using a modified viscosity. This modified viscosity is the sum of the ordinary viscosity and the turbulent viscosity - this turbulent viscosity is calculated at every location in the mesh by a turbulence model. The turbulence models calculate some variables which quantify the turbulence at every location, and use those variables to calculate the turbulent viscosity. Some common turbulence models include (along with some broad pros & cons):

  • K-epsilon model: best for use when the cells are quite large ("coarse") near the walls, and for preliminary simulations.
  • K-omega model: requires smaller (finer) cells near the walls in order to work best. Generally more accurate than K-epsilon, but more computationally expensive, and requires a finer mesh. Is also more sensitive to boundary conditions than K-epsilon.
  • K-omega SST (Shear Stress Transport) model: a combination of the K-epsilon and K-omega models, with an added shear stress transport model. Generally more accurate than K-omega and K-epsilon, but requires a finer mesh and is more computationally expensive than K-epsilon.

Optional: Stability tips

  • CFD simulations can sometimes have variables blow up to infinity or show unphysical oscillations in the flow, if care is not taken with the setup process. The stability of a simulation (i.e. how likely it is to avoid blowing up) can be influenced by choice of numerical scheme, the quality of the mesh, and the selection of boundary conditions.
  • In OpenFOAM, each operation has its own numerical scheme selection - for example:
    • Interpolation (we used the linear scheme for this operation in the example above),
    • Gradient,
    • Divergence - this one is particularly influential to stability,
    • Laplacian.
  • Details on the selection of which scheme to use is beyond the scope of this section, but a good rule-of-thumb is that less accurate schemes (e.g. upwind interpolation) tend to be more stable, and vice versa. Gradient limiters can also be used to add stability at the expense of accuracy.
  • Meshing is not as straightforward as shown here, especially with more complicated geometries. Sometimes, poor quality meshes can cause problems with the simulation. There are numbers that meshing software can calculate to represent the quality of a mesh (called a mesh quality metric).
    • For example, meshes with orthogonal cells perform better than those with non-orthogonal cells. In our example, all the cells were perfectly orthogonal - that is, the faces were perpendicular to the lines between the adjacent cell centroids.
    • How close the cells are to having perpendicular faces is called the orthogonal quality in ANSYS.
    • In OpenFOAM, the opposite of this is given as the non-orthogonality angle.
    • Other examples include cell skewness (lower is better) and cell aspect ratio (closer to 1 is better). It is important to ensure these quality metrics are not too bad, otherwise the stability and accuracy of the simulation can be affected.

Optional: Common tripping points

The solution vector [math]\displaystyle{ \boldsymbol{U} }[/math] does not store the components of velocity, i.e. [math]\displaystyle{ U_x, U_y, U_z }[/math], but instead stores the velocities at each cell centroid, i.e. [math]\displaystyle{ U_1, U_2, U_3, U_4 }[/math]. Indeed, the solution vector [math]\displaystyle{ \boldsymbol{U} }[/math] summarises the x-velocity at every point in our domain, while hypothetically the scalars [math]\displaystyle{ U_x, U_y, U_z }[/math] would only specify the 3D velocity at one point. Since we are currently only concerned with x-velocity, this might be more clear, but in higher dimensions it is important not to get confused.

The term [math]\displaystyle{ p/\rho }[/math] is sometimes called the kinematic pressure. When an incompressible fluid is used, it is sometimes what the CFD software (e.g. OpenFOAM) will confusingly call p. So when finding the pressure drop, you may have to multiply by density. The best way to check whether you need to multiply by density, is to look at the units of p that the CFD software gives you. Ordinary pressure has units [math]\displaystyle{ kgm^{-1}s^{-2} }[/math], while kinematic pressure has units [math]\displaystyle{ m^{2}s^{-2} }[/math].

Optional: Further reading

Hopefully this has been an understandable / intuitive explanation of roughly what most CFD solvers do. If you want to explore the topics mentioned here in more detail, here are some resources that might help:

Fluid-structure interaction (FSI)

FSI installation instructions here

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, OpenFOAM installation instructions here) - 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) - 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) - 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 hope that these instructions will be easily understood by anyone, this is not a beginner's project! If you want a better introduction to CFD, please read the crash-course section above, or any of the linked resources.

This section has moved to a separate page.