**2. 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.

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), and revealed a general good agreement with

> Rapeseed Methyl Ester (RME)

687.8 Kg/m3 880 Kg/m3 846 Kg/m3 790 Kg/m3

371.4ºK 613ºK 536.4ºK 351ºK

371,8KJ/Kg 254KJ/Kg 254KJ/Kg 904KJ/Kg

Molecular mass, MF 100.16 300 198 46

Table 1. Characteristic properties of biofuels compared with n-Heptane and Diesel Fuel.

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

This section describes the mathematical model for turbulent particle dispersion and vaporization assuming that the particles are sufficiently dispersed so that particle-particle

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

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.

Diesel Fuel

(DF-2) Ethanol

experiments.

Fuel density at 288.6ºK, F288.6K

Latent heat of vaporization at atmospheric pressure, LTbn

Boiling temperature at atmospheric pressure, Tbn

fuel used was gasoline.

**2. Mathematical model** 

interaction is negligible.

and the material properties of the phases are constant.

n-Heptane

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

$$\frac{\partial \left(\rho \mathcal{U} J\_i \phi\right)}{\partial X\_i} = \frac{\partial}{\partial X\_i} \left[\Gamma\_\phi \frac{\partial \phi}{\partial X\_i}\right] + S\_\phi \tag{1}$$

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:

$$S\_{\phi} = S\_{\phi\chi} + S\_{\phi\eta} \tag{2}$$

where *S <sup>g</sup>* , specifies the source term of the gas and *S<sup>p</sup>* , specifies the source term of the particle.


Table 2. Terms in the general form of the differential equation.

Numerical Simulation of Biofuels Injection 107

\* ln 1 *vap p T*

*g g D* can be replaced with *K C vap pvap* , assuming a Lewis number of unity.

*vap s*

*Fs*

1

0.38

*Sh Sh F* \* 2 2/ <sup>0</sup> *<sup>M</sup>* (14)

*Nu Nu F* \* 2 2/ <sup>0</sup> *<sup>T</sup>* (15)

*Y*

*L s M*

> 1 *Fs F*

1 1 *<sup>A</sup>*

For any given value of surface temperature, the vapor pressure is readily estimated from the

exp <sup>43</sup> *Fs s*

*cr s s tbn*

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-dimensional parameters *Nu*<sup>0</sup> and *Sh*0 in the classical formulas substituted by \* *Nu* and \* *Sh* respectively. Where are

*T T LT L*

*<sup>b</sup> P a*

*P M <sup>Y</sup>*

*Fs F*

*P M* 

*T* 

*cr bn*

*T T* 

*Y Y <sup>B</sup>*

*M*

where *YFs* is the fuel mass fraction on the droplet surface and defined as:

*Fs*

The latent heat of vaporization is given by Watson (1931) as

*C TT Q m L T B* 

*ggp* \* ln 1 *<sup>M</sup>* (7)

(8)

(9)

(10)

(11)

(12)

(13)

*m D D Sh B*

*m D Nu B*

*pvap*

*K*

*<sup>C</sup>* 

The heat penetrating into the droplet can be expressed by:

The mass transfer number *BM* as defined as

Clausius-Claperyon equation as

expressed as

where a and b are constants of the fuel.

 

and

the quantity

The source terms of the gas phase, *S <sup>g</sup>* and the effective diffusion coefficient , are summarized in table 2 for different depended variables. G is the usual turbulence energy production term defined as:

$$\mathbf{G} = \mu\_t \left[ \frac{\partial \overline{\mathbf{U}\_i}}{\partial \mathbf{X}\_j} + \frac{\partial \overline{\mathbf{U}\_j}}{\partial \mathbf{X}\_i} \right] \frac{\partial \overline{\mathbf{U}\_i}}{\partial \mathbf{X}\_j} \tag{3}$$

and

$$
\mu\_t = \mathcal{C}\_{\mu} \rho \frac{k^2}{\varepsilon} \tag{4}
$$

The turbulence model constants that are used are those indicated by Launder and Spalding (1974) that have given good results for a large number of flows, and are summarized in the next table.


Table 3. Turbulence model constants.

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.

$$\frac{dD\_p}{dt} = -\frac{2\dot{m}}{\pi \rho\_F D\_p^2} \tag{5}$$

$$\frac{dT\_p}{dt} = \frac{6Q\_L}{\pi \rho\_F C\_{P\_p} D\_p^3} \tag{6}$$

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* :

