Utveckling av en Multi-Zonsmodell för NOx Bildning i Diesel Motorer

Development of a Multi-Zone Model for NO

Formation in Diesel Engines

FABRIZIO DIOTALLEVI

Master of Science Thesis Stockholm, Sweden 2007

Development of a Multi-Zone Model for NO _x Formation in Diesel Engines

Fabrizio Diotallevi

Master of Science Thesis MMK 2007:55 MFM109 KTH CICERO

Machine Design SE-100 44 STOCKHOLM

Examensarbete MMK 2007:55 MFM109

Utveckling av en Multi-Zonsmodell för NOx Bildning i Diesel Motorer

Fabrizio Diotallevi

Godkänt

2007-08-20

Examinator

Hans-Erik Ångström

Handledare

Anders Westlund

Uppdragsgivare

KTH CICERO

Kontaktperson

Anders Westlund

Sammanfattning

Detta examensarbete ingår i civilingenjörsutbildning i maskinteknik på KTH, Stockholm.

Målet med projektet var att utveckla en multi-zonsmodell för NOx

bildning i Diesel motorer. Eftersom emissionskraven blir allt hårdare, satsas mycket arbete på att minska skadliga utsläpp samt bränsleförbrukning. Datorsimulationer spelar en viktig roll i denna utveckling eftersom de ersätter dyra och tidskrävande motortester.

Modellen är baserad på en multi-zon approach och använder Zeldovichs mekanism som ger möjligheten att beräkna formationshastighet av NO med hjälp av olika ämneskoncentrationer som är involverade i processen.

Alla de viktigaste fenomenen som sker under Diesel förbränning är modellerade. Jämnviktskoncentrationer har beräknats med stor noggrannhet eftersom de är väldigt viktiga för alla följande steg i modellen.

En Simulink modell har också utvecklats. Den versionen är baserad på den första modellen och den är tänkt för att förenkla kommunikationer med andra mjukvaror, GT-power i detta fall. I stället för iterativa funktioner, använder Simulink look-up tabeller som gör att modellen blir snabbare och lättare att hantera.

Modellen har validerats med mätningar från en encylinder motor baserad på en SCANIA D12 6-cylinder Diesel.

Arbetet har resulterat i en modell som verkar vara tillräckligt noga för att användas som ett verktyg vid simulering av Dieselmotors NO_x- emissioner för den tänkta applikationen, optimering av transient.

Development of a Multi-Zone Model for NOx Formation in Diesel Engines

Fabrizio Diotallevi

Approved

2007-08-20

Examiner

Hans-Erik Ångström

Supervisor

Anders Westlund

Commissioner

KTH CICERO

Contact person

Anders Westlund

Abstract

This thesis work is the last part of the Master of Science education in mechanical engineering at KTH, Stockholm.

The aim of this project was the development of a Multi-zone model for NOx formation in Diesel engines. Because of the stringent emission legislations, great effort is made to decrease the fuel consumption and the harmful emissions of internal combustion engines. Computer simulations play a decisive role in this context because they substitute the expensive and time-consuming laboratory tests.

The model is based on a multi-zone approach and uses the well- known Extended Zeldovich Mechanism which gives a relation to calculate the NO formation rate using the concentration of different gas species involved in the process.

All the most important phenomena involved in the Diesel combustion have been modelled. A special attention is then paid to the equilibrium concentration calculation; this is a significant part of the model which strongly influences the following steps.

A Simulink model has also been developed. This second version is based on the previous model and it is made especially to allow the connection with other external softwares, GT-power in this case. Simulink gives the possibility to create a model which, using look-up tables instead of iterative functions, can be faster and easier to handle if compared to the MATLAB code.

The model has been validated with test sessions on a single cylinder engine based on a heavy duty SCANIA D12 6-cylinder Diesel engine.

This work has finally given a model which seems accurate and precise enough to be used as a tool in the simulation of Diesel engine NOx

emissions.

Notations and Acronyms iii

1 Introduction 1

1.1 The Diesel Engine . . . . 2

1.2 The Diesel Engine Combustion . . . . 4

1.3 The Soot-NO_x Trade-off . . . . 6

1.4 Historical Background . . . . 7

1.5 Present Situation . . . . 8

2 Modelling The Diesel Engine 11 2.1 Heat Release . . . . 11

2.1.1 Gamma . . . . 13

2.1.2 Heat Transfer . . . . 14

2.1.2.1 Convective Heat Losses . . . . 15

2.1.2.2 Radiative Heat Losses . . . . 17

2.1.2.3 Total Heat Losses . . . . 18

2.2 Adiabatic Flame Temperature . . . . 18

2.3 Equilibrium Concentrations . . . . 20

2.3.1 Reaction Rate Constants . . . . 22

2.3.2 Equilibrium Concentration Plots . . . . 25

2.4 Ignition Delay . . . . 29

