Computational Fluid Dynamics

Computational fluid dynamics or CFD is a branch of fluid mechanics that uses numerical analysis and data structures to analyze and solve problems that involve fluid flows. Computers are used to perform the calculations required to simulate the free-stream flow of the fluid, and the interaction of the fluid (liquids and gases) with surfaces defined by boundary conditions. With high-speed supercomputers, better solutions can be achieved, and are often required to solve the largest and most complex problems. Ongoing research yields software that improves the accuracy and speed of complex simulation scenarios such as transonic or turbulent flows. Initial validation of such software is typically performed using experimental apparatus such as wind tunnels. In addition, previously performed analytical or empirical analysis of a particular problem can be used for comparison. A final validation is often performed using full-scale testing, such as flight tests.

CFD is applied to a wide range of research and engineering problems in many fields of study and industries, including aerodynamics and aerospace analysis, weather simulation, natural science and environmental engineering, industrial system design and analysis, biological engineering, fluid flows and heat transfer, and engine and combustion analysis.

What is CFD | Computational Fluid Dynamics?

When an engineer is tasked with designing a new product, e.g. a winning race car for the next season, aerodynamics plays an important role in the engineering process. However, aerodynamic processes are not easily quantifiable during the concept phase. Usually the only way for the engineer to optimize his designs is to conduct physical tests on product prototypes. With the rise of computers and ever-growing computational power (thanks to Moore’s law!), the field of Computational Fluid Dynamics became a commonly applied tool for generating solutions for fluid flows with or without solid interaction. In a CFD analysis, the examination of fluid flow in accordance with its physical properties such as velocity, pressure, temperature, density and viscosity is conducted. To virtually generate a solution for a physical phenomenon associated with fluid flow, without compromise on accuracy, those properties have to be considered simultaneously.

A mathematical model of the physical case and a numerical method are used in a software tool to analyze the fluid flow. For instance, the Navier-Stokes equations are specified as the mathematical model of the physical case. This describes changes on all those physical properties for both fluid flow and heat transfer. The mathematical model varies in accordance with the content of the problem such as heat transfer, mass transfer, phase change, chemical reaction, etc. Moreover, the reliability of a CFD analysis highly depends on the whole structure of the process. The verification of the mathematical model is extremely important to create an accurate case for solving the problem. Besides, the determination of proper numerical methods to generate a path through the solution is as important as a mathematical model. The software, which the analysis is conducted with is one of the key elements in generating a sustainable product development process, as the number of physical prototypes can be reduced drastically.

A computer simulation of high velocity air flow around the Space Shuttle during re-entry.
A computer simulation of high velocity air flow around the Space Shuttle during re-entry.

Background and History

The fundamental basis of almost all CFD problems is the Navier–Stokes equations, which define many single-phase (gas or liquid, but not both) fluid flows. These equations can be simplified by removing terms describing viscous actions to yield the Euler equations. Further simplification, by removing terms describing vorticity yields the full potential equations. Finally, for small perturbations in subsonic and supersonic flows (not transonic or hypersonic) these equations can be linearized to yield the linearized potential equations.

Historically, methods were first developed to solve the linearized potential equations. Two-dimensional (2D) methods, using conformal transformations of the flow about a cylinder to the flow about an airfoil were developed in the 1930s.

Claude-Louis Navier (left) and Sir George Gabriel Stokes
Claude-Louis Navier (left) and Sir George Gabriel Stokes

One of the earliest type of calculations resembling modern CFD are those by Lewis Fry Richardson, in the sense that these calculations used finite differences and divided the physical space in cells. Although they failed dramatically, these calculations, together with Richardson’s book “Weather prediction by numerical process”, set the basis for modern CFD and numerical meteorology. In fact, early CFD calculations during the 1940s using ENIAC used methods close to those in Richardson’s 1922 book.

The computer power available paced development of three-dimensional methods. Probably the first work using computers to model fluid flow, as governed by the Navier-Stokes equations, was performed at Los Alamos National Lab, in the T3 group. This group was led by Francis H. Harlow, who is widely considered as one of the pioneers of CFD. From 1957 to late 1960s, this group developed a variety of numerical methods to simulate transient two-dimensional fluid flows, such as Particle-in-cell method (Harlow, 1957), Fluid-in-cell method (Gentry, Martin and Daly, 1966), Vorticity stream function method (Jake Fromm, 1963), and Marker-and-cell method (Harlow and Welch, 1965). Fromm’s vorticity-stream-function method for 2D, transient, incompressible flow was the first treatment of strongly contorting incompressible flows in the world.

A simulation of the Hyper-X scramjet vehicle in operation at Mach-7
A simulation of the Hyper-X scramjet vehicle in operation at Mach-7

