• CFD, Fluid Flow, FEA, Heat/Mass Transfer
  • +91 987-11-19-383
  • amod.cfd@gmail.com

Solver Setting

Type of Solvers and Solution Control Parameters

This section deals with solution controls for solvers including topics like CFL Number, Time-step for Transient Simulations, Psuedo-time Marching, Parallel Computing, Nodes and Cluster, HPC - High Performance Computing, Threading, Partitioning, MPI - Message Passing Interface and Scalability.

Any program does only one thing: "Solves exactly what one tells it to solve and how to solve, nothing else!"

Solver setting

This activity of a CFD simulatin process involves selection of solver type (pressure-based or density-based, stead-state or transient), selection of [turbulence model, spatial discretization scheme (upwind, central difference), temporal discretization scheme (forward Euler, Backward Euler, Crank-Nicolson scheme, explicit leapfrog method), relaxation factor, PV-coupling method, transient time-step for solver (outer loop iteration), number of inner loop iterations for transient simulations, intermediate output data back-up frequency], gravity ON/OFF, creation of monitor points and monitor planes... The purpose of all these settings is to get a numerical scheme which is consistent, stable and convergent. Here stability refers to the property when errors caused by small perturbations in numerical method remains bound.

Note the the field variables solved during CFD simulation are velocity vector, pressure, temperature (enthalpy). The density at each cell is calculated from equation of state (ideal gas law or real gas laws like Redlich-Kwong and IAPWS equation of state for steam). Note that the "compressible flows" and "compressible gas" are not the same concepts. Compressible flows are characterized by Mach number of the flow whereas compressible behaviour of fluid refers to equation of state such as ideal gas law, essential change in density due to change in pressure and temperature.


Compressibility of Fluids

One of the key consideration is the compressible or incompressible behaviour of the fluid. For liquids under lower pressure conditions (except cases such as hydraulics, deep bore drilling in oil exploration...), they can be assumed incompressible without loss of any accuracy. The compressibility chek for liquids need to be decided based on bulk modulus, operating pressures and coefficient of volume expansion.. The later varies between 0.001 ~ 0.010 [K-1]. The bulk modulus is defined as β = -V × ΔP/ΔV where negative sign refers to the fact that volume decreases as pressure increases. Thus, in terms of density, β = -ρ/Δρ × ΔP. β for water is 2.2 [GPa]. Hence, at pressure of 100 [MPa], the change in water density as compared to density at 1 [atm] would be 1000 [kg/m3] × 0.1 [GPa] / 2.2 [GPa] = 45.5 [kg/m3] which is 4.55% of density of water at 1 [atm].

Bulk Modulus Formula

While the change in density of liquids is more pronouced with change in temperature, the situation is different in cases of gases where density is a strong function of both temperature and pressure changes. Note that the density is calculated based on static conditions (static temperature and static pressure) and not on total conditions (total temperature and total pressure). On the other hand, the stagnation conditions are the true measure of pressure energy and internal energy. There are two modes of reduction in stagnation (or total) pressure: the first mode is due to conversion of pressure energy into velocity (kinetic) energy and the second mode is loss of pressure energy due to wall friction and minor losses (sudden expansion and contraction). The general equations to estimate stagnation and static conditions for all values of flow velocities are described below.

Density as a function of Mach number

Let's consider isothermal flow through a rough pipe of uniform cross-sectional area. The expression for frictional pressure drop is also described. Here, V denotes area-average velocity at the inlet and the outlet.

Isothermal frictional pressure drop in a Pipe

From conservation of mass:

Mass conservation for flow in a Pipe

Static pressure drop taking into account compressibiity of the fluid:

Isothermal  pressure drop in a Pipe including compressibility

Since ρ1 always > ρ2 for isothermal flows, the effect of compressibility of the gas is to increase static pressure drop. Since wall friction is also function of velocity gradient near wall and the free-stream velocity increases as density of gas decrease due to frictional loss, the effect of compressibility is to increase wall friction and subsequently viscous losses. But how much will the compressibility affect static pressure drop? For ρ1 / ρ2 = 1.11 (10% drop in density), V1 = 5 [m/s] and ρ1 = 1.20 [kg/m3], additional pressure drop due to compressibility = 0.5 × 1.20 × 25 × (1.11 - 1.00) = 0.66 [Pa]. However, a 10% reduction in density will call a 10% drop in absolute static pressure. At mean value of 1 [bar], this turns out to be 100,000 x 0.1 = 10,000 [Pa]. This magnitude of pressure drop can be observed in a pipe of diameter 100 [mm] and length of approximately 2.2 [km]!

The operating pressure is used to specify the mean pressure in the system. Though it can be adjusted to convert the gauge pressure into absolute pressure. That is, the operating pressure can be set to zero to get the pressure field directly in absolute values. This is the recommended method for cases involving phase change such as simulation of cavitation in pumps. Note that cavitation occurs when local absolute static pressure in the flow field falls below saturation pressure at given temperature. However, once the operating pressure is set to 0 [Pa], the boundary conditions should also be specified in absolute values, e.g. 1,01,325 [Pa] instead of 0 [Pa].

Operating Pressure in CFD Simulations

