Computational Fluid Dynamics (CFD) - Seminar Report

 Computational Fluid Dynamics (CFD)
This is a quick-and-dirty introduction to the basic concepts underlying CFD. The concepts are illustrated by applying them to simple 1D model problems.  We’ll invoke these concepts while performing "case studies" in FLUENT. Happily for us, these model-problem concepts extend to the more general situations in the case studies in most instances. Since we’ll keep returning to these concepts while performing the FLUENT case studies, it’s worth your time to understand and digest these concepts. 

The Need for CFD 
Applying the fundamental laws of mechanics to a fluid gives the governing equations for a fluid. The conservation of mass equation is 
∂t + ∇ · (ρV ) = 0 
and the conservation of momentum equation is 

∂t + ρ(V · ∇)V = −∇p + ρg + ∇ · τij 
These equations along with the conservation of energy equation form a set of coupled, nonlinear partial differential equations.  It is not possible to solve these equations analytically for most engineering problems. 
However, it is possible to obtain approximate computer-based solutions to the governing equations for a variety of engineering problems. This is the subject matter of Computational Fluid Dynamics (CFD). 

Applications of CFD 
CFD is useful in a wide variety of applications and here we note a few to give you an idea of 
its use in industry. The simulations shown below have been performed using the FLUENT 
CFD can be used to simulate the flow over a vehicle. For instance, it can be used to study 
the interaction of propellers or rotors with the aircraft fuselage The following figure shows 
the prediction of the pressure field induced by the interaction of the rotor with a helicopter 
fuselage in forward flight. Rotors and propellers can be represented with models of varying 

The temperature distribution obtained from a CFD analysis of a mixing manifold is shown below. This mixing manifold is part of the passenger cabin ventilation system on the Boeing 767.  The CFD analysis showed the effectiveness of a simpler manifold design without the need for field testing. 

Bio-medical engineering is a rapidly growing field and uses CFD to study the circulatory and 
respiratory systems. The following figure shows pressure contours and a cutaway view that 
reveals velocity vectors in a blood pump that assumes the role of heart in open-heart surgery. 

CFD is attractive to industry since it is more cost-effective than physical testing. However, one must note that complex flow simulations are challenging and error-prone and it takes a lot of engineering expertise to obtain validated solutions. 

The Strategy of CFD 
Broadly, the strategy of CFD is to replace the continuous problem domain with a discrete domain using a grid. In the continuous domain, each flow variable is defined at every point in the domain. For instance, the pressure p in the continuous 1D domain shown in the figure below would be given as 
p=p(x), 0<x<1 
In the discrete domain, each flow variable is defined only at the grid points.  So, in the discrete domain shown below, the pressure would be defined only at the N grid points. 
pi = p(xi),  i = 1, 2, . . . , N 

Continuous Domain Discrete Domain
0 £ x £ 1 x = x
1, x2, ¼,xN
x=1 x x x
Coupled PDEs + boundary 
conditions in continuous 
1 i N
Grid point
Coupled algebraic eqs. in discrete variables 

In a CFD solution, one would directly solve for the relevant flow variables only at the grid 
points. The values at other locations are determined by interpolating the values at the grid 
The governing partial differential equations and boundary conditions are defined in terms 
of the continuous variables p, V  etc. One can approximate these in the discrete domain in 
terms of the discrete variables pi, Vi etc.  The discrete system is a large set of coupled, 
algebraic equations in the discrete variables. Setting up the discrete system and solving it (which is a matrix inversion problem) involves a very large number of repetitive calculations, a task we humans palm over to the digital computer. 

This idea can be extended to any general problem domain. The following figure shows the grid used for solving the flow over an airfoil. We’ll take a closer look at this airfoil grid soon while discussing the finite-volume method. 

Discretization Using the Finite-Difference Method 
To keep the details simple, we will illustrate the fundamental ideas underlying CFD by applying them to the following simple 1D equation:
dx + um = 0;