The first paper with three-dimensional model was published by John Hess and A.M.O. Smith of Douglas Aircraft in 1967. This method discretized the surface of the geometry with panels, giving rise to this class of programs being called Panel Methods. An intermediate step between Panel Codes and Full Potential codes were codes that used the Transonic Small Disturbance equations. In particular, the three-dimensional WIBCO code, developed by Charlie Boppe of Grumman Aircraft in the early 1980s has seen heavy use.

The next step was the Euler equations, which promised to provide more accurate solutions of transonic flows. The methodology used by Jameson in his three-dimensional FLO57 code (1981) was used by others to produce such programs as Lockheed’s TEAM program and IAI/Analytical Methods’ MGAERO program. MGAERO is unique in being a structured cartesian mesh code, while most other such codes use structured body-fitted grids (with the exception of NASA’s highly successful CART3D code, Lockheed’s SPLITFLOW code and Georgia Tech’s NASCART-GT). Antony Jameson also developed the three-dimensional AIRPLANE code which made use of unstructured tetrahedral grids. In the two-dimensional realm, Mark Drela and Michael Giles, then graduate students at MIT, developed the ISES Euler program (actually a suite of programs) for airfoil design and analysis. This code first became available in 1986 and has been further developed to design, analyze and optimize single or multi-element airfoils, as the MSES program. MSES sees wide use throughout the world. A derivative of MSES, for the design and analysis of airfoils in a cascade, is MISES, developed by Harold “Guppy” Youngren while he was a graduate student at MIT.

Steamlines around an F1 car
Steamlines around an F1 car

The Navier–Stokes equations were the ultimate target of development. Two-dimensional codes, such as NASA Ames’ ARC2D code first emerged. A number of three-dimensional codes were developed (ARC3D, OVERFLOW, CFL3D are three successful NASA contributions), leading to numerous commercial packages.

Governing Equations

The main structure of thermo-fluids examinations is directed by governing equations that are based on the conservation law of fluid’s physical properties. The basic equations are the three physics laws of conservation:

  1. Conservation of Mass: Continuity Equation
  2. Conservation of Momentum: Momentum Equation of Newton’s Second Law
  3. Conservation of Energy: First Law of Thermodynamics or Energy Equation

These principles state that mass, momentum and energy are stable constants within a closed system.

Basically: What comes in, must also go out somewhere else.

The investigation of fluid flow with thermal changes relies on certain physical properties. The three unknowns which must be obtained simultaneously from these three basic conservation equations are the velocity v, pressure p and the absolute temperature T. Yet p and T are considered the two required independent thermodynamics variables. The final form of the conservation equations also contains four other thermodynamics variables; density ρ, the enthalpy h as well as viscosity μ and thermal conductivity k; the last two are also transport properties. Since p and T are considered two required independent thermodynamics variables, these four properties are uniquely determined by the value of p and T. Fluid flow should be analyzed to know vecv, p and T throughout every point of the flow regime. This is most important before designing any product which involves fluid flow.

Observation of fluid motion with the methods of Lagrange and Euler
Observation of fluid motion with the methods of Lagrange and Euler

Furthermore, the method of fluid flow observation based on kinematic properties is a fundamental issue. Movement of fluid can be investigated with either Lagrangian or Eulerian methods. Lagrangian description of fluid motion is based on the theory to follow a fluid particle which is large enough to detect properties. Initial coordinates at time t0 and coordinates of the same particle at time t1 have to be examined. To follow millions of separate particles through the path is almost impossible. In the Eulerian method, any specific particle across the path is not followed, instead, the velocity field as a function of time and position is examined. This missile example precisely fits to emphasize the methods.

Langarian: We take up every point at the beginning of the domain and trace its path till it reaches the end.

Eulerian: We consider a window (Control Volume) within the fluid and analyse the particle flow within this Volume.

Lagrangian formulation of motion is always time-dependent. As a, b and c are the initial coordinates of a particle; x, y, and z are coordinates of the same particle at time t. Description of motion for Lagrangian:

$ \displaystyle x=x(a,b,c,t) ; y=y(a,b,c,t) ; z=z(a,b,c,t)$

In the Eulerian method, u, v and w are the components of velocity at the point (x,y,z) at the time t. Thus, u, v and w are the unknowns which are functions of the independent variables x,y,z and t. Description of motion for Eulerian for any particular value of t:

$ \displaystyle u=u(x,y,z,t) ; v=v(x,y,z,t) ; w=w(x,y,z,t)$

Conservation of Mass is specified as equation below:

$ \displaystyle \frac{{D\rho }}{{Dt}}+\rho (\nabla \cdot \vec{v})=0$

where ρ is the density, v the velocity and ∇ the gradient operator.

$ \displaystyle \vec{\nabla }=\vec{i}\frac{\partial }{{\partial x}}+\vec{j}\frac{\partial }{{\partial y}}+\vec{k}\frac{\partial }{{\partial z}}$

If the density is constant, the flow is assumed to be incompressible and then continuity reduces it to:

$ \displaystyle \frac{{D\rho }}{{Dt}}=0\to \nabla \cdot \vec{v}=\frac{{\partial u}}{{\partial x}}+\frac{{\partial v}}{{\partial y}}+\frac{{\partial w}}{{\partial z}}=0$

Conservation of Momentum which can be specified as Navier-Stokes Equation:

$ \displaystyle \overset{I}{\mathop{{\frac{\partial }{{\partial t}}(\rho \vec{v})}}}\,+\overset{{II}}{\mathop{{\nabla \cdot (\rho \vec{v}\vec{v})}}}\,=\overset{{III}}{\mathop{{-\nabla p}}}\,+\overset{{IV}}{\mathop{{\nabla \cdot \left( {\overline{{\bar{\tau }}}} \right)}}}\,+\overset{V}{\mathop{{\rho \vec{g}}}}\,$

where static pressure p, viscous stress tensor τ̅ and gravitational force ρg.

  • I: Local change with time
  • II: Momentum convection
  • III: Surface force
  • IV: Diffusion term
  • V: Mass force

Viscous stress tensor τ̅ can be specified as below in accordance with Stoke’s Hypothesis:

$ \displaystyle {{\tau }_{{ij}}}=\mu \frac{{\partial {{v}_{i}}}}{{\partial {{x}_{j}}}}+\frac{{\partial {{v}_{j}}}}{{\partial {{x}_{i}}}}-\frac{2}{3}(\nabla \cdot \vec{v}){{\delta }_{{ij}}}$

If the fluid is assumed to be of constant density (cannot be compressed), the equations are greatly simplified that viscosity coefficient μ is assumed constant. Therefore, many terms vanished through equation results in a much simpler Navier-Stokes Equation:

$ \displaystyle \rho \frac{{D\vec{v}}}{{Dt}}=-\nabla p+\mu {{\nabla }^{2}}\vec{v}+\rho \vec{g}$

Conservation of Energy is the first law of thermodynamics which states that the sum of the work and heat added to the system will result in the increase of the energy in the system:

$ \displaystyle \rho \frac{{D\vec{v}}}{{Dt}}=-\nabla p+\mu {{\nabla }^{2}}\vec{v}+\rho \vec{g}$

where dQ is the heat added to the system, dW is the work done on the system and dEt is the increment in the total energy of the system. One of the common types of an energy equation is:

$ \displaystyle \rho \left[ {\overset{I}{\mathop{{\frac{{\partial h}}{{\partial t}}}}}\,+\overset{{II}}{\mathop{{\nabla \cdot (h\vec{v})}}}\,} \right]=\overset{{III}}{\mathop{{-\frac{{\partial p}}{{\partial t}}}}}\,+\overset{{IV}}{\mathop{{\nabla \cdot (k\nabla T)}}}\,+\overset{V}{\mathop{\phi }}\,$

I: Local change with time, II: Convective term, III: Pressure work, IV: Heat flux, V: Source term

– Partial Differential Equations (PDEs)

The Mathematical model merely gives us interrelation between the transport parameters which are involved in the whole process, either directly or indirectly. Even though every single term in those equations has relative effect on the physical phenomenon, changes in parameters should be considered simultaneously through the numerical solution which comprises differential equations, vector and tensor notations. A PDE comprises more than one variable and their derivation which is specified with “∂” instead of “d”. If the derivation of the equation is conducted with “d”, these equations are called as Ordinary Differential Equations (ODE) that contains a single variable and its derivation. The PDEs are implicated to transform the differential operator (∂) into an algebraic operator in order to get a solution. Heat transfer, fluid dynamics, acoustic, electronics and quantum mechanics are the fields that PDEs are highly used to generate solutions.

Example of ODE:

$ \displaystyle \frac{{{{d}^{2}}x}}{{d{{t}^{2}}}}=x\to x(t)$ ; Where T is the single variable

Example of PDE:

$ \displaystyle \frac{{\partial f}}{{\partial x}}+\frac{{\partial f}}{{\partial y}}=5\to f(x,y)$ ; Where both x and y are the variables

What is the significance of PDEs to seeking a solution on governing equations? To answer this question, we initially examine the basic structure of some PDEs as to create connotation. For instance:

$ \displaystyle \frac{{{{\partial }^{2}}f}}{{\partial {{x}^{2}}}}+\frac{{{{\partial }^{2}}f}}{{\partial {{y}^{2}}}}=0\to f(x,y)\to \text{Laplace Equation}$

Comparison between the last equation and equation of Mass Conservation specifies the Laplace part of the continuity equation. What is the next step? What does this Laplace analogy mean? To start solving these enormous equations, the next step comes through discretization to ignite the numerical solution process. The numerical solution is a discretization-based method used in order to obtain approximate solutions to complex problems which cannot be solved with analytic methods because of complexity and ambiguities. As seen in the previous figure, solution processes without discretization merely give you an analytic solution which is exact but simple. Moreover, the accuracy of the numerical solution highly depends on the quality of the discretization. Broadly used discretization methods might be specified such as finite difference, finite volume, finite element, spectral (element) methods and boundary element.

