Baumgarten

134 4 Modeling Spray and Mixture Formation

Fig. 4.31. Measured and predicted penetration. a case 1, b case 2 [125]

Table 4.3. Injection parameters for sprays shown in Figs. 4.30 and 4.31, data from [125]

Injection parameter	Quantity
	Case 1	Case 2
Spray half-cone angle Τ [°]	23	23
Dispersion angle Τ [°]	10	10
Fuel delivery per injection	56.8	56.8
[mm3]
Nozzle Diameter [µm]	560	560
Injection pressure [MPa]	4.76	6.12
Injection duration [ms]	3.86	3.4
Liquid density [g/cm3]	0.77	0.77

Senecal et al. [129] have used the LISA-TAB model to predict the spray formation in the case of outwardly opening nozzles. Due to the pintle, these injectors do not produce a pre-spray. Furthermore, one less equation is needed to initialize the calculation, since the initial thickness of the sheet at the nozzle orifice is determined by the size of the annular gap. The injection velocity is calculated from the conservation of mass.

Stiesch et al. [135] also combined the LISA sheet atomization model with the TAB model in order to simulate the non-evaporating spray development in a pressure bomb and to validate the model at ambient as well as elevated gas pressures and densities. The authors used a pressure-swirl injector with a half-cone angle of 27 degrees, a dispersion angle of 10 degrees, and a nozzle diameter of 0.458 mm. The injection pressure was 4.93 MPa. Figs. 4.32 and 4.33 show a comparison of computational and experimental results for three different timings after start of injection and for two different backpressures. While the photographs are side views of the complete spray, the calculated spray images represent spray slices containing the spray axis. Both series of pictures suggest that the hybrid LISA-TAB model is capable of predicting the spray behavior very well. Especially the recirculating vortex, which starts to form at the spray edges, is accurately predicted,

4.3 Combined Models 135

Fig. 4.32. Measured and calculated spray images, gas pressure: 101 kPa, [135]

Fig. 4.33. Measured and calculated spray images, gas pressure: 366 kPa, [135]

and also the influence of an increase of backpressure (reduced penetration, narrower cone angle and a more distinct vortex at the spray edge) is predicted correctly.

4.3.5 LISA-DDB Model (Hollow-Cone Sprays)

Instead of using the TAB model, Park et al. [104] have combined the LISA sheet atomization model with the Droplet Deformation and Break-up (DDB) model in

136 4 Modeling Spray and Mixture Formation

order to calculate the disintegration process of a hollow-cone spray. The authors have validated the model against experimental data concerning spray penetration, droplet sizes, spray shape, and mean velocity distribution in the case of atmospheric sprays and injection pressures of 5 and 7 MPa. Because the break-up mechanism of the DDB model is very similar to that of the TAB model, there are no significant differences regarding the prediction accuracy of both hybrid models.

4.4 Droplet Drag Modeling

4.4.1 Spherical Drops

The relative velocity between gas and droplet results in a deceleration of the liquid and in an acceleration of the gas phase due to the exchange of momentum. The equation of motion of a spherical drop with radius r moving with a velocity urel relative to the gas is

	Ug	2		4		3	d 2 x		(4.129)
Fdrag			Ul 3 Sr				dt2 ,
	2 urel cD Af

where cD is the drag coefficient, Af = Σr2 is the frontal area of the spherical drop, and x is the coordinate along the droplet trajectory. The drag coefficient is given by that of a rigid sphere [6],

	­	24	§
	°		1
	cD,sphere ®Reg		¨
			©

°¯0.424

where Reg = 2r·urel Υg / µg.

Reg2 / 3 ·

6 ¸¹

½

Re d 000°		,	(4.130)
	¾
Re ! 000	°
	¿

4.4.2 Dynamic Drag Modeling

When a liquid drop enters a gas stream with a sufficiently large Weber number, it deforms and is no longer spherical as it interacts with the gas. This has been observed experimentally by many researchers, e.g. [80, 111]. Hence, the drag coefficient should be a function of its Reynolds number as accounted for in Eq. 4.130 and of its oscillation amplitude as well. Based on these observations, Liu et al. [79] use the TAB model, Fig. 4.34, in order to predict the droplet distortion y (Eq. 4.90) and then modify the drag coefficient by relating it empirically to the magnitude of the drop deformation. Since the drag coefficient of a distorting drop should lie between that of a rigid sphere, Eq. 4.130 (lower limit), and that of a disk, which is about 3.6 times higher (upper limit), a simple linear expression is used for the dynamic drag coefficient:

4.4 Droplet Drag Modeling	137

Fig. 4.34. The dynamic drag model accounts for drop distortion

cD cD,sphere 1 2.632 y .

(4.131)

Hence, the drag coefficients of a spherical drop and a disk are recovered in the undeformed (y = 0) and the fully deformed limits (y = 1).

Compared to the use of the standard rigid sphere drag coefficients, the consideration of dynamic drag results in increased overall drag coefficients and a more realistic calculation of droplet deceleration. Compared to the standard model, the increased droplet deceleration results in reduced penetration and larger droplets due to the shorter time span with high relative velocity between gas and liquid.

It must be pointed out that the TAB model is only used for the calculation of drop distortion and does not influence the break-up modeling. If, for example, the Kelvin-Helmholtz model is used in order to predict drop break-up, then y > 1 does not result in any break-up.

Liu et al. [81] have found that their TAB drop drag model [79] significantly underestimates drop drag effects for high-speed drops and have presented a modification based on the DDB model [59]. According to the DDB model, the drops are deformed due to a pure extensional flow from an initial sphere shape to an oblate spheroid of ellipsoidal cross section with major semi axis a and minor semi axis b, Fig. 4.35. The dynamic drag coefficient is again calculated using Eqs. 4.130 and 4.131, but the dimensionless deformation of the droplet surface is given by

­	§ a		·½
y min®1,	¨		1¸¾	,	(4.132)
¯	© r		¹¿

such that the drag coefficients of a spherical drop and a disk were also recovered in the undeformed and the fully deformed limits (a = 2r or y = 1, [59]), respectively. The drops undergo significant flattening, which changes the frontal area exposed to the airflow. Because the flattening occurs prior to a significant mass loss due to break-up [81], the calculation of the frontal area is also modified and given as Af = Σa2 instead of Af = Σr2.

The authors have extended their model in order to describe the dynamic drag after wall impingement. Before wall impingement, the drop’s frontal area normal to the airflow is the largest cross-sectional area of the deforming drop (Af = Σa2).

138 4 Modeling Spray and Mixture Formation