The solver-setting process encompasses following aspects of numerical solutions.
  • Discretization scheme for momentum, pressure, energy and turbulence parameters: The PRESTO (Pressure Staggering Option) is recommended for pressure because in flows with high swirl numbers.
  • PV-Coupling such as SIMPLE, SIMPLER, PISO
  • Gradient schemes: Green-Gauss node based or Green-Gauss cell-based.

    Discretization Schemes

    On irregular (skewed and distorted) unstructured meshes, the accuracy of the least-squares gradient method is comparable to that of the node-based gradient (and both are much more superior compared to the cell-based gradient).
  • Conservation target and residual levels for convergence criteria. In addition to the "residual norm", it is strongly recommended to set "mass-conservation" or "energy-conservation" as convergence check.
  • Fluid time-scales
  • Solid time-scales (in case of conjugate heat transfer or pure conduction problems). The recommended practice to have solid time-scale set to an order of magnitude (10 times) higher than fluid time-scale.
  • One should ensure that the solution is mesh-independent and use mesh adaption to modify the mesh or create additional meshes for the mesh-independence study. A mesh-sensitive result confirms presence of "False Diffusion".
  • The node-based averaging scheme is known to be more accurate than the default cell-based scheme for unstructured meshes, most notably for triangular and tetrahedral meshes.
  • Note that for coupled solvers which solves psuedo-transient equations even for steady state problems (such as CFX), relaxation factors are not required to be set. The solver applies a false timestep as the convergence process is iterated towards final solution.
  • QUICK or Geo-Reconstruct scheme for the volume fraction in FLUENT

SIMPLE Algorithm: Analogy between analytical solution and numerical simulation

SIMPLE - Analogy

Differences between collocated (non-staggered) and staggered grid layout

  • For collocated solvers (such as ANSYS CFX), control volumes are identical for all transport equations, continuity as well as momentum.
  • CFX is a node-based solver and constructs control volumes around the nodes from element sectors, thus the number of control volumes is equal to the number of nodes and not the number of elements. The input of CFX may be tetrahedron, prism and hexhedron in physical form, the solver internally generates a polyhedral mesh. In a cell centered code, such as FLUENT or STAR-CD / STAR-CCM+, number of elements are the same as the number of control volumes (as the control volumes are same as physical elements) so these are often used interchangeably.
  • All field (or unknown or solution) variables as well as fluid properties are stored at the nodes (the vertices on the mesh).
  • The distinction between collocated and staggered approach is nicely explained in section 4.6 by S. V. Patankar in his book "Numerical Heat Transfer and Fluid Flow" and section 7.2 in the book Computational Methods for Fluid Dynamics by J. H. Ferziger and M. Peric.
  • For staggered solvers like ANSYS FLUENT, values of scalar field variables as well as material properties are stored as cell centres.
  • Due to the difference in the way field variables are stored, the simulation with same mesh, material properties and boundary conditions, the y+ value reported in ANSYS CFX will be approximately twice that reported in FLUENT.

Excerpts from user manuals of commercial tools

  • Smaller physical timesteps are more robust than larger ones.
  • An 'isothermal' or cold-flow simulation is more robust than modeling heat transfer. The Thermal Energy model is more robust than the Total Energy Model.
  • Velocity or mass specified boundary conditions are more robust than pressure specified boundary conditions. Static pressure boundary is more robust than a total pressure boundary.
  • If the characteristic time scale is not simply the advection time of the problem, there may be transient effects holding up convergence. Heat transfer or combustion processes may take a long time to convect through the domain or settle out. There may also be vortices caused by the initial guess, which take longer to move through the entire solution domain. In these cases, a larger timestep may be needed to push things through initially, followed by a smaller timestep to ensure convergence on the small time scale physics. If the large timestep results in solver instability, then a small time scale should be used and more iterations may be required.
  • Sometimes the levels of turbulence in the domain can affect convergence. If the level of turbulence is non-physically too low, then the flow might be "too thin" and transient flow effects may be dominating. Conversely if the level of turbulence is non-physically too high then the flow might be "too thick" and cause unrealistic pressure changes in the domain. It is wise to look at the Eddy Viscosity and compare it to the dynamic (molecular) viscosity. Typically the Eddy Viscosity is of the order of 1000 times the dynamic viscosity, for a fully turbulent flow.
  • The 2nd Order High Resolution advection scheme has the desirable property of giving 2nd order accurate gradient resolution while keeping solution variables physically bounded. However, may cause convergence problems for some cases due to the nonlinearity of the Beta value. If you are running High Res and are having convergence difficulty, try reducing your timestep. If you still have problems converging, try switching to a Specified Blend Factor of 0.75 and gradually increasing the Blend Factor to as close to 1.0 as possible.

Solver Setting: Physics

The flow regime and classification of laminar, transition and turbulent zones are unique to every flow. One of the most widely used flow condition is flow over a cylinder in cross flow. The following summary is an interesting description of how the flow evolve:
  • Re ≤ 5, flow is laminar throughout and no vortices are present downstream of the body.
  • 5 < Re < 40, a vortex pair can be observed in the wake region though flow is still laminar in the boundary layer.
  • 40 ≤ Re ≤ 200, a von Karman vortex street is formed as where pairs of counter-rotating vortices are shed downstream of the body. The vortex shedding frequency is given by a non-dimensional Strouhal number: St = f × D / U.
  • 200 ≤ Re ≤ 300, the flow behind a cylinder undergoes transition to a turbulent wake.
  • 300 ≤ Re ≤ 3×105, the sub-critical Reynolds number regime wake is turbulent while the boundary layer is still laminar.
  • Re > 3×105, a transition from laminar to turbulent boundary layer separation occurs.
In any practical application, the flow phenomena may not be laminar throughout or turbulent throughout. Hence, solver has to deal with different length scales and the result should be validated carefully.

Sover Setting in STAR-CCM+

Solver Setting for Transient Simulation: CFL Number

