• CFD, Fluid Flow, FEA, Heat/Mass Transfer: Feedbacks/Queries - fb@cfdyna.com

CFD Simulation Guidelines

CFD stands for Computational Fluid Dynamics - the branch of mechanical engineering dealing with Fluid Flow and Heat Transfer. Whether you are a fresh graduate with engineering degree or a seasoned professional with industry and design background, you will find plenty of technical stuffs on CFD analyses and CFD software (both commercial such as FLUENT, CFX, STAR-CCM+ and COMSOL as well as open-source Gmsh, OpenFOAM, Salome and ParaView) online. This page is an attempt to provide all relevant information in one place and organized into logical steps. However, one must note that the term CFD is also used in commodity trading where it stands for Contract for Differences.


How many of these names have you heard of? Click on the names to know more about them!

Open-source Companies
Open-source and Post-processors

Online Forum on CFD simulations, Service Providers, Freelancing Websites...

CFD Simulation related Websites

CFD, FEA, Turbulence

For flow and thermal simulation jobs, consultancy, training, automation and longer duration projects, reach me at fb@cfdyna.com. Our bill rate is INR 1,200/hour (INR 9,600/day) with minimum billing of 1 day. All simulations shall be carried out using Open-source Tools including documentations.

Step-by-step guide to CFD, Mesh Generation and OpenFOAM!

Realm of Numerical Simulations!

CFD simulations methods are explained with features of OpenFOAM and ANSYS CFX / FLUENT, STAR-CCM+ and COMSOL. Applications related to turbo-machines, multi-phase flows, aero-acoustics ... have been described along with turbulence modeling and verification of results.

Werner Heisenberg: When I meet God, I am going to ask him two questions: why relativity? And why turbulence? I really believe he will have answer for the first! […]

Aerodynamics

The flow over bluff bodies has applications ranging from electrical high tension wires to cars to aeroplanes.


Flavours of Numerical Methods

Topics in CFD and Structural Analyses


blog icon Visit the blog page to view topics on mathematics, human resource, textbook solutions, data analytics in stock market trading, online quizzes and much more ...


↹ Flow and Thermal Simulation Process and Methods

How do we approach flow and thermal simulations? What steps one need to follow to make first time right (FTR) solution? These are some of the questions which have been answered in following sections with references to multiple pages covering each topic such as Mesh Generation, Boundary Conditions and Turbulence Modeling. For beginners to experienced users, the checklists are needed to ensure correct results and few sample checklists have been created.

It is imperative to have good understanding of mathematical methods adopted to solve linear system of equations, role of relaxation factors in iterative methods and rate of convergence. This has been explained clearly for a simple equation which is valid of any complex set of equations.

Scope Definition

The table below summarize few key pointers to define a sound and quantifiable scope of CFD simulations. These are some of the topics which will lead the concerned stakeholders towards a well-thought scope. However, it may look a bit intimidating and intent is not to insist on all the information right at the beginning of the simulations.

Topic Description of requirement Examples Response
Objective What is the primary aim? Performance evaluation, Capacity upgrade, Design validation, Field failure, Safety, Emission ?
Problem statement What is expected output from CFD simulation? Δp, Temperature, Flow uniformity, Heat flux, Volume fraction, Air-borne noise, Stagnant zone ?
Prior investigation What information is available from past or ongoing study? Test report, Performance values, Field data, Hand calculations ?
Geometrical details What geometrical data is available for area of interest? 2D drawing, 3D CAD data, Site images, Concept sketches ?
Working fluid Water, Air, Non-Newtonian, Mixture (Liquid-Solid, Liquid-Gas, Gas-Gas) Homogeneous mixture, Two-phase flow, Emulsification, Stirred tank ?
Time Domain Steady (not dependent on time) or Transient (p, v, T vary with time) Warm-up, cool-down, moving walls, varying flow rate ?
System or Component Is this a component level or system level simulation? [Comp. level]: pump, fan in ducts, PCB
[Sys. level]: Engine cooling system, Data centre
?

As mentioned earlier, it may look like a daunting task for a person not familiar with finer details of numerical simulations. Hence, a generic list of ingredients needed for CFD simulations can be used to start the discussion. This is summarized below.

  1. Brief description: what is the expected output from simulations - in other words what are design validation requirements?
  2. CAD geometry: STEP, IGES, Parasolid format [data in Creo, Solidworks, CATIA... can also be directly imported into modern pre-processors]
  3. Operating conditions: what physical conditions are to be used for simulations - e.g. how many designs (and modifications), fan speeds, fluid temperature, flow rates, with or without heat transfer...
  4. Is there any calibration or benchmarking with results of particular test conditions and set-up needed?
  5. Any specific software and version to be used?
Step-by-Step Guide to CFD Simulation

CFD Simulation Steps

A typical division of a CFD simulation or as a matter of fact any numerical simulation is in 3 categories namely Pre-Processing, Solution and Post-Processing. However, for beginners into this field, it does not convey much insight. The complete simulation process is explained here in a step-by-step sequence to make it easy to understand.

  • Step-01: Problem Definition - Before you jump on the simulation methods, define the problem statement and expected output, accuracy and turn-around time.

    Problem Definition

    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?

    Logos of Companies

  • Step-02: 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 pre-processor 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.

    CAD Geometry

  • Step-03: Geometry Simplification for Simulation - This step is not required if geometry has been created especially 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 shrink-fit or CAD modeling error... The purpose it to get a water-tight 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.
  • Step-04: Mesh Generation - this is a process of dividing the overall volume into smaller pieces or smaller volumes. A given computational domain may have infinite degree-of-freedom (DOF) or locations where solution need to be found out. However, numerical methods need to reduce the DOF to 'some' finite number of DOF and hence this process is also called Finite-Element Method or Finite-Volume Method. Additional information can be found here for mesh generation methods and recommendations.

    CD Nozzle

  • Step-05: 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 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.
  • Step-06: Solver Setting - The set-up is now ready for solution. All the numerical simulations get reduced to the format [A]{x} = {b} where [A] is large matrix and {x} is the unknowns we are looking to estimate.

    Solver Setting Step

    In this step, the method to model turbulence, wall function, discretization scheme for differential operators such as du/dy and d2u/dy2 is specified. Refer to this page for methods and recommendations. A detailed description of turbulence modeling can be found here.
  • Step-07: 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.

    Solution Troubleshooting

    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. Recommended approach for STEADY STATE simulations is:
    1. Start with first order discretization schemes for all variables, energy OFF, constant material properties (not a function of temperature and pressure)
    2. Change discretization scheme to second order for momentum, keeps others as first order
    3. Reduce under-relaxation factor for momentum to 0.3 and 0.2 for turbulence parameters
    4. Turn energy ON (solver for temperature)
    5. Change material properties to function of temperature and pressure as the case may be
    6. Change discretization scheme for turbulence parameters to second order
  • Step-08: Post-processing - 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 post-processing.

    Post-processing and report preparation

    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 Back-Step" and can be found here! If you have access to ANSYS FLUENT, a step-by-step 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 open-source 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:

    Zones of a Boundary Layer

    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:

    MRF

    Which methods are available to model rotating components and turbo-machines? 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?

    injection Points Python Spiral

    Sometimes mixing of two fluids need to be modeled. These phenomena are called multi-phase 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 multi-phase 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 = 5x2 + 25 when y = 1000. We all know a direct method:

x = [(y - 25)/5]0.5

  Let's try an iterative method using algorithm xi+1 = xi × (1 - Δyi). Here, i refers to the iteration number and Δy refers to deviation of the calculated value of y from desired value of y = 1000, Res-Norm-1 refers to the absolute value where the residual % is normalized (divided) by the residual % obtained in 0th iteration based on guess value, Res-Norm-3 refers to the absolute value when the residual % is normalized (divided) by the residual % obtained in 0th ~ 2nd iterations.

