Flavours of Numerical Methods
Topics in CFD and Structural Analyses
Visit the blog page to view topics on mathematics, human resource, textbook solutions, data analytics in stock market trading, online quizzes and much more ...
StepbyStep Guide to CFD Simulation
A typical division of a CFD simulation or as a matter of fact any numerical simulation is in 3 categories namely PreProcessing, Solution and PostProcessing. However, for beginners into this field, it does not convey much insight. The complete simulation process is explained here in a stepbystep sequence to make it easy to understand.
 Step01: Problem Definition  Before you jump on the simulation methods, define the problem statement and expected output, accuracy and turnaround time. This requires asking questions like
 What would be size of the computational domain (the volume in 3D or area in 2D)? Do I need to model solid or only fluid?
 Where the fluid is entering into the chosen computational domain and how? Is it fully developed or uniform? What would be the turbulence level?
 What is the information available at inlet: pressure, velocity, temperature, turbulence level?
 Where can I put the location where fluid exits from computational domain? What are the known parameter there: pressure, velocity...? Will it be prone to reverse flow? Can I extend the location without affecting the result?
 What are the conditions at the walls: known temperature or heat flux? Do I need to consider the surface roughness?
 Do I need to define the full size or use a reduced scale?
 Step02: Geometry Creation for Flow Domain  Once you are clear with expected goal and outcome from your simulation, prepare the computational domain in a CAD program such as FreeCAD, Solidworks, CATIA or inside the preprocessor directly such as GMSH, Salome, OpenFOAM ... The CAD data stored in format of the program used to create it is called Native format such *.FCStd for FreeCAD. There are neutral file formats which are used to exchange data such as IGES (surface data), STEP (solid data), STL (triangulated surface data) and Parasolid to name few of them. Additional information can be found here for considerations required to create the computational domain.
 Step03: Geometry Simplification for Simulation  This step is not requires if geometry has been created specially for CFD simulation. This step is required when you use CAD data created by design engineers which may contain additional information not needed for numerical simulations. These may be small fillets at corners, bolts and screws, gaps for shrinkfit or CAD modeling error... The purpose it to get a watertight geometry where all the faces are connected with single edges. This step is called geometry simplification and defeaturing (removing features like chamfers, fillets, taper, bolt holes...). Additional information can be found here for defeaturing methods and recommendations.
 Step04: Mesh Generation :  this is a process of dividing the overall volume into smaller pieces or smaller volumes. A given computational domain may have infinite degreeoffreedom (DOF) or locations where solution need to be found out. However, numerical method need to reduce the DOF to 'some' finite number of DOF and hence this process is also called FiniteElement Method or FiniteVolume Method. Additional information can be found here for mesh generation methods and recommendations.
 Step05: Material Property and Boundary Condition Setting  most of the numerical methods are known as boundary value problem. That is they need to satisfy some known conditions or constraints on a predefined locations or boundaries. For fluids, walls, inlet and outlet are some of the basic boundary types. Refer to this page for type of boundaries and recommendations.
 Step06: Solver Setting  The setup is now ready for solution. All the numerical simulations gets reduced to the format [A]{x} = {b} where [A] is large matrix and {x} is the unknowns we are looking to estimate. In this step, the method to model turbulence, wall function, discretization scheme for differential operators such as du/dy and d^{2}u/dy^{2} is specified. Refer to this page for methods and recommendations. A detailed description of turbulence modeling can be found here.
 Step07: Initialization, solution monitoring and troubleshooting  this step the user needs to provide a good guess of pressure, velocity and turbulence then start the simulation run. Then the user need to monitor the progress while program solves equation [A]{x} = {b} iteratively. The program may not result in a solution (known as divergence of solution process) or may lead to a solution but the residuals may not reach to low level recommended by the developers of the program. There might be sharp fluctuations in the calculated values and one will need to adjust the solver parameters to eliminate abnormality. This may be achieved by changing relaxation factors, discretization scheme to first order or even refining and improving the mesh quality.
 Step08: Postprocessing  Once solution is achieved (known as convergence of iterative process), the results need to be presented in visual form using graphs, contour plots and tabulated data. This step is called postprocessing. However, this step is not limited to present the result but share the insight gained through simulation process in the form of design recommendations, justification for accuracy of results and effect of assumptions made during the start of the simulation activity. Refer to this page for methods and recommendations. More information about validation and checking errors associated with simulation can be found here.