2.5 Mixing . . . . 30

3 Formation of Nitrogen Oxides 35 3.1 NO_x Production Mechanisms . . . . 35

3.1.1 Thermal NO_x . . . . 35

3.1.2 Fuel NO_x . . . . 36

3.1.3 Prompt NO_x . . . . 36

3.2 Formation of NO . . . . 36

3.3 Formation of NO₂ . . . . 38 i

4.2 Conditions at Inlet Valve Closing . . . . 44

4.3 Pre-combustion Calculations . . . . 45

4.4 Heat Release . . . . 46

4.5 Zone Creation . . . . 46

4.6 NO Calculation . . . . 49

4.7 Post-combustion Calculations . . . . 49

4.8 Model Outputs . . . . 49

5 Simulink Model 51 5.1 Differences with the MATLAB Model . . . . 55

5.2 Advantages and Disadvantages . . . . 56

6 Experimental Setup 57 6.1 Engine Characteristics . . . . 57

6.2 NOx Measurement . . . . 59

7 Results 61 7.1 Heat Release . . . . 62

7.2 Zone Temperature . . . . 64

7.3 Species Concentration . . . . 66

7.4 NOx . . . . 68

7.5 Comparison with Measured Data . . . . 70

7.5.1 Test 1 . . . . 70

7.5.2 Test 2 . . . . 72

8 Summary and Conclusions 75

9 Future Work 77

10 Acknowledgments 79

Bibliography 81

A Ideal Gas Law 85

B Table of Coefficient Sets for NASA Polynomials 87

ii

Symbols

A Area

B Cylinder bore c Specific heat

c_p Specific heat at constant pressure c_v Specific heat at constant volume d Diameter

e Specific energy E_a Activation energy G Gibbs free energy h_c Heat transfer coefficient h Specific enthalpy

H Enthalpy

k Thermal conductivity k_i⁺ Forward rate constant k_i⁻ Backward rate constant

K_c Equilibrium constant expressed in concentrations K_p Equilibrium constant expressed in partial pressures

m Mass

M Molecular weight n Number of moles

N Crankshaft rotational speed Nu Nusselt number

p Pressure Q Heat transfer Q˙ Heat-transfer rate

r Radius

r_c Compression ratio rpm Revolutions per minute

iii

s Specific entropy

S Entropy

S_p Piston speed

t Time

T Temperature

T_exh Exhaust gas temperature T_wall Cylinder wall temperature u Specific internal energy U Internal energy

v Specific volume

V Volume

V_c Clearance volume

V_d Displaced cylinder volume W Work transfer

α Angle

γ Specific heat ratio η Efficiency

η_c Combustion efficiency θ Crank angle

λ_c Combustion air/fuel ratio λ_gl Global air/fuel ratio µ Dynamic viscosity

ν_i Stoichiometric coefficient of species i

ρ Density

σ Stefan-Boltzmann constant τ_id Ignition delay time

φ_c Combustion fuel/air ratio φ_gl Global fuel/air ratio

iv

BDC Bottom dead center CAD Crank angle degree EGR Exhaust gas recycle EOC End of combustion EOI End of injection EVC Exhaust valve closing EVO Exhaust valve opening IVC Inlet valve closing IVO Inlet valve opening LHV Lower heating value PPM Part per million

SCR Selective catalytic reduction SOC Start of combustion

SOI Start of injection TDC Top dead center

v

Introduction

One of the main causes of the actual environmental pollution is the diffuse use of combustion engines, both in the Otto and Diesel versions. Their negative effects on health, climate and environment are well known and accepted by the whole scientific world. It naturally follows that the emission legislations are getting more demanding in every country. At the same time the source of energy for these engines, i.e. crude oil, is running out at high rates. These are the fundamental reasons that explain why the design and development of more efficient and less polluting engines is extremely important nowadays.

This thesis work is focused on the Diesel engine which is, thanks to its high efficiency, one of the best choices toward diminished fuel consumption. This issue is fundamental for the engine manufacturers as an increasing number of Diesel engines is used every year for personal and commercial vehicles. One of the drawbacks of this kind of engine is the production of relatively high amounts of particulates and oxides of nitrogen (NO_x) which have a negative impact on environment and people, the former being harmful for health and the latter causing acid rains and lung illnesses.

Important technological progresses in injection systems, exhaust aftertreatment, EGR (Ex- haust Gas Recirculation) etc. have been recently made but these are not enough to meet the new demanding emission laws that will be applied in the next years. A more detailed and accurate comprehension of the complex combustion process is needed to act directly on the engine-out emissions and try to reduce them as much as possible. Since the control parameters are increasing sharply and the working conditions can vary in a wide range, computer simula- tions play a significant role as a substitute for the real engine tests which are requiring when it comes to time and money.