Iteration, i xi yiΔyi Res-Norm-1Res-Norm-3
Initial Guess10.00525.0+475.0+47.5%--
114.751112.8-112.8-11.3%0.2375-
213.09881.2+118.8+11.9%0.2501-
314.641096.7-96.70-9.67%0.20360.6032
413.22899.45+100.6+10.1%0.21170.6272
514.551084.1-84.14-8.41%0.17710.5248
613.33913.40+86.60+8.66%0.18230.5401
714.481073.9-73.93-7.39%0.15560.4611
813.41924.57+75.43+7.54% 0.15880.4705
914.421065.4-65.40-6.54%0.13770.4079
1013.48933.77+66.23+6.62%0.13940.4131
1114.371058.1-58.14-5.81%0.12240.3626
1213.54941.50+58.50+5.85%0.12320.3649
1314.331051.9-51.86-5.19%0.10920.3235
1413.59948.11+51.89+5.19%0.10920.3236
1514.291046.4-46.39-4.64%0.09770.2894
1613.63953.82+46.18+4.62%0.09720.2881
1714.261041.6-41.59-4.16%0.08760.2594
1813.67958.79+41.21+4.12%0.08680.2570
1914.231037.3-37.34-3.73%0.07860.2329
2013.70963.15+36.85+3.68%0.07760.2298
2114.201033.6-33.56-3.36%0.07070.2094
2213.73967.00+33.00+3.30%0.06950.2059
2314.181030.2-30.20-3.02%0.06360.1884
2413.75970.40+29.60+2.96%0.06230.1846
2514.161027.2-27.19-2.72%0.05720.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 speed-up 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 xi+1 is calculated from xi. 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 under-relaxation factor as the value is < 1.0).

xi+1 = [xi × (1 - Δyi)] * 0.8 + xi * 0.2

Iteration, i xi yi Δyi Res-Norm-1 Res-Norm-3
Initial Guess10.00525.0475.047.5%--
114.751112.8-112.8-11.3%0.2375-
213.42925.30+74.70+7.50%0.1572-
314.221036.1-36.11-3.61%0.07600.2480
413.81978.54+21.50+2.10%0.04520.1474
514.051011.6-11.56-1.16%0.02430.0794
613.92993.39+6.61+0.66%0.01390.0454
713.991003.7-3.66-0.37%0.00770.0251
813.95997.94+2.06+0.21%0.00430.0141
913.971001.1-1.15-0.11%0.00240.0079
1013.96999.35+0.65+0.06%0.00140.0044
1113.971000.4-0.36-0.04%0.00080.0025
1213.96999.80+0.20+0.02%0.00040.0014
1313.971000.1-0.11-0.01%0.00020.0008
1413.96999.94+0.06+0.01%0.00010.0004
1513.961000.0-0.04+0.00%0.00010.0002
1613.96999.98+0.02+0.00%0.00000.0001
1713.961000.0-0.01+0.00%0.00000.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 column-4 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*xP + B*xN + C*xS + D*xE + E*xW = FP 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*xPn + B*xNn + C*xSn + D*xEn + E*xWn - FPn] where n is the iteration number.

  • Scaled Residuals: The values reported in column-5 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 well-defined 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 column-6 and column-7 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, area-averaged 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 especially 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 bluff-bodies, 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.

A checklist and comparison sheet should be used when a new design is worked out or major upgrades in designs are under considerations.

S. No.Performance Indicator Baseline Upgrade
1Pressure drop≤ 2.0 [kPa] ≤ 1.5 [kPa]
2Temperature rise≤ 15 [K]≤ 12 [K]
3Wetted surface area ≥ 1.0 [m2]≥ 1.5 [m2]
4Capacity or volumetric size ≤ 5.0 [cm3]≥ 7.5 [cm3]
5Power consumption ≤ 50 [W]≤ 40 [W]

Application of CFD in terms of physics, application and industry: sophisticated computer-based analysis technique to evaluate various design alternatives, identify flow and temperature problems, develop solutions and formulate physical validation strategies. Troubleshoot existing design, generate smarter solutions and more optimized design with shorter development cycle to meet ever-increasing demands of working performance and environmental concerns.

S. No.Simulation CategoryPhysics ApplicationIndustry
01Basic
  1. Incompressible Flows
  2. Flow with Heat Transfer
  1. Car external aerodynamics
  2. HVAC: Automotive, Plants, Warehouse
  3. Chargers, Converters, Inverters, LEDs, Battery Packs, Cold Plate
  4. Intake Manifolds, Intake Valves
Automotive, FMCG, Aerospace, Recreational Vehicles, Data Centres, Valves
02Rotating Devices
  1. Incompressible Flows
  2. Swirling Flows
  1. Pumps, fans and blower
  2. Compressors
  3. Passenger car HVAC Module
  4. Cyclone Separators
  5. Turbo-chargers, Motors, Turbo-Generators
Automotive, Home Appliances, Whitegoods, Pharmaceuticals, FMCG, Dairy Industry
03 Multi-phase Flows
  1. Free surface flows
  2. Cavitation
  1. Filling-Emptying
  2. Liquid Sloshing
  3. Erosion Modeling
  4. Defrosting, Defogging, Deicing
  5. Free Shear Flows, Air Curtains
  6. Fume (Smoke) Dispersion
FMCG, Buildings and Real-Estate, Oil and Gas Industry, Medical Devices, Pharmaceuticals
04Advanced Methods
  1. Compressible Flows
  2. Supersonic Flows
  3. Discrete Phase Modeling
  4. Jets and Sprays
  5. Radiative Heat Transfer
  6. Buoyancy-driven Flows
  7. Non-Newtonian Flows
  1. External aerodynamics
  2. Noise and Acoustics
  3. Fuel Injectors
  4. Engine Cooling System
  5. Hot Soak-down, Cold Warm-up
  6. Thermal Comfort Modeling
  7. Self-Priming Pumps
  8. Aerosol Dispersion - Virus Spread
Aerospace, Combustion Engines, Turbines, Automotive, Off-Highway Vehicles, Recreational Vehicles, Heat Recovery Devices, Car Parkings, Road Tunnels
05 Multiphysics
  1. Mechanical - Electrical
  2. Fluid - Thermal - Structural
  3. Chemical - Thermal - Structural
  1. Fire Modeling
  2. FSI: Fluid-Structural Interaction
  3. Radiator Thermal Shock
  4. Engine Thermo-Structural
  5. Combustion and Emissions (NOX, SOX)
Medical Devices, Elevators, Aerospace, Automotive, Burners / Combustors
06 Electronic Cooling
  1. Printed Circuit Boards
  2. Chips and MOSFETs
  3. Thermo-Electric Devices
  4. Cabinets and Enclosures
  5. LED
  1. Natural Convection
  2. Radiation Modeling
  3. Orthotropic Properties
Home Appliances, Whitegoods, ICEPAK, FloTHERM, Sigma Technologies

This list is no way unique and can be arranged in too many ways to suit specific purpose. The attempt is to cover entire gamut of application in a short yet detailed manner.


Allied Methods


DEM

Discrete Element Method: The well-known FEM and FVM method is a continuum-based 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.

    Particle separation DEM


ROM

ROM stands for Reduced Ordered Modeling which is an attempt to expedite the simulation speed with some compromise in accuracy. This is also known as MOR (Model Order Reduction). At the very basic level, ROMs can be the simplification of full 3D simulations to smaller equivalent domains such as 2D section, symmetric and periodic domains. One of the high fidelity ROM is simulation of flow through pipe as a rectangular plane than can be rotated about one of its longer sides to create the full cylinder. Though not all industrial applications can be simplified to this extent!

Reduced order modeling is a fairly generic term and the semantics changes from field to field (e.g. computer science, mechanical engineering, testing...) and includes all methodologies and algorithms that aim to reduce the complexity of a problem.