0≤x≤1;  u(0)=1 (1) 
We’ll first consider the case where m = 1 when the equation is linear. We’ll later consider the m = 2 case when the equation is nonlinear. 
We’ll derive a discrete representation of the above equation with m = 1 on the following 

x x x x
1=0 2=1/3 3=2/3 4=1
This grid has four equally-spaced grid points with Δx being the spacing between successive points. Since the governing equation is valid at any grid point, we have
( du 
+ ui = 0 (2)
where the subscript i represents the value at grid point xi. In order to get an expression for (du/dx)i in terms of u at the grid points, we expand ui−1 in a Taylor’s series:

ui−1 = ui − Δx

Rearranging gives )
( du 
+ O(Δx2)
( du 
dx i
= ui − ui−1 + O(Δx) (3)

The error in (du/dx)i due to the neglected terms in the Taylor’s series is called the truncation error. Since the truncation error above is O(Δx), this discrete representation is termed firstorder accurate. 
Using (3) in (2) and excluding higher-order terms in the Taylor’s series, we get the following discrete equation:
ui − ui−1
+ ui = 0 (4) 
Note that we have gone from a differential equation to an algebraic equation! 
This method of deriving the discrete equation using Taylor’s series expansions is called 
the finite-difference method. However, most commercial CFD codes use the finite-volume or 
finite-element methods which are better suited for modeling flow past complex geometries. 
For example, the FLUENT code uses the finite-volume method whereas ANSYS uses the 
finite-element method.  We’ll briefly indicate the philosophy of the finite-volume method 
next but will keep using the finite-difference approach to illustrate the underlying concepts 
which are very similar between the different approaches with the finite-difference method 
being easiest to understand. 

Discretization Using The Finite-Volume Method 
If you look closely at the airfoil grid shown earlier, you’ll see that it consists of quadrilaterals. 
In the finite-volume method, such a quadrilateral is commonly referred to as a "cell" and a 
grid point as a "node". In 2D, one could also have triangular cells. In 3D, cells are usually 
hexahedrals, tetrahedrals, or prisms. In the finite-volume approach, the integral form of the 
conservation equations are applied to the control volume defined by a cell to get the discrete 
equations for the cell. The integral form of the continuity equation for steady, incompressible 
flow is
V ·ndS =0 (5)
The integration is over the surface S of the control volume and n is the outward normal at the surface.  Physically, this equation means that the net volume flow into the control volume is zero. 
Consider the rectangular cell shown below. 

face 4

face 1 
face 3 
1,v1) 3,v3)
Cell center
y face 2

The velocity at face i is taken to be Vi = ui ˆi+vij.  Applying the mass conservation 
equation (5) to the control volume defined by the cell gives 
−u1Δy − v2Δx + u3Δy + v4Δx = 0 
This is the discrete form of the continuity equation for the cell. It is equivalent to summing 
up the net mass flow into the control volume and setting it to zero. So it ensures that the 
net mass flow into the cell is zero i.e. that mass is conserved for the cell. Usually, though not 
always, the values at the cell centers are solved for directly by inverting the discrete system. 
The face values u1, v2, etc. are obtained by suitably interpolating the cell-center values at 
adjacent cells. 
Similarly, one can obtain discrete equations for the conservation of momentum and energy for the cell. One can readily extend these ideas to any general cell shape in 2D or 3D and any conservation equation. Take a few minutes to contrast the discretization in the finite-volume approach to that in the finite-difference method discussed earlier. 
Look back at the airfoil grid. When you are using FLUENT, it’s useful to remind yourself 
that the code is finding a solution such that mass, momentum, energy and other relevant 
quantities are being conserved for each cell.  Also, the code directly solves for values of 
the flow variables at the cell centers; values at other locations are obtained by suitable 
Assembly of Discrete System and Application of Boundary Condi-
Recall that the discrete equation that we obtained using the finite-difference method was
ui − ui−1
Δx Rearranging, we get
+ ui = 0 
−ui−1 + (1 + Δx)ui = 0
Applying this equation to the 1D grid shown earlier at grid points i = 2, 3, 4 gives
−u1 + (1 + Δx) u2 = 0 (i = 2) (6)
−u2 + (1 + Δx) u3 = 0 (i = 3) (7)
−u3 + (1 + Δx) u4 = 0 (i = 4) (8)
The discrete equation cannot be applied at the left boundary (i=1) since ui−1 is not defined
here. Instead, we use the boundary condition to get
u1 = 1 (9)