The mathematical models can simulate and predict the results of many different phenom- ena if the main mechanisms and processes are correctly modelled and implemented and they represent a fundamental tool for the optimization of internal combustion engines. Several sorts of models have been developed to improve the accuracy of the results but since the combustion process is too complex to be modelled in detail, some simplifications need to be done in every mathematical simulation.

1

The model developed by the author is based on Egnell’s basic ideas (see [1]) and is conceived as a fast algorithm that can predict the NO_x production in a Diesel engine with accuracy and precision. The model is founded a multi-zone approach and takes into account all the significant engine parameters and local conditions that can influence the emissions. The core mechanism of the nitrogen oxides formation is the so called Zeldovich mechanism (see [2]) which has been thoroughly investigated by many researchers all over the world, and is considered to be the main source of NO_x emissions.

The model has been implemented in two different versions with different characteristics.

The first is written in MATLAB code while the second is done using the MATLAB toolbox Simulink. The second version is optimized in order to be connected and to interact with other softwares, e.g. GT-power, to perform a complete cycle simulation which includes a prediction of the exhaust emissions.

This thesis work begins with an introductive part where the basic features of Diesel engines and combustion are explained together with a historical overview. The second chapter discusses in general the Diesel engine modelling pointing out the most important details taken into ac- count in this project. It follows an overview of the basic NO_x formation process. The following chapters present the features of the model itself (chapters 4 and 5), the measurement setup (chapter 6), some results (chapter 7) and a short summary (chapter 8). At the end of this work, after acknowledgments and bibliography, two appendixes are included, with an explanation of the ideal gas law (appendix A) and a list of the NASA coefficients (appendix B), both used in the simulations.

1.1 The Diesel Engine

The Diesel engine is a reciprocating internal combustion engine where the fuel is ignited spon- taneously due to the high temperature and pressure in the combustion chamber, phenomenon which gives the often used name compression-ignition engine (see [20] and [4] for a more detailed description). The principal components of such an engine are:

• the piston, which converts the gases pressure into an oscillating movement,

• the cylinder, which delimits the combustion chamber and allows the movement of the piston,

• the connecting rod, which transmits the power from the piston to the crankshaft,

• the injector, which injects the fuel in the combustion chamber,

• the crankshaft, which permits the conversion of the oscillating movement into a rotation.

The air is inducted in the combustion chamber and then compressed with high compression ratio (12-24) by the piston movement while the fuel is injected just before the top dead center

(TDC) as a high pressure spray. Different injection systems and injection pressures are used depending on the size and power of the engine (see [33]). The most common are:

• unit injectors, located directly on the engine top and driven by a camshaft,

• common rail systems, equipped with injectors that receive pressurized fuel from a common fuel rail.

In normal conditions, a Diesel engine works with a relative air/fuel ratio (λ) much higher than the stoichiometric value. This means that the amount of air contained in the combustion chamber is more than what is needed to burn completely the fuel injected. The excess air is necessary to allow a complete combustion of the fuel.

Diesel engines are available both as four-stroke (see figure 1.1) and two-stroke and a in a large number of sizes and configurations. Big engines are usually used for stationary or marine application and have mostly two-stroke cycles while small and medium engines are used on road vehicles typically with four-stroke cycles.

Figure 1.1: Four-stroke Diesel cycle.

1.2 The Diesel Engine Combustion

The Diesel combustion is a heterogeneous and complex process that has been thoroughly studied for many years. The details are still difficult to understand completely but the principal features of the process are well known since the early years of research (see [2] and [30]).

The fuel is injected in the combustion chamber just before the start of combustion (SOC).

As soon as the fuel drops meet the high-temperature and -pressure air in the cylinder they vaporize, mix with the air and ignite causing the temperature to increase in the chamber. The following process of ignition and combustion of the fuel injected later is then accelerated by the heat released.

This process can schematically be divided into three phases as shown in the heat release curve on figure 1.2:

1. The ignition delay, which is the phase between the beginning of the injection and the start of combustion when the fuel is injected in the combustion chamber and mixes with the compressed air.

2. The premixed combustion, when the fuel injected during the previous phase and already mixed with air ignites and burns rapidly.

3. The mixing controlled combustion, which is the last phase that is controlled by the rate at which fuel and air mix.

A conceptual model of combustion in Diesel engines based on optical measurements has been proposed by John Dec in 1997. A schematic view of a fully developed burning fuel spray is shown in figure 1.3. It shows the different areas of the spray where soot and NO_x are formed.

Soot particles start to form in the fuel rich and premixed part of the flame (light blue and grey parts on the figure) and then flow toward the soot oxidation zone where they oxidize and become bigger.

According to the model, NO_x form on the lean side of the spray external layer where the temperature is high enough and oxygen is available (green line on the figure). Since the cylinder pressure rises during the combustion process, the gases burned early are compressed and their temperature increases to a level which is much higher than the combustion temperature. This causes a significant production of NO which is comparable with the amount of NO produced in the burning flame front.