The panoramic structure of a CFD project and its stages
The panoramic structure of a CFD project and its stages


In all of these approaches the same basic procedure is followed.

  • During pre-processing:
    • The geometry and physical bounds of the problem can be defined using computer aided design (CAD). From there, data can be suitably processed (cleaned-up) and the fluid volume (or fluid domain) is extracted.
    • The volume occupied by the fluid is divided into discrete cells (the mesh). The mesh may be uniform or non-uniform, structured or unstructured, consisting of a combination of hexahedral, tetrahedral, prismatic, pyramidal or polyhedral elements.
    • The physical modeling is defined – for example, the equations of fluid motion + enthalpy + radiation + species conservation.
    • Boundary conditions are defined. This involves specifying the fluid behavior and properties at all bounding surfaces of the fluid domain. For transient problems, the initial conditions are also defined.

  • The simulation is started and the equations are solved iteratively as a steady-state or transient.
  • Finally, a postprocessor is used for the analysis and visualization of the resulting solution.

– Discretization methods

The stability of the selected discretization is generally established numerically rather than analytically as with simple linear problems. Special care must also be taken to ensure that the discretization handles discontinuous solutions gracefully. The Euler equations and Navier–Stokes equations both admit shocks, and contact surfaces.

Some of the discretization methods being used are:

* Finite volume method

The finite volume method (FVM) is a common approach used in CFD codes, as it has an advantage in memory usage and solution speed, especially for large problems, high Reynolds number turbulent flows, and source term dominated flows (like combustion). In the finite volume method, the governing partial differential equations (typically the Navier-Stokes equations, the mass and energy conservation equations, and the turbulence equations) are recast in a conservative form, and then solved over discrete control volumes. This discretization guarantees the conservation of fluxes through a particular control volume. The finite volume equation yields governing equations in the form,

$ \displaystyle \frac{\partial }{{\partial t}}\iiint{Q}dV+\iint{F}d\mathbf{A}=0,$

where Q is the vector of conserved variables, F is the vector of fluxes, V is the volume of the control volume element, and A is the surface area of the control volume element.

* Finite element method

The finite element method (FEM) is used in structural analysis of solids, but is also applicable to fluids. However, the FEM formulation requires special care to ensure a conservative solution. The FEM formulation has been adapted for use with fluid dynamics governing equations. Although FEM must be carefully formulated to be conservative, it is much more stable than the finite volume approach. However, FEM can require more memory and has slower solution times than the FVM.

In this method, a weighted residual equation is formed:

$ \displaystyle {{R}_{i}}=\iiint{{{{W}_{i}}}}Qd{{V}^{e}}$

Where Ri is the equation residual at an element vertex i, Q is the conservation equation expressed on an element basis, Wi is the weight factor, and Ve is the volume of the element.

* Finite difference method

The finite difference method (FDM) has historical importance and is simple to program. It is currently only used in few specialized codes, which handle complex geometry with high accuracy and efficiency by using embedded boundaries or overlapping grids (with the solution interpolated across each grid).

$ \displaystyle \frac{{\partial Q}}{{\partial t}}+\frac{{\partial F}}{{\partial x}}+\frac{{\partial G}}{{\partial y}}+\frac{{\partial H}}{{\partial z}}=0$

where Q is the vector of conserved variables, and F, G, and H are the fluxes in the x, y, and z directions respectively.

* Spectral element method

Spectral element method is a finite element type method. It requires the mathematical problem (the partial differential equation) to be cast in a weak formulation. This is typically done by multiplying the differential equation by an arbitrary test function and integrating over the whole domain. Purely mathematically, the test functions are completely arbitrary – they belong to an infinite-dimensional function space. Clearly an infinite-dimensional function space cannot be represented on a discrete spectral element mesh; this is where the spectral element discretization begins. The most crucial thing is the choice of interpolating and testing functions. In a standard, low order FEM in 2D, for quadrilateral elements the most typical choice is the bilinear test or interpolating function of the form v(x,y)=ax+by+cxy+d. In a spectral element method however, the interpolating and test functions are chosen to be polynomials of a very high order (typically e.g. of the 10th order in CFD applications). This guarantees the rapid convergence of the method. Furthermore, very efficient integration procedures must be used, since the number of integrations to be performed in numerical codes is big. Thus, high order Gauss integration quadrature is employed, since they achieve the highest accuracy with the smallest number of computations to be carried out. At the time there are some academic CFD codes based on the spectral element method and some more are currently under development, since the new time-stepping schemes arise in the scientific world.

* Boundary element method

In the boundary element method, the boundary occupied by the fluid is divided into a surface mesh.

* High-resolution discretization schemes

High-resolution schemes are used where shocks or discontinuities are present. Capturing sharp changes in the solution requires the use of second or higher-order numerical schemes that do not introduce spurious oscillations. This usually necessitates the application of flux limiters to ensure that the solution is total variation diminishing.