After wall impingement, drops with Wein < 80 (Wein: Weber number of incoming drops before wall impingement, Wein = 2rΥlvn /ς, vn: drop velocity component normal to the wall) bounce off the wall, and it is assumed that they have an initial rotation which may be maintained or enhanced by the strong turbulence in the combustion chamber. In this case, the largest cross-sectional area is not always normal to the direction of drop motion, and it is assumed that the drops rotate stochastically. A random number [ from the interval [0, 1] is used to modify the drop’s frontal area

Af * Sa b [ a b	(4.133)
and drag coefficient
c*D cD,sphere 1 2.632 [ y	(4.134)

after they have impacted and rebounded from the wall. After wall impingement, the drops are assumed to start with an initial deformation y = 1, which can decrease again if the droplet velocity slows down and reduces the aerodynamic forces. Because the droplet is rotating, the value of [ cannot be constant for the rest of its lifetime. Hence, [ is only calculated for a time interval equal to that of the droplet’s characteristic oscillation time

Ω rot B	Υl r 3	(4.135)
	,
	ς

where r is the undisturbed droplet radius and B is a constant. Thereafter, a new value of [ is sampled. The drop rotation effect after wall impingement is most important for large drops with high drop-gas relative velocities, but the Weber numbers (Wein) of which are less than 80.

Fig. 4.35. Schematic diagram of the deforming half drop (DDB model)

4.5 Evaporation	139

4.5 Evaporation

In addition to the break-up of the spray and the mixing processes of air and fuel droplets, the evaporation of liquid droplets also has a significant influence on ignition, combustion, and formation of pollutants. The formation of fuel vapor due to evaporation is a prerequisite for the subsequent chemical reactions. The evaporation process determines the spatial distribution of the equivalence ratio, and thus strongly affects the timing and location of ignition. The energy for evaporation is transferred from the combustion chamber gas to the colder droplet due to conductive, convective, and radiative heat transfer, resulting in diffusive and convective mass transfer of fuel vapor from the boundary layer at the drop surface into the gas, Fig. 4.36. This again affects temperature, velocity, and vapor concentration in the gas phase. Hence, there is a strong linking of evaporation rate and gas conditions, and, for this reason, there must always be a combined calculation of heat and mass transfer processes.

In order to describe the evaporation process mathematically, the following assumptions are usually made: the radiative heat transfer is neglected because it is small compared to the convective one. Because it is not feasible to resolve the flow field around all the droplets of a spray, the evaporation modeling is based on averaged flow conditions and average transfer coefficients around the droplets. The droplets are usually assumed to be of spherical shape. Deformation, break-up, collisions, and other interactions of droplets are neglected during the calculation of evaporation. Further on, the droplet’s interior is usually assumed to be well mixed. For this reason, there are no spatial gradients of the relevant quantities like liquid temperature, concentration of fuel components, boiling temperatures, and critical temperatures, heat of evaporation etc. inside the droplet, and only a dependence on time is possible. Furthermore, the solubility of the surrounding gas in the liquid and the effect of surface tension on the vapor pressure are neglected.

Fig. 4.36. Schematic view of drop vaporization

140 4 Modeling Spray and Mixture Formation

In order to determine the transport processes at the gas/liquid interface (mass and energy fluxes), phase equilibrium is assumed. It is presumed that the phase transition (liquid to vapor) is much faster than the vapor transport from the surface into the surrounding gas. Further on it is assumed that even if the conditions in the gas phase or inside the droplet change (e.g. temperature rise), phase equilibrium is always immediately reached.

The concentration of fuel vapor and thus also the properties of the gas mixture in the boundary layer are strongly dependent on radius, Fig. 4.36. In order to get representative values for the calculation of the diffusive mass transport, simplified vapor concentration curves are customarily used.

4.5.1 Evaporation of Single-Component Droplets

Although real fuels consist of a multitude of different components that influence the evaporation process (more volatile but less ignitable components evaporate first, components with higher molecular weight evaporate later), the standard approach today is to use a single-component model fuel. Usually tetradecane (n- C14H30) is used in order to represent the relevant properties of diesel, and octane is used for gasoline.

The temperature change of the liquid droplet can be obtained from an energy balance. The total heat flux

Q	Q	Q			,		(4.136)
drop	heating		evap
Q	m c	p,l	dTdrop			,	(4.137)
heating	drop		dt
Q	'h	dmevap			,		(4.138)
evap	evap		dt

transferred from the hot gas to the liquid droplet results in an increase of droplet

temperature (heating) and in evaporation. In Eqs. 4.136–4.138, mdrop and Tdrop are the droplet mass and temperature, cp,l is the specific heat capacity of the liquid

fuel, hevap is the enthalpy of evaporation, and mevap is the mass that evaporates in the time interval dt. Using Eqs. 4.137 and 4.138 in Eq. 4.136 and solving for the

temperature change yields

dTdrop	1			§ dQdrop		'hevap	dmevap ·
				¨				¸.	(4.139)
dt		m c			dt		dt
			p,l ©					¹
		drop

In order to solve Eq. 4.139, expressions for the quantities dQdrop /dt and dmevap /dt must be derived. A part of the total convective heat transferred from the gas to the

liquid due to the temperature difference T = Tdrop – Tφ (Tφ: temperature of the surrounding gas outside the boundary layer) is needed in order to heat up the

evaporated mass of fuel that is transported from the drop surface (Tdrop = TR) into

4.5 Evaporation	141

the gas (Tφ), see Fig. 4.37. A differential element of the gas atmosphere around the droplet (somewhere between r = R and r = φ) has to be regarded, and an energy balance considering both kinds of heat fluxes, Fig. 4.37, must be written down and solved. After several mathematical integrations and transformations, an equation for the heat transferred to the drop is obtained:

Q

Ο

Σ d

T

T

9

Nu ,

(4.140)

e9 1

drop

g

drop

φ

R

where

Nu

ddrop

(4.141)

Οg

and

9

mevapcp,vap

.

(4.142)

NuΟgΣ ddrop

These equations are derived and discussed in [32] and [33] and are finally utilized in the droplet evaporation model of Borman and Johnson [16], which again is used in most spray and combustion models. The variable 9 is a dimensionless correction factor taking account of the reduced heat transfer due to the simultaneous mass transfer from the drop into the gas atmosphere. In Eqs. 4.140–4.142, Οg is the thermal conductivity of the gas mixture at the drop surface, cp,vap is the specific heat capacity of the fuel vapor (both calculated using the 1/3rd rule of Hubbard et

al. [54], see also Eqs. 4.190 and 4.191, and ddrop is the droplet diameter. The appropriate Nusselt number that includes the effect of a relative velocity between

droplet and gas has been proposed by Ranz and Marshall [112],

Nu	2.0 0.6 Re1/ 2 Pr1/ 3 ,	(4.143)
Re	Ufurel ddrop / µg ,	(4.144)
	Pr µg cp / Og .	(4.145)

Fig. 4.37. Definition of coordinate system and heat fluxes

142 4 Modeling Spray and Mixture Formation

The properties of the gas phase inside the boundary layer are again calculated using the 1/3rd rule.

Next, an expression for the fuel vapor mass flow due to evaporation in Eq. 4.139 must be derived. For the case of a one-dimensional binary diffusion, Fick’s law can be used on a mass basis [142],

mA	Y	§	m	·	Υ D		dYA	,	(4.146)
4Σ r2		A ¨	4Σ r2	¸		AB	dr
		©		¹

where mA is the mass flow of species A, m = mA + mB is the total mixture flow rate, YA is the mass fraction, and DAB is the binary diffusivity. The first term on the right hand side is the mass flow of species A associated with bulk flow, and the second one is the mass flow of species A associated with molecular diffusion. If it is assumed that species B (gas around the droplet) is insoluble in the liquid A, and that outside the boundary layer the concentration of A is constant, there will be no

net transport of B in the boundary layer gas, and							mB = 0. Thus, Eq. 4.146 can be
expressed as	mA dr				dYA .
							(4.147)
	4ΣΥ DAB		r2		1 YA

Assuming that ΥDAB = const, and integrating from r = R to r = f and from YA,R to YA,φ gives

		4Σ RΥ D		§1		YA,φ	·	4Σ RΥ D		ln 1 B ,
				ln ¨			¸				(4.148)
m	A		AB		1	Y			AB
				¨			¸
				©		A,R	¹

where B is the Spalding transfer number,

	B	§	YA,R YA,φ	·
		¨		¸ .	(4.149)
			1 Y
		¨		¸
		©	A,R	¹

Using Eq. 4.148, the evaporating mass flow in Eq. 4.139 can be expressed as

	§1		Yf ,φ		·
	2ddropΣΥ D ln ¨				¸	,	(4.150)
mevap		1	Y
	¨			f ,R	¸
	©				¹

where ddrop = 2r, D is the diffusivity of fuel vapor in air, Yf,φ is the fuel vapor mass fraction outside the boundary layer, and Yf,R is the fuel vapor mass fraction at the

droplet surface. The effect of an increased mass transport due to a relative velocity between droplet and surrounding gas is expressed by the Sherwood number:

	§1		Yf ,φ		·
	ddropΣΥ D ln ¨				¸Sh .	(4.151)
mevap		1	Y
	¨			f ,R	¸
	©				¹

4.5 Evaporation	143

Fig. 4.38. Temperature and mass histories of an evaporating decane droplet

The Sherwood number has a value of Sh = 2.0 if no relative velocity is present. The appropriate Sherwood number, which includes the effect of a relative velocity between droplet and gas, has been proposed by Ranz and Marshall [112],

Sh 2.0 0.6 Re1/ 2 Sc1/ 3 ,	(4.152)
where
Sc µg / Ug D .	(4.143)

The properties of the gas phase inside the boundary layer are calculated using the 1/3rd rule. Assuming equilibrium and using Raoult’s law (e.g. [142]), the fuel vapor fraction in the boundary layer can be calculated as

Yf ,R	pvap Tl		MWf	,	(4.154)
	pcyl		MWmix,R

where pvap(Tl) is the saturated vapor pressure belonging to the droplet temperature,

and MWf and MWmix,R are the molecular masses of fuel and gas mixture at the droplet surface. The vapor pressure can be determined from the Clausius-

Clapeyron equation,

pvap Tl	ª hfg		§	T	·º
	p0 exp «		¨1	0	¸»	,	(4.155)
				Tl
	«RT0		©		¹»
	¬				¼

where (p0, To) is a known point on the vapor pressure curve.

As an example, Fig. 4.38 shows the temperature and mass histories of an evaporating decane droplet as predicted by the above model. The initial and boundary conditions are given in the figure. It can be seen that the droplet temperature increases at first due to the convective heat transfer from the hot gas to the liquid, while the evaporated mass is small. This causes an increased saturated vapor pressure at the droplet surface, which again results in increased diffusive