Some conclusions regarding the NO_x production in the combustion process can be then drawn from Dec’s model (see [22]):

• no NO_x is formed during the premixed phase of the combustion,

• NO_x is firstly found in a thin layer that surround the diffusion-burning flame,

• about one third of the total NO_x production takes place later in the post-combustion process.

Figure 1.2: Combustion phases.

Figure 1.3: Spray model.

Some aspects of Dec’s model are still not completely clear and proved but this model is anyway widely used as a theoretical base for every Diesel combustion model.

1.3 The Soot-NO

Trade-off

A big problem when trying to reduce the undesired emissions in a Diesel engine is the soot-NO_x trade-off. As qualitative shown in figure 1.4, when soot emissions are reduced, NO_x increase and vice versa.

This phenomenon is very important in the development of the modern engines and it is due to the fact that the same factors that facilitate the oxidation and consequent reduction of soot, cause also an increase in the amount of NO_x. The principal factors are:

• high temperatures,

• available oxygen.

In order to form, NO_x needs in fact sufficiently high temperatures and oxygen molecules which can combine with the nitrogen as thoroughly explained in chapter 3.

Figure 1.4: Typical soot-NO_x trade-off.

Some techniques have been developed to overcome this problem that is still a great challenge in the reduction of Diesel engine emissions.

1.4 Historical Background

The history of the Diesel engine started in 1892 when a German engineer, Rudolf Diesel (1858- 1913, figure 1.5), patented a new kind of combustion engine. He was born and grew up in France and already as a kid was very smart; he spoke three languages and received at the age of 12 a medal for good school results in 1870. Later on he moved to Germany where he received a master degree at M¨unchen University in 1880.

Figure 1.5: Rudolf Diesel.

After several years working in the industry he presented the results of his studies about this innovative engine to a company in Augsburg which decided to finance the project and to build a prototype. In august 1893 the first prototype was ready: it was fueled by powdered coal injected with compressed air. Unfortunately a strongly unstable combustion did not allow the engine to run continuously.

It took 4 years to build a more stable engine that reached an efficiency of about 31% which was more than double if compared to the other engines available at that time. A new combus- tion chamber design improved the effective efficiency to a level of about 36% already in 1901, increasing the interest in this engine all over the world.

In 1904 Rudolf Diesel patented the engine in 37 countries and became well known and successful but not happy: several bad business deals brought him to a deep depression. He disappeared during a crossing of the English Channel in September 1913: it was suspected to be a suicide.

The first Diesel engines produced commercially were large low-speed engines used mostly in stationary and marine propulsion machines. At the end of the 1920s the demand for light and fast engines increased and most automotive manufacturers started the production of Diesel

engines that were then widely used in the Second World War.

One of the biggest problems toward the development of small car engines was the injection system mainly because the systems based on air assistance that were used for large engines were not suitable for smaller sizes. Developments in this field were made especially in the USA and in Germany where Bosch became the greatest injection system supplier.

The petroleum crises of the 1970s pushed the development of the engines toward higher efficiency and later lower emissions. Several new devices were then introduced, e.g. new injection systems, turbocharging, intercooling etc. The new techniques increased the power output and the efficiency which reached absolutely high levels.

The most recent innovations introduced the electrical systems to control the injection and other engine features. These aspects are still improving and will have an increasing importance even in the future.

1.5 Present Situation

The most important issues of the actual situation are the emission of harmful substances. The emission legislations in every country, especially Europe, USA and Japan are becoming more requiring and the allowed emission limits are strongly decreasing year by year; therefore the engine manufacturers are forced to put much effort and money to deal with this problem.

Diesel engines have low fuel consumption if compared to the other engines and this is the reason why they’re more suitable to reduce the greenhouse gas emissions. The goal is to maintain a low fuel consumption while meeting the tightening demands for low emissions of NO_x and particles (PM).

Road transports give the largest contribution to NO_x and PM emissions. Actually the air pollution causes thousands of deaths per year in Europe as well as other problems including respiratory and cardiovascular diseases, asthma, acute respiratory symptoms and lung cancer. In particular the particles emitted by the Diesel vehicles are amongst the most dangerous because of their small size and because of their chemical composition.

The European Parliament has introduced gradually tighter emission limits for cars, light trucks and heavy trucks in the last years, concerning in particular nitrogen oxides and particles.

Figure 1.6 shows the evolution of the emission limits for heavy-duty Diesel trucks.

The new standard, i.e. Euro 5, will apply from October 2008, substantially decreasing the emission limits of NO_x for new trucks produced for the EU market. Further reductions, i.e.

Euro 6, are planned to be applied in the future, but the new limits have been established just for passenger cars.

Concerning NO_x the most common and advanced techniques available today to reduce the emissions are:

• EGR (Exhaust Gas Recycle), which tries to reduce the NOx production already in the combustion process. It uses the exhaust gases which come from the exhaust pipe to diminish the combustion temperature.

Figure 1.6: EU emission limits.

• HCCI (Homogeneous Charge Compression Ignition), which combines homogeneous charge, typical of the spark-ignition engines, with compression ignition as in normal Diesel engines.

• SCR (Selective Catalytic Reduction), an aftertreatment method which uses urea solution injected into the exhaust system to reduce NOx forming nitrogen and water.

The recent emission standards require a combination of the engine-out reduction technology with advanced exhaust aftertreatments and a big effort is needed in both directions to obtain the best possible results in terms of low emissions.

Modelling The Diesel Engine

The modelling process of the Diesel combustion and emission formation is rather complicated and it requires a deep understanding of the several phenomena involved. On the other hand the results of this combustion process can be considered rather repeatable and the cycle-to-cycle variations relatively small. These are great advantages when trying to simulate numerically this process because they give the opportunity to obtain reliable results if the model is accurate enough.

Different kinds of models have been developed, most of which are zero-dimensional. This means that they don’t provide any information on the spatial configuration of the combustion process inside the cylinder but they provide only global and averaged values. The content of the cylinder is usually divided in various zones which are not real spatial entities but only group the areas with common properties. Depending on the number of zones it is possible to obtain for example two-zone or multi-zone models.

In short in two-zone models the combustion products are grouped in a common zone usually called burned zone, while the fresh charge is contained in the unburned zone. If the burned zone is instead divided in many different parts, one for each time step, then a multi-zone model is obtained, as done in the present work as well as in Egnell’s.

In this chapter all the theoretical steps for the modelling of the principal phenomena in the Diesel combustion, for example the heat release and the chemical species equilibrium, are explained. Some comparisons between different methods have been done to choose the most suitable equations for the model.

2.1 Heat Release

The first step in the creation of a combustion model is the analysis of the rate of heat release which is then used for the calculation of the burned fuel fraction. This is done writing the classical first law of thermodynamics for a system where there is no mass flow as in the period between inlet valve closing (IVC) and exhaust valve opening (EVO), if blow-by, valve leakage

11

and fuel injection are neglected. The heat release is then obtained as:

dQ_hr = dW + dU + dQ_w (2.1)

where:

dW is the work due to piston movement, dU is the change of internal energy,

dQ_w is the heat transfer to the cylinder wall.

Further changes are needed to obtain a new expression of the heat release. The work per- formed can be written as:

dW = p dV (2.2)

By differentiating the ideal gas law, equation 2.3 can then be derived:

m · dT = 1

R(p dV + V dp) (2.3)

The variation of internal energy depends on the temperature according to:

dU = m · c_v· dT (2.4)

Combining equations 2.3 and 2.4 gives:

dU = c_v

R(p dV + V dp) (2.5)

Substituting equations 2.5 and 2.2 yields:

dQ_hr = c_v

R (p dV + V dp) + p dV + dQ_w (2.6)

The first law heat release equation is then obtained by substituting R = c_p− c_v and γ = ^c_c^p

v: dQhr = γ

γ − 1p dV + 1

γ − 1V dp + dQw (2.7)

where the cylinder volume (V) is given by the geometry of the engine and the pressure (p) is measured and is one of the model inputs.

From the heat release trace it is possible to obtain the amount of fuel needed to release such energy, using the fuel lower heating value (LHV).

2.1.1 Gamma

The ratio of specific heats γ depends on temperature and gas composition and it is a definitely important parameter when determining the heat release in the combustion chamber. Therefore an accurate estimation of its value is crucial to obtain a reliable starting point for the model which in turn will affect all the following steps.

Three different expressions of γ have been investigated to find out which one is the most suitable:

1.

γ = 1.35 − 6.0 · T + 1.0 · T² (2.8)

Equation 2.8 is an approximate relationship taken from [21] which gives γ = 1.3 and γ = 1.27 at temperatures of 1000 K and 2000 K respectively. The gamma values calculated with this equation should be within circa ±0.02 of the correct values for Diesel fuel under a wide range of operating conditions if a correct enough temperature is used.

2.

γ = γ₀− k₁exp −k2

T

!

(2.9) where:

γ₀ = 1.38, k₁ = 0.2, k₂ = 900.

Equation 2.9 is taken from [1] and should give quite an accurate result over a temperature range of 1000 - 3000 K with those specific constants (k₁ and k₂) and reference value (γ₀).

3.

γ = gamma(G, T ) (2.10)

Equation 2.10 is a MATLAB function personally made which determines γ for a given mixture of gases, the input vector G, at the temperature T. It uses gas properties from the JANAF tables and performs a weighted average which depends on the gas composition and temperature (see [24]).