– Turbulence models

In computational modeling of turbulent flows, one common objective is to obtain a model that can predict quantities of interest, such as fluid velocity, for use in engineering designs of the system being modeled. For turbulent flows, the range of length scales and complexity of phenomena involved in turbulence make most modeling approaches prohibitively expensive; the resolution required to resolve all scales involved in turbulence is beyond what is computationally possible. The primary approach in such cases is to create numerical models to approximate unresolved phenomena. This section lists some commonly used computational models for turbulent flows.

Turbulence models can be classified based on computational expense, which corresponds to the range of scales that are modeled versus resolved (the more turbulent scales that are resolved, the finer the resolution of the simulation, and therefore the higher the computational cost). If a majority or all of the turbulent scales are not modeled, the computational cost is very low, but the tradeoff comes in the form of decreased accuracy.

In addition to the wide range of length and time scales and the associated computational cost, the governing equations of fluid dynamics contain a non-linear convection term and a non-linear and non-local pressure gradient term. These nonlinear equations must be solved numerically with the appropriate boundary and initial conditions.

* Reynolds-averaged Navier–Stokes

Reynolds-averaged Navier–Stokes (RANS) equations are the oldest approach to turbulence modeling. An ensemble version of the governing equations is solved, which introduces new apparent stresses known as Reynolds stresses. This adds a second order tensor of unknowns for which various models can provide different levels of closure. It is a common misconception that the RANS equations do not apply to flows with a time-varying mean flow because these equations are ‘time-averaged’. In fact, statistically unsteady (or non-stationary) flows can equally be treated. This is sometimes referred to as URANS. There is nothing inherent in Reynolds averaging to preclude this, but the turbulence models used to close the equations are valid only as long as the time over which these changes in the mean occur is large compared to the time scales of the turbulent motion containing most of the energy.

RANS models can be divided into two broad approaches:

* Boussinesq hypothesis

This method involves using an algebraic equation for the Reynolds stresses which include determining the turbulent viscosity, and depending on the level of sophistication of the model, solving transport equations for determining the turbulent kinetic energy and dissipation. Models include k-ε (Launder and Spalding), Mixing Length Model (Prandtl), and Zero Equation Model (Cebeci and Smith). The models available in this approach are often referred to by the number of transport equations associated with the method. For example, the Mixing Length model is a “Zero Equation” model because no transport equations are solved; the k-ɛ is a “Two Equation” model because two transport equations (one for k and one for ɛ) are solved.

* Reynolds stress model (RSM)

This approach attempts to actually solve transport equations for the Reynolds stresses. This means introduction of several transport equations for all the Reynolds stresses and hence this approach is much more costly in CPU effort.

* Large eddy simulation

Large eddy simulation (LES) is a technique in which the smallest scales of the flow are removed through a filtering operation, and their effect modeled using sub-grid scale models. This allows the largest and most important scales of the turbulence to be resolved, while greatly reducing the computational cost incurred by the smallest scales. This method requires greater computational resources than RANS methods, but is far cheaper than DNS.

Volume rendering of a non-premixed swirl flame as simulated by LES.
Volume rendering of a non-premixed swirl flame as simulated by LES.
* Detached eddy simulation

Detached eddy simulations (DES) is a modification of a RANS model in which the model switches to a sub-grid scale formulation in regions fine enough for LES calculations. Regions near solid boundaries and where the turbulent length scale is less than the maximum grid dimension are assigned the RANS mode of solution. As the turbulent length scale exceeds the grid dimension, the regions are solved using the LES mode. Therefore, the grid resolution for DES is not as demanding as pure LES, thereby considerably cutting down the cost of the computation. Though DES was initially formulated for the Spalart-Allmaras model (Spalart et al., 1997), it can be implemented with other RANS models (Strelets, 2001), by appropriately modifying the length scale which is explicitly or implicitly involved in the RANS model. So, while Spalart–Allmaras model-based DES acts as LES with a wall model, DES based on other models (like two equation models) behave as a hybrid RANS-LES model. Grid generation is more complicated than for a simple RANS or LES case due to the RANS-LES switch. DES is a non-zonal approach and provides a single smooth velocity field across the RANS and the LES regions of the solutions.

* Direct numerical simulation

Direct numerical simulation (DNS) resolves the entire range of turbulent length scales. This marginalizes the effect of models, but is extremely expensive. The computational cost is proportional to Re3. DNS is intractable for flows with complex geometries or flow configurations.

* Coherent vortex simulation