In a general situation, one would apply the discrete equations to the grid points (or cells 
in the finite-volume method) in the interior of the domain. For grid points (or cells) at or 
near the boundary, one would apply a combination of the discrete equations and boundary 
conditions. In the end, one would obtain a system of simultaneous algebraic equations with 
the number of equations being equal to the number of independent discrete variables. The 
process is essentially the same as for the model equation above with the details being much 
more complex. 
FLUENT, like other commercial CFD codes, offers a variety of boundary condition options such as velocity inlet, pressure inlet, pressure outlet, etc.  It is very important that you specify the proper boundary conditions in order to have a well-defined problem. Also, read through the documentation for a boundary condition option to understand what it does before you use it (it might not be doing what you expect). A single wrong boundary condition can give you a totally wrong result. 

Solution of Discrete System 
The discrete system (10) for our own humble 1D example can be easily inverted to obtain 
the unknowns at the grid points. Solving for u1, u2, u3 and u4 in turn and using Δx = 1/3, 
we get 
u1 = 1 u2 = 3/4 u3 = 9/16 u4 = 27/64
The exact solution for the 1D example is easily calculated to be 
uexact = exp(−x) 
The figure below shows the comparison of the discrete solution obtained on the four-point 
grid with the exact solution. The error is largest at the right boundary where it is equal to 

In a practical CFD application, one would have thousands to millions of unknowns in 
the discrete system and if one uses, say, a Gaussian elimination procedure naively to invert 
the matrix, you’d graduate before the computer finishes the calculation.  So a lot of work 
goes into optimizing the matrix inversion in order to minimize the CPU time and memory 
required.  The matrix to be inverted is sparse i.e. most of the entries in it are zeros since the discrete equation at a grid point or cell will contain only quantities at the neighboring points or cells; verify that this is indeed the case for our matrix system (10). A CFD code would store only the non-zero values to minimize memory usage. It would also generally use an iterative procedure to invert the matrix; the longer one iterates, the closer one gets to the true solution for the matrix inversion. 

Grid Convergence 
While developing the finite-difference approximation for the 1D example, we saw that the truncation error in our discrete system is O(Δx). So one expects that as the number of grid points is increased and Δx is reduced, the error in the numerical solution would decrease and the agreement between the numerical and exact solutions would get better. 
Let’s consider the effect of increasing the number of grid points N on the numerical solution of the 1D problem.  We’ll consider N = 8 and N = 16 in addition to the N = 4 case solved previously.  We repeat the above assembly and solution steps on each of these additional grids.  The resulting discrete system was solved using MATLAB. The following figure compares the results obtained on the three grids with the exact solution. As expected, the numerical error decreases as the number of grid points is increased. 

When the numerical solutions obtained on different grids agree to within a level of tolerance 
specified by the user, they are referred to as "grid converged" solutions.  The concept of 
grid convergence applies to the finite-volume approach also where the numerical solution, if 
correct, becomes independent of the grid as the cell size is reduced.  It is very important 
that you investigate the effect of grid resolution on the solution in every CFD problem you 
solve.  Never trust a CFD solution unless you have convinced yourself that the solution is 
grid converged to an acceptance level of tolerance (which would be problem dependent). 