Excerpt from web.stanford.edu/group/frg/active_research_themes/reducedmodel.html: "Reduced-order models (ROMs) are usually thought of as computationally inexpensive mathematical representations that offer the potential for near real-time analysis. While most ROMs can operate in near real-time, their construction can however be computationally expensive as it requires accumulating a large number of system responses to input excitations."

As per mathworks.com/discovery/reduced-order-modeling.htm: Many techniques such as singular value decomposition (SVD), eigenvalue decomposition (EVD), principal component analysis (PCA), lookup tables, interpolation and curve fitting are also useful as part of model order reduction. As per in.mathworks.com/discovery/reduced-order-modeling.html - "Reduced order modeling (ROM) and model order reduction (MOR) are techniques for reducing the computational complexity or storage requirement of a computer model, while preserving the expected fidelity within a controlled error. Working with surrogate models can simplify analysis and control design. Scientists and engineers use ROM techniques to create system-level simulations, design control systems, optimize product designs, and build digital twin applications."

Excerpts from "An Overview of Reduced Order Modeling Techniques for Safety Applications" by D. Mandelli et. al. The safety modeling of a nuclear system is not only a thermal-hydraulic problem but several other models are required: neutron transport, thermo-mechanics, chemistry, fracture propagation... Note how the overall problem is not only multi-physics but also multi-scale both in the spatial scale but also in the temporal scale. The concept of Reduced Order Modeling can be categorised into three main categories:

  • Reduced physics: use of simulator codes that employ simplified physics problems. An example is the use of diffusion codes to solve neutronic problems instead of transport codes. In this category we also include the possibility to use in the same simulation run high and low fidelity models depending on the boundary conditions of the simulation.
  • Reduced dimensionality: a simulation run can be seen as a trajectory in the phase space and a single point in the input space. The dimensionality of these spaces can be very high for the complex analyses. This category includes all methods than aim to reduce the dimensionality of these spaces and project the original problem into the reduced space.
  • Surrogate model: surrogate models are mathematical objects that emulate the behavior of a code by learning its input/output relations and reconstructing such relations through a regression/interpolation based approach.

Machine Learning (ML) in Numerical Simulations: It is proven that Monte Carlo method can be used to solve Conduction Equation without mesh generation and matrix inversion. However, same does not apply to FEA or CFD. Any ML and/or ANN (Artificial Neural Network) does not solve the underlying physics - they act as what is called "surrogate models" to work on the data generated by FEA or CFD or any numerical simulations, play the role of a decision maker. In case of Fluid Dynamics, ANN cannot be even a surrogate model to turbulence modeling (please note the word 'modeling' after turbulence) because if turbulence can be predicted, numerically computed or interpolated, the methods of CFD simulations would not remain worthy of a ML or NN.


FSI - VIV - FIV

Fluid-Structure Interaction - this is a multi-physics 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 two-way coupling where CFD and FEM method need to share the information to define time-dependent 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 a-priori. Application includes dynamic response of reed valves in a hermetically sealed reciprocating compressor, motion of swing check valve in water pipelines, equilibrium position of check valves...

Excerpts from Fluid-Structure Interaction (FSI) case study of a cantilever using OpenFOAM and DEAL.II with application to VIV, Johan Lorentzon, June 17, 2009: The Fluid-Structure Interaction (FSI) appears as a physical phenomenon in engineering such as static load, drag, and Flow-Induced Vibrations (FIV) or Flow-Induced Pulsations (FIP). FIV is further classified into flutter, galloping and vortex-induced vibrations (VIV). The static load on the structure originates mainly from the total pressure difference between front and wake side of the structure. The dynamic pressure can, due to instability in the structure where the damping is less than the energy transfered by the vortex-induced wake, create a resonance with the structure with a negative damping known as flutter. Typical examples of applications in which flutter is considered are aircraft wings, tall buildings, rotating blades and long-span bridges. The VIV concerns bridges, chimneys, heat exchanger tubes, power lines and undersea constructions such as sea cables/risers or biomechanical applications as blood vessels.

Reference: Flow Induced Pulsations Caused By Split Flow In A T-Branch Connection by Stefan P.C. Belfroid et al. ASME 93884 - "By using a vortex blob method it has been demonstrated that cross flows even up to 25% of the main flow lead to significant FIP source strength. With increasing cross flow merging from the side branch into the main flow, the maximum source strength shifts to slightly lower Strouhal number due to the increased convective velocity of the vortex structures. The source strength decreases rapidly for increasing cross flow. As a rule of thumb, the source strength becomes negligible for cross flows larger than 25% of the total flow."

FIP Tee Junction

When the mean flow is diverted completely into the side branch (case 'b') no pulsating flow is encountered. However, for cases 'a' and 'c', strong pulsations have been measured. For calculating acoustic resonances of a piping system, one dimensional simulation tools are quick and reasonable accurate. A full acoustic solution including travelling and standing waves can be carried out in the flow system.


CAA

Computational Aero-Acoustics - this is another multi-physics 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 post-processing activity where the fluctuating component of pressure is converted into noise using density-based 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 available 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.

MHD

Magneto-Hydro-Dynamics (MHD) or Hydro-Magnetics is the physical-mathematical macroscopic theory that deals with dynamics of magnetic fields in electrically conducting fluids [ plasmas, liquid metals, salt water and electrolytes]. It does not apply to effect of magentic field on fluids containing magnetic particles. The presence of magnetic fields leads to forces that in turn act on the fluid thereby potentially altering the topology and strength of the magnetic fields themselves.

The set of equations that describe MHD are a combination of the Navier–Stokes equations of fluid dynamics and Maxwell’s equations of electromagnetism. These differential equations have to be solved simultaneously, either analytically or numerically.

The study of ferromagnetic particle suspensions is relevant for better understanding of wet-state separation processes in mineral processing and recycling. Stokesian dynamics is a particle-based, mesh-free method developed for nm to mm size particles in static or flowing liquid suspensions. LB-LES-DEM approach (Lattice-BoltzmannLarge Eddy Simulation-Discrete Element Method) has been investigated to evaluate the combined effect of liquid phase flow and magnetic effects though limited number of particles can be used due to the high computational cost.


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 (Navier-Stokes equation).

PowerFLOW from DASSAULT SYSTÈMES

Excerpts: "High fidelity transient Lattice Boltzmann based solution, accurate across most flow regimes (laminar to transonic) to solve the most complex CFD design problems in Transportation & Mobility and Aerospace & Defense."


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 × Q2 where Q is flow rate. A network is created with nodes and branches where nodes represent the connection of two or more branches (the flow resistances). Similar to Kirchhoff'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. Hardy-Cross method is one of the widely used method to solve flow network models. 1D model for a single cylinder engine in GT-Power described below is an example of FNM.

1D: single cylinder GT-Power


Modelica

Excerpts - "The Modelica Language is a non-proprietary, object-oriented, equation based language to conveniently model complex physical systems containing, e.g. mechanical, electrical, electronic, hydraulic, thermal, control, electric power or process-oriented subcomponents."

DYMOLA (DYnamic MOdelling LAboratory) is a user interface and Modelica language compiler owned and developed solely by Dassault Systemes. DYMOLA enables the user to write, compile and simulate Modelica based models. There are other Commercial Modelica Simulation Environments such as ANSYS Simplorer, SimulationX from ESI ITI GmbH and Simcenter Amesim from Siemens. There is an open-source and free simulation environment tool for Modelica named OpenModelica.


SU2

SU2 stands for Stanford University Unstructured. This is an open-source solver having its own mesh format. The other way to interact with SU2 is to use of CGNS for mesh and ParaView / Tecplot for post-processing. The work-flow of SU2 is as follows. It is a text-based solver similar to OpenFOAM where a single configuration file is equivalent to the 0, controlDict, transportProperties and fvSolution dictionaries. A converter from CGNS to the SU2 format has been built into SU2 and hence the mesh does not need to be explicitly converted into SU2 format.

