Thursday, 15 October 2015

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.

Friday, 9 October 2015

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.