Numerical Simulation of Biofuels Injection

The more important biofuels currently under investigation are the bio-alchohols and their derived ethers, and the vegetable oils and their derived esters. Methyl esters of rapeseed oils or soybean oils have been tested in Diesel engines, and in spite of the strong dispersion of the published results, there are indications that their use is a promising solution to the problems originated with the raw vegetable oil due to their higher viscosity, boiling temperature, final temperature of distillation and point of obstruction of cold filter (Tinaut, 2005).

dynamics of a single droplet.Abramzon and Sirignano (1989) and Berlemont et al. (1995) presented droplet vaporization models for the case of a spray in stagnant surroundings, and showed that the convective effects were most relevant.The same type of configuration was studied by Chen and Pereira (1992), and the predictions were found to follow satisfactorily the measurements presented.More recently, Sommerfeld (1998) presented a study on stationary turbulent sprays, using a droplet evaporation model based on the model proposed by Abramzon and Sirignano (1989) If special attention is dedicated to the biofuels injection and evaporation, then practically no numerical or experimental studies can be found.Recently Bai et al. (2002) presented a most relevant numerical study of a spray in wind tunnel using the Arcoumanis et al. (1997) experiments, but concentrated on the development of the spray impingement model and the fuel used was gasoline.

Mathematical model
This section describes the mathematical model for turbulent particle dispersion and vaporization assuming that the particles are sufficiently dispersed so that particle-particle interaction is negligible.
The particle phase is described using a Lagrangian approach while an Eulerian frame is used to describe the effects of both interphase slip and turbulence on particle motion using random-sampling techniques (Monte Carlo).It is also assumed that the mean flow is steady and the material properties of the phases are constant.
When vaporizing droplets are involved in the simulations, two-way coupling must be accounted for since the phase change modifies the characteristics of the fluid phase.The vapour produced by the droplets is a mass source for the fluid; moreover the vaporization process generates modifications in the momentum and energy balances between both phases.Fluid phase equations then contain many extra-source terms.It is assumed that the vapour production does not significantly modify the fluid phase density.
The method to solve the continuous phase is based on the solution of the conservation equations for momentum and mass.Turbulence is modelled with the "k-" turbulence model of Launder and Spalding (1974), which is widely and thoroughly tested, and was found to predict reasonably well the mean flow (Barata, 1998).In order to reduce the numerical errors to an acceptable level, the higher-order QUICK scheme of Leonard (1979) is used to evaluate the convection terms.A similar method has been used for threedimensional (Barata, 1998) or axisymmetric flows (Shuen et al., 1985;Lilley, 1976;Lockwood & Naguib, 1975) and only the main features are summarized here.
The governing equations (continuity, momentum, turbulent kinetic energy, dissipation, enthalpy, and vapour mass fraction) constitute a set of coupled partial differential equations that can be reduced to a single convective-diffusive conservation equation of the form where   is the effective diffusion coefficient for quantity  .The term on the left-hand side is the convection term, whilst the first and the second terms on the right-hand side are the diffusion term and the source term, respectively.
The source term S  as divided into two parts, which yields the following expression: where g S  , specifies the source term of the gas and p S  , specifies the source term of the particle.
Vaporization phenomena are described in the present study assuming spherical symmetry for heat and mass transfers between the droplet and the surrounding fluid, and convection effects are taken into account by introducing empirical correlation laws.
The main assumptions of the models are: spherical symmetry; quasi-steady gas film around the droplet; uniform physical properties of the surrounding fluid; uniform pressure around the droplet; and liquid/vapor thermal equilibrium on the droplet surface.
The effect of the convective transport caused by the droplet motion relative to the gas was a accounted for by the so called "film theory", which results in modified correlations for the Nusselt and Sherwood numbers.For rapid evaporation (i.e boiling effects) additional corrections were applied.The infinite droplet conductivity model was used to describe the liquid side heat transfer taking into account droplet heat-up.Hence, the following differential equations for the temporal changes of droplet size and temperature have to be solved.
Under the assumption of steady state conditions in the gas film and assuming a spherical control surface around the droplet, the total mass flow through this surface will be equal to the evaporation rate m  : The heat penetrating into the droplet can be expressed by: The mass transfer number M B as defined as where Fs Y is the fuel mass fraction on the droplet surface and defined as: For any given value of surface temperature, the vapor pressure is readily estimated from the Clausius-Claperyon equation as exp 43 where a and b are constants of the fuel.
The latent heat of vaporization is given by Watson (1931) as The parameters * Nu and * Sh are the "modified" Nusselt and Sherwood numbers, and tend to 0 Nu and 0 Sh , respectively, as T F and M F tend to the unity.
In the case of an isothermal surface and constant physical properties of the fluid, the problem has a self-similar solution and the correction factors M F and T F do not depend on the local Reynolds number.It was found that the values M F and T F are practically insensitive to the Schmidt and Prandtl numbers and the wedge angle variations, and can be approximated as where   FB is given by Nu and 0 Sh are evaluated by the Frossling correlations: The evaporation rate m  with convection is: The Schmidt number and the Prandtl number are equal assuming a Lewis number of unity.Equation 20 has the advantage that it applies under all conditions, including the transient state of droplet heat-up, whereas Eq. ( 31) can only be used for steady-state evaporation.
Finally the evaporation rate m  is: And the equations for the temporal changes of droplet size and temperature are: Of the air/vapor mixture in the boundary layer near the droplet surface according to Hubbard et al. (1973), the best results are obtained using the one-third role (Sparrow & Gregg, 1958), where average properties are evaluated at the following reference temperature and composition: For example, the reference specific heat at constant pressure is obtained as The dispersed phase was treated using the Lagrangian reference frame.Particle trajectories were obtained by solving the particle momentum equation through the Eulerian fluid velocity field, for a sufficiently high number of trajectories to provide a representative statistics.
The equations used to calculate the position and velocity of each particle were obtained considering the usual simplification for dilute particle-laden flows (Shuen et al., 1985).Static pressure gradients are small, particles can be assumed spherical and particle collisions can be neglected.Since , the effects of Basset, virtual mass, Magnus, Saffman and buoyancy forces are negligible (Arcoumanis et al., 1997;Lockwood & Naguib, 1975).In dilute flows of engineering interest, the steady-state drag term is the most important force acting on the particle.Under these conditions the simplified particle momentum equation is: The mathematical expression for the relaxation time,  p , is The critical issues are to determine the instantaneous fluid velocity and the evaluation of the time, t, of interaction of a particle with a particular eddy.
The time step is obviously the eddy-particle interaction time, which is the minimum of the eddy lifetime, FL  , and the eddy transit time, t c .The eddy lifetime is estimated assuming that the characteristic size of an eddy is the dissipation length scale in isotropic flow: where A and B are two dependent constants (Shirolkar et al., 1996).
The transit time, t c , is the minimum time a particle would take to cross an eddy with characteristic dimension, l e , and is given by where d v  is the relative velocity between the particle and the fluid (drift velocity).
A different expression for the transit time is also recommended in the literature (Shitolkar et al., 1996;Shuen et al., 1983;Gosman & Ioannides, 1981), and was used in the present work: where the drift velocity is also estimated at the beginning of a new iteration.