One of the concern is to chose a right time-step value in transient simulation to have the best compromize between computational time and accuracy. Can a very very small time step be unconditionally stable? CFL (Courant–Friedrichs–Lewy) number is mathematically represented as Co = a.Δt/Δx. Here, a is the characteristic wave speed of the system. Why wave speed, there is no wave phenomena in the flow? It may not be apparent though the velocity field varies both in space and time domain - which is typical characteristics of wave propagation phenomena. The wave characterics is the reason Courant number is a measure of how much information 'u' traverses a cell or element of computational mesh (Δx) in a given time-step (Δt). For most use cases, 'a' can be assumed equal to the velocity scale while estimating Δt for a given Δx.

In an explicit Eulerian method, the Courant number cannot be > 1 otherwise it will lead a condition where information is propagating through more than one grid cell at each time step. This is explained in the physical domain of dependence and numerical domain of dependence in following diagram. For a stable solution, numerical domain of dependence must cover the physical (real) domain of dependence.

TIme discretization scheme

In multi-dimensional problem, depending upon the spatial and temporal discretization scheme, the CFL condition varies as shown below.

CFL criteria for multi-dimensional problems

Thus, Courant number can be reduced either by reducing the time step or by increasing the mesh size (i.e. coarsening of the mesh). The later option may appear counter-intuitive though mathematically that remains a fact. It is difficult for most of the simulation engineers to believe the fact that coarsening the mesh help them achieve a better convergence.

  • The CFL number scales the time-step sizes that are used for the time-marching scheme of the flow solver.
  • A higher value leads to faster convergence but can lead to divergence and unstable simulations. The (stiff) CFL condition can be relaxed or softened by using implicit Euler such as Crank–Nicolson which is unconditionally stable.
  • The inverse of this is also true where choosing a smaller value in an unstable simulation improves convergence.
  • ANSYS FLUENT uses implicit (first order and second order) formulation for pressure-based solver. As per the user guide: "It is best to select the Coupled pressure-velocity coupling scheme if you are using large time steps to solve your transient flow, or if you have a poor quality mesh."
  • In general, the solution at the end of (inner iterations in) current time-step is used as initial value for iterations at next time-step. ANSYS FLUENT provides option to extrapolate the solution variable values for the next time step using a Taylor series expansion. This can improve the convergence of the transient calculations by enabling and the absolute residual levels are lowered during inner iterations.
So how should one chose the best time step, Δt? There is not unique option and need few trials at the early stage of transient phenomena. A recommended method is to calculate starting value of Δt = Δx/V/10 where Δx is the minimum grid spacing in computational grid. Observe the number of inner iterations chosen solve needs to converge at each time step. The ideal number of iterations per time step is 5 to 10. If the solver needs substantially more number of time steps, the Δt is too large to capture the transient effect. If the solver takes only a few iterations per time step, Δt is unnecessarily small and can be increased without loss of accuracy or likelihood of convergence issues.

To verify that your choice for time-step was proper, check contours of the Courant number within the domain. In FLUENT, select Cell Courant Number under Velocity... For a stable and accurate calculation, the Courant number should not exceed a value of 20-40 in most sensitive transient regions of the domain.

Note that CFX reports two types of Courant number:

Courant Number inn CFX

Material Properties

Sometime underlying physics and type of solver chosen governs the material definition. For example, a density based solver requires density to be function of pressure and temperature. Similarly, simulation of Joule-Thomson phenomena (cooling effect due to throttling) requires material properties to be modeled as real-gas. An ideal gas will never cause this effect. The simulation related to steams required careful understanding of state of steam: wet-steam vs. dry steam. An incompressible flow may required only density and viscosity of a fluid whereas phase change and/or reacting mixture phenomena will require definition of reference enthalpy and entropy.

Material Properties - Reacting and Multiphase flows

Wet Steam: The steam with droplets of water is known as wet steam and it characterized by dryness fraction. Most of the boilers do not generate dry steam as steam bubbles breaking through the heated surface will pull tiny water droplets along with it. In general, when a substance exists partly as liquid and partly as vapor (gas) phase, its quality is defined as the ratio of the mass of vapor to the total mass. On the hand, mass fraction of the condensed phase is known as wetness factor.

Dry Steam: Saturated or superheated steam are known as dry steam.

For water, the triple point T = 0.01 [°C] and P = 0.6113 [kPa], is selected as the reference state where the internal energy, enthalpy and entropy of saturated liquid are assigned a zero value. Similarly, saturated refrigerant R-134a is specified enthalpy of 0 [J/kg] at temperature -40 [°C] and pressure of 51.25 [kPa]. Thus, standard state enthalpy of liquid water is 0 and reference temperature is 273.15 [K]. The standard state enthalpy of steam is heat of vaporization at T = 0.01 [°C] and P = 0.6113 [kPa] and equals 2501 [kJ/kg]. The saturation pressure of water as per ASHRAE: p(T in [K]) in [Pa] = exp[c1 / T + c2 + c3 * T + c4 * T2 + c5 * T3 + c6 * Ln(T)] where c1 = -5800.22, c2 = +1.3914993E+00, c3 = -4.8640239E−02, c4 = +4.1764768E−05, c5 = −1.4452093E−08, c6 = +6.54597.

Porous Domains

Flow geometry such as heat exchangers with closely spaced fins, honeycomb flow passages in a catalytic converters, screens or perforated plates used as protection cover at the from of a tractor engine ... are too complicated to model as it is. They are simplified with equivalent performance characteristic, knowns as Δp-Q curve. These curves are either generated using empirical correlations from textbooks or using a CFD simulation for smalled, periodic / symmetric flow arrangement. The simplified computational domain is known as "porous zone" in case it is represented as a 3D volume or pressure or porous jump in case it is represented as a plane of zero thickness. In a similar fashion, the performance data of a fan can be specified including the swirl component.