Simulation ActivityMethod
Mesh GenerationSU2 can be used to generate for simple shapes. For most practical operations,mesh has to be generated in other pre-processors such as OpenFOAM, GMSH, Pointwise, ICEM CFD...
Convert mesh into SU2 formatGMSH and CENTAUR can save the mesh directly in SU2 format. The Python script written by Laurent Berdoulat can be found here to convert OpenFOAM mesh into SU2 format. The input file template for 2D meshes is here and the input file template for 3D meshes is here.
Convert mesh into SU2 formatFor commercial programs, name boundaries (used as marker tags in SU2) appropriately before exporting the mesh into CGNS format. Use MESH_FORMAT= CGNS* in configuration file.
Solver SettingDefine material, apply boundary conditions, select solver type, set iteration controls and relaxation factors... as per pre-defined format in a configuration file.
Post-processing settingsSet the output format for results: OUTPUT_FORMAT= TECPLOT or OUTPUT_FORMAT= PARAVIEW, OUTPUT_FILES = PARAVIEW_ASCII
Make the simulation runF:\SU2\SU2_CFD.exe wedgeShock.cfg, for parallel computing: mpirun -n 8 SU2_CFD.exe wedgeShock.cfg or use Parallel Computation Script - parallel_computation.py
Result Filesflow.dat or flow.vtk - full volume flow solution
surface_flow.dat or surface_flow.vtk - flow solution along the wall surfaces
surface_flow.csv - file containing values along the wall surfaces
restart_flow.dat - restart** file in an internal format for restarting this simulation in SU2
history.dat or history.csv - file containing the convergence history information
SU2_CFD (Computational Fluid Dynamics Code): Solves direct, adjoint and linearized problems for the Euler, Navier-Stokes, and Reynolds-Averaged Navier-Stokes (RANS) equation sets, among many others.
SU2_DOT (Gradient Projection Code): Uses the surface sensitivity, the flow solution and the definition of the geometrical variable to evaluate the derivative of a particular functional (e.g. drag, lift...).
SU2_DEF (Mesh Deformation Code): Computes the geometrical deformation of an aerodynamic surface and the surrounding volumetric grid. Once the type of deformation is defined, SU2_DEF performs the grid deformation by solving the linear elasticity equations on the volume grid.
SU2_MSH (Mesh Adaptation Code): Performs grid adaptation using various techniques based on an analysis of a converged flow solution, adjoint solution and linearized problem to strategically refine the mesh about key flow features.
SU2_SOL (Solution Export Code): Generates the volumetric and surface solution files from SU2 restart files.
SU2_GEO (Geometry Definition Code): Geometry preprocessing and definition code. In particular, this module performs the calculation of geometric constraints for shape optimization.

*SU2 will not use any specific boundary conditions that are embedded within the CGNS mesh. However, it will use the names given to each boundary as the marker tags.

**The restart files can be used as input to generate the visualization files using SU2_SOL module. This will be done automatically if the parallel_computation.py script is used, but it can also be done manually at any time by calling SU2_SOL.


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. These phenomena is governed by Joukowsky equation: Δp = [± ρ c Δv] where ρ is the density of fluid, c is (water-hammer) wave speed and δv is change in velocity. Similarly, a pump tripping or shut-down is equivalent to valve closure which also causes water hammer - which is equivalent to opening of a control valve to increase the flow rate in the pipe. Note:

  • Pressure rise in maximum for valve closing time ≤ 2L/c where L = length of the pipe
  • On the other hand, maximum pressure rise cannot exceed 2 × static head in the pipe
  • The velocity of pressure wave when fluid is contained in a pipe or conduit is < the wave speed in unconfined media, due to elasticity of pipe walls.

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 steady state flow field. In industry, there are many 1D tools which are based on Methods of Characteristics (MOC) to solve water hammer phenomena. Some of the tools are WANDA, HAMMER, AFT Impulse...


Erosion

Erosion Zone

[Image Reference: Springer] Loss of material due to impact of solid particles. For example, sand particles hitting at high velocity of a solid wall causing loss of material by "cutting or chipping mechanism" and/or by "deformation mechanism". Erosion may also occur by abrasive process where particle slip along the wall causing rubbing action between particles and the wall. Erosion can also occur due to cavitation and liquid particle drop falling over a wall. The erosion rate is function of impact velocity, impact angle, particle properties (sharpness ratio), hardness of the wall and hardness of the solid particle. Experiments have shown that there is exponential correlation between particle impact velocity and rate of erosion.

The erosion process is not a mathematically well-defined problem and the estimation of erosion location as well as erosion rate relies heavily on empirical correlations such as erosion rate E = C × FS × Vn × F(θ) where C and n are constants, FS is a particle shape factor (~1.0 for particles with sharp corners say diamond, cone and square shapes, ~0.5 for semi rounded shapes such as cylinders, disks and ellipsoids and ~0.2 for fully rounded shapes such as spheres). V is impact velocity and F(θ) is impact angle function. The "paint erosion test" is a quick and reliable approach to find "erosion pattern". It is analogous to measurement of contact area for bolted joints.

Erosion modeling using CFD: This is a 3-step process described below.

  1. Flow modeling: this step is similar to any typical CFD simulation assuming no solid particle in the flow. This is example of one way coupling where the presence of (discrete phase) solid particles is assumed not to affect the continuous phase.
  2. Lagrangian particle tracking - after steady flow is reached (Eulerian flow field is established), particles are injected into the flow and tracked based on drag, rebound and collision phenomena.
  3. Estimation of rate of erosion, loss of material per unit surface area per unit time (kg/m2/s).

Erosion function ductile and brittle material


CFD for Pneumatic Conveying Process: The viscous and pressure forces cause the particles to move against friction and gravity. The minimum velocity required to make the particles move is called "particle pick-up velocity". Saltation velocity is the value at which majority of the solids become entrained in the airstream. Air volume flow rate optimization is the critical process to achieve smooth conveying without slugging and plugging of the conduits. The conveying process also has to deal with "dilute phase" (air velocity > saltation velocity, suspension flows) and "dense phase" (air velocity < saltation velocity) conveying. The dilute and dense phases are also distinguished by solid loading (volume or mass fraction of the solids). The assessment of conveying capacity is similar to multiphase flow simulation where the solid phase is treated as continuum. Applications includes conveying of coffee beans, grains, granular sugar, cement powder, lime, spices...

air-assisted pneumatic Conveyorg

Choking Velocity - This equals the actual gas velocity in a vertical pipeline at which particles in a homogeneous mixture with the conveying gas settles out of the gas stream. Bulk Density (Fluidized) - it is the apparent (equivalent) density of a material in its fluidized or suspension state. It is generally lower than either the packed or loose bulk density due to the air absorbed into the voids. Angle of Repose - The angle of repose of a material is the angle formed between the horizontal and sloping surface of a piled material, which has been allowed to form naturally without any conditioning. This helps estimate the slope of the channel for gravity assisted conveying. Terminal Gas Velocity - The terminal gas velocity in a pneumatic conveying system is the velocity of the gas as it exits the system. It is also known as the ending gas velocity and conveying line exit velocity. Floodability - describs the tendency of a material to aerate and act as a fluid, squeeze material quickly in your fist - if it squirts through fingers, then it is floodable. Floodable materials are difficult to restrain in controlled feeding applications. Pneumatic Conveying Design Guide by David Mills is another good reference to start with.


Atomization and Spray Jets

Nozzle Spray Atomization

The near- and far-field break-up and atomization of a liquid jet (aerodynamic break-up of liquids) in multi-phase fluid systems fits in the broader category of multi-phase fluid flow. Jet and spray formation has many applications such as inhalers in medical applications, fire-extinguishers, Diesel engines, water chillers... Primary atomization of the liquid jet occurs near the outlet of the nozzle. Turbulent fluctuations of the liquid jet induce perturbations on the jet surface, which grow and break the jet into droplets. The length scale of turbulence is the dominant length scale of atomization, which also determines the resulting droplet size (Chryssakis et al.)1.