$$\dot{m} = \pi \overline{\rho\_{\mathcal{S}} D\_{\mathcal{g}}} D\_{\mathcal{g}} Sh \, ^\* \ln \left( 1 + B\_M \right) \tag{7}$$

and

106 Fuel Injection in Automotive Engineering

summarized in table 2 for different depended variables. G is the usual turbulence energy

*<sup>j</sup> i i <sup>t</sup> j i j U U U*

 

*X XX*

2

*k*

 *C* 

The turbulence model constants that are used are those indicated by Launder and Spalding (1974) that have given good results for a large number of flows, and are summarized in the

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

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 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

*<sup>p</sup>* 2

*dD m*

*dt*

*dt*

2

3

*F p*

*F*

*F P p*

6

*p L*

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

*dT Q*

*<sup>g</sup>* and the effective diffusion coefficient

, are

(3)

(4)

<sup>3</sup> Pr*<sup>t</sup> Sct* Pr *Sc*

*<sup>D</sup>* (5)

*C D* (6)

*C k <sup>P</sup>* / *<sup>g</sup>*

 / *D*

*t*

 *C*

0.09 1.44 1.92 1.0 1.3 1.1 0.6 0.85

effects are taken into account by introducing empirical correlation laws.

the droplet; and liquid/vapor thermal equilibrium on the droplet surface.

*G*

The source terms of the gas phase, *S*

production term defined as:

<sup>1</sup> *C*

2

Table 3. Turbulence model constants.

*k*

and

next table.

*CC*

solved.

the evaporation rate *m* :

$$\dot{m} = \pi \frac{\overline{K\_{\text{vap}}}}{\overline{C\_{\text{prop}}}} D\_p N \nu^\* \ln\left(1 + B\_T\right) \tag{8}$$

the quantity *g g D* can be replaced with *K C vap pvap* , assuming a Lewis number of unity. The heat penetrating into the droplet can be expressed by:

$$Q\_L = \dot{m} \left( \frac{\overline{C\_{vap}} \left( T\_{\infty} - T\_s \right)}{B\_M} - L \left( T\_s \right) \right) \tag{9}$$

The mass transfer number *BM* as defined as

$$B\_M = \frac{Y\_{\rm Fs} - Y\_{\rm Fos}}{1 - Y\_{\rm Fs}} \tag{10}$$

where *YFs* is the fuel mass fraction on the droplet surface and defined as:

$$Y\_{Fs} = \left[1 + \left(\frac{P}{P\_{Fs}} - 1\right) \frac{M\_A}{M\_F}\right]^{-1} \tag{11}$$

For any given value of surface temperature, the vapor pressure is readily estimated from the Clausius-Claperyon equation as

$$P\_{Fs} = \exp\left(a - \frac{b}{T\_s - 4\Im}\right) \tag{12}$$

where a and b are constants of the fuel.

The latent heat of vaporization is given by Watson (1931) as