All the porous media formulation take the form: Δp = -L × (A.v + B.v2) where v is the 'superficial' flow velocity and negative sign refers to the fact that pressure decreases along the flow direction. The 'superficial velocity' is calculated assuming there is no blockage of the flow. L is the thickness of the porous domain in the direction of the flow. Here, A and B are coefficients of viscous and inertial resistances.

In FLUENT, the equation used is: Δp/L = -(μ/α.v + C2.ρ/2.v2) where α is known as 'permeability' and μ is the dynamic viscosity of fluid flowing through the porous domain. This is a measure of flow resistance and has unit of [m2]. Other unit of measurement is the darcy [1 darcy = 0.987 μm2], named after the French scientist who discovered the phenomenon.

STAR-CCM+ uses the expression Δp/L = -(Pv.v + Pi.v2) for a porous domain.

The pressure drop is usually specified as Δp = ζ/2·ρ·v2 where ζ is 'equivalent loss coefficient' and is dimensionless. Darcy expressed the pressure gradient in the porous media as v = -[K/μ]·dP/dL where 'K' is the permeability and 'v' is the superficial velocity or the apparent velocity determined by dividing the flow rate by the cross-sectional area across which fluid is flowing.

Steps to find out viscous and inertial resistances:

  1. Calculate the pressure drop vs. flow velocity data [Δp-v] from empirical correlations or wind-tunnel test or simplified CFD simulations.
  2. Divide the pressure drop value with thickness of the porous domain. Let's name it as [Δp'-v curve].
  3. Calculate the quadratic polynomial curve fit coefficients [A, B] from the curve Δp'-v. Ensure that the intercept to the y-axis is zero.
  4. In STAR-CCM+ these coefficients 'A' and 'B' can be directly used as Pv and Pi which are viscous and inertial resistance coefficients respectively.
  5. Divide 'A' by dynamics viscosity of the fluid to get inverse of permeability that is 1/α to be used in ANSYS FLUENT as "viscous resistance coefficient".
  6. Divide 'B' by [0.5ρ] where ρ is the density of fluid, to get C2 to be used in ANSYS FLUENT as "inertial resistance coefficient".
  7. The method needs to be repeated for the other 3 directions. If the flow is primarily one-directional, the resistances in other two directions need to be set to a very high value, typically 3 order of magnitude higher.
The GUI to set porous domain in ANSYS FLUENT is as shown below. For 3D cases, direction vectors for any two principal axes need to be specified, the third direction is automatically calculated by FLUENT. However, one must be consistent in specification of direction vector and resistance coefficients.

Porous Media Setting in ANSYS FLUENT

In case porous domain is not aligned to any coordinate direction, the direction of unit vector along the flow and across the flow directions can be estimated from following Javascript program. Note that empty field is considered as 0.0. There is no check if a text value is specified in the input fields and the calculator will result in an error.

First point - X1:  
First point - Y1:  
First point - Z1:  
Second point - X2:  
Second point - Y2:  
Second point - Z2:  

Temperature Drop across a Porous Domain

If flow through the porous domain is assumed (a) incompressible, one-dimension, steady state and (b) pressure drop varies linearly with the flow direction x:

temperature rise in a porous domain

Initial Condition

The field variables (u, v, w, p, T, k , ε...) need to be specified certain starting value. It can all be set to 'ZERO' in their respective units or some other physically realistic values. In CFD simulation process, it is called 'initialization' which can be domain (volume specific) or global (applicable to entire computational nodes). For steady-state simulation, the initial condition is a 'wise guess' of the converged steady-state solution. In transient simulations, the initial condition provides the (acutal) spatial distribution at t = 0 and hence temporal starting point for the simulation. Thus, for transient simulations, the intermediate result will be strongly influenced the initial condition.

In case of multi-phase flow, an initial distribution of volume fraction of the constituents is important. The method to set initial conditions in ANSYS FLUENT is called patching which requires regions to be marked as shown below. In OpenFOAM setFields (requires setFieldsDict) utility is used to initialize field variables in various zones.

Region Adaptation in ANSYS FLUENT

Convergence History

The residuals of continuity, momentum (u, v and w), energy (T), k and ε is used to check the convergence of iterative solution process. Convergence refers to the situation where the calculated value is assumed to tend to correct inversion of equation [A]{x}={b}.

Sometimes k (designated as Tke in STAR CCM+) and and ε (referred to as Tdr in STAR CCM+) may be as high as few thousand and the results are reasonable. The scale can be reduced by turning off the normalization for k and ε. Normalization takes the very first iteration's values to normalize all subsequent residuals. When solution already start with a good solution or guess value for k and ε, the first value is already very low, and residuals may not reduce further. It then sometimes can happen that the residual increases instead of decreasing even though the absolute (without normalization) value might still be very low.

In ANSYS FLUENT, the denominator to calculate the scaled residual for the continuity equation of a pressure-based solver is the largest absolute value of the continuity residual in the first 5 iterations. When iteration starts with a good solution, the largest residual of first 5 iterations may be low. Hence in subsequent iterations, significant reduction in order of magnitude of the residuals may not be observed even though the convergence can be deemed to be of good accuracy. Note 'scaling' and 'normalization' of residuals are different and interpretation should be different as well.