2In the plain-orifice atomizer model, the operation mode of the nozzle is determined based on the Reynolds number, the cavitation number and the critical values for the inception of cavitation and flipping (ANSYS, 2016). The Reynolds number based on hydraulic head is defined as ReH = (D ρL / μ) [2(p1 − p2)/ρL]0.5 where D is the nozzle diameter, ρL is the liquid density and μ is the liquid viscosity. The upstream and downstream pressures are denoted by p1 and p2, respectively. The cavitation number is defined as K = (p1 − pv)/(p1 − pv) where pv is the vapour pressure of liquid at operating temperature. For short, sharp-edged nozzles, the inception of cavitation occurs approximately at cavitation number K ~ 1.9.

1Chryssakis, C.A., Assanis, D.N. and Tanner, F.X., 2011. Atomization Models, in Handbook of Atomization and Sprays, Ed. N. Ashgriz, Springer, New York, USA.

2ANSYS Fluent Theory Guide, section Atomizer Model Theory

Additional information on Jet Break-up can be found here.


DOE in CFD: When CFD simulations are used to check impact of various parameters such as mesh size, turbulence model / wall function, inlet boundary condition, discretization scheme and P-V coupling, a Design-of-Experiments can be used to make the evaluation more scientific and robust. Here k = number of factors = 5. If each of these parameters considered are at two levels (two types or values), there are total 25 = 32 runs to be completed. It have been observed that the interactions involving > 3 factors are insignificant and a fractional factorial can be used.

One of the key concept in both CFD and DOE is DOF (Degrees-of-freedom). In order to understand a classroom with capacity N. The first student has N choices to take a seat, the second student has N-1 options and so on. The last student will have no choice. Thus, in case of N dataset, the DOF is N-1.

A full quadratic model contains all main effects, all quadratic effects and all possible two-way interactions. For a set of k-factors, the number of terms in a full quadratic model is N = (k+1)(k+2)/2. Thus, for k = 4, N = 15 ad for k = 5, N = 21. Central Composite Design (CCD) and Box-Behnken Design (BBD) are used to reduce the number of terms. This is part of the family known as Response Surface Model (RSM) which looks for quadratic and/or higher order trends as well as assumes that all variables are significant.

The interactions are designated as main effects, 2-way interactions, 3-way interactions and so on. The number of runs for main effect = kC1 = k, number of runs for 2-way interactions = kC2 = k/2 × (k-1), number of runs for 3-way interactions = kC3 = k/6 × (k-1) × (k-2) ... This is also known as degrees-of-freedom where each main effect will have two degrees of freedom when each factor has three levels. Each two-factor interaction has (3−1) × (3−1) = 4 degrees of freedom. The ABC interaction has (3−1) ×(3−1)×(3−1) = 8 degrees of freedom.

When there is a curvilinear (non-linear) relation between the response and a quantitative factor, it is not possible to detect such a curvature effect with two levels. A mix-level design refers to experiments where factors have unequal levels. For example, factor A has two levels, factor B has 3 levels and factor C has 4 levels.

Fractional Factorial Example

This table is rich in information and following explanation help decode them.
  • The 23-1 design is formed by selecting only those treatment combinations that have a + in the ABC column. Thus, ABC is called the GENERATOR of this particular fractions. Sometimes the generator such as ABC is also called WORD.
  • Furthermore, the identity column I is always +, so we may call I = ABC, the defining relation for the fractional design.
  • In general, the defining relation for a fractional factorial will always be the set of all columns that are equal to the identity column I.
  • We observe that values in columns are: A = BC, B = AC, C = AB. Hence, it is impossible to differentiate between A and BC, B and AC and C and AB. In fact when we estimate A, B, and C we are really estimating A+BC, B+AC, C+AB. Two or more effects that have this property are called ALIASES. In the above table, A and BC are aliases, B and AC are aliases and so on. When it is said that A is alias to BC, it means they are exactly equal.
  • Rules of finding aliases:
    • Step-1: Find the generating relation, I. Here, I = ABC.
    • Step-2: Multiply both sides of the generating relation by an effect. E.g. B × I = B × ABC.
    • Step-3: Note B × I = B. Delete all square terms arising out of this multiplication. That is, B = AB2C results in B = AC.
    • Thus, for a 4 factor factorial, alias of AB can be calculated as: AB × I = AB × ABCD = A2B2CD. That is: AB = CD.

Two most important properties of any factorial designs are 'balance' and 'orthogonality'. Balance requires that each level of a factor is run the same number of times in an experiment. For example, a design with 4 factors, first with 2 levels, second with 3 levels, third with 4 levels and the last with 5 levels. To generate a balanced design, at least 60 runs are needed. Orthogonality refers to pairwise linearly independent columns, useful for assessing factor significance. In vector mechanics, two vectors are orthogonal if the dot product of the vector and hence the sum of the products of their corresponding elements is 0.

Taguchi Experiment

Let's consider an experiment that has 4 variables each at 3 levels. A full factorial experiment would require 34 = 81 experiments. A Taguchi experiment with a L9(34-2) orthogonal array will require only 9 runs. Thus, Taguchi's orthogonal arrays are highly fractional factorial designs. Here L9 means array requires 9 runs. All Taguchi's arrays are orthogonal but highly unbalanced.

RunFactor-AFactor-BFactor-CFactor-D
11111
21222
31333
42123
52231
62312
73132
83213
93321

Adjoint Solvers

Adjoint: In linear algebra, Adjoint or Adjugate refers to conjugate transpose of a matrix. For a square matrix where inverse exists, AH ≡ adj[A] := Inverse[A] Det[A] An analogous concept applied to an operator instead of a matrix is also known as the Hermitian conjugate.The adjoint operator is widely used in both Sturm-Liouville theory and quantum mechanics.

adjoint Matrix Example

Optimization of shape and size in CFD simulations require incremental change in geometry and re-meshing. The CFD programs are now incroporating adjoint solvers which runs in addition to the fluid solver to find an optimum geometry to minimize (e.g. pressure drop) or maximize (e.g. flow uniformity) specified field variables. Both the Single Object Optimization (SOO) and Multiple Objectives Optimization (MOO) are possible.

Similar concepts Full Order Modeling (FOM) and Reduced Order Modeling (ROM) which are based on mathematical concepts such as MDO (Multidisciplinary Design Optimization), Higher Dimension Model (HDM), Proper Orthogonal Decomposition (POD) or Singular Value Decomposition (SVD).

There are 4 types of optimization techniques:

Optimization Method Salient Features
0D | Empirical
  1. Optimization by hand based on engineering knowledge and experience
  2. Iterative and sequential where design concepts are refined based on results from previous iteration.

zero-Dimensional Optimization

1D | Parameter optimization
  1. Engineer has to define parameters such as diameters, bend radii, wall thickness, cone angle...
  2. The levels or values of parameters are changed during the optimization. Many combination of parameters are used to make CFD runs and are often solved using Design of Experiments
  3. Optimized solution will only lie within the defined parameter space

One-Dimensional Optimization

2D | Surface or Shape Optimization
  1. Surfaces (a proposed shape) have to be identified, points on the surface are selected as a design variable and optimizer moves nodes to change geometry
  2. Often morphing techniques are used to modify and move the local surface which may also stretch the attached volume mesh
  3. These problems are also solved either by (a)many CFD runs and DoE or (b)adjoint solvers which may calculate the sensitivity of the flow depending on so-called cost functions like pressure drop
  4. In case of Adjoint solvers, CFD mesh is morphed in a further step based on the sensitivity result. Due to calculation of gradients, adjoint solvers need high numerical accuracy and good convergence.

adjoint Optimization

