The numerical simulation of fluid flow and heat/mass-transfer phenomena requires the numerical solution of the Navier-Stokes and energy-conservation equations coupled with the continuity equation. Numerical or false diffusion is the phenomenon of inserting errors in the calculations that compromise the accuracy of the computational solution. The Taylor series analysis that reveals the truncation/discretization errors of the differential equations terms should not be termed as false diffusion. Numerical diffusion appears in multi-dimensional flows when the differencing scheme fails to account for the true direction of the flow. Numerical errors associated with false diffusion are investigated via two- and three-dimensional problems. A numerical scheme must satisfy necessary criteria for the successful solution of the convection-diffusion formulations. The common practice of approximating the diffusion terms via the central-difference approximation is satisfactory. Attention is directed to the convection terms since these approximations induce false diffusion. The equations of all the conservation equations in this study are discretized by the finite volume method.
- false diffusion
- numerical dispersion
- computational errors
- finite volume method
Numerical diffusion is a significant source of error in numerical solution of conservation equations and can be separated into two components namely cross-stream and streamwise numerical diffusion. The former occurs when gradients in a convected quantity exist perpendicular to the flow and the direction of the flow is oblique to the grid lines, i.e. due to the multi-dimensional nature of the flow. The latter happens when gradients in a convected quantity exist parallel to the flow even in one-dimensional situations [1, 2, 3].
According to Vahl Davis et al.  a theoretical approximation for the false-diffusion term in two-dimensional geometries is:
where is the resultant velocity and is the angle, which lies between 0 and 90°, made by the velocity vector with the x-direction. Expression (1) reveals that there is no false-diffusion present for θ = 0 and 90°, and that it is at its maximum at θ = 45°.
The first-order accurate upwind-difference scheme is more stable than the second-order accurate central-difference scheme or higher-order accurate schemes, e.g., QUICK at high grid Peclet numbers. Non-linear schemes, e.g., van LEER scheme appear highly stable but they do not deal with the flow direction inclination at all [5, 6, 7, 8].
One way to overcome diffusion errors is to use an upwind approximation which essentially follows the streamlines. This approach, originally derived by Raithby , is formally called the skew-upwind differencing scheme. The CUPID (Corner UPwInDing) scheme retains the general objectives of the Raithby approach but uses an entirely different 2D formulation that eliminates the shortcomings of the original scheme. The SUCCA (Skew Upwind Corner Convection Algorithm) scheme is based on a simplified formulation of the CUPID scheme [10, 11].
In this chapter the performance of various numerical schemes is compared as far as the accuracy of the numerical flow prediction is concerned. Various cases of numerical problems that include one- and two-phase flows, heat and mass transfer, geometry of two and three-dimensions, grid of Cartesian and curvilinear coordinate system, straight and inclined inflows are investigated and the conclusions are presented. Among the numerical schemes that are evaluated is the SUPER (Skew Upwind and cornER algorithm) numerical scheme that takes into account the 3D flow orientation and follows all the rules for the successful solution of the equations .
Numerical dispersion comes about when the convective scheme used is unstable and large gradients are present . The use of the upwind or hybrid numerical scheme ensures the stability of the calculations but the first-order accuracy makes them prone to streamwise numerical diffusion errors. Higher-order schemes involve more neighbour points and reduce the streamwise false-diffusion by bringing in a wider influence. These schemes are based on one-dimension formulations that do not take into account the flow direction which makes them unstable at high grid Peclet numbers . Furthermore numerical schemes that take into account the flow direction become unstable at high grid Peclet numbers in regions of high convection discontinuities. The phenomenon of the numerical dispersion causes spatial oscillations that lead to many numerical errors and non-physical results [1, 2].
Generally the numerical schemes should follow some basic rules in order to be able to calculate a successful solution:
The formulation of the scheme should ensure the flux consistency at the cell faces in adjacent control volumes satisfying the conservative property.
Conservation principle should be retained for all the variables and in the whole domain.
The scheme should recognize the relationship between the magnitude of the Peclet number and the directionality of influencing known as transportiveness property.
The influence coefficients in the finite volume equations should be always positive.
The neighbouring coefficients for the solution control volume of central node P should obey the following relation: , so as to be consistent with the differential transport equation.
Satisfying the Scarborough criterion ensures that the resulting matrix of coefficients is diagonally dominant which is a sufficient condition for a converged iterative method.
Linearizing the source term with a negative slope, so that the incidence of unbounded solutions is reduced.
The formulation of the scheme should be easy to implement without increasing the computational cost.
The method employed in this study follows the so called “control volume” approach which is developed by formulating the governing equations on the basis that mass, heat and momentum fluxes are balanced over control volumes. If there are two phases the control volume can be regarded as containing a volume fraction of each phase that obeys the following relation: . Each phase is treated as a continuum in the control volume under consideration. The phases share the control volume and they may, as they move within it, interpenetrate [14, 15].
1.1. Mathematical modeling
where the dependent variable of each phase ; , the volume fraction of each phase, (m3/m3); , the density of each phase, (kg/m3); , the velocity vector of each phase, (m/s); , within-phase diffusion coefficient; and , within-phase volumetric sources, (kg/m3 s).
The dependent variables solved for are:
Pressure that is shared by the two phases, (Pa)
Volume fractions of each phase, (m3/m3)
Three components of velocity for each phase, , (m/s)
Turbulence kinetic energy and dissipation rate of turbulence for the first phase, and (m2/s2)
Temperature (K), enthalpy (J/kg)
The phase volume fraction equation is obtained from the continuity equation:
where = phase volume fraction, (m3/m3); = phase density, (kg/m3); = phase velocity vector, (m/s); and = net rate of mass entering phase i from phase j, (kg/m3 s), if there is phase change.
The discretized by the finite volume method form of all the conservation equations is solved by the SIMPLEST and IPSA algorithms embodied in the Computational Fluid Dynamics (CFD) code PHOENICS .
2. Numerical problems
2.1. Inclined flow of a scalar quantity
The physical problem presented in this section describes the transport of a scalar quantity (C1) in a two-dimensional geometry. The flow of C1 is inclined to the grid lines (θ = 45°) and the natural diffusion is assumed equal to zero. This physical problem is considered as a benchmark for comparing the performance of various numerical schemes.
The boundary conditions applied to the grid are the following:
Air enters the domain from the west side with stable mass flow rate , uniform values of the two components of velocity and uniform value of the convected quantity (C1) .
Air enters the domain also from the south side with stable mass flow rate , uniform values of the two components of velocity and uniform value of the convected quantity (C1) .
At the two outlets (north and east side) of the domain the air is supposed to exhaust at an environment of fixed uniform pressure.
The dimensions of the domain are 1 × 1 m. The numerical results presented are based on the independent grid of 33 × 33 cells. For the discretization of the convected term in the transport equation of the scalar quantity C1 the numerical schemes: (a) UPWIND, (b) van LEER, (c) SUCCA, (d) SUPER are applied [6, 8, 11, 12].
In Figure 1 the velocity vector distribution predicted by the SUPER numerical scheme is presented.
In Figure 2 the vertical (concentration) distribution of the scalar quantity C1 in the middle of the domain is presented.
The vertical distribution of the scalar quantity C1 that is transferred by air with inclined direction is predicted more abrupt by the numerical schemes (SUCCA, SUPER) that take into account the flow direction than the conventional numerical scheme (UPWIND) and the non-linear numerical scheme (van LEER). It is concluded that the numerical diffusion errors are minimized when the discretization schemes (SUCCA, SUPER) are applied to predict the transport of the scalar quantity C1 when the air-flow direction appears the major inclination angle 45° to the grid lines. The performance of the van LEER numerical scheme is also satisfactory in predicting the scalar concentration distribution. The conventional UPWIND numerical scheme that does not take into account the phenomenon of false-diffusion presents the poorest accuracy.
2.2. Tubular flow of a scalar quantity and heat conduction in the radial direction
A physical problem used in this study for assessing the numerical schemes is the numerical prediction of the energy transfer in a triple tube heat exchanger . The heat exchanger has a three-dimensional curvilinear geometry similar to a circular cylinder. Three fluids are being considered to flow inside the domain that is separated by solid. Chilled water (10°C) flows in the inner tube, hot water (70°C) flows in the inner annulus and normal tap water (18°C) flows in the outer annulus. The three fluids follow a co-current flow. The numerical model is first validated and then used to compare the performance of various numerical schemes.
The boundary conditions applied for solving the above physical problem are the following:
Water with uniform properties (temperature, density, thermal capacity, kinematic viscosity) enters the three different inlets of the heat exchanger. The heat exchanger is separated to three components that are filled with water of different properties and the solid parts that prevent the water vertical flow. Heat is exchanged between the fluids through conduction. The solid is steel at 27°C of uniform properties (density, viscosity, specific heat, thermal conductivity, thermal expansion coefficient, compressibility). At the three outlets of the heat exchanger the boundary condition of zero mass flow is applied.
In Figure 3 the numerical prediction of the temperature distribution along the heat exchanger is validated against experimental data. The average% relative error of the numerical prediction does not exceed the value of 20%.
In Figure 4 the numerical results of the enthalpy distribution at the lateral plane in the middle of the heat exchanger applying three different numerical schemes (HYBRID, van LEER, SUPER) [6, 8, 12] is presented.
Comparing the performance of the numerical schemes in predicting the transfer of the scalar quantity H1 it is concluded that the numerical prediction is equally satisfactory for the three discretization schemes.
2.3. Airflow in a backward facing step
A physical problem appropriate for assessing the accuracy of various numerical schemes is the prediction of the airflow velocity vectors distribution in a backward facing step.
The separation of the airflow, the recirculation zone formed forward the step and the reattachment length are the main characteristics of the flow in this geometry. The numerical solution of the airflow field in a backward facing step is of major concern due to the difficulty in accurate prediction of these complex phenomena. The flow field may become more complex in the presence of particles.
In this study the two-dimensional backward facing step of Benavides and Wachem  is investigated.
According to the numerical results that agree well with the experimental data the performance of the numerical schemes (UPWIND, HYBRID, van LEER, SUPER) is equally satisfactory. The same conclusion is derived for the case of air-particles flow in a backward facing step.
2.4. Inclined air-particles flow
A physical problem that has been used for assessing the performance of various numerical schemes is the dispersed air-particles flow in a model room geometry (Figure 7) . A two-phase Euler-Euler mathematical model is applied to calculate the oblique inflow in the interior of the internal space. The accuracy of the four numerical schemes: (a) the conventional first-order UPWIND numerical scheme, (b) the second-order HYBRID scheme, (c) the non-linear van LEER and (d) the flow-oriented SUPER scheme [6, 8, 12] are compared in the case of inclined inflow (θ = 45°).
2.4.1. Boundary conditions
At the inlet the mass flow rate of each phase is multiplied by its volume fraction. The dispersed turbulent air-particles flow enters the room with uniform velocity (0.225 m/s) in an inclined direction. Turbulence is modeled by the RNG k-ε model . Particles are assumed to be transported and dispersed due to turbulence of the carrier fluid (air).
The turbulence kinetic energy of the air-phase that is applied at the inlet is defined as [13, 21] where . The mean air inlet velocity and the turbulence intensity, considered as 6%. The dissipation rate is given by , where the turbulence length scale is taken as = 0.07 (hydraulic diameter of the duct) and = 0.0845 an empirical constant of the turbulence model. The mathematical model uses the logarithmic “wall functions” near the solid surfaces.
At the outlet both air phases are supposed to exhaust at an environment of fixed uniform pressure. At the walls the no-slip and no-penetration condition is applied for both phases.
As far as the accuracy of the four numerical schemes in the case of inclined inflow (θ = 45°) is concerned the higher-order and non-linear schemes (hybrid and van Leer numerical scheme) present a similar performance with the first-order upwind numerical scheme and a different performance than the flow-oriented scheme (SUPER scheme). The vertical w1 velocity distribution predicted by the SUPER scheme presents a more abrupt and accurate profile due to the successful minimizing of the false-diffusion errors.
2.5. Heat and mass transfer
A physical problem that has been used to evaluate the performance of the numerical accuracy is the water-vapor condensation of humid air in the three dimensional geometry of a real-scale indoor space. A two-phase flow Euler-Euler mathematical model has been developed, wherein the humid air and water droplets are being treated as separate phases. The two phases exchange momentum and energy and, as the temperature drops below the dew point of humid air, mass transfer and phase change of water vapor to liquid takes place. The flow of humid air inside the room is buoyancy driven in the temperature range of 290–303 K. The properties of humid air (enthalpy, relative humidity, concentration of water vapor, saturation vapor pressure) vary with the temperature . The dimensions of the domain are: width (X) × height (Y) × length (Z) = 4.0 × 4.0 × 8.0 m.
2.5.1. Boundary and initially conditions
The interior of the domain is filled with humid air of 303 K and no liquid phase. The temperature on the surface of the walls is 303 K and on the surface of the cold bottom is 290 K. The density of humid air is 1.16 kg/m3 and the water droplets have a constant density (996 kg/m3) and specific heat (4190.0 ) at the dew point temperature. The initial pressure inside the domain is equal to the atmospheric pressure. The initial relative humidity condition tested is 90%. The overall water vapor mass of humid air at the initial temperature 303 K is taken from the psychrometric chart .
In Figure 9 the vertical temperature distribution in the middle of the domain at the height (0–4 m) is presented.
The temperature of humid air near the floor is below the dew point (301.2 K) and water phase humidity covers the whole surface. The hot air close to the floor that comes into contact with the cold surface, reduces its temperature and flows down due to the gravity. The remained hot air flows up to the roof.
Temperature profile in the region of major gradient near the floor surface is predicted more abrupt by the SUPER scheme.
In Figure 11 the vertical absolute humidity ratio (kg H2O/kg of dry air) predicted by the three different numerical schemes (HYBRID, van LEER, SUPER) is presented at time 360 s.
Heat convection is accompanied by mass transfer and phase change of humid air. The larger gradient of temperature profile predicted by the SUPER scheme  leads to the formation of larger amount of water phase. Comparing the performance of the discretization schemes a more accurate solution of the condensation procedure is observed when applying the SUPER scheme.
Numerical diffusion is the main source of errors in computational calculations. In this study the performance of various numerical schemes is evaluated in order to find the limitations of the convection formulations. The first-order UPWIND scheme presents more stability and overcomes the streamwise diffusion but such the higher-order schemes do not deal at all with the cross-stream diffusion. The second-order hybrid numerical schemes take into account the influence of more neighboring points but appears unstable at high Peclet numbers. The non-linear van LEER numerical scheme is highly stable and reduces satisfactory the false diffusion errors. The flow oriented SUPER scheme overcomes the phenomenon of the numerical diffusion in most of the cases investigated without increasing the computational cost.
The first author (D. P. Karadimou) gratefully acknowledges the financial support from the State Scholarships Foundation of Greece through the “IKY Fellowships of Excellence for Postgraduate studies in Greece - SIEMENS” Program.