Under Relaxation Factors (URF): As per ANSYS documentation and training material, Pressure and Momentum URFs MUST sum to 1.00 and flow stability is best with smaller pressure URF (0.03 is not uncommon).

Thin Walls

Thin Wall Modeling in FLUENT

Mesh Reordering - Node Renumbering

Both elements and nodes are numbered where elements are described as a set of nodes forming its vertices. The more compact is the arrangement of elements and nodes, lesser will be the memory requirements. Some terms associated with elements and node arrangements in a mesh are as follows.

  • Bandwidth: It is the maximum difference between neighboring cells in a zone i.e. if each cell in the zone is numbered in increasing order sequentially, bandwidth is the maximum differences between these indices.
  • Excerpt from user manual - FLUENT: Since most of the computational loops are over faces, you would like the two cells in memory cache at the same time to reduce cache and/or disk swapping i.e. you want the cells near each other in memory to reduce the cost of memory access.
  • In general, the faces and cells are reordered so that neighboring cells are near each other in the zone and in memory resulting in a more diagonal matrix that is non-zero elements are closer to diagonal. Refer to "banded-matrices" in context with numerical simulations.

Solver Setting for Multi-Phase Flows

Multi-phase flows have wide applications in process, automotive, power generation and metal industries including phenomena like mixing, particle-laden flows, CSTR - Contunuously Stirred Tanks Reactor, Water Gas Shift Reaction (WGSR), fluidized bed, fuel injection in engines, bubble columns, mixer vessels, Lagrangian Particle Tracking (LPT). Some of the general characteristics and categories of multi-phase flow are described below along with setting parameters in ANSYS FLUENT

 Multiphase flow regimes are typically grouped into five categories: gas-liquid (which are naturally immiscible) flows and (immiscible) liquid-liquid flows, gas-solid flows, liquid-solid flows, three or more phase flows. As can be seen, the immiscibility is a important criteria. In a multi-phase flow, one of the phase is usually continuous and the other phase(s) are dispersed in it. The adjective "Lagrangian" indicates that it relates to the phenomena of tracking a moving points ("fluid particles") - such as tracking a moving vehicle on a road. On the other hand, the adjective "Eulerian" is used to describe correlations between two fixed points in a fixed frame of reference - such as counting the type of vehicles and their speeds while passing through a fixed point on the road.

Gas-liquid flows are further grouped into many categories depending upon the distribution and shape of gas parcels. Three such types are described below.

Bubbly Flow: it represents a flow of discrete gaseous or fluid bubbles in a continuous phase.

Bubbly Flow - Gas Liquid

Slug Flow: This is characterized by flow of large gas bubbles in liquid.

Annular Flow - Gas Liquid

Annular Flow: Here one of the phase if confined to area near the wall forming an annular section.

Annular Flow - Gas Liquid

Some other types of flows are particle-laden flow such as air carrying dust particles, slurry flow where particles are transported in a liquid, hydrotransport which describes densely-distributed solid particles in a continuous liquid such as cement concrete mix. Gas assisted mixing of solid such as fluidized-bed and settling tank where particles tend to sediment near the bottom of the tank forming thick sludge are some other examples of multi-phase flows.

The two dominant method of multi-phase simulations are listed below.

multiphase Methods

In the Euler-Euler approach, the different phases are treated mathematically as interpenetrating continua where each phase coexist as continuous function in space and time.
  • A phase cannot be occupied by the other phases in space, the concept of "phasic volume fraction" is used and assumed to be continuous functions of space and time. Distribution of phases are plotted by volume fraction whereas in DPM (Euler-Lagrange approach) the distribution of phases are plotted by particle track.
  • Globally, the sum of this phasic volume fraction equals to one.
  • Conservation equations for each phase are derived as a set of equations having similar structure for all phases, analogous of momentum equations.
  • These equations are closed by constitutive relations obtained from empirical information.

multi-phase flow: VOF and Eulerian in ANSYS FLUENT

In transient VOF simulations, ANSYS FLUENT solves the volume fraction equation(s) once for each time step (outer loop). To solve the volume fraction equation(s) at every (inner) iteration within a time step, use TUI: /define/models/multiphase/volume-fraction-parameters [solve vof every iteration?] yes. This will update convective flux coefficients in the other transport equations based on the updated volume fractions at each iteration, though this a computationally more expensive and less stable option. While using sliding meshes or dynamic meshes with layering and/or remeshing, using the solve vof every iteration? option will yield more accurate results, although at a greater computational cost.

multi-phase flow: boiling in ANSYS FLUENT

Note: in ANSYS FLUENT, in steady-state discrete phase modeling, particles do not interact with each other and are tracked one at a time in the domain.

particle Injection Setting

Note that if the field "Number Of Continuous Phase Iterations Per DPM Iteration field" set to 0, FLUENT will not perform any discrete phase iterations. Excerpts from ANSYS FLUENT: rule of thumb to follow when setting the parameters above is that if you want the particles to advance through a domain consisting of N mesh cells into the main flow direction, the Step Length Factor × N ≈ Max. Number of Steps.

DPM Trajectory Equation

Newton Second Law for Particle Tracking

Here subscript 'p' denotes 'particle'. Various empirical correlations of drag laws for viscous forces over solid particles are:

Drag Laws Solid Particles

ANSYS FLUENT calculates drag coefficient of non-spherical particle using shape factor defined as [s/S] where 's' is the surface area of a hypothetical sphere having the same volume as the particle and 'S' is the actual surface area of the particle.