The coherent vortex simulation approach decomposes the turbulent flow field into a coherent part, consisting of organized vortical motion, and the incoherent part, which is the random background flow. This decomposition is done using wavelet filtering. The approach has much in common with LES, since it uses decomposition and resolves only the filtered portion, but different in that it does not use a linear, low-pass filter. Instead, the filtering operation is based on wavelets, and the filter can be adapted as the flow field evolves. Farge and Schneider tested the CVS method with two flow configurations and showed that the coherent portion of the flow exhibited the $ \displaystyle -\frac{{40}}{{39}}$ energy spectrum exhibited by the total flow, and corresponded to coherent structures (vortex tubes), while the incoherent parts of the flow composed homogeneous background noise, which exhibited no organized structures. Goldstein and Vasilyev applied the FDV model to large eddy simulation, but did not assume that the wavelet filter completely eliminated all coherent motions from the sub filter scales. By employing both LES and CVS filtering, they showed that the SFS dissipation was dominated by the SFS flow field’s coherent portion.

* PDF methods

Probability density function (PDF) methods for turbulence, first introduced by Lundgren, are based on tracking the one-point PDF of the velocity, fV(v;x,t) which gives the probability of the velocity at point x being between v and v+dv. This approach is analogous to the kinetic theory of gases, in which the macroscopic properties of a gas are described by a large number of particles. PDF methods are unique in that they can be applied in the framework of a number of different turbulence models; the main differences occur in the form of the PDF transport equation. For example, in the context of large eddy simulation, the PDF becomes the filtered PDF. PDF methods can also be used to describe chemical reactions, and are particularly useful for simulating chemically reacting flows because the chemical source term is closed and does not require a model. The PDF is commonly tracked by using Lagrangian particle methods; when combined with large eddy simulation, this leads to a Langevin equation for sub filter particle evolution.

* Vortex method

The vortex method is a grid-free technique for the simulation of turbulent flows. It uses vortices as the computational elements, mimicking the physical structures in turbulence. Vortex methods were developed as a grid-free methodology that would not be limited by the fundamental smoothing effects associated with grid-based methods. To be practical, however, vortex methods require means for rapidly computing velocities from the vortex elements – in other words they require the solution to a particular form of the N-body problem (in which the motion of N objects is tied to their mutual influences). A breakthrough came in the late 1980s with the development of the fast multipole method (FMM), an algorithm by V. Rokhlin (Yale) and L. Greengard (Courant Institute). This breakthrough paved the way to practical computation of the velocities from the vortex elements and is the basis of successful algorithms. They are especially well-suited to simulating filamentary motion, such as wisps of smoke, in real-time simulations such as video games, because of the fine detail achieved using minimal computation.

Software based on the vortex method offer a new means for solving tough fluid dynamics problems with minimal user intervention. All that is required is specification of problem geometry and setting of boundary and initial conditions. Among the significant advantages of this modern technology;

  • It is practically grid-free, thus eliminating numerous iterations associated with RANS and LES.
  • All problems are treated identically. No modeling or calibration inputs are required.
  • Time-series simulations, which are crucial for correct analysis of acoustics, are possible.
  • The small scale and large scale are accurately simulated at the same time.
* Vorticity confinement method

The vorticity confinement (VC) method is a Eulerian technique used in the simulation of turbulent wakes. It uses a solitary-wave like approach to produce a stable solution with no numerical spreading. VC can capture the small-scale features to within as few as 2 grid cells. Within these features, a nonlinear difference equation is solved as opposed to the finite difference equation. VC is similar to shock capturing methods, where conservation laws are satisfied, so that the essential integral quantities are accurately computed.

* Linear eddy model

The Linear eddy model is a technique used to simulate the convective mixing that takes place in turbulent flow. Specifically, it provides a mathematical way to describe the interactions of a scalar variable within the vector flow field. It is primarily used in one-dimensional representations of turbulent flow, since it can be applied across a wide range of length scales and Reynolds numbers. This model is generally used as a building block for more complicated flow representations, as it provides high resolution predictions that hold across a large range of flow conditions.

– Two-phase flow

The modeling of two-phase flow is still under development. Different methods have been proposed, including the Volume of fluid method, the level-set method and front tracking. These methods often involve a tradeoff between maintaining a sharp interface or conserving mass. This is crucial since the evaluation of the density, viscosity and surface tension is based on the values averaged over the interface. Lagrangian multiphase models, which are used for dispersed media, are based on solving the Lagrangian equation of motion for the dispersed phase.

Simulation of bubble horde using volume of fluid method
Simulation of bubble horde using volume of fluid method

– Solution algorithms

Discretization in the space produces a system of ordinary differential equations for unsteady problems and algebraic equations for steady problems. Implicit or semi-implicit methods are generally used to integrate the ordinary differential equations, producing a system of (usually) nonlinear algebraic equations. Applying a Newton or Picard iteration produces a system of linear equations which is nonsymmetric in the presence of advection and indefinite in the presence of incompressibility. Such systems, particularly in 3D, are frequently too large for direct solvers, so iterative methods are used, either stationary methods such as successive overrelaxation or Krylov subspace methods. Krylov methods such as GMRES, typically used with preconditioning, operate by minimizing the residual over successive subspaces generated by the preconditioned operator.