Figure 2.1 shows the values of γ calculated with the three expressions above as a function of the temperature. As it can be seen the first and the second curves are quite similar even if Egnell’s has a smaller range of variation. It is easy to notice that the third curve, the self-made function, gives higher γ values in the whole considered temperature range.

Figure 2.1: Ratio of specific heats γ determined with three different expressions.

This difference is then reflected in figure 2.2 where the cylinder pressure is calculated with equation 2.11 and compared with the measured data. The third γ values give the best approx- imation of the real cylinder pressure while the others give a slightly lower result.

From this analysis of three different methods to calculate gamma it is possible to conclude that the best approximation is probably the third which uses the information about the com- position of the gas mixture. Since the importance of the specific heat ratio is crucial for several calculations, e.g. the heat release, these γ values will then be used throughout the following steps, even if the computing time will slightly increase due to the higher complexity of the function.

p1

p₂ =

V2

V₁

^γ

(2.11)

2.1.2 Heat Transfer

The heat losses due to heat transfer from the high temperature charge inside the cylinder to the wall may be divided into two main contributions, i.e. the convective and the radiative losses:

dQ_w = dQ_conv+ dQ_rad (2.12)

These terms can reach up to 10% of the total rate of heat release and consequently cannot be neglected if a high accuracy is requested in the calculation of the combustion and post combustion temperatures.

Figure 2.2: Comparison between measured and calculated cylinder pressure.

2.1.2.1 Convective Heat Losses The convective losses are determined with:

Q˙_conv = A · h_c· (T_gas− T_wall) (2.13) where:

A is the area of the surface [m²], h_c the heat transfer coefficient [_{K m}^W2], T_gas the average gas temperature [K],

T_wall is the estimated average wall temperature [K].

Equation 2.13 is a simplification of the convective heat transfer which uses single-zone ap- proach where the contributions of the unburned and burned zones are not determined but only a total value is given. This is an acceptable approximation which has just a small influence on the determination of the zone temperatures.

Two expressions to estimate h_c have been compared in this work. The first is Annand equation (see [21]):

h_cB k

!

= a ρ S_pB µ

!b

(2.14) where:

B is the cylinder bore [m],

k is the thermal conductivity [_{m K}^W ], a and b are constants,

ρ is the density [_m^kg3],

µ is the dynamic viscosity [P a · · · s], S_p is the mean piston speed [^m_s].

The second is Woschni equation (see [21]):

N u = 0.035 Re^m (2.15)

Assuming B as characteristic length and assuming k ∝ T^0.75 and µ ∝ T^0.62 equation 2.15 can be written as:

hc = 3.26 B^−0.2p^0.8T^−0.55w^0.8 (2.16) where:

B is the cylinder bore [m], p is the cylinder pressure [kPa], T is the gas temperature [K],

w is the average gas velocity [^m_s] estimated with the following equation:

w =

"

C₁S_p+ C₂ VdTivc

p_ivcV_ivc(p − p_m)

#

(2.17) where:

C₁ = 2.28 for the combustion and expansion period, S_p is the mean piston speed [^m_s],

C₂ = 3.24·10⁻³ for the combustion and expansion period, V_d is the displaced volume [m³],

Tivc is the temperature at inlet valve closing [K], p_ivc is the pressure at inlet valve closing [Pa], V_ivc is the volume at inlet valve closing [m²],

p is the cylinder pressure [Pa],

p_m is the motored cylinder pressure [Pa].

Figure 2.3 shows the heat transfer rate determined with the two expressions above. The curves have a rather similar trend but the values calculated with equation 2.14 are higher in a certain range. Both the results can be considered a good approximation of the real convective heat transfer. Annand equation has been finally implemented in the model since, as shown in [21], Woschni equation tends to give slightly lower heat release if compared with real data. An overestimation of the heat transfer would give excessive heat release energy and consequently an amount of burned fuel higher than the actual injected one.

Figure 2.3: Convective heat transfer rate calculated with two different expressions.

2.1.2.2 Radiative Heat Losses

The radiative heat losses are rather hard to estimate with a simple expression. The most commonly accepted and widely used formula is the one proposed by Annand (see [2]):

˙

q_rad = β σ (T_gas⁴ − T_wall⁴ ) (2.18) where:

β is a coefficient assumed equal to 0.6 (see [2]), σ is the Stefan-Boltzmann constant,

T_gas the average gas temperature [K],

T_wall is the estimated average wall temperature [K].

2.1.2.3 Total Heat Losses

Figure 2.4 shows the contribution of the convective (Annand equation is used) and radiative losses to the overall loss due to heat transfer through the combustion chamber wall. The total heat transfer reaches a maximum value of about 20 J/CAD. This figure can be compared with the heat release graphs shown in section 7.1 to understand the relative importance of the heat transfer in the modelling process.

Figure 2.4: Convective and radiative heat transfer rate.

2.2 Adiabatic Flame Temperature