As explained above, the underlying physics of DPM is described by ordinary differential equations (ODE) as opposed to the Navier-Stokes equation for continuous flow which is expressed in the form of partial differential equations (PDE). Therefore, DPM can have its own numerical mechanisms and discretization schemes, which can be completely different from other numerics (such as discretization schemes) to solve the governing equations. The Numerics tab gives users control over the numerical schemes for particle tracking as well as solutions of heat and mass equations related to particles - continuous phase interactions.

DPM Numerics

Type of coupling: continuous (carrier) phase and discrete phase

One-way coupling between continuous and discrete phase

The approach differs in ways injected particles (number, volume fraction, injection speed) are treated or the way particle-particle collisions are treated. For example: one particle can represent either one or a group of particles (sometimes known as parcels). In Discrete Particle Method (DPM), particle particle interaction can be treated in following 3-ways:
  • no collisions at all (used for dilute suspensions): one-way or two-way coupling, Newton's equation of motion for dispersed phase.
  • 3-way coupling - soft-sphere approach (the discrete element method, DEM): the instantaneous inter-particle contact forces due to oblique collisions (the normal, damping and sliding forces) are computed using equivalent mechanical elements springs, dashpots and sliders respectively. This approach requires a very small time step which can be as low as < 1 [μs].
  • 3-way coupling - hard-sphere approach: interactions between particles are pair-wise (binary collision), instantaneous, momentum-conserving. Faster than soft-sphere method due to feasibility to use longer time steps.
For situations where volume concentration of discrete phase is high (> 10%), DDPM or Dense DPM (described in ANSYS FLUENT) has been developed which is based on statistical treatment of the particle-particle interactions and can be classified as CPFD (Computation Particle Fluid Dynamics).

Two-way coupling between continuous and discrete phase

In two-way coupling, dispersed phase modifies carrier phase turbulence where bigger particles usually enhance turbulence and smaller particles damp turbulence. This leads to situation where carrier phase turbulence becomes anisotropic. On the other hand, carrier phase turbulence enhances dispersion (fluctuating flow field) and clustering of the dispersed phase (stagnant flow field). When the flow is turbulent, it imparts a random motion on the particles based in instantaneous components of fluctuating velocity. However, RANS flow solution does not resolve all the small–scale turbulent eddies in the flow. The effect of eddies on particle trajectories is resolved using stochastic tracking. Here, a number as specified by "number of tries", particle streams are released from the same point. Each one is given a random "kick" in each grid cell based on the calculated turbulent intensity. This will indicate how turbulence will modify the trajectories.

Type of injection

Injection Types

In ANSYS FLUENT, for surface injection, mass flow rates of the particles from each face (each element of the boundary zone) can be scaled by the ratio of the area of the face particles are released from to the total area of the surface using option "Scale Flow Rate By Face Area" under Point Properties. Note that surfaces invariably have non-uniform distributions of points. In order to generate a uniform spatial distribution of particle streams released from a surface in 3D, one will need to create a bounded plane surface with a uniform distribution using the Plane Surface panel.

In ANSYS FLUENT, particle initial conditions (position, velocity, diameter, temperature and mass flow rate) can also be read from an external file to describe injection distribution. The file has the following form with all of the parameters in SI units. All the parentheses are required, name is optional.

(( x y z u v w diameter temperature mass-flow) name)

Mass flow rate: Only for coupled calculations (discrete phase affecting flow field of continuous phase), the mass of particles per unit time that follows the trajectory can be defined in the injection setting option. Also note that the number of streams (say 100) is not the same as the number of particles. The DPM uses a parcel method where many particles are tracked in each parcel and ANSYS FLUENT allocates the number of particles in each parcel based on the particle diameter, density and mass flow rate. Numper of particles in each parcel NP = mPS × Δt / mP where mPS = mass flow rate in particle stream, mP = mass of each particle and Δt is the time step. Under steady state operations, it is assumed that any particle released from the same location with the same conditions will follow the same trajectory. P_FLOW_RATE macro can be used in a UDF to return the flow rate of a parcel in [kg/s].

Example: particle release from a surface - assuming there are 200 faces in the surface.

  • The particle injected is sand with density of 2000 [kg/m3] and diameter 250 [μm]. Thus, mass of each particle is 8.18 [μg]
  • Total mass flow rate of sand particles = 0.1 [kg/s].
  • Total number of particles = 0.10 [kg/s] / 8.18×10-9 [kg] = 1.22×10-9 [particles/s]
  • Thus, the 200 trajectories releases from 200 faces of the surface represents 1.22×10-9 actual [particles/s].
  • If one defines a range of diameters by specifying "minimum diameter", "maximum diameter" and "number of diameters", the total number of trajectories computed will increase to 200 × "number of diameters".
  • If stochastic tracking is switched on to account for Turbulent Dispersion and number of tries specified is 20, now the total number of trajectories computed will be 200 × "number of diameters" × "number of tries".

Eulerian Wall Film
A wall film is a thin layer of liquid over a solid surface such as water dripping due to condensation. The surface tension is one of the dominant parameter which controls the thickness and flow rate of such liquid films.

multi-phase flow: Eulerian Wall Film in ANSYS FLUENT

Setting for Evaporation-Condensation Modeling

Interface Tracking

interface Tracking schemes in Mulitphase

The geometric reconstruction scheme represents the interface between fluids using a piecewise-linear approach. It assumes that the interface between two fluids has a linear slope within each cell and uses this linear shape for calculation of the advection of fluid through the cell faces.In ANSYS FLUENT this scheme has better accuracy as compared to other schemes (such as donor-acceptor scheme) and is applicable for general unstructured meshes.

Advanaced Solver Setting: AMG