Dealing with Nonlinearity 
The momentum conservation equation for a fluid is nonlinear due to the convection term 
(V ·∇)V. Phenomena such as turbulence and chemical reaction introduce additional non-
linearities.  The highly nonlinear nature of the governing equations for a fluid makes it 
challenging to obtain accurate numerical solutions for complex flows of practical interest. 
We will demonstrate the effect of nonlinearity by setting m = 2 in our simple 1D exam-ple (1):
dx + u2 = 0;
0≤x≤1;  u(0)=1 
A first-order finite-difference approximation to this equation, analogous to that in (4) for m=1,is
ui − ui−1
+ ui = 0 (11) 
This is a nonlinear algebraic equation with the ui term being the source of the nonlinearity. 
The strategy that is adopted to deal with nonlinearity is to linearize the equations about a guess value of the solution and to iterate until the guess agrees with the solution to a 
specified tolerance level. We’ll illustrate this on the above example. Let ugi  be the guess for ui. Define 
Δui = ui − ugi 
Rearranging and squaring this equation gives 

ui = ugi + 2ugi Δui + (Δui)2 

Assuming that Δui ≪ ugi , we can neglect the Δui term to get
ui ≃ ug

+ 2ugi Δui = ugi  + 2ugi (ui − ugi )

ui ≃ 2ug i ui − ugi 
The finite-difference approximation (11) after linearization becomes 
ui − ui−1 
Δx + 2ugi ui − ugi  = 0 (12)
Since the error due to linearization is O(Δu2), it tends to zero as ug → u. 
In order to calculate the finite-difference approximation (12), we need guess values ug  at the grid points. We start with an initial guess value in the first iteration. For each subsequent iteration, the u value obtained in the previous iteration is used as the guess value. 
Iteration 1:  ug1) = Initial guess 
Iteration 2:  ug2) = u(1) 

Iteration l:   ugl) = u(l−1) 
The superscript indicates the iteration level. We continue the iterations until they converge. 
We’ll defer the discussion on how to evaluate convergence until a little later. 
This is essentially the process used in CFD codes to linearize the nonlinear terms in the conservations equations, with the details varying depending on the code.  The important points to remember are that the linearization is performed about a guess and that it is necessary to iterate through successive approximations until the iterations converge. 

Direct and Iterative Solvers 
We saw that we need to perform iterations to deal with the nonlinear terms in the governing equations. We next discuss another factor that makes it necessary to carry out iterations in practical CFD problems. 

Δx ug4 
In a practical problem, one would usually have thousands to millions of grid points or cells 
so that each dimension of the above matrix would be of the order of a million (with most of 
the elements being zeros). Inverting such a matrix directly would take a prohibitively large 
amount of memory. So instead, the matrix is inverted using an iterative scheme as discussed 
Rearrange the finite-difference approximation (12) at grid point i so that ui is expressed in terms of the values at the neighboring grid points and the guess values: 

ui = ui−1 + Δx ugi 
If a neighboring value at the current iteration level is not available, we use the guess value 
for it. Let’s say that we sweep from right to left on our grid i.e. we update u4, then u3 and 
finally u2 in each iteration. In the mth iteration, ui−)1 is not available while updating ui and
so we use the guess value ugl)i−1
for it instead:
+ Δx ugl)2
Since we are using the guess values at neighboring points, we are effectively obtaining only 
an approximate solution for the matrix inversion in (13) during each iteration but in the 
process have greatly reduced the memory required for the inversion.  This tradeoff is good 
strategy since it doesn’t make sense to expend a great deal of resources to do an exact matrix 
inversion when the matrix elements depend on guess values which are continuously being 
refined. In an act of cleverness, we have combined the iteration to handle nonlinear terms 
with the iteration for matrix inversion into a single iteration process. Most importantly, as 
the iterations converge and ug → u, the approximate solution for the matrix inversion tends 
towards the exact solution for the inversion since the error introduced by using ug  instead 
of u in (14) also tends to zero. 
Thus, iteration serves two purposes: 
1. It allows for efficient matrix inversion with greatly reduced memory requirements. 
2. It is necessary to solve nonlinear equations. 
In steady problems, a common and effective strategy used in CFD codes is to solve the 
unsteady form of the governing equations and "march" the solution in time until the solution 
converges to a steady value. In this case, each time step is effectively an iteration, with the 
the guess value at any time level being given by the solution at the previous time level. 

10 uil) = ugl)