All these steps have been explained through an example simulation for "Flow Over a BackStep" and can be found here! If you have access to ANSYS FLUENT, a stepbystep guide with GUI screenshots can be found in this PDF file. Similarly, if you have experience in using any commercial tool say CFX and want to try opensource product such as OpenFOAM, a comparison of methods is summarized in this page.
In addition to the steps, as a thermal and fluid flow analyst, you may need to make many other considerations as summarized below.
 Boundary layer height: How do I know well in advance what should be the first layer height to resolve boundary layer properly? While no unique method exist, most of the flow conditions can be approximated to either of the 3 conditions  flow inside a duct, flow over a bluff body and flow over a flat plate. This page provides an initial estimate of first layer height for required Y^{+} requirements.
 Rotating parts: Which methods are available to model rotating components and turbomachines? A short summary on simulation methods of rotating walls and turbomachines has been summarized here.
 Steady or transient? Steady state refers to the condition when field variables (velocity, pressure, temperature ...) are not function of time. Sometimes simulation does not converge when a steady state solution approach is used and one of the reason can also be inherent transient nature of the flow phenomena.
 Porous zones: The simulation parts may involve components with not only very thin passage but also large number of them. Examples include heat exchangers like automotive radiators, HVAC condensers, honeycomb geometry in catalytic converters, perforated grilles... The need to be modeled in a special way known as porous domain (3D) and porous jump (2D).
 More than one phases? Sometimes mixing of two fluids need to be modeled. These phenomena are called multiphase flows such as flows in emulsification, mixing tanks, continuously stirred tank reactors (CSTR), filling of a tank, sloshing in a automotive fuel tank... A basic description of multiphase flows and methods to simulate such phenomena in OpenFOAM is described here.
 Buoyancy and Natural Convection: This is one area which had wide applications in cooling of electronic components and HVAC. However, such simulations are prone to divergence and special considerations are required. One of them is to chose which assumption is best applicable to model the buoyancy effect.