This equation has no solution when
;; e pf i p i lu u   , that is, when the linearized stopping distance of the particle is smaller than the eddy size.In such a case, the particle can be assumed to be trapped by the eddy, and the interaction time will be the eddy lifetime.
The instantaneous velocity at the start of a particle-eddy interaction is obtained by random sampling from an isotropic Gaussian pdf having standard deviations of 2/3k and zero mean values.
The above isotropic model was extended in the present work to account for crosscorrelation's and anisotropy.To obtain the fluctuating velocities u' f and v' f , two fluctuating velocities u' 1 and u' 2 are sampled independently, and then are correlated using the correlation coefficient R uv : The interaction between the continuous and dispersed phase is introduced by treating particles as sources of mass, momentum and energy to the gaseous phase.The source terms due to the particles are calculated for each Eulerian cell of the continuous phase and are summarized in Table 4, and can be divided into two parts, which yields the following expression: where i S  specifies the source term due to inter-phase transport and m S  takes into consideration the transfer caused by evaporation.
To represent the temporal changes of droplet size and temperature Chen and Pereira (1992) used the following equations.In the last equation is assumed that the prevailing mode of heat transfer is forced convection, no evaporation occurs during the preheating period and the temperature is uniform across the droplet radius.For the forced convection the Ranz and Marshall (1952) correlation has taken the place of the Nusselt Number.
The solution of the governing equations was obtained using a finite-difference method that used discretized algebraic equations deduced from the exact differential equations that they represent.In order to reduce the numerical diffusion errors to an acceptable level, the quadratic upstream-weighted interpolation scheme was used (Leonard, 1979).Nevertheless, the usual grid independence tests were performed.The computational domain (see Fig. 1) has six boundaries where dependent values are specified: an inlet plane and outlet planes, a symmetry plane, and three solid walls at the top, bottom and side of the channel.At the inlet boundary, uniform profiles of all dependent variables are set, while at the outflow boundaries, the gradients of dependent variables in the perpendicular direction are set to zero.On the symmetry plane, the normal velocity vanishes, and the normal derivates of the other variables are zero.At the solid surfaces, the wall function method described in detail by Launder and Spalding (1974) is used to prescribe the boundary conditions for the velocity and turbulence quantities, assuming that the turbulence is in state of local equilibrium.
The cross section of the computational domain is 0.05 x 0.05 m, whilst the channel length is 0.273 m.The droplets injection is perpendicular to the crossflow and the location of the injection point is 0.023 m far from the inlet plane (Z in /H = 0.46).
The monosize array of droplets of 230m of diameter is injected with an initial velocity V p =-1m/s and a temperature of 293K or 443K through a crossflow with W c =10m/s.The wall temperatures are 800K.