The Algebraic Multi Grid solver is one among the family of multi-grid solvers. The adjective multi-grid refers to the concept that the solution process is implemented more than one grid (levels) even though the user has created only one grid.

Settings of a multi-grid solver in ANSYS FLUENT

In AMG solver, a 'virtual' coarse mesh equations are created by merging the original control volumes. The merged coarse control volume meshes are in general very irregular in shape. The coarse mesh equations thus impose conservation reduce the error components at longer wavelengths. The following image is taken from CFX user manual:

AMG control volume merging approach

Here, the results from coarse-mesh at each level is interpolated from/on the nearest finer mesh and vice-versa.

FSI: Fluid-Structure Interaction in ANSYS

FSI has wide application to check the effect of fluid forces on structural components such tank sloshing, flapping of wings and helicopter blades, pouring of liquid from a non-rigid container, pressure-limiting valves, leakage in fuel injectors operating at high pressure (such as Diesel engines) and check valves (deformable flow control devices) to name few.

FSI simulations are classified under a broader category of "Multi-physics Simulations" and not limited only to structural simulations. Fluid-thermal interaction, Fluid-thermal-electric interaction, Fluid-piezoelectric interaction, Structural-magentic interaction ... are other types of multi-physics simulations.

FSI setting in ANSYS workbench

Geometry import setting

Geometry import options

FSI Coupled flow and structure simulation

FSI System Coupling

In addition to the settings described above, following steps are required for a FSI simulation:
  • Solution update can be done through System Coupling component only.
  • System Coupling component overrides few settings in participant solver so that the time duration and time step settings are consistent across all participant solvers.
  • The surface of solid boundary in contact with fluids need to be set as "Fluid-Solid Interface"in Transient/Static Structural Model.
  • In fluid solver such s FLUENT, the "System Coupling" option needs to be activated under "Dynamic Mesh" tab.
  • For FLUENT, all regions of type 'Wall' are shown in System Coupling setup.
  • For (Transient/Static) Mechanical, all regions of type "Fluid-Solid Interface" appear in System Coupling component.
  • Create solid-fluid coupling: select a fluid (FLUENT) and solid (Mechanical) region pair (use Ctrl key for multiple selection) and select "Create Data Transfer" from right click pop-up menu.
  • "Data Transfer" defines the Source (participant solver, deforming region and variable being transfered e.g. force, displacement), Target (participant solver, deforming region and variable being transfered e.g. force, incremental displacement) and Data Transfer Controls (data transfer frequency, under-relaxation factor and convergence).
  • Post-processing a FSI result:
    1. Use 'Result' cell in Transient/Static Structural or Fluid Flow (FLUENT) for solver-specific post-processing
    2. Connect structural 'Solution' cell directly to FLUENT system 'Results' cell or (c)Add a 'Results' System (ANSYS CFD-Post) for unified post-processing of structural and fluid results.

Parallel Processing

In order to reduce the simulation run time, parallel processing methods have been developed where the arithmetic involved with matrices are broken into segments and assigned to different processors. Some of the the terms associated with this technology are described below. The image from "Optimising the Parallelisation of OpenFOAM Simulations" by Shannon Keough outlines the layout of various components in a cluster.

HPC Cluster Layout

  1. HPC: High Performance Computing - A generic term used to describe the infrastructure (both hardware and software) for parallel processing.
  2. Cluster: A collection of workstations or nodes connected with each other with high speed network such as 1 GB/s Ethernet network or InfiniBand.
  3. Node: Each independent component of a cluster or HPC set-up. Each node has following configurations:
    • Chassis: For example - HP Z820
    • CPU: For example - Intel XEON-E5 2687W
    • Cores: For example - 8 cores per CPU @ 3.1GHz
    • RAM: For example - 8GB DDR3-1600 Registered ECC memory
    • Storage: For example - 1TB SATA HDD plus network file system (NFS) server
    • Operating System: For example - CentOS 6.5
    • Application Program: For example - CFD Software: OpenFOAM or FLUENT V17.2
  4. Hyper-threading: It is a method which allows each core on the CPU to present itself to the operating system as two cores: one real and one virtual. The operating system can then assign jobs to the virtual cores and these jobs are run when the real core would otherwise be idle (such as during memory read/write) theoretically maximising the utilisation of the CPU.
    • Excerpts from STAR-CCM+ User Manual: "Generally, for best performance, it is recommended that you turn off any features that artificially increase the number of cores on a processor, such as hyper-threading. In situations where hyper-threading is turned on to benefit other applications, you should generally avoid loading more processes than there are physical cores.
       In addition to turning off hyper-threading, it is recommended that you set your BIOS settings to favor performance rather than power saving. Almost all heavily parallelized applications suffer performance problems if the CPU frequencies are spun up and down. In most parallel applications and STAR-CCM+ in particular, when one core is operating at a reduced frequency, all other cores are running at lower frequencies as well.
  5. Memory Intensive Applications: This refers to the programs whose performance is more influenced by RAM than the clock-speed of the CPU. For example - OpenFOAM
  6. Partitioninng: This is nothing but a "work distribution" among CPUs. However, tracking of nodes and elements are important. Excerpts from FLUENT user's manual: "Balancing the partitions (equalizing the number of cells) ensures that each processor has an equal load and that the partitions will be ready to communicate at about the same time. Since communication between partitions can be a relatively time-consuming process, minimizing the number of interfaces can reduce the time associated with this data interchange."

    Mesh Partitioning and Load Balancing in ANSYS FLUENT

  7. MPI: Message Passing Interface - Each compute node is connected to every other compute node though the network and relies on inter-process communication to perform functions such as sending and receiving arrays, summations over all cells, assembling global matrix - remember A.x = b. Inter-process communication is managed by a message-passing library by an appropriate implementation of the Message Passing Interface (MPI) standard such as OpenMPI, Intel MPI, Cary MPI, IBM Platform MPI, Microsoft MPI, MPICH2 or vendor's own version of MPU.
  8. Scalability: It refers to decrease in turnaround time for solutions as the number of compute nodes increases. However, beyond a certain point the ratio of network communication time (refer to definition of Cluster above) to computational time increases, leading to reduced parallel efficiency, so optimal system sizing is important for any simulations - determined from a ratio of the time to compute and the time that is taken to exchange data between cluster nodes.