Convergence and relaxation factors
This is one of the most talked about (and sometimes dreaded) term used in numerical simulation field. The term 'convergence' refers to the process where the iterative method moves towards "expected final result". It is best explained using a problem from algebra where we need find out value of x in equation y = 5x^{2} + 25 when y = 1000. We all know a direct method:
x = [(y  25)/5]^{0.5}.
Let's try an iterative method using algorithm x_{i+1} = x_{i} × (1  Δy_{i}). Here, i refers to the iteration number and Δy refers to deviation of the calculated value of y from desired value of y = 1000, ResNorm1 refers to the absolute value where the residual % is normalized (divided) by the residual % obtained in 0^{th} iteration based on guess value, ResNorm3 refers to the absolute value when the residual % is normalized (divided) by the residual % obtained in 0^{th} ~ 2^{nd} iterations.
Iteration, i  x_{i}  y_{i}  Δy_{i}  ResNorm1  ResNorm3 
Initial Guess  10.00  525.0  +475.0  +47.5%     
1  14.75  1112.8  112.8  11.3%  0.2375   
2  13.09  881.2  +118.8  +11.9%  0.2501   
3  14.64  1096.7  96.70  9.67%  0.2036  0.6032 
4  13.22  899.45  +100.6  +10.1%  0.2117  0.6272 
5  14.55  1084.1  84.14  8.41%  0.1771  0.5248 
6  13.33  913.40  +86.60  +8.66%  0.1823  0.5401 
7  14.48  1073.9  73.93  7.39%  0.1556  0.4611 
8  13.41  924.57  +75.43  +7.54%  0.1588  0.4705 
9  14.42  1065.4  65.40  6.54%  0.1377  0.4079 
10  13.48  933.77  +66.23  +6.62%  0.1394  0.4131 
11  14.37  1058.1  58.14  5.81%  0.1224  0.3626 
12  13.54  941.50  +58.50  +5.85%  0.1232  0.3649 
13  14.33  1051.9  51.86  5.19%  0.1092  0.3235 
14  13.59  948.11  +51.89  +5.19%  0.1092  0.3236 
15  14.29  1046.4  46.39  4.64%  0.0977  0.2894 
16  13.63  953.82  +46.18  +4.62%  0.0972  0.2881 
17  14.26  1041.6  41.59  4.16%  0.0876  0.2594 
18  13.67  958.79  +41.21  +4.12%  0.0868  0.2570 
19  14.23  1037.3  37.34  3.73%  0.0786  0.2329 
20  13.70  963.15  +36.85  +3.68%  0.0776  0.2298 
21  14.20  1033.6  33.56  3.36%  0.0707  0.2094 
22  13.73  967.00  +33.00  +3.30%  0.0695  0.2059 
23  14.18  1030.2  30.20  3.02%  0.0636  0.1884 
24  13.75  970.40  +29.60  +2.96%  0.0623  0.1846 
25  14.16  1027.2  27.19  2.72%  0.0572  0.1696 
Look how the calculated value of x still makes y off the final value by 2.72% even after 25 iterations. What are the ways to speedup this process  that is how can we reach to the correct value with lesser number of iterations? One option is to change the way the value of x_{i+1} is calculated from x_{i}. This method is known as use of "Relaxation Factors" in numerical simulation technology. Let's try the approach: 0.8 in the following equation is the relaxation factor (also called underrelaxation factor as the value is < 1.0).
x_{i+1} = [x_{i} × (1  Δy_{i})] * 0.8 + x_{i} * 0.2
Iteration, i  x_{i}  y_{i}  Δy_{i}  ResNorm1  ResNorm3 
Initial Guess  10.00  525.0  475.0  47.5%     
1  14.75  1112.8  112.8  11.3%  0.2375   
2  13.42  925.30  +74.70  +7.50%  0.1572   
3  14.22  1036.1  36.11  3.61%  0.0760  0.2480 
4  13.81  978.54  +21.50  +2.10%  0.0452  0.1474 
5  14.05  1011.6  11.56  1.16%  0.0243  0.0794 
6  13.92  993.39  +6.61  +0.66%  0.0139  0.0454 
7  13.99  1003.7  3.66  0.37%  0.0077  0.0251 
8  13.95  997.94  +2.06  +0.21%  0.0043  0.0141 
9  13.97  1001.1  1.15  0.11%  0.0024  0.0079 
10  13.96  999.35  +0.65  +0.06%  0.0014  0.0044 
11  13.97  1000.4  0.36  0.04%  0.0008  0.0025 
12  13.96  999.80  +0.20  +0.02%  0.0004  0.0014 
13  13.97  1000.1  0.11  0.01%  0.0002  0.0008 
14  13.96  999.94  +0.06  +0.01%  0.0001  0.0004 
15  13.96  1000.0  0.04  +0.00%  0.0001  0.0002 
16  13.96  999.98  +0.02  +0.00%  0.0000  0.0001 
17  13.96  1000.0  0.01  +0.00%  0.0000  0.0001 
This is an excellent example of relaxation factor improving the speed of convergence. Readers are advised to try out the calculation with different values of initial guess (say 12, 15 instead of 10) and relaxation factor (say 0.2, 0.5 instead of 0.8) to see the effect on rate of convergence (number of iterations required to find result x = 13.96).
 Absolute Residuals: The values reported in column4 of the tables above are absolute residuals, "in this example" estimated as difference of exact solution and value calculated during iterations. In most of the simulations, the result (or target value) is not known in advance. Hence, for a system of linear equation [A]{x} ={b} the absolute residual is sum of residuals at each node. In 2D case:
A*x_{P} + B*x_{N} + C*x_{S} + D*x_{E} + E*x_{W} = F_{P} where P denotes the calculation node and N, S, E, W denotes the north, south, east and west nodes respectively. Here, absolute residual is
ε^{n} = [A*x_{P}^{n} + B*x_{N}^{n} + C*x_{S}^{n} + D*x_{E}^{n} + E*x_{W}^{n}  F_{P}^{n}] where n is the iteration number.
 Scaled Residuals: The values reported in column5 of the tables above are scaled residuals. They have been divided by the calculated value of dependent variable and hence indicate relative change with respect to expected result. In comparison to "absolute residuals" which has no bounds and bounds themselves change from problem to problem, scaled residual has welldefined bounds on a scale of 0~1 or 0~100 and easy to interpret for all types of iterations. ANSYS FLUENT reports this residual by default. It also has option to report normalized residual as explained below.
 Normalized Residuals: The values reported in column6 and colum7 of the tables above are normalized residuals. Note that the residuals are further divided by another residual hence there cannot be an upper bound defined. What if the residual is very low in first iteration (occurs when the initial guess is very close to exact result)? The subsequent normalized residuals will be high even if the solution approaches (converges) to a correct and accurate result. In order to avoid this ambiguity, one of the following two approaches can be adopted.
 Normalize by average of first 'few' iterations
 Normalize by maximum residual obtained in first few iterations. ANSYS FLUENT uses this method with default value of first '5' iterations and can be set by user.
 Convergence of Residuals:Typically the recommended practice it to have residuals (scales or normalized) of the order of 10^{4} or lower (this threshold varies across the industries and desired accuracy  for example an accuracy in lift coefficient or drag coefficient prediction of 0.001 for airfoils). However, even in cases where such low residuals are not achieved, it does not necessarily mean that the simulation results are poor. In the absence of good residual convergence, more emphasis should be placed on the convergence of integral values (known as monitor points) such as lift or force coefficients, areaaveraged wall shear stress ... Yet what constitutes good convergence of a such integral value needs to be evolved for each field of applications (as it may not always be asymptotic especialy in cases of large region of separation and flow instabilities in wake regions).
 Oscillating Residuals: In many applications such as turbomachines, flow around ships and bluffbodies, there is always some amount of separation and unsteadiness in wake regions. In such cases, residuals show a "steady oscillation" pattern where such oscillations tend to be periodic and reaches a constant amplitude about a constant mean value. The outcomes of such simulation can be considered fairly accurate and meaningful with average value used in final calculations of field variables.
Allied Methods
 DEM: Discrete Element Method: The wellknown FEM and FVM method is a continuumbased calculation techniques where the media is assumed to conform to continuum hypothesis (where the properties of the matter are considered as continuous functions of space variables x, y an z. Though all matters are composed of atoms and molecules, the concept of continuum assumes a continuous distribution of mass within the matter with no void or empty space).
 DEM is application similar to CFD which applies to flow (rather movement) of granular media such as grains, stone chips, sand.
 In DEM, the masses are numerically represented as a group of points which is subsequently is associated with a physical particle in the mass.

