Baumgarten

194 4 Modeling Spray and Mixture Formation

Fig. 4.71. Momentum cell

The second assumption is that the mass flux, the tangential momentum, and the dynamic pressure, which are added to the liquid film due to impinging drops, are averaged over the cell area. Further on, it is assumed that the velocity profile in the cross-film direction, which is needed when integrating across the film thickness, is either laminar or turbulent.

After integrating in the cross-film direction, the continuity equation yields

		ωΓ	Nside	&	&
	Υl Awall		¦	( v film n )iΓili Υl Awall mimp mevap ,				(4.297)
		ωt
			i 1

and after some rearrangement it results in

ωΓ

1

Nside

&

&

mimp

mevap

¦( v film n )iΓili

,

(4.298)

ωt

A

Υ

A

Υ

A

wall

i 1

l

wall

l

wall

where Awall is the&wall cell area, li and Γi are the length of side i and the film thickness of side i, v film is the mass-averaged film velocity relative to the wall, and

mimp is the difference of incoming mass flux due to impinging drops and outgoing

mass flux due to splashing. The quantity mevap is the rate of fuel vaporization and must be predicted by a film vaporization model, Sect. 4.5.4.

The momentum equation is derived in the same manner. Fig. 4.71 shows a typical momentum cell containing gas (upper part) and the liquid film at the bottom. The momentum equation yields

	ω Γ v&film		Nside &	& &
Υl Awall			¦v film ( v filmn )i ΥlΓiliΙi
		ωt
			i 1		(4.299)
Nside				Nedge &
	&		&		&
¦ p n		Γili Υl AwallΓ g Mtang ¦ Ω			Ai Υl AwallΓ a ,
i 1				i 1

and after some rearrangement the following form can be obtained:

4.8 Wall Impingement 195

ω Γ v&film

1

Nside

&

&

&

¦v film ( v filmn )i ΓiliΙi

ωt

A

wall i 1

Nside

&

Nedge

&

A

(4.300)

¦

p n Γ l

Mtang

¦

Ω

i i

&

i

i

&

i 1

Γ g

i 1

Γ a .

Υl Awall

Υl Awall

Υl Awall

The left-hand side of Eq. 4.300 is the material derivate of the film momentum. The convective momentum term (second term) is an approximation of the integration of the non-linear convective term in the cross-film direction (due to the velocity profile), where

Ι

1

4t

1

Γ

t

Γ

♣

Γt ∙

2

(4.301)

1

Γ

♦1

÷

Γ

♥

≠

is used in order to compensate the effect of using the mass-averaged film velocity instead of a momentum-averaged velocity derived from the exact velocity profile. In order to calculate the displacement thickness Γt and the momentum thickness 4t, a velocity profile has to be specified in the cross-film direction.

The right-hand side of Eq. 4.300 is the sum of all relevant forces. The first term describes the pressure force. The pressure

p pamb pimp

(4.302)

does not vary in the cross-film direction and includes the ambient pressure pamb and the dynamic pressure pimp due to drop impingement and splashing effects, which can be calculated using the momentum equation for one-dimensional flows,

Ndrop	A	Nsplash		Aj
pimp Υl ¦ vn,k2	k	Υl ¦ vn,2	j		.	(4.303)
	A			A
k 1	wall	j 1		wall

In Eq. 4.303, vn,k is the velocity component of the incident drop normal to the surface, and vn,j is the normal velocity of the droplets resulting from splashing. Ndrop

and Nsplash are the total numbers of incident and splashed droplets, and Ak and Aj are the corresponding droplet areas.

The second term on the right hand side describes the gravity effect, which is important for wall films on vertical surfaces. The third term is the change of tangential momentum (index Ω) due to spray impingement and splashing and is given as

Ndrop	Nsplash	mj v&Ω , j .
Mtang ¦ miv&Ω ,i ¦			(4.304)
i 1	j 1

196 4 Modeling Spray and Mixture Formation

In Eq. 4.304, mi and mj are the masses of the incident drops and the droplets resulting from splashing.