3D | Topology optimization
  1. This method does not need any parameterization of the geometry - begins without a predesigned shape
  2. Design space (possible flow region) is defined and meshed. Topology optimization is generally formulated as a material distribution problem - during the optimization run, the finite volumes is turned ON or OFF off
  3. Tosca fluid is one of the first commercial software systems for CFD topology optimization.

Excerpts from "Topology Optimisation of Fluids Through the Continuous Adjoint Approach in OpenFOAM" by Prof. Hakan Nilsson - "The Topology Optimization Method (TOM) consists in determine the material distribution in a design domain to maximize or minimize an objective function subject to certain constraints. To maximize/minimize the objective function at flow devices is done by adding the velocity field u multiplied by a scalar field α to the momentum equations, so regions with a high value of α determine a low permeability area (solid portion) and regions with a low value of α determine a high permeability area (fluid portion)." This paper AIAA-2007-3947 describes the implementation of topology optimization in OpenFOAM.

For incompressible steady state Navier-Stokes equations, the problem can be written as:

  minimize cost function J = J(U, p, α).

   such that R ≡ (U . ∇)U + ∇p − ∇.[2μD(U)] + αU = 0, strain rate tensor D(ν) = ∇U + (∇U)T

   ∇.U = 0

Note:
  1. The cost function can be split into contributions from the domain Ω and boundary Γ i.e. J = JΩ + JΓ
  2. The cost function is only depending on values calculated within the domain, any contribution from the flux through the boundaries is 0
  3. To calculate the magnitude of the cost function within any domain, it needs to be integrated over the volume of the domain.
  4. Adjoint velocity [ua] and adjoint pressure [pa] are not velocities or pressure in the physical meaning, but a representation of these quantities for the Lagrangian [L] minimization problem: L := J + ∫Ω(Ua, pa)R
  5. R is set of 3 components of the momentum equation and one continuity equation
  6. Using the "frozen turbulence" assumption the variation of the eddy viscosity μt is neglected

Here, α is the porosity. In OpenFOAM, adjointShapeOptimizationFoam solves both the primal flow (U) and the adjoint flow (Ua) and at the same time optimizes the geometry for minimized cost function say pressure loss. The solver is built around a case of optimization of a duct shape by applying blockage in regions causing pressure loss which are estimated using the adjoint method. The solver is also programmed for shape optimization with respect to pressure loss. The redundant material is shown by high value of α (closer to alphaMax specified in transportProperties dictionary). A detailed description of adjointShapeOptimizationFoam solver can be found here. The pitzDaily case (flow over a back step with pinched outlet) gives following output:

Initial Velocity

adjointShapeOptimizationFoam - Initial Velocity

Optimized Porosity

adjointShapeOptimizationFoam - Porosity Alpha

Optimized Velocity

adjointShapeOptimizationFoam - Optimized Velocity


Excerpts from blogs.sw.siemens.com/simcenter/topology-optimization-cfd-creating-designs-like-nature:- The topology optimization method is based on a level set approach. Unlike the traditional porosity method this approach is not just blocking the flow from the cells with source terms. It is rather defining a moving interface around the predicted flow path. And, this has various advantages! Getting a cleaner interface between solid and fluid which results in designs with less kinks and folds. Eliminating the spare pockets of flow that the porosity method creates. This significantly reduces the clean-up and CAD reproduction time and ultimately reducing the cost of production.

Excerpts from NAFEMS paper "Adjoint-based Topology Optimization - Maximizing Heat Transfer of a Brake Cooling Duct" - Additive manufacturing has made it possible to push the limits of designs beyond the skills and imagination of traditional design engineers. Topology optimization closes the gap between the possibilities and the automation in obtaining those organic optimized shapes. A key class of engineering problems, fluid and thermal applications, are benefiting from the advances on those two design and manufacturing technologies.


Shape Optimzation using Mesh Deformation: The shape optimization process uses the morphing procedure to deform the geometry of interest. Specified control points are displaced by amounts calculated using mesh sensitivities that the adjoint solver provides. The displaced positions of the specified control points are set to maximize/minimize the cost function of interest. At each boundary of computational domain, the conditions of how should be morphed needs to be specified. E.g. the simulation of an airfoil within a wind tunnel, only the airfoil is subject to adjoint shape optimization with the vertices at the boundary being morphed. The surrounding wind tunnel boundaries remain fixed.

Advantage of Adjoint Methods

The adjoint method is an efficient means to predict the influence of many design parameters (inputs) on sensitivity of the objective (output). The advantages of the adjoint method are:
  1. The computational cost for obtaining the sensitivities of an objective does not increase with an increasing number of design variables
  2. The computational cost is essentially independent of the number of design variables because the adjoint method requires only a single flow solution and a single adjoint solution for any number of design variables.
Examples of the types of problems to which the adjoint method is applicable are:
  • What effect does the shape of a duct (input variable) have on the pressure drop (objective)?
  • What is the influence of inlet conditions or location of bend (input variable) on flow uniformity at the outlet (objective)?
  • What areas of the airfoil surface (input variable) have the biggest impact on lift and drag (objectives)?

SIMULIA Tosca Fluid: Shape Optimization for Fluid Flow

  1. TOSCA is a powerful topology optimization tool for producing design concepts that often can't be manufactured using traditional techniques
  2. Such designs are driven primarily by the functional requirements of the part and not the manfacturability
  3. Additive Manufacturing such as 3D printing makes topology optimization results accessible for production
Two options are available for optimization with SIMULIA Tosca
Isight: Parametric OptimizationTosca: Non-parametric Optimization
Size, shape and location of the cutout is unknown and it is constrained by parametersSize, shape and location of the cutout is unknown and thus is non-parametric optimization

optimization Simulia Tosca

Excerpts from "CFD Topology and Shape Optimization of Ford Applications using Tosca Fluid" by Dr. Anselm Hopf, Ford Motor Company, Aachen, Germany: "The roots of Tosca Fluid have been developed at Daimler research in the early 2000s under the name AutoDuct (Klimetzek, 2006), (Moos, 2004)". One of its inventors, Dr. Klimetzek, explains the methodology of CFD topology optimization by looking to the nature of a river delta. Zones of backflow and recirculation have been eliminated by accumulation of sand, called sedimentation. Tosca Fluid is using the same bionic idea and thus called bionic optimization method. CFD topology optimization has been also implemented using mathematical adjoint methods (Othmer, 2006). Tosca Fluid has been implemented for STAR-CCM+ and ANSYS Fluent."


Examples of Shape and Topology Optimizations
shape Optimization Wings
Ref: J. Reuther and A. Jameson. Aerodynamic shape optimization of wing and wingbody configurations using control theory
shape Topology Optimization WingFrame
Ref: L. JiaQi, X. JunTao, and L. Feng. Aerodynamic design optimization by using a continuous adjoint method

Morphing: In morphing technique, a new design variant is created as a complex combination of two or more objects which are geometrically different but topologically similar. Mesh morphing can be used for various purposes such as move onto a new known shape by just updating nodal positions, move onto a shape predicted by the CAE solution (shape optimization using adjoint solver) or move CFD wall boundry according to erosion/deposition.


Fire Modeling

CFD Fire Modeling

Fire modeling is partly modeling of a combustion phenomena. However, in contrast to a well-designed combustor, fire is random in nature. Still CFD simulations can be used to an accurate prediction for the spread of smoke and other toxic substances in the course of a fire.

Zone Model

Mathematical models for the simulation of fires depends on the use of zone models. These divide the compartment into two volumes in which the average properties are calculated as a function of time. With a zone model the fire room is divided into two distinct regions of gases, a hot layer towards the ceiling (which results from the hot buoyant gases arising from the fire forming a plume) and a cooler layer near the floor. When the plume reaches the ceiling it propagates along the ceiling, then slowly descends filling up the space. The height of this layer is dependent on the presence of openings into the compartment. The two regions of the zone model are treated as internally homogeneous control volumes for which the average properties are calculated by applying fundamental laws and correlations, as a function of time. For most situations that a fire engineer encounters, a zone model would provide an acceptable description of what would happen in the event of a fire.