Further, interaction among the particles (friction, collision, deformation) are modelled by springs and dashpots that are applied in normal and shear directions (at the point of contact or collision) respectively.
 In civil engineering, this method is being applied for failure of slopes, calculation of angle of repose.
 Some more applications include design of deflectors and particle chains used to slow down moving particles, design of transfer chutes, spillage from conveyors, placement of rockboxes...
 Ease of applicability to moving boundary problems is one of the key characteristics of this method. Recently, it is also being used to study the fracture of brittle materials.
 One such program is Newton DEM software from Advanced Conveyor Technologies Inc widely used in mining applications.
 Another program dealing with DEM is Rocky DEM  example of particle separation using vibrating screen in Rocky DEM
 FSI: FluidStructure Interaction  this is a multiphysics application of CFD methods where the effect of fluid forces cause deflection of the structure large enough to affect the flow field. Thus, it is a twoway coupling where CFD and FEM method need to share the information to define timedependent displacement of walls and fluid velocity/pressure fields. As evident, this is necessarily an unsteady (transient) phenomena and selection of time step is crucial and most of the times not known apriori.
 CAA: Computational AeroAcoustics  this is another multiphysics application of CFD methods to predict the noise sources and solve noise propagation using other solvers such as ACTRAN. In CFD simulation environment, acoustics is just a postprocessing activity where the fluctuating componenet of pressure is converted into noise using densitybased analogies like (Lighthill’s analogy, Powell’s analogy, the Ffowcs Williams–Hawkings [FWH] analogy and Curle’s analogy) or few other analogies like Phillips’ analogy and Lilley’s analogy, Howe’s analogy and Doak’s analogy. The noise source can be calculated both for steady and transient turbulent flows, however the option avaialable will vary for steady state and transient cases.
 BEM: Boundary Element Method  also known as the "boundary integral equation method" or "boundary integral method": As the name suggest, this method discretizes only the boundaries of the computational domain. Widely used in field problems with unbounded exterior domain but bounded boundaries such as magnetic field calculation and acoustics  this method does not need to model the entire domain required in FEM/CFD.
 Mathematical definition: "boundary element method" (BEM) refers to any method for the approximate numerical solution of a "boundary value problem" (BVP) of "partial differential equations" (PDE).
 This method is limited to only for linear partial differential equations with constant coefficients whereas problems with nonlinear differential equations require coupling between FEM and BEM.
 LBM: Lattice Boltzmann Method  a discrete computational method based on Boltzmann's equation. This is an alternate method to make fluid dynamics simulations based on lattice gas methods. The method allows particles to move on a discrete lattice (analogous to mesh) and collide such that local collisions conserve mass (continuity) and momentum (NavierStokes equation).
 FNM: Fluid Network Method  a flow and pressure drop calculation approach based on analogy with electric resistances and circuits. 3D flow paths are condensed into equivalent flow resistance 'R' defined as Δp = ρ × R × Q^{2} where Q is flow rate. A network is created with nodes and branches where nodes represents the connection of two or more branches (the flow resistances). Similar to Kirchoff's law for electrical networks, the flow and pressure drop across each hydraulic resistance can be calculated using mass balance at each node and zero pressure drops across all the flow loops. HardyCross method is one of the widely used method to solve flow network models. 1D model for a single cylinder engine in GTPower described below is an example of FNM.
Water Hammer: This is a special type of transient phenomena which calculates pressure generated due to sudden change in flow velocity and subsequent calculation of pressure wave propagation. Hence, computations related to water hammer involves both CFD and wave speed calculation techniques. Suppose water is flowing into a straight long pipe at constant flow rate and suddenly a control valve is use to reduce the flow rate by 50% or even stopped the flow completely, a pressure wave will get generated which will travel upstream at a speed closer to speed to sound in water. This phenomena is governed by Joukowsky equation: Δp = [± ρ c Δv] where ρ is the density of fluid, c is (waterhammer) wave speed and δv is change in velocity.
This is known as fundamental equation of water hammer. Thus, head developed due to water hammer, Δh = [c Δv / g] where g is acceleration due to gravity. This is inherently a transient phenomena and CFD methods require simulation for steady state condition and then implementation of a transient phenomena over a converged steasy state flow field. In industry, there are many 1D tools which are based on Methods of Charactaristics (MOC) to solve water hammer phenomena. Some of the tools are WANDA, HAMMER, AFT Impulse ...