The effect of viscous shear forces is given by the fourth term on the right-hand side,

Nedge	Nside
¦ Ω&i Ai	¦(Ω&edgeΓ l )j Ω&wall Awall Ω&liq / air Awall ,	(4.305)
i 1	j 1

where Ω&edge&is the viscous shear along the edges of the film cell, Ω&wall is the wall shear, and Ω liq/air is the shear at the liquid-gas interface.

The last term on the right-hand side describes the effect of an acceleration of the solid surface, which is only relevant for fuel films on piston or valves.

The film model has been combined with a film evaporation model and has been successfully validated against experimental data under diesel engine conditions [133]. The relevant engine parameters and operating conditions are summarized in Table 4.8. As an example, Fig. 4.72 shows a comparison of measured and simulated film thickness as a function of crank angle at different positions on the piston

Table 4.8. Engine parameters and operating conditions [133]

Bore [mm]	150
Stroke [mm]	225
Compression ratio [/]	14
Wall temperature [°C]	400
Swirl ratio [/]	0.0
Fuel injected per hole [mm3/cycle]	25.17
Orifice diameter [mm]	0.35
Start of injection [deg. CA]	-19
Injection duration [deg. CA]	34

Fig. 4.72. Comparison of film thickness. a simulation, b measurement, data from [133]

4.9 Ignition	197

bowl surface for the full load condition. Despite some small discrepancies between measurement and experiment, a good overall qualitative and quantitative agreement is obtained. A high level of agreement is also attained in the case of reduced loads and different injection timings, emphasizing the suitability of the model for the simulation of film development and film evaporation.

4.9 Ignition

4.9.1 Auto-Ignition

The auto-ignition of hydrocarbons in diesel engines is a chain-branching process including the four reaction classes of chain initiation, chain propagation, chain branching, and chain termination. After the start of injection, ignition occurs after a certain induction time, the ignition delay. During this time delay, fuel evaporates until a first region of ignitable mixture with an air-fuel ratio of 0.5 < Ο < 0.7 is formed. Furthermore, the chemical reactions in this region have to produce enough fuel radicals in order to start the combustion process. The chain initiation produces these first radicals from stable fuel molecules. This reaction proceeds slowly, because stable molecules are involved in the process. Then, if a certain radical concentration is reached, the chain propagation and the chain branching reactions form additional radicals. The chain propagation reactions change the nature of the radicals but not their number. Some of the chain propagation reactions produce radicals, which then take part in the chain branching reactions that increase the number of radicals and result in a considerable acceleration of the reactions, leading finally to explosion. The ignition delay is strongly temperaturedependent, a rise of temperature decreases this time.

The multi-stage ignition process can be divided into three temperature regions, the low temperature reactions (cool flame regime), the intermediate temperature region, and the high temperature oxidation. Because details concerning the relevant reactions in the three temperature regimes are given in Sect. 6.4.2, only a brief description of the relevant mechanisms will be given in the following. The cool flame regime typically occurs at gas temperatures between 600 and 800 K. Reactions proceed slowly with only a small temperature rise. However, for increasing temperatures, the production of radicals by the cool flame reactions is reduced because the reverse reactions become faster (degenerate chain branching). Hence, this intermediate temperature region is characterized by the so-called negative temperature coefficient (NTC), which represents an increased ignition delay for increased temperatures. As soon as the temperature, which is increased by the heat release of the cool flame reactions and the further compression of the cylinder charge, is reached, the high temperature chain branching reactions (T > 1000 K) lead to explosion.

Auto-ignition takes place on time scales that are relatively long compared to the relevant hydrodynamic time scales [134]. Hence, the influence of convective and

198 4 Modeling Spray and Mixture Formation

diffusive species transport must be taken into account. This is usually done by solving mass conservation equations for a certain radical indicator species, also including source terms for the chemical reactions. These source terms are expressed by Arrhenius-type reaction rates. Ignition timing is predicted by defining a certain threshold of the indicator species that has to be reached. Then, the ignition model is switched off and the calculation continues with a combustion model. Details concerning appropriate combustion models are given in [134], for example.