Fire modeling using CFD methods (where the compartment or fire-room is split into many small cells or control volumes, within which conservation of mass, energy and momentum are enforced by applying the Navier-Stokes equations) is also called field modeling. One of the well known and thoroughly investigated phenomena in fire modeling is the "trench effect". It describes now a well established mechanism for fire propagation on enclosed slopes such as escalators or stairwells in building structures.

As observed during Kings Cross (London) Underground fire disaster in the 1987, studies into structural fires have revealed that the trench effect can produce extreme fire behaviour and rapid fire spread. Similarly, a wildfire will spread more rapidly up a slope than it will on flat terrain. The slope of the terrain brings the ground and hence the fuel upon it into closer proximity with the flames downstream the fire. This in turn extends the preheating range and allows for faster rates of spread.

Reference: Modelling Smoke Flow Using Computational Fluid Dynamics by TN Kardos, Fire Engineering Research Report 96/4, December 1996: "The phenomenon arises from a sudden entry of fresh air via a gravity current, into a closed space in which a fire has been burning. This space has accumulated an excess amount of fuel vapour, and with an ignition source results in a deflagration. the fuel can not be burned and will start to accumulate in the gaseous layer. This unburned fuel is known as excess pyrolyzates. If a new ventilation source is introduced, such as a fire fighter opening a door or a window breaks due to thermal stress, the hot gases within the room flow out of the top of the vent and the cooler ambient air would flow in at the bottom.

The flow of cold ambient air into the room is known as the gravity current. The gravity current travels towards the ignition source through the bottom of the vent while hot gases containing excess pyrolyzates exit at the top. In the interface between the gravity current and the layer of hot gases from the leading edge or nose of the gravity current, and the region immediately behind the nose, mixing of the ambient air and excess pyrolyzates occurs.

As a result of mixing oxygen and the unburnt fuel, this region will tend to be within the flammable limits. When the gravity current reaches an ignition source the flammable mixture ignites and travels out along the interface of the gravity current and the layer of hot gases. The flow of gases behind the flame front is sufficiently turbulent to increase mixing in this area, this results in combustion in which the flame front travels below the speed of sound, termed turbulent deflagration, within the compartment. This combination of events helps to drive the excess pyrolyzates out of the compartment in an external fireball."

gravity Current Fire


CGNS

CGNS = CFD General Notation System which provides a general, portable and extensible standard for the storage and retrieval of computational fluid dynamics (CFD) analysis data. In other words, it is a standard for recording and recovering computer data associated with the numerical solution of the equations of fluid dynamics. CGNS facilitates the exchange of CFD data between sites, applications codes (software vendors) and across computing platforms and to stabilize the archiving of CFD data.

CGNS database comprise of:
  • grid coordinates and elements
  • flow solution data, including nondimensional parameters
  • multizone interface connectivity, including abutting and overset
  • boundary conditions
  • flow equation descriptions
  • time-dependent flow
  • reference states
  • dimensional units and nondimensionalization information associated with data
  • convergence history
  • association to geometry definition
  • topologically based hierarchical structures

CFD Simulations for FE Engineers

Most of the numerical simulations such as Finite Element Analysis (FEA) and Fluid Mechanics simulations such as CFD and CAA have lot of things in common. It is easy to participate in cross-application simulations activities if not to become a practicing analyst in both these disciplines. Following sections try to summarize the similarities and differences in FEA vs. CFD) approaches and activities.

Can you describe what are the differences between two representations of the same 3D space?

CAD vs. FE representation

Differences between FEA and CFD simulations

S. No.FEACFD
01We model what we see! Let's call this as concept of modeling the "positive volume". We usually model what we do NOT see. Let's call this approach as modeling the "negative volume". Note that there is nothing called a 'negative' volume.
02Geometry Examples: CHT Domain:

CHT in a pipe


Q: what would be the simulation domain to model flow over this motor-bike with a rider?

Flow over a motor-bike

FEA Domain:

FEA of a hollow pipe

CFD Domain:

CFD pipe flow

03There is no concept of 'walls', 'inlet' and 'outlet' boundry conditions The 'walls', 'inlet' and 'outlet' are the most basic ingredients of CFD simulations, there is no fluid without a 'wall'
04FEA works on force and moment balance, calculation of displacement field which translates into strain and stress CFD works on mass and energy balance, calculation of velocity and pressure fields which translate into pressure drops, viscous/shear forces and pressure forces on the walls
05FEA cannot use hanging elements, it implies crack CFD uses and even recommended sometimes to use hanging elements
06FEA requires fine mesh near high stress gradient such as bends CFD requires fine mesh near high velocity gradient such as walls
07Mesh Examples: the activity involves creating a closed volume - both for FEA and CFD mesh. Mesh of the FEA Domain:

FEA Mesh of a hollow pipe

CAD geometry works on concept of edge, surface and volume. FEA / CFD pre-processor works on concept of topology and material points. Mesh of the CFD Domain:

CFD Mesh pipe flow

08FEA puts limit on mesh aspect ratio < 10 in most of the cases CFD accepts aspect ratio up to 50 in cells or elements near the walls (called boundary layer or prism layer)

Example of computational domain creation: the building planner is interested in knowing the flow pattern around the high-rise buildings. The information from simulation can be used for many purposes such as locating chimneys, gap between two adjacent towers to avoid humming noise...

Input: the 3D geometry of buildings

Flow around the building geometry

Computational domain - the space around the buildings

Flow around the building - computational domain

Computational domain - the solid model of buildings subtracted from the cuboid

Flow around building negative volume

Similarities between FEA and CFD simulations

S. No. Description Remark
01 Both applications uses all types of elements: tri, quad, tetra, wedge / prism, hexahedrons There is no concept of higher order elements in CFD
02 Same pre-processor can be used to generate mesh for FEA and CFD. Boundary layer mesh generation is nothing but a finer mesh near the walls. Why wedge (prism) and hexahendron elements are used in boundary layer and not tetrahedrons? The answer lies in definition of aspect ratio and tetra-collapse?
03 Both conformal and non-conformal contacts can be used in FEA and CFD. In CFD, it is known as 'interfaces'. Since the number of interfaces can be very high in CFD, it is normally recommended to keep less number of interfaces or avoid altogether.
04 Periodicity and symmetry / anti-symmetry FEA uses symmetry but does not use concept of periodicity. CFD simulations do not use concept of anti-symmetry.

Symmetry Anti-symmetry B.C.

05 Most FEA and CFD programs can solve thermal conduction phenomena without modeling the fluid flow if temperature and other heat transfer boundary conditions are known. In FEA, it is used to find thermal stress. In CFD, it is used to make a quick calculation of hot spot temperatures.

Tetra-Collapse: The height of the Tetra element is measured from each of the four nodes to its opposite face and then divided by the square root of the area of opposite face. The minimum of four resulting values is then normalized by dividing it by 1.24. As the size of tetrahedron collapses, the value approaches to 0. While a perfect tetrahedron has a value of 1.0, tetra-collapse values < 0.2 are generally not acceptable in most of the simulations.

Stretched and sliver tetrahedrons

The exchange of information between FEA and CFD engineers can be in terms of topology checks and mesh quality parameters required for FLUENT Mesher and STAR-CCM+