$$L\left(T\_s\right) = L\_{thn} \left(\frac{T\_{cr} - T\_s}{T\_{cr} - T\_{bn}}\right)^{-0.38} \tag{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-dimensional parameters *Nu*<sup>0</sup> and *Sh*0 in the classical formulas substituted by \* *Nu* and \* *Sh* respectively. Where are expressed as

$$\text{Sls}^\* = \text{2} + \left(\text{Sls}\_0 - \text{2}\right) / \text{F}\_{\text{M}} \tag{14}$$

$$\mathcal{N}\Lambda\iota^\* = \mathcal{Z} + \left(\mathcal{N}\iota\_0 - \mathcal{Z}\right) / F\_T \tag{15}$$

Numerical Simulation of Biofuels Injection 109

*<sup>p</sup> <sup>g</sup> <sup>M</sup> L T*

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

> 3 *<sup>s</sup> <sup>r</sup> <sup>s</sup>*

> > 3 *F Fs*

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

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

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:

> ; ; ;

*u u <sup>g</sup> <sup>t</sup>*

<sup>2</sup> 24 18 Re *p p*

Re *fp f p*

*f D p D C* 

*f*

 *V VD* 

*p*

*f i p i i*

*<sup>p</sup>*, is

*p i* 1

*p*

*p*

*u*

field, for a sufficiently high number of trajectories to provide a representative statistics.

0.276Re Pr <sup>1</sup>

*dt C D pvap F* 

*Fr Fs*

For example, the reference specific heat at constant pressure is obtained as

*dD K B p vap M*

 

*Cv D C F*

*F*

*K B*

12 ln 1 

2

*dt dT*

and composition:

be neglected. Since

 *p* 

The mathematical expression for the relaxation time,

where Re*p* is the particle Reynolds number,

*ap F P P M*

 

<sup>1</sup> <sup>1</sup> 4 ln 1 <sup>2</sup> <sup>3</sup> 0.276Re Pr <sup>1</sup>

*F P M*

*M vap s*

*<sup>T</sup> <sup>T</sup> <sup>T</sup> <sup>T</sup>* (25)

*<sup>Y</sup> <sup>Y</sup> <sup>Y</sup> <sup>Y</sup>* (26)

*<sup>p</sup> Ar <sup>p</sup> <sup>r</sup> Fr <sup>p</sup> Tr C Y C at T Y C at vap Ar <sup>F</sup>* (27)

*<sup>f</sup>* 200 , the effects of Basset, virtual mass, Magnus, Saffman and

(28)

(29)

(30)

*B C T T*

 

 

3 1 2 1

(23)

(24)

 

*s*

The parameters \* *Nu* and \* *Sh* are the "modified" Nusselt and Sherwood numbers, and tend to *Nu*0 and *Sh*<sup>0</sup> , respectively, as *FT* and *FM* 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 *FM* and *FT* do not depend on the local Reynolds number. It was found that the values *FM* and *FT* are practically insensitive to the Schmidt and Prandtl numbers and the wedge angle variations, and can be approximated as

$$F\_M = F\left(B\_M\right)\_\prime \; F\_T = F\left(B\_T\right) \tag{16}$$

where *F B* is given by

$$F(B) = \left(1 + B\right)^{0.7} \frac{\ln\left(1 + B\right)}{B} \tag{17}$$

*Nu*0 and *Sh*0 are evaluated by the Frossling correlations:

$$Nu\_0 = 2 + 0.552 \,\mathrm{Re}^{\prime \text{J}} \,\mathrm{Pr}^{\prime \text{J}} \tag{18}$$

$$\text{Sh}\_0 = 2 + 0.552 \,\text{Re}^2 \,\text{Sc}^2 \,\text{Sc}^3 \tag{19}$$

The evaporation rate *m* with convection is:

$$\dot{m} = \pi \overline{\rho\_g D\_g} D\_p \ln\left(1 + B\_M\right) \left(2 + \frac{0.552 \operatorname{Re}^{\lambda\_2} \mathcal{S} c^{\lambda\_3}}{F\_M}\right) \tag{20}$$

and

$$\dot{m} = \pi \frac{\overline{K\_{vap}}}{\overline{C\_{p\,vap}}} D\_p \ln(1 + B\_T) \left( 2 + \frac{0.552 \,\mathrm{Re}^{\prime\prime} \mathrm{Pr}^{\prime\prime} \mathcal{N}^{\prime}}{F\_T} \right) \tag{21}$$

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:

$$\dot{m} = 2\pi \frac{\overline{K\_{vap}}}{\overline{C\_{pvap}}} D\_p \ln\left(1 + B\_M\right) \left(1 + \frac{0.276 \,\mathrm{Re}^{\frac{1}{2}} \sqrt{2} \,\mathrm{Pr}^{\frac{1}{3}}}{F\_M}\right) \tag{22}$$

And the equations for the temporal changes of droplet size and temperature are:

The parameters \* *Nu* and \* *Sh* are the "modified" Nusselt and Sherwood numbers, and

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 *FM* and *FT* do not depend on the local Reynolds number. It was found that the values *FM* and *FT* are practically insensitive to the Schmidt and Prandtl numbers and the wedge angle variations, and can be

0.7 ln 1 <sup>1</sup>

*C F*

And the equations for the temporal changes of droplet size and temperature are:

<sup>2</sup> <sup>3</sup> 0.276Re Pr <sup>2</sup> ln 1 1

*p M FM*

*pvap T*

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.

*p T*

*Sc m DD B*

<sup>2</sup> <sup>3</sup> 0.552Re ln 1 2 *ggp M*

<sup>2</sup> <sup>3</sup> 0.552Re Pr ln 1 2 *vap*

*FB B*

*Nu*0 and *Sh*0 are evaluated by the Frossling correlations:

 

*K m DB*

*Kvap m DB Cpvap*

The evaporation rate *m* with convection is:

Finally the evaporation rate *m* is:

*F FB F FB <sup>M</sup> MT T* , (16)

(17)

<sup>0</sup> *Nu* 2 0.552Re Pr (18)

<sup>0</sup> *Sh* 2 0.552Re *Sc* (19)

1 1

1 1

1 1

 

(22)

*M*

*F*

 (20)

 (21)

*B*

*B* 

1 1 2 3

1 1 2 3

tend to *Nu*0 and *Sh*<sup>0</sup> , respectively, as *FT* and *FM* tend to the unity.

approximated as

and

where *F B* is given by

$$\frac{dD\_P}{dt} = -\frac{4K\_{vap}\ln\left(1 + B\_M\right)}{C\_{pvap}\rho\_F D\_P} \left(1 + \frac{0.276\,\text{Re}^2\text{/}2\text{ Pr}^2\text{/}3}{F\_M}\right) \tag{23}$$

$$\frac{dT\_p}{dt} = \frac{12K\_g \ln\left(1 + B\_M\right)}{\text{Cov}\_{ap}\rho\_F D\_P^2 C\_{P\_\varepsilon}} \left(1 + \frac{0.276 \,\text{Re}^{\frac{1}{2}} \,\text{Pr}^{\frac{1}{2}}}{F\_M}\right) \left(\frac{C\_{\text{vap}}\left(T\_w - T\_s\right)}{B\_M} - L\left(T\_s\right)\right) \tag{24}$$

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:

$$T\_r = T\_s + \frac{T\_w - T\_s}{3} \tag{25}$$

$$Y\_{Fr} = Y\_{Fs} + \frac{Y\_{Fs} - Y\_{Fs}}{3} \tag{26}$$

For example, the reference specific heat at constant pressure is obtained as

$$\mathbf{C}\_{p\_{np}} = Y\_{Ar} \begin{pmatrix} \mathbf{C}\_{p\_{Ar}} & at & T\_r \end{pmatrix} + Y\_{Fr} \begin{Bmatrix} \mathbf{C}\_{p\_r} & at & T\_r \end{Bmatrix} \tag{27}$$

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 *p <sup>f</sup>* 200 , 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:

$$\frac{\partial \boldsymbol{u}\_{p;i}}{\partial t} = \frac{1}{\tau\_p} (\boldsymbol{u}\_{f;i} - \boldsymbol{u}\_{p;i}) + \boldsymbol{g}\_i \tag{28}$$

The mathematical expression for the relaxation time, *<sup>p</sup>*, is

$$\tau\_p = \frac{24\rho\_p D\_p^2}{18\mu\_f C\_D \,\mathrm{Re}\_p} \tag{29}$$

where Re*p* is the particle Reynolds number,

$$\text{Re}\_p = \frac{\rho\_f \left| \overrightarrow{V\_p} - \overrightarrow{V\_f} \right| D\_p}{\mu\_f} \tag{30}$$

Numerical Simulation of Biofuels Injection 111

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/3*k* and zero

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

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

> *S SS p*

*<sup>i</sup>* specifies the source term due to inter-phase transport and *S*

To represent the temporal changes of droplet size and temperature Chen and Pereira (1992)

*vap*

<sup>2</sup> 12 \* \* 1 0.3Re \* Pr *F*

ln 1 ( 1 0.23Re ( )

1 1

*T T*

*s*

ln 1 *<sup>e</sup> c p*

*<sup>l</sup> <sup>t</sup>*

where the drift velocity is also estimated at the beginning of a new iteration.

was obtained from the measurements.

This equation has no solution when *<sup>e</sup> p f* ; ; *<sup>i</sup> <sup>p</sup> <sup>i</sup> l uu*

mean values.

where

*uv*

*R*

expression:

where *S*

used the following equations.

correlation coefficient *Ruv*:

2 2 ' '

*f f*

consideration the transfer caused by evaporation.

*dt*

4

*vap p vap p*

*dD K C*

*p s*

*D C*

*FPP*

*dt C L T*

*p s g*

*dT T T K*

' ' *f f*

*u v*

*u v*

; ;

2

(39)

, that is, when the linearized stopping

<sup>1</sup> ' ' *u u <sup>f</sup>* (40)

*i m* (42)

1 2

2 3

(43)

(44)

*<sup>m</sup>* takes into

1 2 ' '1 ' *f uv uv v Ru Ru* (41)

*p f i p i*

*u u*

 

Note that the physical properties of *<sup>f</sup>* and *<sup>f</sup>* should be evaluated at the reference temperature *Tr* and are

$$
\mu\_f = \mathbf{Y}\_{Ar} \left( \mu\_A \text{ at } T\_r \right) + \mathbf{Y}\_{Fr} \left( \mu\_F \text{ at } T\_r \right) \tag{31}
$$

$$
\rho\_{\text{vap}} = \left(\frac{Y\_{Ar}}{\rho\_A} + \frac{Y\_{Fr}}{\rho\_F}\right)^{-1} \tag{32}
$$

and *CD* is the drag coefficient (Shirolkar et al., 1996),

$$\mathbf{C}\_{D} = \left(\frac{24}{\mathrm{Re}\_{p}}\right) \left(1 + 0.15 \,\mathrm{Re}\_{p}^{0.687}\right) \tag{33}$$

for Re*p*<103.

The particle momentum equation can be analytically solved over small time steps, *t* , and the particle trajectory is given by

$$\boldsymbol{\mu}\_{p;i}^{\rm NENV} = \boldsymbol{\mu}\_{f;i} + \left(\boldsymbol{\mu}\_{p;i}^{\rm OLD} - \boldsymbol{\mu}\_{f;i}\right) \boldsymbol{e}^{-\Lambda t/\tau\_p} + \boldsymbol{g}\_i \tau\_p \left[\mathbf{1} - \boldsymbol{e}^{-\Lambda t/\tau\_p}\right] \tag{34}$$

$$\mathbf{x}\_{p;i}^{NEW} = \mathbf{x}\_{p;i}^{OLD} + \frac{\Delta t}{2} \left( \mathbf{u}\_{p;i}^{NEW} + \mathbf{u}\_{p;i}^{OLD} \right) \tag{35}$$

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, *tc*. The eddy lifetime is estimated assuming that the characteristic size of an eddy is the dissipation length scale in isotropic flow:

$$l\_{\varepsilon} = B \frac{k^{3/2}}{\varepsilon} \approx C\_{\mu}^{3/4} \frac{k^{3/2}}{\varepsilon} \tag{36}$$

$$
\tau\_{\rm FL} = A \frac{k}{\varepsilon} \approx 0.2 \frac{k}{\varepsilon} \tag{37}
$$

where *A* and *B* are two dependent constants (Shirolkar et al., 1996).

The transit time, *tc*, is the minimum time a particle would take to cross an eddy with characteristic dimension, *le*, and is given by

$$t\_c = \frac{I\_c}{\left| \overline{\overline{v\_d}} \right|} \tag{38}$$

where *<sup>d</sup> 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:

*Ar Fr*

 

*A F Y Y*

 

<sup>24</sup> 0.687 1 0.15Re

/ /

, and the eddy transit time, *tc*. The eddy lifetime is estimated assuming

(34)

Re *D p p*

The particle momentum equation can be analytically solved over small time steps, *t* , and

; ;; ; <sup>1</sup> *p p NEW OLD t t u u u ue g e pi f i pi f i i p*

;; ;; <sup>2</sup> *NEW OLD NEW OLD pi pi pi pi t xx uu* 

The critical issues are to determine the instantaneous fluid velocity and the evaluation of the

The time step is obviously the eddy-particle interaction time, which is the minimum of the

3/2 3/2 3/4

that the characteristic size of an eddy is the dissipation length scale in isotropic flow:

*k k lB C*

0.2 *FL k k*

The transit time, *tc*, is the minimum time a particle would take to cross an eddy with

*<sup>e</sup> <sup>c</sup> d*

*<sup>l</sup> <sup>t</sup>*

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:

*e*

where *A* and *B* are two dependent constants (Shirolkar et al., 1996).

*A*  

1

*<sup>f</sup> Y at T Y at T Ar A r Fr F r* (31)

*<sup>f</sup>* should be evaluated at the reference

(32)

(33)

 

(35)

(36)

(37)

*<sup>v</sup>* (38)

*<sup>f</sup>* and

*vap*

 

*C*

time, *t*, of interaction of a particle with a particular eddy.

Note that the physical properties of

the particle trajectory is given by

and *CD* is the drag coefficient (Shirolkar et al., 1996),

temperature *Tr* and are

for Re*p*<103.

eddy lifetime, *FL*

where *<sup>d</sup> v*

characteristic dimension, *le*, and is given by

$$t\_c = -\tau\_p \ln \left( 1 - \frac{l\_e}{\tau\_p \left| u\_{f;i} - u\_{p;i} \right|} \right) \tag{39}$$

where the drift velocity is also estimated at the beginning of a new iteration.

This equation has no solution when *<sup>e</sup> p f* ; ; *<sup>i</sup> <sup>p</sup> <sup>i</sup> l uu* , 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/3*k* 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 *Ruv*:

$$\mathfrak{u}'\_f = \mathfrak{u}'\_1 \tag{40}$$

$$
\boldsymbol{\upsilon}^{\prime}\_{f} = \boldsymbol{R}\_{\boldsymbol{uv}} \boldsymbol{\mu}^{\prime}\_{1} + \sqrt{1 - \boldsymbol{R}\_{\boldsymbol{uv}}^{2}} \boldsymbol{\mu}^{\prime}\_{2} \tag{41}
$$

where 2 2 ' ' ' ' *f f uv f f u v R u v* was obtained from the measurements.

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:

$$S\_{\phi p} = S\_{\phi i} + S\_{\phi m} \tag{42}$$

where *S<sup>i</sup>* specifies the source term due to inter-phase transport and *S<sup>m</sup>* 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.

$$\frac{dD\_p}{dt} = -\frac{4K\_{vap}}{C\_{p\_{vap}}} \ln\left(1 + \frac{C\_{p\_{vap}}}{L(T\_s)}(T\_o - T\_s)\left(1 + 0.23\,\text{Re}^2\right)\right) \tag{43}$$

$$\frac{dT\_p}{dt} = 12K\_g \ast \left(\frac{(T\_\infty - T\_s)}{\rho\_F D\_p^2 C\_{p\_\mathbb{F}}}\right) \ast \left(1 + 0.3 \text{Re}^{\mathcal{Y}\_2} \text{Pr}^{\mathcal{Y}\_3}\right) \tag{44}$$

Numerical Simulation of Biofuels Injection 113

The monosize array of droplets of 230m of diameter is injected with an initial velocity Vp=- 1m/s and a temperature of 293K or 443K through a crossflow with Wc=10m/s. The wall

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

temperatures are 800K.

Fig. 1. Flow configuration

temperature distribution.

**3. Results** 

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.


Table 4. Dispersed phase source terms.

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 (Zin/H = 0.46).

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)

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,

*<sup>i</sup> S*

*C S* <sup>3</sup> *km <sup>k</sup>*

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 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

,

*<sup>p</sup> i j <sup>P</sup>*

*V C*

1 2 *UUS j j <sup>m</sup>*

*m*

*p p <sup>p</sup> <sup>i</sup> <sup>j</sup> m N <sup>V</sup>*

, *p p ia <sup>p</sup> i j mNu <sup>V</sup>*

\* *vap*

*p p <sup>p</sup> <sup>i</sup> <sup>j</sup> m N <sup>V</sup>*

1 2 *US U S U U S <sup>j</sup> Ujm <sup>j</sup> Ujm j j <sup>m</sup>*

*p p P s rs*

*m N C T TT*

*A*

 

correlation has taken the place of the Nusselt Number.

the usual grid independence tests were performed.

*p p tt t jp jp i <sup>p</sup> i j*

, *<sup>A</sup> p tbn p L*

*SY*1,*<sup>p</sup>* 0 0

*SY*2,*<sup>p</sup>* 0 ,

*<sup>p</sup> ij P N Lm Q V C* 

,*<sup>p</sup>* 0 ,

*u u g t*

*<sup>p</sup> S*

, *Ui S <sup>p</sup>* , , ,

*Sk*,*<sup>p</sup> US U S <sup>j</sup> Uji <sup>j</sup> Uji*

,*<sup>p</sup> C S* <sup>3</sup> *ki <sup>k</sup>*

Table 4. Dispersed phase source terms.

the turbulence is in state of local equilibrium.

injection point is 0.023 m far from the inlet plane (Zin/H = 0.46).

*m N*

*V*

*S*

*S*

*ST* ,*<sup>p</sup>*

*S* The monosize array of droplets of 230m of diameter is injected with an initial velocity Vp=- 1m/s and a temperature of 293K or 443K through a crossflow with Wc=10m/s. The wall temperatures are 800K.

Fig. 1. Flow configuration