Results
To assess the computational method two sources of experimental data were used for the case of a polydisperse spray with a co-axial flow at atmospheric pressure (Heitor & Moreira, 1994) or elevated pressures (Barros, 1997).The method yielded reasonable results and revealed capabilities to improve the knowledge of the particle dispersion phenomena in more complex configurations.An example of the results obtained is shown in Figure 2, and a more complete analysis can be found in previous publications (Barata et al., 1999;Barata & Silva, 2000).The method was then extended to the case of an evaporating spray in a crossflow, and the evaporation models used by Chen and Pereira (1992) and Sommerfeld (1998) were tested and compared with the present model (see Barata, 2005 for details).
Figure 3 presents a parallel projection of the droplet trajectories in the vertical plane of symmetry (X=0) for two volatile fuels: n-Heptane and Ethanol.The former is used to define the zero limit of the anti-knock (resistance to pre-ignition) quality of fuels, while the other can be used to increase the octane number of gasoline.The higher volatility level of Ethanol can be inferred from Fig. 3b) by the more uniform distribution at the right side of the domain and the trajectories in the direction of the top wall (at Y=0.05m).Due to the high volatility level of both fuels the droplets are injected and start almost immediately to evaporate, which gives rise to smaller droplets that follow quite closely the gaseous flow.Further downstream of the injection point, the trajectories of the droplets of Ethanol are more directed downwards than those of n-Heptane due to the higher fuel density and higher latent heat of vaporization.As a consequence, although a colder region near the injector is observed with Ethanol, the domain shows in general a much more uniform temperature distribution.Figure 4 shows the droplet temperature and diameter variation time for the different fuels and test conditions used in the present work that are summarized in Table 5.The horizontal part of the line of the temperature variation with time (Fig. 4a) reveals the equilibrium of the evaporation process that corresponds to the horizontal part of the droplet diameter variation with time.It should be pointed out that since the droplet is moving in the direction of the Z coordinate (with the crossflow), the ambient temperature may not be constant, and the evolution of the droplet diameter with time is also influenced by its velocity.The Rapeseed Methyl Ester (RME) and the Diesel Fuel (DF-2) have a higher boiling temperature and do not attain the equilibrium temperature in the first 50miliseconds (Fig. 5a).As a consequence of the "pre-heating period" (about 20miliseconds in Fig. 4b) the droplets diameters remain approximately constant and the evaporation starts later (Fig. 4b).
Figure 5 shows parallel projections of the droplet trajectories in the vertical plane of symmetry, and confirms the main evaporation characteristics of the DF-2 and RME described in the previous paragraph.The pattern is similar for DF-2 and RME, although in the latter case there is a higher concentration of droplets in the core of the deflected monosize spray.This result is consistent with the slightly poorer evaporation characteristics of the RME deducted from Fig. 5, and taking into account the average time that a droplet takes to reach the right hand side of the domain (at Z=0.3m) it is expected that its diameter would be (in average) at the exit of the channel about 92.2% of the initial diameter.So, in general the diameters of the droplets will be larger with DF-2 and RME, and the dispersion will be more difficult, because the slip between the gaseous phase and the dispersed phase will be more pronounced.Some collisions with the bottom wall are observed, but were not taken into account in the present study, although this phenomena has been investigated and reported elsewhere (see Barata & Silva, 2005).
Increasing the injection temperature of RME improves the evaporation of the droplets, and a more uniform distribution is obtained (Fig. 5c).As shown in Fig. 4, to obtain the equilibrium stage of evaporation near the injection point a pre-heating of 150K is required, which will be particularly difficult to implement in most of the practical situations.The evaporative characteristics of RME can be further analysed with the help of Fig. 6 that shows a three-dimensional perspective of the mass fraction distribution with and without additional pre-heating.The results obtained with DF-2 and RME without additional preheating (Figs.6a and b) show that the mass fraction of fuel is always less than 0.04.For the DF-2 there is a larger evaporation near the injector, but further downstream the RME gives the higher values.When the additional pre-heating of 150K is used with RME, the domain shows a large region with a concentration of fuel vapour greater than 0.06, and the resulting pattern is quite similar to those obtained with more volatile fuels such as n-Heptane or Ethanol.

Conclusion
An Eulerian/Lagragian approach has been presented to calculate evaporating sprays through a crossflow.A method developed to study isothermal turbulent two-and threedimensional dispersion was extended to the case of an array of evaporating biofuel droplets.
The droplet diameter, temperature and mass fraction distributions were found to be strongly dependent on the fuels used.Rapeseed Methyl Esters exhibit similar evaporating characteristics to DF-2, which indicates that it can be successfully used as an alternative fuel www.intechopen.com Numerical Simulation of Biofuels Injection 143 in many applications that utilize diffusion flames.The use of RME in homogeneous combustion systems may require a prohibitive level of pre-heating, and the use of Ethanol (obtained from sugar or starch crops) may be a better alternative.

Aknowledgements
The present work has been performed in the scope of the activities of the AFTUR project.
The financial support of the European Commission under Contract nº ENK5-CT-2002-00662 is gratefully acknowledged.
13)Equations 7 and 8 for m  are similar to the expressions for the droplet vaporization rate predicted by the classical model, with the values of the non-

Table 1 .
Characteristic properties of biofuels compared with n-Heptane and Diesel Fuel.
Terms in the general form of the differential equation.

Table 4 .
Dispersed phase source terms.

Table 5 .
Summary of test conditions.