The most widely used auto-ignition model is the so-called Shell model, which was developed by Halstead et al. [47]. The name of the model stems from the affiliation of the authors. The model was originally developed to predict knock in spark ignition engines and was later adjusted to predict auto-ignition in diesel engines (e.g. Kong et al. [67], Sazhina et al. [124]). Because it is not possible to model all of the several hundreds of relevant reactions during ignition, the model is based on a class chemistry concept and includes only eight reaction steps between five species. It represents a virtual mechanism between generic species and is formulated to reflect the multistage ignition behavior of hydrocarbon air mixtures including a degenerate chain branching mechanism. The eight reaction steps are given as

RH O2

kq o2R*

(chain initiation)

(4.306)

R

*

kp

oR

*

P heat

(chain propagation)

(4.307)

R*

f1kp oR* B

(chain propagation forming B)

(4.308)

R

*

f4 kp

oR

*

Q

(chain propagation forming Q)

(4.309)

R

*

Q

f2 kp

oR

*

B

(chain propagation forming B)

(4.310)

B

kb o2R*

(degenerate branching)

(4.311)

R

*

f3kp

otermination

linear termination

(4.312)

2R*

kt otermination

quadratic termination.

(4.313)

In Eqs. 4.306–4.313, R* represents the radical, RH is the fuel, Q is an unstable intermediate agent, B is the branching agent, and P represents oxidized products. The concentrations of the different species can be calculated solving the differential equations for their change rates, which are given as

	d[R* ]	2kq [RH][O2 ] 2kb [B] f3kp [R* ] - kt [R* ]2 ,									(4.314)
	dt
		d[B]	f k	p	[R* ] f	2	k	p	[R* ][Q] k	[B] ,	(4.315)
		dt	1						b

4.9 Ignition 199

d[Q]

f4kp

[R* ] - f2kp[R* ][Q] ,

(4.316)

dt

d[O2 ]

pkp [R* ] ,

(4.317)

dt

d[RH] [O2 ] - [O2

](t 0)

[RH]

0)

.

(4.318)

dt

p m