Iterative Convergence 
Recall that as ug  → u, the linearization and matrix inversion errors tends to zero.  So we continue the iteration process until some selected measure of the difference between ug  and u, refered to as the residual, is "small enough". We could, for instance, define the residual R as the RMS value of the difference between u and ug  on the grid: 
R≡ i=1
(ui − ugi )2 
It’s useful to scale this residual with the average value of u in the domain.  An unscaled residual of, say, 0.01 would be relatively small if the average value of u in the domain is 5000 but would be relatively large if the average value is 0.1. Scaling ensures that the residual is a relative rather than an absolute measure.  Scaling the above residual by dividing by the average value of u gives
For the nonlinear 1D example, we’ll take the initial guess at all grid points to be equal 
to the value at the left boundary i.e. ug1) = 1.  In each iteration, we update ug , sweep 
from right to left on the grid updating, in turn, u4, u3  and u2  using (14) and calculate the residual using (15).  We’ll terminate the iterations when the residual falls below 10−9 (which is referred to as the convergence criterion).  Take a few minutes to implement this procedure in MATLAB which will help you gain some familiarity with the mechanics of the implementation.  The variation of the residual with iterations obtained from MATLAB is shown below.  Note that logarithmic scale is used for the ordinate.  The iterative process converges to a level smaller than 10−9 in just 6 iterations. In more complex problems, a lot more iterations would be necessary for achieving convergence. 

1 2 3 4 5 6
Iteration number 

11 N

The solution after 2,4 and 6 iterations and the exact solution are shown below.  It can easily be verified that the exact solution is given by

uexact =
The solutions for iterations 4 and 6 are indistinguishable on the graph.  This is another 
indication that the solution has converged. The converged solution doesn’t agree well with 
the exact solution because we are using a coarse grid for which the truncation error is 
relatively large.  The iterative convergence error, which is of order 10−9, is swamped out 
by the truncation error of order 10−1.  So driving the residual down to 10−9  when the 
truncation error is of order 10−1  is a waste of computing resources. In a good calculation, 
both errors would be of comparable level and less than a tolerance level chosen by the user. 
The agreement between the numerical and exact solutions should get much better on refining 
the grid as was the case for m = 1. 

Iteration 2 
Iteration 4 
0.9 Iteration 6




0 0.2 0.4 0.6 0.8 1

Some points to note: 
1. Different codes use slightly different definitions for the residual. Read the documenta-
tion to understand how the residual is calculated. 
2. In the FLUENT code, residuals are reported for each conservation equation. A discrete 
conservation equation at any cell can be written in the form LHS = 0. For any iteration, 
if one uses the current solution to compute the LHS, it won’t be exactly equal to 
zero, with the deviation from zero being a mesaure of how far one is from achieving 
convergence. So FLUENT calculates the residual as the (scaled) mean of the absolute 
value of the LHS over all cells. 
3. The convergence criterion you choose for each conservation equation is problem- and 
code-dependent. It’s a good idea to start with the default values in the code. One may 
then have to tweak these values. 

Numerical Stability 
In our previous 1D example, the iterations converged very rapidly with the residual falling 
below the convergence criterion of 10−9 in just 6 iterations. In more complex problems, the 
iterations converge more slowly and in some instances, may even diverge.  One would like 
to know a priori the conditions under which a given numerical scheme converges.  This is 
determined by performing a stability analysis of the numerical scheme. A numerical method 
is referred to as being stable when the iterative process converges and as being unstable 
when it diverges. It is not possible to carry out an exact stability analysis for the Euler or 
Navier-Stokes equations. But a stability analysis of simpler, model equations provides useful 
insight and approximate conditions for stability.  A common strategy used in CFD codes 
for steady problems is to solve the unsteady equations and march the solution in time until 
it converges to a steady state.  A stability analysis is usually performed in the context of 
While using time-marching to a steady state, we are only interested in accurately obtaining the asymptotic behavior at large times.  So we would like to take as large a time-step Δt as possible to reach the steady state in the least number of time-steps. There is usually a maximum allowable time-step Δtmax  beyond which the numerical scheme is unstable.  If Δt > Δtmax, the numerical errors will grow exponentially in time causing the solution to diverge from the steady-state result.  The value of Δtmax  depends on the numerical discretization scheme used.  There are two classes of numerical shemes, explicit and implicit, with very different stability characteristics as we’ll briefly discuss next. 