Multigrid has the advantage of asymptotically optimal performance on many problems. Traditional solvers and preconditioners are effective at reducing high-frequency components of the residual, but low-frequency components typically require many iterations to reduce. By operating on multiple scales, multigrid reduces all components of the residual by similar factors, leading to a mesh-independent number of iterations.

For indefinite systems, preconditioners such as incomplete LU factorization, additive Schwarz, and multigrid perform poorly or fail entirely, so the problem structure must be used for effective preconditioning. Methods commonly used in CFD are the SIMPLE and Uzawa algorithms which exhibit mesh-dependent convergence rates, but recent advances based on block LU factorization combined with multigrid for the resulting definite systems have led to preconditioners that deliver mesh-independent convergence rates.

– Unsteady aerodynamics

CFD made a major breakthrough in late 70s with the introduction of LTRAN2, a 2-D code to model oscillating airfoils based on transonic small perturbation theory by Ballhaus and associates. It uses a Murman-Cole switch algorithm for modeling the moving shock-waves. Later it was extended to 3-D with use of a rotated difference scheme by AFWAL/Boeing that resulted in LTRAN3.

– Biomedical engineering

CFD investigations are used to clarify the characteristics of aortic flow in detail that are otherwise invisible to experimental measurements. To analyze these conditions, CAD models of the human vascular system are extracted employing modern imaging techniques. A 3D model is reconstructed from this data and the fluid flow can be computed. Blood properties like Non-Newtonian behavior and realistic boundary conditions (e.g. systemic pressure) have to be taken into consideration. Therefore, making it possible to analyze and optimize the flow in the cardiovascular system for different applications.

Simulation of blood flow in a human aorta
Simulation of blood flow in a human aorta

– CPU versus GPU

Traditionally, CFD simulations are performed on CPUs. In a more recent trend, simulations are also performed on GPUs. These typically contain slower but more processors. For CFD algorithms that feature good parallelism performance (i.e. good speed-up by adding more cores) this can greatly reduce simulation times. Lattice-Boltzmann methods are a typical example of codes that scale well on GPUs.

Mesh Convergence

Multitasking is one of the plagues of the century that generally ends up with procrastination or failure. Therefore, having planned, segmented and sequenced tasks is much more appropriate to achieving goals: this has also been working for CFD. In order to conduct an analysis, the solution domain is split into multiple sub-domains which are called cells. The combination of these cells in the computational structure is named mesh.

Mesh of a Formula 1 car
Mesh of a Formula 1 car

The mesh as simplification of the domain is needed, because it is only possible to solve the mathematical model under the assumption of linearity. This means that we need to ensure that the behavior of the variables we want to solve for can assumed to be linear within each cell. This requirement also implies that a finer mesh (generated via mesh refinement steps) is needed for areas in the domain where the physical properties to be predicted are suspected to be highly volatile. Errors based on mesh structure are an often-encountered issue which results in the failure of the simulation.

This might happen because the mesh is too coarse and does not cover all effects that happen in this single element one by one, but rather cover multiple effects that then change as the mesh gets finer. Therefore, a study of independency needs to be carried out. The accuracy of the solution enormously depends on the mesh structure. To conduct accurate solutions and obtaining reliable results, the analyst has to be extremely careful on the type of cell, the number of cell and the computation time. The optimization of those restrictions is defined as mesh convergence which might be sorted as below:

  1. Generate a mesh structure that has a quite low number of elements and carries out analysis. Before, assure that the mesh quality and coverage of CAD model is reasonable to examination.
  2. Regenerate mesh structures with a higher number of elements. Conduct analysis again. Compare results in accordance with properties of examined case. For instance, if a case is an examination of internal flow through a channel, pressure drop at critical regions might be used as comparison.
  3. Keep on ascending to a number of elements where results converge satisfactorily with previous one(s).

Therefore, errors, based on mesh structure, can be eliminated and optimum value for number of elements might virtually be achieved as to optimize calculation time and necessary computation resources. An illustration is shown in Figure 4 that looks into static pressure change at imaginary region X through increase in number of elements. According to Figure 4, around 1,000,000 elements would have been sufficient to conduct a reliable study.

Example of mesh convergence analysis
Example of mesh convergence analysis

Convergence in Computational Fluid Dynamics

Creating a sculpture requires a highly talented artist with the ability to imagine the final product from the beginning. Yet a sculpture can be, for example, a simple piece of rock in the beginning, but might become an exceptional artwork in the end. A completely gradual processing throughout carving is an important issue to obtain the desired unique shape. Keep in mind that in every single process, some of the elements, such as stone particles, leftovers, are thrown away from the object. CFD also has a similar structure that relies on gradual processing during the analysis. In regions that are highly critical to the simulation results (for example a spoiler on a Formula 1 car) the mesh is refined into smaller elements to make the simulation more accurate.

Mesh refinement example
Mesh refinement example