S. No.Description Checkpoint Actual value in the mesh
01Mesh Type Surface mesh, all tri-elements-
02Minimum size 2.0 [mm]-
03Maximum size 16 [mm]-
04Total 2D cell count2,00,000-
05Topology: number of free edges 0-
06Topology: number of self-intersections0-
07Topology: number of penetrating elements0-
08Topology: number of duplicate elements0-
09Topology: number of over-lapping elements0-
10Quality: skewness [0 is ideal, 1 is degenerate]≤0.60-
11Quality: aspect ratio ≤5.0-
12Quality: surface deviation*≤0.1 [mm]-
13PID / Zone naming as per naming convention agreed upon-
* for external flows where drag prediction is important, the suface deviation should be ≤ 0.005 [mm]
Oursourcing of CFD Simulations: Why and When?

The engineering services has grown to the extent that a new category of companies ESP (Engineer Service Providers) has come into vogue. These companies supplement the resources (work-forces, hardwares, software) to brand owners known as OEM (Original Equipment Manufacturers).

Why companies outsource the simulation tasks?
  1. Peak load scenario / short term requirement: The in-house capacity is typically sized for long-term average condition. Field failure investigations, data conversion to digital format (e.g. 2D to 3D) fall under this category.
  2. Cost and flexibility: Companies not only get (employment) cost savings, they do not need to resort to lay-off in weak market conditions. Substitution of capital for labor is the classic method of obtaining operating leverage in both product- and service-oriented businesses. When the tasks cannot be automated because human judgment must be exercised, cheap labor can often be substituted for expensive labor as a means of obtaining operating leverage.
  3. Uncertainty factors: Long term forecast of the utility of additional human resource and infrastructure (hardware, software, test labs) is not easy and reliable. A transition from Time and Material to Fixed Price mode falls under this category.
  4. "Make or Buy" decisions: If you notice how the OEMs work, they make and/or assemble strategically important parts in-house (say engines) and buy specialized components say valve seats. Similarly, the decision to outsource any design or simulation activity is strongly influenced by the dependency on other process and services.
  5. Franchising Model: Sometimes bigger firms don't want to do everything on their own, they consider going the franchising route. Many enterprises consider ESP an extension of their engineering organization.
  6. Track Records : Most of the outsourcing decisions are based on the track records of the service providers. Hence, it is a vicious cycle.
  7. Professional network and marketing outreach: In the engineering industry, few firms dominate each niche, making it extremely difficult for small startups to break into specialized fields and capture their slice of the total market. They have big pool of executives who keep in touch with anybody their firm has connected with in past. They develop relationships with new clients while ensuring that queries and complaints their current clients have get answered satisfactorily.
  8. Generalist vs. Specialists: While there is no rule, many decisions are based on complexity (and associated risk) of the work.

Few myths on benefits of outsourcing:
  1. Benefits of different time-zones: Any particular project or task does not get expedited just because there are engineers working in different time zones. None of the projects are executed in a mode where the tasks are handed over from one time zone to other. Most of the time, the tasks are in totality executed by engineers of same time zone.
  2. Higher value at lower cost: Except few exceptions (such as brand premium) cost and quality are closely linked. An Apple iPhone costs far higher than mobile phones from other brands having similar functionality. Though it is difficult to quantify if an iPhone owner gets higher value for money as far as utility of a "smart phone" is concerned.
  3. Engineering service providers can compete against larger engineering firms and OEMs: No OEM or product engineering firm outsource more that the work equivalent to 5% of their total workforce. Hence, claim that service providers make a large impact on top-line and bottom-line of an OEM is a bit of exaggeration!

Evaluation and Comparison of CFD Simulation Tools

The pricing structure of modern CFD tools are getting modular. There are three distinct category of licenses: Pre-Post, Solver and Parallel operation (known as HPC: High Performance Computing). A type of license which packs all 3 features in one is sometimes call a 'Bundled' license. They are designated differently such as Perpetual, Leased or Term License, Student Version, Academic, Professional, Power, Premium, Network License, Batch License, Floating License (WAN-type, LAN-type, Global), On-Demand, Enterprise, Node-locked or CPU-Locked, User-locked or Named Single-user License, Credits-based... Most of commerical tools have mix of licenses: For example, STAR-CCM+ Lite license does not allow to set-up a case with thermal radiation. ANSA uses credit based license system where with 100 credits one user per computer can open simultaneously as many ANSA sessions as they want. ANSYS FLUENT requires HPC Packs to make simulations in parallel with number of cores more than 4. However, one thing which is often missed is the tools required to clean-up and defeature the simulation geometry. The comparison of two CFD codes needs to be done on mainly 5 factors:

Geometry Clean-upPre-Processing (Meshing)Pre-Processing (Solve Set-up)Separate HPC? Post-Processing
ANSAANSAFLUENTY FLUENT
ANSAFLUENT FLUENTYTecplot
STAR-CCM+STAR-CCM+STAR-CCM+NSTAR-CCM+
SpaceClaimANSAFLUENTYFLUENT
SpaceClaimHyperMeshSTAR-CCM+NSTAR-CCM+

Similar way, the features of software also need to be compared such as:

  • Which all neutral CAD formats (STEP, Parasolid, IGES, STL) be imported?
  • What native CAD formats (Creo, CATIA, Inventor...) be imported?
  • Can the pre-processor read tesselated geometry (STL) or the surface mesh generated in other programs?
  • Can the solver read volume mesh generated in other programs?
  • What is the licensing split between pre-processing and full High Performance Computing (HPC) for solver?
  • Does the program have capability to export data to post-process in open source post-processors such as ParaView?
  • Which programming language is used to record the session and replay, write automation scripts and user functions?
  • Can the software interpolate data from previous results, make re-start simulations and exchange data with other solvers?

Note that as on V2022, ANSYS ICEPAK needs the mesh to be generated inside ICEPAK environment. STAR-CCM+ has a Power license option which provides unlimited parallel operations. ANSYS FLUENT provides modular HPC Packs where you add a module as you ramp-up the parallel operations. All these options point to the fact that the cost-effectiveness of any particular set of tools as tabulated above shall depend on the number of users. For small number of users, modular licensing may not help. On the other hand, as the number of users increases significantly, modular structure of licenses may turn out to be more effective. Other factor which needs to be considered is the available HPC hardware. If the solver runs are limited to fewer number of cores, a modular licensing may be beneficial.


Problem-solving and issue / obstacle / hindrance-resolution

A precise definition of any "technical issue" is 50% of its resolution. An important aspect to remember is ability to distinguish 'cause' and 'effect' where one tends to swap the two.

Effect: Water started boiling, Cause: Water was heated and the molecules started moving quickly.

Effect: The ocean have tides, Cause: The moon creates gravitational pull.

Effect: Sales is in downtrend, Cause: There is no new order and/or delay in ongoing projects.

The cause-and-effect discussion can be organized in one of these two primary ways:

  • Start with the cause and then talk about the effect
  • Start with the effect and then talk about the cause

CFD Simulations in Various Phases of Product Development Cycle

Preliminary Design Review (PDR) is performed during the Design and Engineering Phase when the System Developers create a baseline high-level system architecture with a reasonable expectation of satisfying the requirements. This phase of development does not require any simulations activities and are mostly based on past experiences and domain knowledge. A PDR is conducted before the start of detailed design work.

Preliminary Design Review (PDR) is conducted to ensure new technologies are mature enough to be integrated into a product subsystem to form its allocated baseline. A Critical Design Review (CDR) is focused on determining if a system can meet its stated performance requirements within cost, schedule, and risk. CDR is the phase when numerical simulations may need to be performed to get some baseline insight into the hydraulic and thermal performances.

Contact us
Disclaimers and Policies

Step-by-step approach to CFD simulations using commercial tool ANSYS CFX / FLUENT or open-source product like OpenFOAM, textbook solutions to Fluid Mechanics problem, quiz on CFD and FEA are covered on this site. Various physics like simulations of multiphase flows, rotating machines are covered with examples. The content on CFDyna.com is being constantly refined and improvised with on-the-job experience, testing and training. Examples might be simplified to improve insight into the physics and basic understanding. Linked pages, articles, references, and examples are constantly reviewed to reduce errors, but we cannot warrant full correctness of all content.