Explicit and Implicit Schemes 
The difference between explicit and implicit schemes can be most easily illustrated by applying them to the wave equation 
∂t + c ∂x = 0 
where c is the wavespeed. One possible way to discretize this equation at grid point i and time-level n is
ui − ui−1
+cui−1 −ui
= O(Δt,Δx) (16) 
The crucial thing to note here is that the spatial derivative is evaluated at the n−1 time-level.
Solving for ui gives
ui = 1−
( cΔt 
)] )
( cΔt
ui−1 + ui 1 (17)
This is an explicit expression i.e. the value of ui at any grid point can be calculated directly 
from this expression without the need for any matrix inversion. The scheme in (16) is known as an explicit scheme.  Since ui  at each grid point can be updated independently, these schemes are easy to implement on the computer.  On the downside, it turns out that this scheme is stable only when 
C≡ cΔt 
Δx ≤ 1 
where C is called the Courant number. This condition is refered to as the Courant-Friedrichs-
Lewy or CFL condition. While a detailed derivation of the CFL condition through stability 
analysis is outside the scope of the current discussion, it can seen that the coefficient of ui−1 
in (17) changes sign depending on whether C > 1 or C < 1 leading to very different behavior in the two cases. The CFL condition places a rather severe limitation on Δtmax. 
In an implicit scheme, the spatial derivative term is evaluated at the n time-level:
ui − ui−1

+cui −ui−1  = O(Δt,Δx)
In this case, we can’t update ui at each grid point independently. We instead need to solve a system of algebraic equations in order to calculate the values at all grid points simultaneously. It can be shown that this scheme is unconditionally stable for the wave equation so that the numerical errors will be damped out irrespective of how large the time-step is. 
The stability limits discussed above apply specifically to the wave equation. In general, 
explicit schemes applied to the Euler or Navier-Stokes equations have the same restriction 
that the Courant number needs to be less than or equal to one.  Implicit schemes are not 
unconditonally stable for the Euler or Navier-Stokes equations since the nonlinearities in 
the governing equations often limit stability.  However, they allow a much larger Courant 
number than explicit schemes. The specific value of the maximum allowable Courant number 
is problem dependent. 
Some points to note: 
1. CFD codes will allow you to set the Courant number (which is also referred to as 
the CFL number) when using time-stepping. Taking larger time-steps leads to faster 
convergence to the steady state, so it is advantageous to set the Courant number as 
large as possible, within the limits of stability, for steady problems. 
2. You may find that a lower Courant number is required during startup when changes 
in the solution are highly nonlinear but it can be increased as the solution progresses. 

Turbulence Modeling 
There are two radically different states of flows that are easily identified and distinguished: laminar flow and turbulent flow. Laminar flows are characterized by smoothly varying velocity fields in space and time in which individual "laminae" (sheets) move past one another without generating cross currents. These flows arise when the fluid viscosity is sufficiently large to damp out any perturbations to the flow that may occur due to boundary imperfections or other irregularities.  These flows occur at low-to-moderate values of the Reynolds number. In contrast, turbulent flows are characterized by large, nearly random fluctuations in velocity and pressure in both space and time. These fluctuations arise from instabilities that grow until nonlinear interactions cause them to break down into finer and finer whirls that eventually are dissipated (into heat) by the action of viscosity. Turbulent flows occur in the opposite limit of high Reynolds numbers. 
A typical time history of the flow variable u at a fixed point in space is shown in Fig. 1(a). 
The dashed line through the curve indicates the "average" velocity. We can define three types 
of averages: 
1. Time average
2. Volume average
3. Ensemble average

No comments:

Post a Comment

leave your opinion