The adiabatic flame temperature has a key role in the determination of the NO_x emissions since the process which governs the nitrogen oxides production is deeply influenced by the temperature in the combustion chamber. In the present model the flame temperature has been calculated with an iterative function which compares the enthalpy of the reactants before the combustion with that of the products and keeps this value constant (see [6] and [7]).

The first step is to guess a combustion temperature depending on the reactants temperature, on the lambda value and on other parameters:

T_products = T_reactants+ ∆T (2.19) The equilibrium concentration at this guessed temperature is then determined solving the equi- librium equation system. By using the equilibrium concentration it is then possible to determine the enthalpy of the products with the JANAF data. This value is finally compared with the enthalpy of the reactants calculated with the same procedure.

∆H = Hproducts− Hreactants (2.20)

The error in the enthalpy balance is used to correct the adiabatic flame temperature starting an iterative and time consuming procedure which stops when an accurate enough value is obtained:

Tproducts,new = Tproducts,old− ∆H

c_v (2.21)

The values of temperature calculated with this iterative function have been compared with other results to verify their accuracy. One of the programs used for comparison is the one developed by Depcik (see [17]).

The same process is followed if EGR is present, taking in account the EGR percentage and composition during the enthalpy calculations. An example is shown in figure 2.5 where 3 cases with different global lambda values (1, 4 and 7 respectively) are analyzed.

Figure 2.5: Adiabatic flame temperature with different EGR ratios.

It is possible to notice how the adiabatic flame temperature decreases when the amount of EGR increases and the different influence of the global lambda ratios. The highest variation in

the flame temperature is obtained with the stoichiometric value while the difference is smaller with larger ratios due to the higher amount of O₂ and N₂ in the EGR gases.

The values shown in figure 2.5 have not been compared to other results or validated with experimental data but give only an overall view of the EGR effects on the adiabatic flame temperature.

2.3 Equilibrium Concentrations

One of the most critical aspects for the calculation of the thermodynamic properties of the system is the knowledge of the precise composition of the gas mixture inside the combustion chamber at every time step. It is then necessary to develop a chemical equilibrium scheme which includes a certain number of species, reactants and products of the combustion process.

Many approaches have been used by other authors who created algorithms considering a different number of chemical species. The most simple models deals with the problem using only six species as products of the combustion process, namely CO2, H2O, N2, O2, CO and H₂ (see [1]). The model developed by NASA instead considers a very high number of chemical species (see [15]). This model though is used rather rarely in combustion engine simulation due to its complexity and long computational times.

Some species are present in significant quantities only for particular engine conditions and operation parameters which are not found so often in practice. Therefore, for the calculations in this work, the number of important species is reduced to 10, i.e. the six mentioned above plus O, H, OH and NO. Several models consider also N and in some cases Ar (see [25] and [26]) but an analysis of the results has showed that these species highly increase the complexity of the model giving negligible improvements in the results, at least at the requested accuracy levels.

To calculate the equilibrium concentration of 10 chemical species, a system with an equal number of equations is needed. The chemical reactions which include the considered species are:

CO + H₂O ⇔ H₂+ CO₂ (2.22)

CO + 1

2O₂ ⇔ CO₂ (2.23)

1

2O₂ ⇔ O (2.24)

1

2H₂ ⇔ H (2.25)

OH + O ⇔ O₂+ H (2.26)

N₂+ O₂ ⇔ N O + N O (2.27)

The typical reaction times have been calculated and verified to be shorter than the time step used. Reaction 2.27 is actually the slowest and could be an error source if considered at the

equilibrium. However, in the typical high pressure and high temperature conditions of a normal engine, this error may be considered small enough and consequently neglected.

The chemical equation which describes the combustion is expressed as follows:

C_nH_m+



n +m 4



λ (O₂+ 3.773 N₂) ⇔ CO₂+ CO + H₂O + H₂+ O₂+ N₂+ H + O + N O + OH (2.28) From now on, the 10 chemical species considered in the model will be denoted with a number as follows:

CO₂ = 1 CO = 2 H₂O = 3 H₂ = 4 O₂ = 5 N₂ = 6 H = 7 O = 8 N O = 9 OH = 10

Denoting the reaction coefficients of each species as a_j with j=1, 2,...10 and writing the C-H-O-N atom balances to the equations 2.22 to 2.27, the first 4 equations of the system can be obtained:

• Carbon balance:

n = a₁+ a₂ (2.29)

• Hydrogen balance:

m = 2 (a₃+ a₄) + a₇ + a₁₀ (2.30)

• Oxygen balance:

2



n + m 4



λ = 2 a₁+ a₂+ a₃+ 2 a₅+ a₈+ a₁₀ (2.31)

• Nitrogen balance:

2 · 3.773



n + m 4



λ = 2 a6+ a9 (2.32)

The remaining 6 system equations are obtained writing the dynamic equilibrium of the reactions above, 2.22 to 2.27:

KP 1 = P_CO2P_H2

P_COP_H2O (2.33)

KP 2 = P_CO

P_CO^qP_H2O (2.34)

K_{P 3} = P_O2

√P_O (2.35)

KP 4 = P_H2

√P_H (2.36)

KP 5 = P_COP_O

P_OP_OH (2.37)

K_{P 6} = P_{N O}²

P_N2P_O2 (2.38)

(2.39) where P_i is the partial pressure of the species i defined as follows:

Pi = a_iP

P10

j=1 a_j (2.40)

A thorough explanation of the meaning of the equilibrium constant Kp and of the calculation procedure is given in section 2.3.1. Substituting the stoichiometric coefficients a_i with i=1, 2...10, the following equations can be written:

K_{P 1} = a₁a₄

a₂a₃ (2.41)

K_{P 2} =

v u u t

a₁ a₂

2P10 i=1 a_i

a₅P (2.42)

K_{P 3} = a₇

P10 i=1 a_i ·

P10

i=1 a_i· (R T )^0.5

a₄ (2.43)

K_{P 4} = a₈

P10 i=1 a_i ·

P10

i=1 a_i· (R T )^0.5

a₅ (2.44)

K_{P 5} = a₅a₇ a10 a8

(2.45) K_{P 6} = a9

P10 i=1 a_i ·

P₁₀

i=1 ai· (R T )²

a₅a₆ (2.46)

The previous 4 atom balances (2.29 to 2.32) added to the 6 equations above (2.41 to 2.46) give the system of 10 equations needed to find the reaction coefficients of each species, i.e. a₁ to a₁₀. This system can be solved using the Newton-Raphson method for non-linear systems (see [32]). An iterative function has been developed by the author and compared with the MATLAB function fsolve. The results showed a really good concordance and similar computational times if the initial guess was chosen carefully.

2.3.1 Reaction Rate Constants

Every chemical reaction can be expressed as follows (see [5] and [8]):

n

X

i=1

ν_RiM_Ri ⇔

m

X

j=1

ν_PiM_Pi (2.47)

where:

R_i are the reaction reactants, P_j are the reaction products,

ν are the stoichiometric coefficients.

All the reactions work in both directions transforming reactants in products and vice versa until an equilibrium situation is reached, when the reaction speed in both directions is equal and the amount of products and reactants does not change anymore. It is then possible to define a forward and a backward reaction rate as follows:

R⁺ = k⁺

n

Y

i=1

[M_Ri]^νRi (2.48)

R⁻ = k⁻

m

Y

j=1

[M_Pi]^νPi (2.49)

where:

k⁺ is the forward reaction rate constant, k⁻ is the backward reaction rate constant.

The reaction rate constants have the Arrhenius form:

k = A · T^be^−EART (2.50)

where:

A is the pre-exponential factor, T is the temperature,

b is a temperature exponent, EA is the activation energy, R is the universal gas constant.

When the equilibrium situation is reached the following expression which defines the equi- librium constant K_c is obtained:

k⁺ k⁻ =

Qm

i=1[M_Pi]^νPi

Qn

i=1[M_Ri]^νRi =

n

Y

i=1

[Mi]^νi = Kc (2.51)

The constant expressed by the equation above is based on the species concentrations usually given in moles per cubic centimeter. Another useful constant K_p based instead on partial pressures can be written:

Kp = Kc(R T )^−∆n (2.52)

If the difference between the number of product moles and the number of reactant moles is 0, then K_p and K_c are equal.

Two different methods have been used in this work to calculate the reaction rate constants.

The rate constants for the reactions 2.22 and 2.23 have been determined using JANAF Thermochemical Tables (see [16]) according to the following:

log₁₀K_p,reaction =^X

i

ν_ilog₁₀K_p,i (2.53)

where:

K_p,reaction is the equilibrium constant for the reaction,

K_p,i is the equilibrium constant for the formation of the species i,

Since the JANAF values are given only at certain temperatures, an interpolation method is needed. For this purpose an equation of the type proposed by Lapuerta (see [10]) has been used:

log10Kp,i = A ln

 T 1000



+ B

T 1000

+ C + D

 T 1000



+ E

 T 1000

²

(2.54) The calculated coefficients are shown in table 2.1

Table 2.1: Coefficients used in equation 2.54

Species A B C D E

CO₂ -0.0149 20.6070 -0.0326 -0.0602 0.0040 CO 0.60385 5.9144 5.1071 -0.61367 0.055736 H2O -1.0438 12.4052 -2.9927 0.6755 -0.0769

An external function has been developed in order to lighten the principal code. This function receives as input the temperature and gives directly as output the equilibrium constants of the two reactions:

K_p,1 = 10^log(^K_p,CO2)^−log(^Kp,CO)^−log(^K_p,H2O) (2.55)