First to Higher Order Blending

This is just a quick post on first-to-higher-order blending. First-to-higher-order blending is used when you want something in between a first order upwind and a second order upwind scheme in your simulation. Say you are trying to solve, and you get good convergence in first order, but when you move to second order you have trouble getting a converged solution - you can use first-to-higher-order blending to make this change gradually.

I've written a scheme journal file that'll do this gradual ramp-up. Scheme can be used just like a normal journal file. I've found the main difference between scheme and journal files is that scheme code gives you more options. Here I'm using scheme to do a loop (which to my knowledge isn't possible in an ordinary journal file). I've found that one downside is that scheme is more difficult to cancel & quit halfway through a calculation. It can also be more difficult to write.  

This scheme file slowly ramps up the first-to-higher-order blending value from 0 to 1, linearly over 1000 iterations (0,0.001,0.002,...). It'll then run the full second (or whatever) order sim. for another 1000 iterations. To use it take your first order case file, then change the "solution methods" from first order upwind to the new discretisation you want, then read the scheme file. You'll first need the scheme code saved as a .scm file in the directory. I do this using notepad & recommend it (copy and paste the code, go to save as, type: all files, add .scm to the end of the file name).

;;------------------------------------------------------------------
;;scheme function for ramping-up to second order
;;ramps up over 1000 iterations
;;JRH 18_12_15 cfdyourself.blogspot.co.uk
;;------------------------------------------------------------------
(do ((i 0 (+ i 1)) ) ((> i 999.1));;do loop for 1000 iterations
(define r (/ i 1000))
(ti-menu-load-string (format #f "/solve/set/numerics no yes no no no ~a " r));;
(ti-menu-load-string "/solve/iterate 1 ");;run the simulation for 1 iteration
);;end of the do loop
(ti-menu-load-string "solve/set/numerics no yes no no no 1 ")
(ti-menu-load-string "solve/iterate 1000 ")

Numerical Diffusion vs. Numerical Dispersion

Numerical diffusion and dispersion are a major error source in CFD studies. In informal discussion the terms may be used interchangeably or along with others (eg: “Numerical Dissipation”, “Numerical viscosity” etc.) but the two terms have quite separate formal meanings.

Numerical diffusion is the tendency for transported variables to diffuse more than they should. Compared to exact solutions with diffusion terms (viscosity, mass diffusivity etc.) that are physically realistic, the numerical diffusion is an error that adds to this diffusion, “smudging” the results. In a grid independence study (GIS) this is one of the errors you are testing for. Say you observe a velocity profile in your GIS, with mesh refinement you will see that peaks will become more defined as this “smudging” effect is reduced.

Figure 1 – The effect of numerical diffusion (Andersson et al., 2009).

Numerical diffusion occurs when 1^st order discretisation are used such as the first order upwind scheme.  It is sometimes referred to as “numerical viscosity” because the error behaves as if there was an addition to the diffusion term (which is viscosity in the case of fluid flow).

Numerical dispersion occurs when a higher order discretisation scheme is used to improve accuracy of the result. Numerical dispersion often takes the form of so-called ‘spurious oscillations’. The difference between numerical diffusion and dispersion is often depicted as an approximation to a step change as shown in Figure 2. The spurious oscillations cause problems in the flow, for example if this step change was in a variable that can only be between 0 and 1 (such as mass fraction for example), spurious oscillations can lead to unphysical values.

Figure 2 – The difference between numerical diffusion and numerical dispersion (Wendt, 2009).

This is due to the truncation error of the discretisation. If the truncation is an odd-order method (such as the first order unwind method) the leading truncation error is even. If the truncation is, say, a second order upwind method, the leading truncation error is odd. Even order derivatives in the truncation error contribute to numerical diffusion, odd order derivatives contribute to numerical dispersion. These represent real and imaginary parts respectively in the Von Neumann amplification factor (more information in Hirsch, (2007)).

Ideally we would like to eliminate both errors from our simulation result. But this is hindered by Godunov’s Theorem: “linear numerical schemes for solving PDE’s having the property of not generating new extrema (monotone scheme) can at most be first order accurate”. First-order upwind schemes are the least diffusive of the first order schemes.

Some more points to finish:

• Shock waves are often resolved using artificial viscosity, deliberately including an extra viscosity term to produce something like the LHS of Figure 2.
• FLUENT contains an option called “first to higher order blending” that will switch to first order
• For higher-order methods a loophole exists, involving the non-linear limiter functions that enforce monotonicity (Wendt, 2009).

That’s all I know about it so far. I’m interested in finding out more regarding how to actually identify dispersive errors in your simulation. That’s what I’m hopefully going to look into in the future.

References:

Andersson, B. et al. (2012). Computational Fluid Dynamics for Engineers. Cambridge University Press.

Hirsch, C. (2007). Numerical Computation of Internal and External Flows, 2^nd ed. Elsevier.

Wendt, JF. (2009). Computational Fluid Dynamics an Introduction. Springer.

Buoyancy effects in heat exchangers & CFD models

Conventionally, two regimes of flow are typically discussed in the context of heat exchange. Forced flow, where the flow is driven by a pump, fan, stirrer or other prime mover, is by far the most common process used industrially. This mode provides efficient, predictable, controllable heat transfer. Free convection (or natural convection) is the opposite case. In this case, the flow is driven entirely by buoyancy – hotter fluid is denser, and so rises in the surrounding fluid. Industrial examples include tank heaters (used in storage tanks containing viscous or high melting point liquids). It’s worth pointing out that both forced flow and free convection can be laminar or turbulent.

In heat exchangers a combination of these two flow regimes can occur, known as mixed convection. In mixed convection, the buoyancy of the heated or cooled fluid influences the forced flow. Gas flows under high temperature gradients and low flowrates are particularly sensitive to mixed convection effects. The contribution of mixed convection changes the velocity profile, and can lead to new flow structures or changes in the boundary layer. So the possibility of mixed convection has consequences in laminar-turbulent transitions, prediction of heat transfer and in computational fluid dynamics modelling of heat exchange surfaces.

Particularly in CFD, when buoyancy forces are in equilibrium with velocity forces, this can make convergence difficult.

You can determine whether mixed convection is important in your modelling work using the Richardson number. The Richardson number determines whether the flow is dominated by buoyancy (Ri above 1) or forced flow (Ri below 1). Richardson numbers between 0.1 and 10 indicate mixed convection. Richardson number is the ratio of the Grashof number to the square of the Reynolds number.

In my simulation work, I’ve used the Richardson number to try and give me an idea of whether buoyancy effects were affecting the convergence. It turned out that the Richardson number in my work is very small, and so this event is unlikely. I also compared the convergence of two runs (with gravity on and with gravity off) and got similar results. Buoyancy relies on gravity to take effect, so this showed me that the problems I am having with convergence are unlikely to be due to buoyancy effects in the fluid. So I am satisfied with this level of detail, but mixed convection heat exchange and buoyancy modelling in CFD is a rabbit’s hole of a subject area, and there are many aspects of the problem I haven’t covered in this blog entry.

Including a diameter-dependent Cunningham slip correction factor in Fluent.

Very small particles (typically < 1 micron) do not obey Stokes law for drag when they are suspended in a gas. The particles are so small that they experience interactions with the gas not as a continuum, but as individual collisions. The main outcome of this is that the assumption of a no-slip boundary on the particle’s surface is no longer true. Small particles have less drag than is calculated by Stokes law, and this can be corrected for using Cunningham slip correction factor.

The Cunningham slip correction factor can be set in Fluent when defining the injection. So, say you have a particle injection of 0.1um diameter, you can set the drag law to “Stokes-Cunningham” and add the Cunningham slip correction factor as a constant (say “11.0”) that is relevant to the size of that particle.

But what if you’re particle is changing in size as it moves through the domain ? What if you particle is growing through condensation, or losing mass through evaporation ? What if you inject a distribution of particle sizes rather than just one size ? In these cases, the particle will have a diameter-dependent Cunningham slip correction factor – it won’t be just one value.

There’s no option to support this in fluent, so I have wrote a quick UDF that deals with it. I have tried to stick to the Fluent defaults as much as possible. To calculate the drag, I’ve used the same source as Fluent for it’s default, “spherical” drag law (Morsi and Alexander, 1972). To calculate the Cunningham Slip correction factor, I’ve used the formula given in the Fluent Theory Guide (page 386, eqn 16.69). it’s worth noting here that the constants used in the Cunningham correction formula here are for air, if you’re looking at a different gas you’ll need different constants (see the Wikipedia article: https://en.wikipedia.org/wiki/Cunningham_correction_factor). I’ve used the mean free path formula for air, taken from Flaggan and Sienfield (2012, p. 296) and added the values that are constants assuming a pressure of 1atm:

equation 1: mean free path of air

Fluent’s receives the drag force in the form of the dimensionless group 18*Cd*Re/24. So as you’ll see below the drag force is converted to this before it is returned to Fluent, and the Cunningham correction is also applied as appropriate.

Here is the resulting UDF, which must be “complied” and must be “hooked” as described in the Fluent UDF manual for Define_dpm_drag. Please feel free to use my code, as long as you reference my blog.

/*the following implements a diameter-dependent Cunningham slip correction factor*/

/*Actual drag coeffs used by Fluent have been referenced-chased (Morsi, SA. Alexander, AJ. (1972). An investigation of particle trajectories in two-phase flow systems)*/

/*CFDYouself.blogspot.co.uk, JRH, 09.10.15 */

#include "udf.h"

#include "math.h"

DEFINE_DPM_DRAG(cunningham_drag_force, Re, p)

{

cell_t c = P_CELL(p);

Thread *t = P_CELL_THREAD(p);

real lambda;

lambda = 0.0002135166*C_MU_L(c, t)*sqrt(C_T(c, t));

real Cc;

Cc = 1 + (2 * lambda / P_DIAM(p))*(1.257 + 0.4*exp(-1.1*P_DIAM(p) / (2 * lambda)));

real Cd, drag_force;

if (Re < 0.1)

{

Cd = 24 / Re;

drag_force = 18*Cd*Re/(24*Cc);

return (drag_force);

}

else if (Re < 1)

{

Cd = 22.73 / Re + 0.0903 / (Re*Re) + 3.69;

drag_force = 18 * Cd*Re / (24*Cc);

return (drag_force);

}

else if (Re < 10)

{

Cd = 29.1667 / Re - 3.8889 / (Re*Re) + 1.222;

drag_force = 18 * Cd*Re / (24*Cc);

return (drag_force);

}

else /*if (R < 100)*/

{

Cd = 46.5 / Re - 116.67 / (Re*Re) + 1.6167;

drag_force = 18 * Cd*Re / (24*Cc);

return (drag_force);

}

/*included incase needed later - not needed here.

else if (Re < 1000)

{

Cd = 98.33 / Re - 2778 / (Re*Re) + 0.3644;

drag_force = 18 * Cd*Re / 24;

return (drag_force);

}

else if (Re < 5000)

{

Cd = 148.62 / Re - 4.75e+4 / (Re*Re) + 0.357;

drag_force = 18 * Cd*Re / 24;

return (drag_force);

}

else if (Re < 10000)

{

Cd =  -490.546/ Re - 5.79e+5 / (Re*Re) + 0.46;

drag_force = 18 * Cd*Re / 24;

return (drag_force);

}

else

{

Cd =  -1662.5/ Re - 5.42e+6 / (Re*Re) + 0.5191;

drag_force = 18 * Cd*Re / 24;

return (drag_force);

} */

}

Morsi, SA. Alexander, AJ. (1972). An investigation of particle trajectories in two-phase flow systems. Journal of Fluid Mechanics, 55(2), pp 193-208.

Ansys, (2013). Ansys fluent theory guide release 15.0. Canonsburg: Ansys inc.

Flaggan, R. and Sienfield, J. (2012) Fundamentals of air pollution engineering. Mineloa: Dover Publications.

Correction for Heterogeneous Nucleation

heterogeneous nucleation

Frequently studies which use classical nucleation theory account for heterogeneous nucleation using a heterogeneous nucleation correction factor. This is normally denoted f(θ) or something similar, and is inserted into the classical (homogeneous) nucleation rate equation to account for the effect of nucleation on a surface. In this blog entry I’m going to cover this factor. The “flat surface” correction is common because in most of the cases commonly studied the foreign nuclei is much larger than the critical cluster that is forming on its surface, so the contact surface between them can be considered flat. A more complicated version of f(θ) can be derived for a curved foreign nucleus, derivation is covered in Vehkamaki (2006) ch. 7. This subject is something of a departure from the normal subjects discussed on this blog, as this correction to my knowledge is not available in default settings of any particle formation laws etc. in Fluent (Fluent is the software I currently use exclusively for CFD). But It’s still of interest to the particle formation stuff I’m working on currently. 

Heterogeneous nucleation has an effect that is similar to catalysis in reactions. It lowers the energy barrier. Heterogeneous nucleation happens faster at lower degrees of supersaturation due to this change in the critical Gibbs free energy.

Figure 1 – A graph showing the effect of heterogeneous nucleation on critical energy (Zang, 2015).

The correction factor was originally derived by Fletcher in 1958, and is sometimes referred to as the “Fletcher factor” (Kalikmanov, VI. 2013). But Equation 1 reduces to the flat-surface correction when the curved surface is very large.

 

Equation 1 – The long-winded version of Fletcher factor, for curved surfaces. Dimensionless versions are common. Here a is the ratio between the foreign nucleus and the cluster; m is the cos of the (cluster) contact angle; w = √1+a2-2ma.

The flat surface correction is:

Equation 2 – The flat-surface correction. υ is the contact angle (between the critical droplet on the surface).

For the heterogeneously nucleating critical droplet, the radius is exactly the same as for the homogeneous case, but the volume of the nucleating particle changes because the particle becomes a spherical cap instead of a sphere. Equation 2 corrects for this reduced volume (which reduces the critical Gibbs free energy, the critical Gibbs free energy is the point where the (-ve) volume free energy change balances the (+) free energy change needed to create a surface). The angle can be found by young’s equation if the surface tensions between the gas, cluster and substrate are known.

Figure 2 - Effect of heterogeneous nucleation on the nucleation rate. shamelessly copied-and-pasted from Coulson 2.

How it is applied in the nucleation rate equation (shown in Equation 3). The biggest difference comes in the form of the effect on ΔG*het. Typically the pre-exponential factor may be known experimentally (typical value is around 1017 cm-2s-1 for heterogeneous vapour condensation). The monomer flux (W*) is the sum of both the molecular-bombardment style growth present in homogeneous nucleation (R* = monomer flux, s* = surface area of critical droplet) but also adatom impingement (W*s) – surface diffusion of adsorbed atoms to the nucleus periphery.

Equation 3 - Difference between heterogeneous and homogeneous nucleation, (Cooper, 2015).

Most sources simply state that the dependency of Jhet on the corrected free energy, and do not discuss the pre-exponential factor in detail. However, Fletcher (1959) himself states that heterogeneous nucleation requires changes to this pre-exponential factor (though he does not cover what these changes are in the source). Talanquer and Oxtoby (1996) give an example of an analytical expression for this pre-exponential factor in their paper. However, overall this is all there is to the classical treatment of heterogeneous nucleation.

Cooper, S. (2015). Crystallisation Kinetics. Durham University. Available: http://community.dur.ac.uk/sharon.cooper/lectures/cryskinetics/handoutsalla.html accessed: 17 Aug 2015.

Fletcher, NH. (1959). On ice-crystal production by aerosol particles. Journal of Meteorology. 16(2). pp. 173-180.

Kalikmanov, VI. (2013). Nucleation theory. Berlin: Springer.

Vehkamaki, H. (2006). Classical nucleation theory in multicomponent systems. Berlin: Springer.

Zang, L. (2015). Lecture 12 Heterogeneous Nucleation: a surface catalysed process. College of engineering University of Utah. Available: http://www.eng.utah.edu/~lzang/ accessed: 12 Aug 2015.