Making FLUENT Runs on Clusters
FLUENT runs on a cluster can be started either in [A] interactive mode (as if it is running on your local machine with all functionality of a GUI) or [B] in a batch mode with GUI OFF and graphics ON or [C] batch mode with GUI OFF and graphics OFF. In case of mode 'C', the runs can only be controlled and monitored either through text files and/or through FLUENT Remote Visualization interface.

Graphics refers to the pictures, videos, charts and animations. The graphics components of a computer control and enhance how these 'graphics' are displayed on computer screen. Usually the graphics components are on a separate card that is plugged into a slot on the motherboard though these can and are built directly into the motherboard as well. Some other names used for the graphics components are video card, video adapter nd display adapter.

Similar to an operating system which deals with how a computer starts and shuts down, the graphics driver is a program that enables computer to use the graphics card and communicate with CPU, HDD and other components. A 'driver' in general is a software component that lets the operating system and a 'device' communicate with each other. E.g. there is USD 'driver' to make a USB mouse work. The graphics card sends the signal to the computer display and makes the data visible.

Note that computer 'monitors' are not same as "graphics cards or graphics drivers". Graphics card is a hardware and graphics driver is software to make it work. AMD Radeon a graphics card used in my laptop manufactured by Lenovo.

components windows 10

Graphics cards Lenovo Laptop

Starting FLUENT in batch mode in Linux

fluent3d -g -ibatch.jou-tN -driver opengl -rx 19.2 -mpi=hp -cnf=hosts.txt>& batch.out
2d -gr -t8 x1119.1intel
2ddp -gu -t32null 19.0ms
3ddp -t96mpich2

Starting FLUENT in batch mode in Windows

fluent3d -g -ibatch.jou-tN -driver opengl -rx 19.2 -mpi=hp -cnf=hosts.txt-wait
2d -gr -t8 x1119.1intel-hidden
2ddp -gu -t32null 19.0ms
3ddp -t96mpich2
  • -g indicates that the program (Cortex) is to be run without the GUI and graphics (Linux) or without graphics and with GUI minimized in the taskbar (Windows).
  • -gr will run Cortex without graphics. This option can be used in conjunction with the -i option to run a job in background mode.
  • -gu will run Cortex without the GUI but with the graphics window(s). (On Windows: it will run FLUENT keeping the window minimized that is the GUI will be available)
  • -i reads the specified journal file say batch.jou
  • -tN where N = number of processors = 2, 3, 4, ... 8, ... 32, ...
  • -driver allows to specify the graphics driver to be used in the solver session
  • -gu -driver null allow the generation and export of graphics files from FLUENT simulation run in batch mode
  • hosts.txt is the file with the names of the computational nodes (or the hosts)
  • Fluent 19.0 onwards, the solver does not differentiate between a serial mode and a parallel mode. Now the serial mode is in fact the general parallel mode with just one parallel task, -t1 option explained above.
  • For a ANSYS FLUENT session that is already in progress, the graphics window display can be disabled using the TUI command: display → set → rendering-options → driver → null
  • For a ANSYS FLUENT session that is already in progress, to re-enable a graphics window display that had been previously disabled, use the TUI command: (close-all-open-windows); display → set → rendering-options → driver → opengl

Reference: https://doku.lrz.de/pages/viewpage.action?pageId=27394531
> WARNING: Rank:0 Machine mpp2-xxxxxxx has 37 % of RAM filled with file buffer caches.
This can cause potential performance issues. Please use -cflush flag to flush the cache.
(In case of any trouble with that, try the TUI/Scheme command '(flush-cache)'.) 
In order to clear the buffer cache, add following lines at the beginning the journal file.
(rpsetvar 'cache-flush/target/reduce-by-mb 12288)
The number 12288 [Mb] = 12 x 1024 [kb] is the specification of the amount of system memory which cannot be flushed out due to boot image of Linux OS files in the memory of the diskless compute nodes. This means that on a compute node with 64 Gb, it tries to flush 64 [Gb] - 12.0 [Gb] = 52 [Gb] of memory. The value of 12288 also varies from cluster-to-cluster and may also change over the time.

Dimensionless Numbers

  • CFL NUMBER: CFL, it is ratio of characteristic flow velocity to the speed with which time-marching solution advances over a mesh.
  • Knudsen number: Kn, defined as the ratio of the mean-free path to the characteristic length scale.
  • Mach number: M, defined as the ratio of the characteristic velocity to the sound speed
  • Reynolds number: Re, defined as the ratio of the product of the characteristic velocity × characteristic length to the kinematic viscosity. This is one of the widely used number required to distinguish laminar and turbulent flow regimes.
  • Strouhal number: Sh, defined as the ratio of the mean-free time to the characteristic time scale. This number is used to express vortex-shedding frequency of wake regions behind bluff-bodies.
Contact us
Disclaimers and Policies

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.