(t

The quantity m depends on the number of hydrogen atoms in the original fuel molecule CnH2m, p = (n(2-ϑ) + m)/2m, and ϑ | 0,67 is the CO/CO2 ratio. The rate coefficients

f1	Af 1 exp E f 1	/ RT [O2	]x1 [RH]y1 ,	(4.319)
	f2 Af 2 exp E f 2 / RT ,			(4.320)
f3	Af 3 exp E f 3	/ RT [O2	]x3 [RH]y3 ,	(4.321)
f4	Af 4 exp E f 4	/ RT [O2	]x4 [RH]y4 ,	(4.322)
ki Ai exp Ei / RT , i = 1, 2, 3, 4, q, b, t ,				(4.323)

are of Arrhenius-type, and

Fig. 4.73. Comparison of simulated and measured ignition delays [47], RON: Research Octane Number

200 4 Modeling Spray and Mixture Formation

ª		1			1					1	º	1
kp «											»	.	(4.324)
	k	[O	2	]		k	2		k	[RH]
¬	1							3			¼

The model includes 26 parameters that have to be adjusted for each fuel of interest. More details are given in Kong et al. [67], for example.

The Shell model is the most-used model in simulations of diesel engine autoignition today. It is capable of describing the negative temperature coefficient effect as well as the influence of temperature, pressure, and fuel-air mixture on ignition delay. Fig. 4.73 shows a comparison between simulated and experimentally obtained ignition delays as a function of temperature for three primary reference fuels of different octane quality.

However, if auto-ignition processes with ignition delays much longer than those of conventional diesel engines are regarded (e.g. HCCI processes), the Shell model is no longer capable of predicting the auto-ignition process with sufficient accuracy. In these cases, more complex chemical models that include a much more detailed description of the low temperature reactions have to be used, see also Sect. 6.4.2.

4.9.2 Spark-Ignition

Gasoline-air mixtures in spark-ignition engines do not auto-ignite due to compression. In these engines, the start of combustion is initiated by the energy of an electrical spark. The spark forms a high-temperature zone between the two electrodes of the spark plug, in which radicals are produced and combustion begins. In this zone, enough heat must be released in order to heat up the neighboring mixture and to initiate a flame front that can propagate into the combustion chamber without any further energy input. During the time the electrical spark exists, a highly ionized plasma with temperatures of about 60000 K [134] forms a channel between the electrodes and then expands to become a spherical ignition kernel, while the temperature decreases to approximately 4000 K, Fig. 4.74. Further details about the kernel growth process are given in [83] and [51], for example. The time span needed for this ignition kernel growth is in the microto millisecond range. The flame front starts from this ignition kernel and propagates into the combustion chamber. At first, it propagates with laminar flame speed, and as the flame surface increases, turbulence increases the flame speed.

Fig. 4.74. Ignition kernel growth

4.9 Ignition	201

In CFD models, a phenomenological sub-model usually describes the growth of the ignition kernel. If the kernel has reached a certain size, the calculation proceeds with a combustion model for turbulent flames.

The modeling of ignition kernel growth and flame speed is usually based on the assumption that the temperature inside the kernel is uniform and equal to the adiabatic flame temperature of the associated fuel-air mixture. Because the time needed to form a spherical ignition kernel is in the same range as the computational time step t0 (about 1 µs), an isothermal ignition kernel of adiabatic flame temperature is assumed to exist at the end of the first numerical time step of ignition. The radius rk,0 of this kernel is determined from an energy balance: the sum of the electrical energy Wsp of the spark discharge and the chemical energy released due to the combustion of the mixture inside the kernel results in a rise of the temperature inside the kernel from the unburned temperature Tu to the adiabatic temperature Tad:

	4	Srk3,0 Uk cp Tad Tu	KspWsp 't0	4	Srk3,0 Uk LHVmix .	(4.325)
	3			3

In Eq. 4.325, LHV is the lower heating value of the air-fuel mixture per gram mixture, and Κsp is the energy transfer efficiency, which can be assumed to be Κsp | 1.0 in this first phase of kernel growth. From Eq. 4.325 the radius

	ª	3K W		't	0	º1/ 3
rk ,0	«	sp	sp			»	(4.326)
	¬«4SUk cp Tad		Tu LHVmix ¼»

of the kernel after the first time step can be derived. During the following time steps the kernel size increases. Because all of the energy of the mixture inside the kernel with radius rk is already released, only the surplus volume, which is the product of the kernel surface and the increase of kernel radius, can be released. Thus, the energy balance gives

2 drk		Uu cp Tad Tu	2 drk
4Srk			KspWsp 4Srk		Uu LHVmix .	(4.327)
	dt			dt

The energy transfer in this later phase of kernel growth is less efficient than in the first phase, the value of Κsp beeing about 30–50% [134]. The quantity drk/dt is called plasma velocity spl. Eq. 4.327 yields

		K W
	spl	sp	sp	.	(4.328)
		4Srk2 Uu cp Tad Tu LHVmix

Because spl is inversely proportional to rk2, it is reduced with increasing kernel radius. However, as the kernel growths, the effect of turbulence on the effective flame velocity increases and leads to an increasing velocity of the flame front propagating into the combustion chamber, see Fig. 4.75.

202 4 Modeling Spray and Mixture Formation

Fig. 4.75. Characteristic behavior of the effective flame speed after ignition [134]

Fan et al. [36] published the discrete particle ignition kernel (DPIK) model, which is based on the ideas of ignition kernel modeling as presented above. In this model, the initial flame kernel is represented by Lagrangian marker particles, allowing to reduce the influence of grid size effects on the ignition process. The initial ignition kernel is assumed to be spherical. The initial diameter dk,0 is assumed to be in the range of the gap size between the two electrodes, usually 1 mm. The radial kernel growth rate is calculated as a function of the laminar flame speed slam and the turbulent kinetic energy k of the flow field,

u	k	drk		Tad	s	2k	.	(4.329)
		dt		Tu	lam	3

The parameter Tad /Tu accounts for the effect of thermal expansion. Tad is the adiabatic flame temperature, and Tu is the local unburned gas temperature, which is determined assuming adiabatic compression from the conditions inside the kernel before (index 1) and after the start of ignition (index 2):

Tu	§ pk ,2 ·				N 1 / N	(4.330)
	¨			¸	.
T		p
k ,1	©		k ,1 ¹

The laminar flame speed slam is given by a relation of Metghalchi and Keck [87], which is

s

1.0

2.1R

ª

26.32 84.72 ) 1.13