Convergence is a major issue for computational analysis. The movement of fluid has a non-linear mathematical model with various complex models such as turbulence, phase change and mass transfer. Apart from the analytical solution, the numerical solution goes through an iterative scheme where results are obtained by the reduction of errors among previous stages. The differences between the last two values specify the error. When the absolute error is descending, the reliability of the result increases, which means that the result converges towards a stable solution.

The criteria for convergence vary with the mathematical models such as turbulence, multi-phase, etc. How do analysts decide when the solution is converged? Convergence should go on and on until a steady-state condition has been obtained, even if the aimed case is transient, which indicates results change through time. Convergence has to be realized for each time-step as if they all are a steady-state process. What are the criteria for convergence? The rate of accuracy (acceptable error), complexity of the case and calculation time have to be considered as major topics to carry out an optimal process. The residuals of equations, like stone leftovers, change over each iteration. As iterations get down to the threshold value, convergence is achieved. For a transient case, those processes have to be achieved for each of the time steps. Furthermore, convergence might be diversified as follows below:

  • Can be accelerated by parameters as initial conditions, under-relaxation and Courant number.
  • Does not always have to be correct, yet solution can converge, preferred mathematical model and mesh would be incorrect or have ambiguities.
  • Can be stabilized within several methods like reasonable mesh quality, mesh refinement, using discretization schemes first- to second-order.
  • Ensure that solution should be repeatable if necessary, as to refrain ambiguity.

Applications of Computational Fluid Dynamics

Where there is fluid, there is CFD. Having mentioned before, the initial stage to conduct a CFD simulation is specifying an appropriate mathematical model of reality. Rapprochements and assumptions give direction through solution processes to examine the case in the computational domain. For instance, fluid flow over a sphere / cylinder is a repetitive issue that has been taught by the lecturer as an example in fluid courses. The same phenomenon is virtually available in the movement of clouds in the atmosphere which is indeed tremendous.

Examples of Karman vortex streets; numerical result (top) and real-life example of clouds (bottom)
Examples of Karman vortex streets; numerical result (top) and real-life example of clouds (bottom)

– Incompressible and Compressible flow

If compressibility becomes a non-negligible factor, this type of analysis helps you to find solutions in a very robust and accurate way. One example would be a Large Eddy Simulation of flow around a cylinder.

– Laminar and Turbulent flow

Different turbulence models play a role in this type of analysis. A lot of computing power is required to solve turbulence simulations and its complex numerical models. The difficulty of turbulence is the simulation of changes over time. The entire domain where the simulation takes place needs to be recalculated after every time step.

The analysis of a ball valve is one possible application of a turbulent flow analysis.

– Mass and Thermal transport

Mass transport simulations include smoke propagation, passive scalar transport or gas distributions.  Heat exchanger simulations are one possible application.

– Different Types of CFD Applications

Computational Fluid Dynamics tools diversify in accordance with mathematical models, numerical methods, computational equipment and post-processing facilities. As a physical phenomenon could be modeled with completely different mathematical approaches, it would also be integrated with unlike numerical methods simultaneously. Thus, a conscious rapprochement is the essential factor on the path to developing CFD tools. There are several license-required commercial software solutions, though there are also open source projects available.

Websites: 1. Wikipedia 2. Simscale 3.
Methink 4. Bakker 5. ScienceDirect
Books and Articles: 1ba. M. Kawaguti, “Numerical Solution of the NS Equations for the Flow Around a Circular Cylinder at Reynolds Number 40”, 1953, Journal of Phy.Soc. Japan, vol. 8, pp. 747-757 / 2ba. Jameson, A., Schmidt, W. and Turkel, E., “Numerical Solution of the Euler Equations by Finite Volume Methods” / 3ba. Runge-Kutta, “Time-Stepping Schemes”, AIAA paper 81-1259, presented at the AIAA 14th Fluid and Plasma Dynamics Conference, Palo Alto California, 1981 / 4ba. Milne-Thomson, L.M. (1973). Theoretical Aerodynamics. Physics of Fluids A. 5. Dover Publications. p. 1023. ISBN 978-0-486-61980-4 / 5ba. Hess, J.L.; A.M.O. Smith (1967). “Calculation of Potential Flow About Arbitrary Bodies”. Progress in Aerospace Sciences. 8: 1–138. Bibcode:1967PrAeS…8….1H / 6ba. Hirt, C.W.; Nichols, B.D. (1981). “Volume of fluid (VOF) method for the dynamics of free boundaries”. Journal of Computational Physics.

R&D Website
R & D is the research and development center of Acamech that develops Mechanical Engineering Problems; including Renewable Energy Technology, Strctural Mechanics, Acoustics, Micro and Nano Scale especially MEMS, Fluid Flow and Thermal Systems, Computational Fluid Dynamics, Finite Element Analysis, Computer-Aided Design, Concepts and Relevant Technology.

Leave a Reply

Your email address will not be published. Required fields are marked *