2

º

¬

¼

lam

Tu

º

(4.331)

ª

º2.18 0.8 ) 1

ª pk ,2

0.16 0.22 ) 1

«

»

«

»

298

¬

K ¼

¬101.3 kPa ¼

In the case of gasoline [135]. In Eq. 4.331 ) is the equivalence ratio in the spark region, and R is the residual mass fraction. Hence, the diameter of the ignition kernel can be calculated as

dk 2uk t tign dk ,0 .

(4.332)

A one-step reaction is chosen during this early stage of combustion,

References	203

C8 H18 12.5O2 8CO2 9H2O ,

(4.333)

and the burn rate is calculated by

d Υ

­

Υ

f

Υ

O

½

i

C

°

°

¦

s

MW C

min ®

,

2

¾

.

(4.334)

dt

MWO

Csto,O

w

°MWf Csto, f

°

l

i

sto,i

¯

2

2

¿

In Eq. 4.334 MWi is the molecular weight of species i, Csto,i are the stoichiometric coefficients, and Cw is a constant used to account for the wrinkling effect of the

kernel surface. For example, Fan et al. [36] use a value of Cw = 80 while Stiesch et al. [135] use Cw = 20. The quantity

¦	N p,cell		Σ dk2		(4.335)
	N	p,total		V
				cell

is the flame surface density within a particular cell, where Np,total is the total number of marker particles, and Np,cell is the number of particles in the specific cell having a volume Vcell.

The combustion is initiated by the ignition model, and if the kernel size reaches the order of the integral length scale of turbulence,

dk t Ck lt ,

(4.336)

the model switches to an appropriate turbulent combustion model. In Eq. 4.336 lt = 0.16 k1.5/Η is the turbulence length scale, and Ck = 3.5.

References

[1]Agmon N, Alhassid Y, Levine RD (1979) An Algorithm for Finding the Distribution of Maximal Entropy. J Comp Phys, vol 30, no 2, pp 250–258

[2]Ahmadi M, Sellens RW (1993) A Simplified Maximum-Entropy-Based Drop Size Distribution. Atomisation and Sprays, vol 3, pp 291–310

[3]Akao F, Araki K, Mori S, Moriyama A (1980) Deformation Behaviours of a Liquid Droplet Impinging onto a Hot Metal Surface. Trans I S I J, vol 20, pp 737–743

[4]Allocca L, Corcione, FE, Fusco A, Papetti F, Succi S (1994) Modelling of Diesel Spray Dynamics and Comparison with Experiments. SAE paper 941895

[5]Amsden AA, O’Rourke PJ, Butler TD (1989) KIVA-II – A Computer Program for Chemically Reactive Flows with Sprays. Los Alamos National Laboratories, LA- 11560-MS

[6]Amsden AA, Ramshaw JD, O’Rourke PJ, Dukowicz JK (1985) KIVA.. A Computer Program for Twoand Three-Dimensional Fluid Flows with Chemical Reactions and Fuel Sprays. Los Alamos National Labs, Rept LA-10245-MS

[7]Arcoumanis C, Gavaises M, French B (1997) Effect of Fuel Injection Process on the Structure of Diesel Sprays. SAE paper 970799

[8]Ashgriz N, Poo JY (1990) Coalescence and Separation in Binary Collisions of Liquid Drops, J. Fluid Mech., vol. 221, pp. 183–204