ISSN 0280-5316

ISRN LUTFD2/TFRT--5887--SE

Modeling of NOx formation

in heavy duty engines

Kenan Muri

Department of Automatic Control

Lund University
September 2011
 Document name
Lund University
MASTER THESIS
Department of Automatic Control Date of issue
Box 118 September 2011
SE-221 00 Lund Sweden Document Number
ISRN LUTFD2/TFRT--5887--SE
Author(s) Supervisor
Kenan Muri Ola Stenls, Scania CV AB Sdertlje, Sweden
Rolf Johansson Automatic Control Lund, Sweden
(Examiner)
Sponsoring organization

Title and subtitle

Modeling of NOx formation in heavy duty engines. (Modellering av NOx bildning i HD-motorer)

Abstract
In order to apply combustion feedback control in diesel engines to control pollution formation,
good models are needed to estimate the states. One of those would be NOx levels which are
regulated by European law. There are conventional sensors on the market that can measure
NOx but these are too slow for incycle feedback control.
The purpose of this thesis was to derive a model to estimate NOx where the engine very well
could have a built in EGR function. The model seems promising, however additional tuning and
modification is needed for it to be fully useful in practice.
The model presented in this thesis is a zerodimensional twozone model. The NOx model
assumes that combustion locally occurs at stoichiometric conditions. Model validation is done
against singlecylinder engine NOx output data received from the Royal Institute of Technology
(KTH). Leastsquares fitting methods are also used to simplify calculations considerably. The
results indicate that the model is applicable in practice.

Keywords

Classification system and/or index terms (if any)

Supplementary bibliographical information

ISSN and key title ISBN

0280-5316
Language Number of pages Recipients notes
English 80
Security classification

http://www.control.lth.se/publications/
Acknowledgements

First off, I would like to thank my supervisor Dr. Ola Stenls for all the
encouragement and support during this master thesis project. I would also like to
thank my examiner Professor Rolf Johansson for valuable inputs on the report.

I am also thankful to Associate Professor Per Tunestl at the combustion engine

department in Lund for inputs on the NOx model. Mikael Lindstrm has been
essential in this thesis by providing me with data from the Royal Institute of
Technology (KTH). Christer Lundberg at NMED also provided me with valuable
data on heat release.

Last but not least, I would like to thank all the good people at NMEB, the engine
control group, for many interesting discussions around the fika table.

The thesis work was done at Scania CV AB in Sdertlje.

Page 4 (80)
Page 5 (80)
Contents
1 Introduction ... 9

1.1 Motivation.....9

1.2 Problem Formulation........10

1.3 Methodology......11

2 Background .....11

2.1 The Exhaust Gas Recirculation......11

2.2 The Selective Catalytic Reduction.....12

2.3 Models for NOx estimation purposes.....13

2.3.1 NOx measurement...14

2.4 A conceptual model of diesel combustion........15

2.5 The FPGA......16

3 Theory ...18

3.1 The rate of heat transfer......18

3.2 The Exhaust Gas Recirculation......20

3.2.1 Background..20

3.2.2 The Zeldovich Mechanism....20

3.3 The two zone combustion.......24

3.4 Stoichiometric calculations......27

3.5 Calculations and estimations of parameters........28

3.5.1 Estimating and cp.....28

3.5.2 Calculation of molar masses and SAFR.....29

3.5.3 The average gas velocity estimate......29

3.5.4 Approximation of equilibrium concentrations.....30

3.5.5 A description of the EGR dynamics.........30

3.5.6 The cylinder geometry.......34

3.5.7 Mass relation between the two zones........35

3.5.8 Heat transfer between zones.......36

3.6 Calculating equilibrium concentrations.........37

Page 6 (80)
 3.7 Adding products to the post combustion zone............39

4 Implementation in MATLAB.......42

4.1 Heat release calculations...............42

4.2 Temperature calculations...............42

4.3 Equilibrium calculations..................43

4.4 Zeldovich Mechanism.................43

4.5 C code implementation...............43

5 Implementation considerations for FPGA code..44

5.1 The division operation................44

5.1.1 Thesis algorithm.....44

5.1.2 Modified thesis algorithm...46

5.2 Tabulation of data................48

5.2.1 The least-squares method.....48

5.2.2 Functions of a single variable...49

5.2.3 Data storage in memory maps..49

6 Results...50

6.1 Influence of the EGR...................50

6.2 Calculation of heat release.................52

6.2.1 Estimation results of and the impact on heat release.55

6.3 Temperature estimations....................58

6.4 Equilibrium concentrations.....................60

6.5 Model NOx output....................63

6.6 The division algorithm.....................66

6.7 Interpolation results.....................69

Page 7 (80)
7 Calculations of time efficiency....72

7.1 Division algorithm........................72

7.2 Specific heat interpolation......................72

8 Discussion.....74

9 Conclusions...76

10 Future work.....77

11 References........ 79

Page 8 (80)
1 Introduction

1.1 Motivation

In todays modern society, where dependence on fossil fuel is an environmental

problem, more and more engine producers are forced to face tightening emission
legislations. Beside the traditional carbon dioxide debate there are also other
pollutants that need to be kept at low levels because of their health hazardous
impact. There are a total of four exhaust pollutants that are regulated by law:
carbon monoxide, hydrocarbons, particulate matter and nitric oxides. The first
two are easily brought down to low emission levels by simple oxidation in a
catalytic converter. On the other hand, much more work needs to be done in the
areas of NOx (nitrogen oxides) and PM (particulate matter).1

In December 2007 the European Commission proposed a new set of regulations in

the heavyduty vehicle industry to lower the amount of pollutions. The new sets
of regulations are known by name as Euro VI. The Euro VI proposal tightens the
emission standards for heavy duty vehicles in the areas of NOx, PM, HC
(hydrocarbons) and ammonia emissions. The motivation is that the road
transport is a major source of pollutants in Europe. The road transports are in
fact the largest contributors to NOx emissions in urban societies. Therefore there
is of course a need for heavy duty vehicle companies to adapt to the new
regulations and to contribute in the overall mission to lower pollution of the
environment.2

Date CO NOx HC PM
(g/kWh) (g/kWh) (g/kWh) (g/kWh)
Euro IV October 05 1.5 3.5 0.46 0.02
Euro V October 08 1.5 2.0 0.46 0.02
Euro VI October 13 1.5 0.4 0.13 0.01
Table 1.1 Regulations set by the European Commission (for the stationary
case). One should pay special notice to the tightening regulations on NOx
emissions. [2]

When talking about NOx one should keep in mind that it consists of two separate
molecular compounds, the nitric oxide and the nitric dioxide. Both compounds are
hazardous to people and animals alike. The nitric dioxide can produce ozone
according to reaction R1.1, which on ground level is a carcinogenic molecule:

R1.1

Inhalation of either nitric oxide or ozone is dangerous and in a worst case

scenario, lethal.

Page 9 (80)
At Scania CV AB, advanced engine technology is developed to fulfill the Euro VI
standards set by the Commission of the European Union. For future engine
development, there is of course an interest in having good models that can
describe the NOx formation in heavy duty engines. The purpose of this thesis was
to derive a theoretical model for NOx estimations and also to implement the
model in C for future real time executions on an FPGA. The idea is that the
model later on can be used in the area of control theory, for effective and fast
calculations of NOx concentrations in diesel feedback control applications.3

1.2 Problem formulation

In this master thesis, the investigated problem is twofold. First, a suitable zero
dimensional model must be derived that will satisfy two necessities. First and
foremost, the model must track the true NOx formation sufficiently accurate.
Secondly, the model has to be simple enough to be implemented on a real time
running FPGA. In the case of combustion in a diesel engine, the time constraints
are hard to satisfy. It will be essential to fully utilize the parallel structure in the
FPGA, as much as possible. This is described in detail in section 2.5. Besides the
parallelism, some considerations are needed when it comes to the
implementation of algorithms. Although the scope of this thesis does not include
a finished FPGA implementation, these aspects will be taken into account
nevertheless. The code has to be kept simple and avoidance of division operations
will be crucial for future practical utilization of the NOx model (explained in
section 2.5).

1.3 Methodology

There are three clear steps in this master thesis, which can be stated as follows:

1. Theoretical derivation of a NOx model.

2. Running and testing the model on a PC in MATLAB.
3. Implementing the model in C.

In the first step, a literature study is done to build up a framework for the
theoretical work which is needed to derive the NOx model. The model is then
derived using an approach that will be possible to implement later on. This does
render some approximations in the theoretical modeling process; however they
should generally be sufficient and theoretically valid. In the second step, the
paper model is implemented, run and tested in MATLAB. When the test results
are satisfying enough, the model will be implemented in C. The resulting C code
should then work as a template for a future implementation that will be done on
an FPGA.

Page 10 (80)
2 Background

There are fundamentally two key approaches in maintaining NOx at sufficiently

low levels, namely EGR and SCR.

2.1 The Exhaust Gas Recirculation

The basic concept of Exhaust Gas Recirculation is that exhaust gases can be
used to control formation rates of NOx. To see why this is the case; one has to
keep John Decs conceptual model4 of combustion in mind (see below, section 2.4).
The NOx formation takes place in oxygen rich areas with sufficiently high
temperatures. This can also be seen by a NOx formation rate approximation5 of
initial NOx formation, where also is a temperature dependent variable:

Eq. (2.1) 2

The EGR, consisting of exhaust gases (that is nitrogen, carbon dioxide, oxygen
and water), will reduce the oxygen concentration somewhat. This is because of
the fact that only a specific quantity of gas can pass during the intake stroke,
where EGR is mixed with air. Thus, the NOx formation rate will be reduced
according to Eq. (2.1).

There is also a major thermal effect because of the and in the EGR, that
will have an impact on the mean specific heat value of the charge. The reason is
the higher heat capacity of the two products, compared to the specific heat of
clean air. The combustion temperature will then be lower than it would be if only
clean air was used. The variables in Eq. (2.1) are described more extensively in
the theory section. This is not the only method to decrease NOx emissions to meet
European regulations.3 However, the EGR will be an essential part in the NOx
modeling procedure while the same cannot be said about the Selective Catalytic
Reduction.

Fig. 2.1 A schematic picture of the EGR function. Some of the exhaust gases are
recirculated in the engine. The EGR is first cooled in the EGR cooler, then
released by an EGR control unit (at a desired reference value) and finally mixed
with fresh air at the intake [6].

Page 11 (80)
2.2 The Selective Catalytic Reduction

The idea of EGR was to reduce the formation of NOx itself; however this does not
reduce NOx that is already present. To get rid of already formed NOx, one needs
to find suitable reaction mechanisms to break down the NOx into molecular
compounds that are less harmful to the environment. One common way is to add
a Selective Catalytic Reduction (SCR) in the after treatment system of exhaust
gases. A reductant, for example urea or ammonia, is used to break down the
harmful nitrogen oxide molecules. In the case of urea, this will lead to the
following reactions:

R2.1 4 4 4 6 (Standard SCR)

R2.2 2 2 3 (Fast SCR)

R2.3 8 6 7 12 (Slow SCR)

The nitrogen oxide is broken down into nitrogen and water. For the reaction to
take place there is also a need of a suitable catalyst. There are different catalysts
that could be used, where one of them is a titanium based ceramic honey comb
catalyst.7 One problem however, is that the Selective Catalytic Reduction
contributes to global warming if its not implemented wisely in the after
treatment system. When there is an overdosing of the reductant in the SCR, a
substantial quantity of laughing gas ( ) is formed. Laughing gas has 300 times
greater impact, per unit weight, on the global warming than carbon dioxide.

Fig. 2.2 A simplified schematic picture of the SCR where harmful NOx is
converted to nitrogen and water. [8]

Page 12 (80)
2.3 Models for NOx estimation purposes

The possibility to produce good models for NOx level estimations has been
researched extensively during the years. The range of complexity of the models
derived is wide. The models span from simple zero dimensional models to a
variety of advance three dimensional models. This does not necessarily imply
that three dimensional models are better to use. In fact, the advance models with
high level of complexity have not been successfully implemented for the purposes
of real time use. The authors opinion is that these complex models should rather
be seen as tools to a better conceptual understanding of NOx formation during
combustion in an engine. The conceptual knowledge gained by this could very
well be used to improve zero dimensional models that are applied in real time
systems.

The main difference between a zero dimensional model (ZDM) and a three
dimensional model (TDM) is that the ZDM does not take spatial properties into
account in describing combustion phenomena, such as NOx formation. The TDM
is a more ambitious approach to the spatial properties and its contribution to
NOx formation. For example, one could take into account the turbulence,
formation of droplets and the flow field in the combustion engine to produce a
high level TDM.9

In this thesis however, a ZDM will be implemented and run. The main reason is
that the ZDM can be sufficiently good at estimations while being a lot faster than
its TDM counterpart. It is possible to evade gruesome and time consuming
partial differential equation solvers by using a ZDM instead. The thesis ZDM
assumes two zones, an unburned and a burned zone. This concept of modeling
will be explained in depth in the theory section (see theory section, section 3).

Fig. 2.3 An example of a 3dimensional simulation model of fuel vapor,

produced with computational fluid dynamics (CFD). The CFD approach is also
used in calculations of NOx formation in a TDM. [10]

Page 13 (80)
2.3.1 NOx measurement

It is in fact possible to measure NOx levels using a NOxsensor. Nevertheless, a

NOx model is still useful because of the fact that it updates the emission levels
considerably faster than the sensor. This property will be essential in crank angle
resolved diesel feedback control. The sensor does also not work as good at cold
starts and it takes some amount of time for it to work properly. It is also a good
idea to have some kind of backup when the sensor cease to work, as it could do,
for example, due to a harsh environment (high pressure etc.). The NOx model can
also fulfill the function of surveying the NOxsensor and detecting when it starts
to fail or drift away.

In the future, with sufficiently accurate and fast NOx models, it would be possible
to remove a NOxsensor in the engine system (in a normal case, there are several
sensors). More specifically, it would be the NOxsensor right before the SCR. This
would of course decrease cost for both the manufacturer and the purchaser of
heavy duty vehicles.

Fig. 2.4 A NOxsensor used to measure NOx levels. [11]

Page 14 (80)
2.4 A conceptual model of diesel combustion

Before embarking on a mission to model the NOx formation one should of course
seek a deeper understanding of how combustion takes place in a diesel engine
and where the NOx is formed. John Dec4 and his research team presented a fairly
thorough conceptual model of the combustion process. The work was done by
using optically accessible diesel engines. In combination with laser based
techniques, they used LaserInduced Incandescence (LII) and elastic scattering
to determine distribution of liquid fuel, soot and particle sizes. To quantify local
lambdas, Rayleigh Scatter was used and a Laser Induced Fluorescence (LIF)
enabled the possibility to measure local concentrations of , hydrocarbons and
radicals. 9

Fig. 2.5 Quasisteady Diesel combustion plume as presented by Dec (1997).

Courtesy Dr. John E. Dec (Sandia National Laboratories).

As can be seen in Fig. 2.5, the NOx production takes place on the lean side of the
diffusive layer. This is because of the fact that the area is rich on oxygen while at
the same time the temperatures are sufficiently high.

Another interesting fact is that the combustion can be divided into two separate
steps. In the first step, the combustion takes place in a fuel rich and premixed
flame where, although temperatures are high, there is not enough oxygen for
NOx formation to occur. Instead of NOx, initial soot particles are formed in this
section. The second step is the combustion at the diffusive layer (orange). The
content after the combustion is directly added to the postcombustion zone,
where NOx is formed on the lean side of the diffusion flame. This fact simplifies
the calculations somewhat in the NOx modeling process.

The conclusion from Decs research is that the NOx formation takes place on the
lean side of the diffusive flame and in the post combustion zones after the end of
combustion.9

Page 15 (80)
2.5 The FPGA

Models and control algorithms in engine development are implemented in

embedded systems. This is of course due to practical reasons; a PC would not be a
suitable option because of its size. The embedded system is essentially an
information processing device suitable for products where the processor isnt the
most important part of the product. The FPGA (Field Programmable Gate Array)
is a silicon chip with a whole net of logic blocks that, depending on the purpose,
can be configured to solve a specific task, while at the same time being
reconfigurable. This means that a programmed task can be erased and replaced
with a completely new and different task. The FPGA also consists of an I/O pad
that is used for communication with the outside process.

There is a fundamental benefit gained by using the FPGA instead of an ordinary

processor. If an input signal is sent to the device the operations, programmed by
the developer, propagate through the FPGA in a parallel manner. The different
logical blocks in the FPGA are interconnected through small wires which enables
the transport of intermediate results and calculations between units.

The FPGA is best suited for tasks that use integer arithmetic and exhibits a high
degree of parallelism. In reality this means that algorithms have to be
transformed from ordinary floating point calculations to a fixed point
implementation. However, this is rarely a trivial task and in the case of models
used in this thesis, a lot of considerations need to be made to maintain a
sufficient degree of model stability and precision.12 As stated, the FPGA is
especially suited for tasks with parallel structures. Thus, to exploit the full
potential of the FPGA circuit some amount of time should be spent on dividing
the task in parallel sub tasks, where it is possible to do so.13

Fig. 2.3 The principal layout of the FPGA. The connect block works as a
switching unit. [13]

Page 16 (80)
One major setback with the Field Programmable Gate Array, in terms of
calculations for practical use, is its insufficient handling of division operations.
When using the FPGA one should avoid unnecessary divisions if possible. In
some cases this is not possible due to the fact that the divisor is an unknown
variable that changes, where an example of such a case is the division dp/p.
There are a couple of rather neat ways of solving this problem. One of them is to
use a Taylor approximation of a suitable function (in this case: ln p) and applying
basic differential calculus, where it will be possible to evade gruesome divisions.3

Fig. 2.4 Example of an Altera Stratix IV GX FPGA [14]

The FPGA can either be programmed by using Low Level Design (LLD) or a High
Level Design (HLD) approach. In the LLD, implementation is done using a
hardware description language. This is not a preferred way for software
developers because of the fact that a lot of the implementation time goes to
telling the FPGA which logic block is used for what task. The focus should
instead be on programming the functionality. For this purpose it is possible to
use the HLD approach instead. This basically means that a high level
programming language is used for behavioral descriptions which then are
translated to LLD language. Luckily, the FPGA is user friendly enough that it is
possible to use a large variety of high level languages, such as C, C++, Java
etcetera. It is also possible to use graphical programming languages such as
LabVIEW.

In this thesis, the C programming language will suffice for the implementation of
a crank angle resolved NOx model. C is often used in implementation and
development of embedded systems, and most engine developers at Scania use C
in engine software implementation.

Page 17 (80)
3 Theory

3.1 The rate of heat transfer

To derive a suitable combustion model for a diesel engine, the first law of
thermodynamics can preferably be applied. First and foremost, the heat release
can be stated in mathematical terms as following9:

Eq. (3.1)

In this case the corresponds to change in the internal energy, the heat
loss to the surrounding environment, for example the cylinder walls or due to
radiation heat transfers. The combustion also performs a desirable physical work
. Rewriting the term as = m c dT, where m equals mass and c is the
specific heat at a constant volume, gives:

Eq. (3.2) m c dT

This equation can be rewritten further by differentiating the ideal gas law, where
change in amount of substance is neglected:

Eq. (3.3)

Inserting the right side of Eq. (3.3) in Eq. (3.2) and realizing that the work done
is , gives:

Eq. (3.4) c pdV

V
Eq. (3.5) c

In an ideal gas the constant , the specific gas constant, can be calculated as
. Eq. (3.5) can then be written as:

V
Eq. (3.6) c

Differentiating Eq. (3.6) with regards to the crank angle, with , results in a
differential equation on the form8:

Eq. (3.7)

Solving the equation will give the heat released by the combustion of fuel.
Additionally, if the equation is solved, the mass of the burned fuel in the
combustion engine can be calculated by integration: .

Page 18 (80)
The is the lower heating value of the fuel. The relation between the
combustion heat and the heat produced by burning the fuel is simply
. The total fuel mass in the cylinder, at a certain time step, is
not equal to the burned mass of fuel. A term has to be added for the unburned
fuel, i.e. introducing a combustion efficiency coefficient.

If Eq. (3.7) ought to be useful pressure must be measured. This is not the case in
most commercial vehicles today. Therefore an alternative heat release equation
can be derived for NOx modeling purposes. This is just a simple rewriting of Eq.
(3.7) on a form where pressure is assumed to be the unknown variable, instead of
the heat release:

What is missing for the equation to be useful is a model of the injection of fuel
and a model of the combustion process. There are several different approaches to
model fuel injection and combustion. A rather simple model of combustion is an
assumption that the heat release is proportional to the difference between the
energy content of the injected fuel and the accumulated heat release as
. Injection can on the other hand be modeled as a dependency of
the injector needle lift. A common approach is to use a top hat model of the
injector needle lift. The approach is suggested in several papers and due to the
fact that todays development of injectors goes toward fast piezoelectric units it
seems fairly reasonable.1 However, in this thesis pressure is measured and used
for modeling of NOx formation. This because of the fact that pressure sensors are
getting commercialized and are already so on the market.15

Fig. 3.1 A top hat model (purpledashed) of the injector needle lift compared to
the real result (blue). [1]

Page 19 (80)
3.2 Calculation of NOx formation rate

3.2.1 Background

Principally, NOx can be formed in three different ways: prompt NOx formation,
thermal NOx formation and fuel NOx formation. In the prompt formation, the
nitrogen in the air reacts with unburned hydrocarbons ( ) where nitrogen
radicals are produced. The free radicals are reactants in further chemical
reactions where NOx is finally formed. The fuel NOx formation can only take
place if there are nitrogen atoms chemically bound in the fuel used for
combustion. When the bounded nitrogen reacts with the oxygen in the air, NOx is
produced. Although there are three different ways for NOx to be formed, the most
contributing mechanism to NOx production is the thermal NOx formation. In this
case the combustion provides sufficiently high temperatures for the stable oxygen
and nitrogen molecules to react and form NOx.16

A common assumption is that the NOx formed is equal to the formed. This is
however not completely true due to the fact that there are several reaction
mechanisms that will form molecules as well. Nevertheless, the assumption
is somewhat valid and practical for further modeling of NOx formation.5

3.2.2 The Zeldovich mechanism

There are many possible reaction paths where NOx could be formed. Luckily
enough Yakov Zeldovich, a Russian physicist, identified two main reaction paths,
that could explain most of the NOx formation, R3.1 and R3.2 as:

R3.1

R3.2

However, this was proven to be too insufficient and therefore a third reaction,
here R3.3, was introduced in the extended Zeldovich mechanism.5

R3.3

For the final step in the theoretical NOx formation modeling, some basic
definitions are needed.

An arbitrary chemical product formation rate can be expressed in the following

equations:

Eq. (3.9)

Eq. (3.10)

Eq. (3.11)

Page 20 (80)
 and corresponds to the products and reactants respectively. The and
are the stoichiometric coefficients of products and reactants, and . The
and are the coefficients of the rates of the forward and backward reactions
respectively. These coefficients can be expressed on an Arrhenius form9 i.e.:

Eq. (3.12) e

Where:

is a temperature dependent constant. However, this temperature

dependency is fairly weak.

is the temperature dependent activation energy.

is the universal gas constant.

Although, there are tabulated values for the different . The concentrations of
reactants and products are noted with the expression. The and
represent the corresponding forward direction and backward direction reaction
rates respectively.9

Finally, using the notations given in Eq. (3.9) Eq. (3.11), the NOx formation rate
can be written as:

Eq. (3.13)

O N N O N OH N NO NO O NO H

In the combustion engine, the formation of N is sufficiently less than the

formation of other species that it is possible to neglect it completely, i.e. 0.

The formation rate of , using the same notation as in Eq. (3.13) is:

Eq. (3.14)

O N N O N OH N NO NO O NO H

With 0 the is given by:

Eq. (3.15) NO O OH N O N NO O NO H

O N NO O NO H
Eq. (3.16) N
NO O OH

Page 21 (80)
Substitution of with the expression above in Eq. (3.13) results in:

Eq. (3.17) 2

The next step is to assume equilibrium concentrations of O, O , OH, H and at

the local temperature and pressure in the post combustion zones . This
9

assumption is due to the fact that most of the formation takes places in the
combustion zones where these conditions are met. This further implies that the
backward and forward reaction rates are equal, where denotes the
equilibrium concentrations:

Eq. (3.18)

Eq. (3.19)

Eq. (3.20)

Thus Eq. (3.17) can be expressed, on an easy form, as:

Eq. (3.21)

As can be seen in Eq. (3.18) Eq. (3.20), the only equilibrium concentrations that
have to be calculated are , and . An interesting case is when
. Then Eq. (3.21) can be simplified as follows:

Eq. (3.22) 2 2

The can be written in the form of , which results in the

expression:

Eq. (3.23) 2

One lesson from Eq. (3.23) is that the NOx formation is dependent on
temperature, oxygen concentrations and the concentration of nitrogen. This is an
expected result from the experimental research conducted by Dec et al.

If only the Original Zeldovich mechanism is used, i.e., the first two reactions, the
corresponding equation to (Eq.3.22) will then be13:

.
(Eq.3.24)

Page 22 (80)
There are reasons that one would prefer the Original Zeldovich. First and
foremost, the third reaction is relevant at temperature under 1400 Kelvin. Also,
it is a bit simpler to implement Eq. (3.24) in C code. After all, high calculation
speed is also a requirement as well as accuracy.

The and can be calculated from5 :

/
Eq. (3.23) 3.6 10

. /
Eq. (3.24) 2.129 10

The equilibrium concentration of , , can be derived by using following

reaction mechanism:

R3.4 2

This results in an equilibrium concentration:

Eq. (3.25) NO K NO O N

Finally, the NOx formation rate is calculable in Eq. (3.23) and the concentration
of is:

Eq. (3.26)

The resolution is most often given in CAD and not time. To get the corresponding
time step at each incremental CAD step, following substitution is done:

Eq. (3.27) t and t

Eq. (3.28)

The function is dependent on time through the rate per minute and the is
the angular velocity. It is now possible to calculate the NOx formation rate and
thus the NOx concentration for each time step taken or rather, for each new
change in crank angle degree.

Page 23 (80)
3.3 The two zone combustion

The theory above is not enough for a good model of the NOx formation rate.
Combustion has to be modeled in a manner that is sufficient for estimations but
also for implementation purposes. In the two zone approach to combustion, the
zero dimensional model assumes two separate zones during the combustion. In
the first zone, Zone 1 (burned zone), combustion takes place and NOx is formed.
The second zone, Zone 2 (unburned zone), is assumed to only consist of unburned
gas where no NOx formation occurs. To maintain a constant , mass flows from
Zone 2 to Zone 1 and except this mass flow, masses are maintained in both zones.
What in fact does transfer is the heat released by the combustion in Zone 1.
Optionally, even these transfers can be neglected as well, due to calculation time
issues, in a more simplified model of the combustion process.

The burned zone is however split up in two additional subzones, as illustrated in

Fig. 3.2. In the first subzone, the combustion zone, fuel and air reacts in a
combustion reaction and forms the products , , , , and . These
products are directly added to the second subzone, the post combustion zone.
Here the newly formed products are added to the products from combustion at
previous time steps. It is in fact in this zone that NOx is assumed to form. The
combustion zone is merely there to describe the heat release and temperature
increase, but the zone that utilize the Zeldovich mechanism is in fact the post
combustion zone. This decoupling of NOx formation and combustion greatly
simplifies the modeling process and also the code implementation procedure.

Fig. 3.2 The concept of the two zone model with an additional subzone, the
post combustion zone.

Page 24 (80)
This is however still not enough for a useful combustion model, and therefore
some other assumptions and statements are added accordingly9:

(1) Gases in the cylinder can be treated according to the ideal gas law.
(2) The sum of the volumes of Zone 1 and Zone 2 has to equal the volume of
the cylinder.
(3) Temperature differences between Zone 1 and Zone 2 depend solely on
ROHR and a difference in heat transfers to the cylinder walls, due to
different temperature gradients.
(4) The heat losses to cylinder walls can be described with the Woschni heat17
transfer equation as stated:

Eq. (3.29)

. . . .
Eq. (3.30) 3.26

= cylinder pressure
= average gas velocity
= mean gas temperature
= heat transfer coefficient
= cylinder bore
= cylinder surface area

(5) The pressure is uniform in the cylinder. This means that Zone 1 and
Zone 2 has exactly the same pressure at a given time.
(6) The temperature is the same in both zones, just before the injection of fuel
and combustion.
(7) Concentrations of residual gases are the same in both zones.

The burned fuel mass is given in Eq. (3.8). The lambda value in Zone 1 is known
and therefore the total mass in the zone can be calculated:

Eq. (3.31) 1 , stoichiometric air-to-fuel ratio

The temperatures in Zone 1 and Zone 2 also changes due to compression and
expansion. Using the ideal gas law and the fact that the product is constant,
an expression for , changes due to compression, and ,changes due to
expansion, is derived.

Page 25 (80)
The temperature change due to compression (or expansion) is given by the
relation13:

Eq. (3.32)

This implies a due to isentropic compression (expansion):

Eq. (3.33) 1

There is also a temperature change because of the combustion itself. The

was previously derived. However, this will not suffice, because one has
to take the heat losses into account. This will give a :

Eq. (3.34)

And the temperature change given by the ROHR is (where the air can be diluted
with EGR):

Eq. (3.35)

The heat transfer from the warmer gas in Zone 1 to the colder Zone 2 contributes
to a temperature gradient from the first zone to the second zone.

Eq. (3.36)

Where

Eq. (3.37)

The variable is the maximum temperature difference between the two zones
and is the pressure of a motored cycle at CAD . SOC and EVO are the
Start of combustion and Exhaust valve opens respectively.9 Although, this
approach is rather gruesome and therefore another path can be taken, if there is
a desire to include heat transfers between zones in the model. The alternative
calculation path is derived in section 3.5.8 (see below).

Finally, the temperature in Zone 1 can be calculated at each , the mass of

burned fuel and air per is now known, and the pressure can be measured. The
volume of Zone 1 is then given by the ideal gas law as:
RT

Page 26 (80)
This volume is used for calculations of gas concentrations in Zone 1, where the
NOx formation is occurring.

3.4 Stoichiometric calculations

The chemical reaction of combustion in a diesel engine is given by18:

R3.5 3.773 1
3.773

In some cases there is also an EGR function in the diesel engine that will have an
impact on the reaction mechanism. However, this is a fairly simple addition to
R3.5:

R3.6
3.773
1 3.773

The amount of substance of the fuel, , is:

Eq. (3.38)

The is the specific molar mass of the fuel and the mass of fuel burnt is
previously known. Every concentration needed in Eq. (3.21) can now be
calculated and used to determine the NOx formation at each CAD step, following
the assumption that is constant in Zone 1. The content of the combustion zone
is directly added to the post combustion zone (see section 2.4). The post
combustion content produced at the current time step is added to the sum of
contents produced in the previous time steps and that will result in the total post
combustion zone at current time. Despite of this, some considerations need to be
made regarding the temperature change when fresh combustion products are
added, and this is due to a difference in temperatures between the contents.
When the temperatures are determined it is possible to do simple calculations of
the different combustion product concentrations needed in the NOx formation
rate formula, assuming homogeneous mixtures in the zones.

Page 27 (80)
3.5 Calculations and estimations of parameters

3.5.1 Estimating and cp

To successfully implement the above theory in algorithms, parameters have to be

known in advance. For example the in the Rate of Heat Releaseequation is not
trivial to determine at all. One way is to omit heat transfers and mass flows
during the adiabatic compression phase.3 The heat release equation can then be
simplified as follows:

Eq. (3.39) 0 and 0

Thus

Eq. (3.40) 0

This gives an estimate of :

Eq. (3.41)

The combustion and expansion stroke will of course change the estimate.
Assuming small deviations, this impact can be neglected.

The parameter , needed in temperature calculations, can be calculated using

NASA Glenn tabulation of specific heats19. Another approach to calculate the
specific heat is to utilize the estimate of . First and foremost the specific ideal
gas law constant can be stated as:

Eq. (3.42)

Where R is the ideal gas constant and is the molar mass for the mixture of
gases. Further, the ratio between and is known through the estimate.

Eq. (3.43)

Substitution of in Eq. (3.42) gives:

Eq. (3.44)

The downside is that the approach assumes a fairly good estimate of and the
mass contained in Zone 1, which generally cannot be guaranteed. Therefore, this
method is neglected and not presented in the result section. It is much more
convenient to use NASA tabulations and interpolations of thermo chemical data.

Page 28 (80)
3.5.2 Calculation of molar masses and SAFR

The is the molar mass of Zone 1. This of course is also a parameter that
has to be calculated on beforehand. One fitting way is as follows:

Eq. (3.45) 1

The mass of Zone 1 is previously known, given by Eq. (3.31) and thus this result
can be used in Eq. (3.44). Additionally, a calibration parameter could be added to
account for EGR. The , and can be calculated with some help from
R3.5 and R3.6 as:

Eq. (3.46)

.
Eq. (3.47)
.

is the atomic mass of atom .

This will give an air molar mass and a :

. .
Eq. (3.48) 28.96 /
.

Eq. (3.49)

3.5.3 The average gas velocity estimate

The average gas velocity , needed in Eq. (3.30) can be estimated by

calculating12:

Eq. (3.50)

The variable is the mean piston speed while is the volume displacement.
Also, the , and is the temperature, pressure and volume at valve inlet
closing respectively. Pressure is the current pressure in the cylinder and is
the motored cylinder pressure, at the same CAD as pressure . and are
empirical constants. To simplify the estimate of , the second term is omitted in
this thesis. Thus the estimate of is:

Eq. (3.51)

Where the constant can be determined emperically9:

6.18 at the gas exchange period

2.28 at compression, expansion and combustion period

Page 29 (80)
3.5.4 Approximation of equilibrium concentrations

When calculating the radical concentrations in Eq. (3.23) and Eq. (3.24), one
assumption has to be made. The equations are essentially derived from the
equilibrium equations:

Eq. (3.52)

Eq. (3.53)

The amount of substances of oxygen and water is known through equilibrium

calculations (see section 3.6) and the volume is known from the ideal gas law. The
concentrations at equilibrium for Eq. (3.52) will be, by applying basic chemistry:

Eq. (3.54)

The following approximation is made, based on the fact that the amount of
substance of the oxygen radical is far less than the amount of substance of the
oxygen:

Eq. (3.55)

The same assumption is made for Eq. (3.53), i.e. the amount of substance of the
radical is less than the same, for oxygen radicals and water molecules.

3.5.5 A description of the EGR dynamics

If the reaction mechanism in R3.6 ought to be useful, the EGR inflow has to be
known. The approach to solve this problem is to establish a model for mass flows.
One way to do so is described in Fig. 3.3.

Fig. 3.3 Schematic picture of gaseous mass flows to the engine. The variables p
and T are directly measurable.

Page 30 (80)
The intake pressure and temperature is measured and the ideal gas law is
applicable as long as most of the in the EGR is not liquefied in too large
quantities. The equation of interest, in this case, is3:

Eq. (3.56)

Where is the density, is the engine displacement, is the volumetric

efficiency degree, the N is the rate per minute and the is the intake mass
flow.

Although, the problem is that the density is unknown due to the fact that air is
mixed with EGR. A suitable trick is simply to rewrite Eq. (3.56) as:

Eq. (3.57)

Now, using the time derivative of the ideal gas law gives (at constant inlet
pressure and temperature):

Eq. (3.57)

And finally,

Eq. (3.58)

Knowing that the air mass flow is measured, the can be calculated:

Eq. (3.59)

The amount of substance of the EGR in the combustion zone can easily be
determined using assumption (7) in the twozone model (see section 3.3).

Eq. (3.60)

The parameter is given by the reaction mechanism in R3.6. It should

be noted that the fractions of carbon dioxide and water in the recirculated
exhaust gas need to be determined. However, this is a trivial task that is possible
to solve by calculations of particular fractions in the exhaust gas in R3.5. This
possibility is based on the fact that the water and carbon dioxide in the EGR
originates from the combustion itself. The combustion reactions are fully
described by the reaction mechanism in R3.5. The only thing that has to be
known is how much air is added during the intake stroke.

Page 31 (80)
It is possible to calculate the amount of air taken in by integration:

Eq. (3.61) 1 %

is the total intake volume at a volumetric efficiency, .

The proportion between oxygen and nitrogen in the air is 1:3.773, as can be seen
in reaction R3.6.

One last step in the EGR is to understand the coefficients , ,

and . A percentage of the combustion products in R3.5 will be in the EGR,
this percentage will be labeled EGR%. Hence, if 20% of the exhaust gas is re
circulated then % 20%. Given this information, it is possible to quantify
the coefficients during different running conditions, where the injection of fuel
may very well vary.

The EGR dynamics are modeled assuming a homogeneous EGRtank, where a

percentage of the exhaust gases enter. The composition of the EGR is fully
explained by the reaction mechanism in R 3.5. The local value will still be
maintained around 1, which is a necessity in this NOx model (see section 3.3).
Hence, the resulting principal equation in the combustion zone (see Fig. 3.1),
during stationary conditions, will be:

Eq. (3.62)

1 % %

It is although rarely the case that the fuel consumption remains unchanged. The
variation in fuel consumption will render different carbon dioxide, water and
oxygen concentrations. However, treating the EGR tank as a well stirred tank of
gas, this does not introduce any calculation problems. The amount of substances
and, assuming a well stirred tank, concentrations can still be calculated, utilizing
a global lambda . Oxygen in the EGR is assumed not to react with fuel. This is
purely a mathematical trick due to a desire for algebraic simplicity.

Eq. (3.63)

1 % /4 /2 3.773 %

1 3.773
4 2 2
1 % 1
2 4 2
3.773 %
4 2
1 3.773
4 2 2

Page 32 (80)
A comment is needed on the local lambda. In a case where no EGR is added, the
lambda is easy to understand. If EGR is added, the assumption in the theoretical
model is altered. As previously stated, the EGR has both a thermal effect and an
effect on the oxygen concentration. In Eq. (3.63) this is clearly seen as a de facto
decrease in the amount of substance ratio between oxygen and the fuel.

This could also be seen as a change in local combustion lambda by a factor

1 % % 1 , which will result in a oxygen concentration
decrease.

The EGR will have an influence on the oxygen concentration and specific heat
value during the combustion phase. The easiest way to account for EGR influence
on the specific heat value is to calculate the specific heat capacities for each
species in the gas mixture. The average specific heat for the gas mixture can then
be calculated as:

Eq. (3.64)

The specific heat of each species is calculated using NASA Glenn interpolations of
thermo dynamical data.

The exhaust gas content, produced each cycle, is added to the EGR tank in
following way, if seen as a function of time:

Eq. (3.65)

The equilibrium concentrations, relevant for NOx formation, can be determined

using the fact that the split up of the EGR, in the two different zones, will be:

Eq. (3.66)

Eq. (3.67)

Now the equilibrium concentrations in the burned zone can be calculated in

diesel engines that utilize an inner EGR function.

What is not seen in Eq. (3.63) is a delay due to the fact that it takes some time for
the recirculated gases to really recirculate. This fact is hard to formulate using
chemical nomenclature, but mathematically and algorithmically this is fairly
easy, assuming a well stirred tank of EGR, at a mean global lambda .

Page 33 (80)
3.5.6 The cylinder geometry

The cylinder geometry will be of major importance in the model. The volume,
volume angle derivative and the surface area are, respectively18:

Eq. (3.68) sin 1 cos sin

Eq. (3.69) sin 1

Eq. (3.70) 2 1 cos sin

is the bore diameter, is the ratio between the conrod length and half the
stroke. and are compression and displacement volume respectively.

The piston position and velocity will be functions of the CAD. These relations are
given, with some approximations18, by following equations:

Eq. (3.71) 1 cos 2 sin

Eq. (3.72) sin 2 sin 2

The and are the crank radius and rod length respectively. The angular
velocity is a function of the rate per minute, where is the frequency and is
the angular velocity in rpm:

Eq. (3.73) 2 2

With the above equations, the geometry of the cylinder and the cylinder motion
are completely described.
-3 -3
x 10 x 10
3 2
dV/dfi (m3/CAD)

1
Volume (m3)

2
0
1
-1

0 -2
-400 -200 0 200 400 -400 -200 0 200 400
CAD CAD

0.08
Surface area (m2)

0.06

0.04

0.02

0
-400 -200 0 200 400
CAD

Fig. 3.4 Geometric calculations with the parameter values:

= 127 mm, = 154 mm, and 255 mm.

Page 34 (80)
3.5.7 Mass relation between the two zones

One interesting issue is the contribution of the unburned zone to the burned zone
in terms of mass flow. First of all, the initial amount of gas in the cylinder before
compression is:

Eq. (3.74)

Mass flows from the unburned zone to the burned zone in a way that keeps the
local constant at 1. The relation between the amount of substances in Zone 1
(burned zone) is:

Eq. (3.75) ,

The equivalent relation in Zone 2, the unburned zone, is:

Eq. (3.76) ,

Rewriting Eq. (3.76) by using Eq. (3.75) gives the amount of air in the unburned
zone, zone 2:

Eq. (3.77)

The EGR in Zone 2 is easy to calculate when statement (7) is applied (see section
3.3 above). The total substance of EGR is known from integration of Eq. (3.59),
thus:

Eq. (3.78)

Where is the volume of the unburned zone and the is the total volume at
a particular CAD. The volume of the burned zone, , was previously calculated
from the ideal gas law. The fact that volume is conserved gives:

Eq. (3.79)

Eq. (3.80) 1

All the variables on the right side of Eq. (3.80) are known through previous
calculations. The mass flow from the unburned zone to the burned zone is now
completely known. The amount of air, relevant for the NOx formation, can thus
be calculated, assuming 21% oxygen in the air content.

A comment on is needed in this case. The regular takes air/fuel mass ratio
into account. However, concerns the equivalent amount of substances instead.

Page 35 (80)
Nevertheless, is simple to calculate from the reaction formula stated in R3.5:

Eq. (3.81)

3.5.8 Heat transfer between zones

One of the last remaining puzzle pieces in the modeling part is the heat transfer
from the burned zone to the unburned zone. The Eq. (3.37) and Eq. (3.38) are not
on a desirable form for an FGPA implementation. The heat transfer can instead
be calculated from the fact that the unburned zone, in a normal case, would get a
temperature rise from compression only. The isentropic compression would result
in following, with a resulting temperature and compression pressure :

Eq. (3.82)

However, the result will be different because of the heat release in the burned
zone and a pressure rise will occur:

Eq. (3.83)

The temperature difference will then be:

Eq. (3.84)

In fact, this equation can be simplified even further by realizing that is

(see Eq. 3.74).

Eq. (3.85)

The amount of heat that is needed to warm up the content in the unburned zone
is thus:

Eq. (3.86) ,

If the EGR contribution is not neglected (i.e. the air is diluted with combustion
products from previous cycles), the mass in the unburned zone is approximated
and the , specific heat at constant pressure, of air is calculated using tabulated
specific heats. If one would want to calibrate the model due to the
approximations, especially of global lambda, an empirical constant could be
added in Eq. (3.86):

Eq. (3.87)

Page 36 (80)
The parameter would have to be empirically tabulated for a successful
implementation, or determined by NASA tabulation of specific heats. Because of
this tuning problem and the fact that the calculation would be delayed one step
due to an unknown unburned mass, a common assumption in NOx modeling is to
omit heat transfers between zones. Instead the unburned gases are assumed to
change temperature solely due to isentropic compression and expansion. This
assumption implies .

3.6 Calculating equilibrium concentrations

The combustion products in R3.5 will, unfortunately, not be the only ones. It is
fairly well known that combustion also produces carbon monoxide in the gas
water shift reaction. The two following reactions have a huge impact on the
concentrations in the post combustion zone9:

R3.7 H

R3.8 CO

The equilibrium, for R3.7 and R3.8, will be at:

Eq. (3.88)

Eq. (3.89)

If the amount of substance, in relation to one mole fuel, of , , , ,

and are called , , , , and , respectively, a system of nonlinear
equations can be derived, using normalized pressure. If Eq. (3.88), Eq. (3.89) and
atom balances are utilized, the following nonlinear equations will hold to be true,
where is the normalized pressure9:

Eq. (3.90)

Eq. (3.91)

Eq. (3.92)

Eq. (3.93) 2

Eq. (3.94) 2 2 2

Eq. (3.95)
.

Page 37 (80)
The equations are simple to solve with for example the fsolve function in Matlab.
However, if the objective is to implement an algorithm in an embedded system,
then the fsolve function is not a suitable approach. In this thesis a modified
NewtonRaphson method was successfully derived. First and foremost Eq. (3.90)
and Eq. (3.91) were slightly rewritten:

Eq. (3.96)

Eq. (3.97)

Five functions were then defined and differentiated symbolically to form the
Jacobi matrix ( , , , , , fuel and ):

Eq. (3.98)

Eq. (3.99)

Eq. (3.100)

Eq. (3.101) 2

Eq. (3.102) 2 2 2

For practical reasons, the resulting Jacobi is omitted in this report (see Matlab
code instead). The resulting algorithm is given as:

Eq. (3.103)

The term, a relaxing square matrix, is needed because of the often large
magnitudes involved, especially in Eq. (3.98). A unit matrix would cause
instability and the values would diverge to infinity, no matter how close the
initial guess is set to the true value. is added to guarantee stability, therefore
two values in were set to be substantially less than 1. This does inflict some
convergence rate issues. However, there are reasons to believe that the above
algorithm can be improved substantially by fitting at certain points of interest.
Another possible approach would be to make dependent on the magnitude of
.

The equilibrium constants, and , are calculated with:

Eq. (3.104) ln

The is the Gibbs free energy for species i, at atmospheric pressure (1 bar).
These values are calculated at each temperature using tabulated interpolations
of entropy and enthalpy, that can be found in NASA Glenn thermo chemical

Page 38 (80)
tables.19 The temperature range of interest, in the burned zone, is around
1500 3000 K.

The Gibbs free energy can then be calculated as:

Eq. (3.105)

is the entropy and is the enthalpy of the species.

In the case of , another interpolation4 can be used:

Eq. (3.106) ln 2.743 1.761 1.611 0.2803

Now the system described by Eq. (3.90)Eq. (3.95) can be solved and the
equilibrium substances determined, when nfuel is known.

3.7 Adding products to the post combustion zone

Each time step the combustion products produced at the current time are added
to the combustion products that were produced at previous time steps. This can
be stated as:

Eq. (3.107)

is the mass of the total amount of combustion products, is the added

contribution at the current time step, t. This will of course have an impact on the
temperature of the content. The new and warmer product mass, , will heat up
the rest of the post combustion zone.

The calculations can be rather complicated if an assumption is not made about

the values in question. For simplicity the is assumed to be the same for the
content given by and the content in the mass . When the heat exchange
in the post combustion zone has been considered, the ideal gas law will be
applicable. The resulting temperature, when the new content is added to the
combustion zone, is:

Eq. (3.108)

The old content will be cooled down due to isentropic expansion during the
combustion and expansion phase. This is of course an approximation of the
physical reality which assumes that there are no heat losses to the unburned
charge.

Page 39 (80)
A good question to ask oneself is if it is reasonable to assume only one burned
zone in the NOx model. What if one would assume a new separate zone at each
additional combustion step taken (CADstep), instead of only one where new
components are simply added to the old ones? Under the assumption that the
pressure is instantaneously uniform in the cylinder and zones ( is equal to
the amount of time steps between start and end of the heat release) are used, a
calculation can be done. Application of the gas law, at a pressure and at time ,
gives:

Eq. (3.109)

is the total volume of the burned zones created at time .

Rewriting Eq. (3.109) gives the equation:

Eq. (3.110)

The mean burned zone temperature at step , , , is then:

Eq. (3.111) ,

Remembering that there is a temperature drop (increase) due to isentropic

expansion (compression), the temperature of an element i combusted at time ,
, , can be expressed:

Eq. (3.112) , ,

Inserted in Eq. (3.111) this renders a , :

Eq. (3.113) , ,

At time a new zone is created, and when the gas law is applied the result
will be a recognizable :

Eq. (3.114) , ,

This expression is in fact very similar to the one used in the single burned zone
model, Eq. (3.108), considering that the initial condition, for the temperature of
the first element combusted, is:

Eq. (3.115) , ,

Page 40 (80)
 Hence, if the pressure is uniform throughout the cylinder and combustion occurs
at a constant or nearconstant local lambda, nothing essential can be gained by
keeping the zones apart, except for a higher resolution of the temperature
distribution. Thus, the single burned zone used in this thesis is fairly relevant
and also easy to implement and keep track of.

This does however not imply that a multizone model is worthless, even if the
pressure is uniform and lambda constant. Because, by studying Eq. (3.26), one
can see that the NOx formation rates have a nonlinear relation to temperature.
A mean temperature will give a different value of NOx formation rates than a
multizone model with distributed temperatures. However, a single burned zone
seems feasible even when temperature nonlinearities are considered. Especially
considering the fact that zones most probably do mix during combustion. The
burned zones do in fact also mix with the air and EGR content in the unburned
zone. This is the so called air dilution. However, this effect is often neglected in
zero dimensional models due to the inherent complexity of describing such a
phenomena in a physically accurate and time efficient way.9

Page 41 (80)
4 Implementation in MATLAB

Theoretical work aside, the model must also be implemented in an efficient and
easy way for future use in practice. As stated in previous sections, the model
should also be useful when ran online. Therefore, even if it works in MATLAB,
the model has to be kept simple with regards to time consumption. The model
was scripted for diagnostic purposes, i.e. a complete data set was run from
beginning to end and results of NOx concentrations were acquired. In this
section, the final MATLAB implementation of the NOx model is briefly described.

4.1 Heat release calculations

The MATLAB code utilize a fairly trivial solution algorithm to the heat release
equation. A function ROHRSolver was implemented for this purpose.
The rate of heat release was calculated with Euler discretization of the
derivatives in Eq. (3.7).

Eq. (4.1)

Volumes and volume derivatives are calculated analytically with Eq. (3.69).
However, the pressure has to be differentiated with an approximation. The
convective heat loss utilizes the Woschni heat transfer equation22, as previously
given. can either be set to a fixed value or approximated with the algorithm
presented in section 3.5.1. To keep things simple, gamma was fixed at 1.3, which
proved to give a fairly good fit in the heat release solution.

4.2 Temperature calculations

A function calcT calculates and returns the temperature in the burned zone. The
calculations assume that the temperatures of the reactants are equal to the
unburned zone temperature. A simplified assumption of the unburned zone
temperature is that it only changes due to isentropic compression and expansion:

Eq. (4.2)

NASA interpolations of specific heats were implemented in the code and used to
calculate the temperature in the burned zone:

Eq. (4.3)

Finally, the new combustion products were then added to the total burned zone.
This will result in a temperature change in the burned zone. The easiest way to
account for this is with the usage of temperature averaging between old content

Page 42 (80)
and the newly added component (see Eq. (3.107)). However, this approach
assumes temperature constant specific heats, as previously stated.

Isentropic expansion will also result in a temperature decrease in the burned

zone that cannot be neglected. Luckily, this is implemented in the same manner
as the temperature change in the unburned zone, as stated in Eq. (4.2). It is in
the temperature calculations that the EGR contribution is taken into account.

4.3 Equilibrium calculations

When temperatures and masses are known, it is fairly straightforward to

calculate the equilibrium concentrations. The equations Eq. (3.98) Eq. (3.102)
are iteratively solved for a given set of , , , . is set to 1, is measured
pressure, while and are calculated according to section 3.6. A MATLAB
function burnedZone executes the procedure.

4.4 Zeldovich mechanism

Lastly, a function ZeldovichMechanism calculates the radical concentrations and

the NOx concentration formation rate. Integration according to Eq. (3.26) returns
the NOx concentration at a given CAD. No additional post processing of data is
done, except for a trivial recalculation of NOx concentration, from mole/cm3 to
ppm.

4.5 C code implementation

The C code implementation follows the same principal layout as the MATLAB
implementation. However, some other considerations are taken into account, as
described in detail below (see section 5).

Solve heat release equation

Calculate zone temperatures

Calculate equilibrium
concentrations

Calculate NOx formation

Fig. 4.1 Flow chart of the MATLAB implemented model.

Page 43 (80)
5 Implementation considerations for FPGA code

5.1 The division operation

It is a well known fact that an implementation in FPGA should avoid division

operations. A study of the equations in the theory section (see section 3) implies
an appalling fact. There are a lot of divisions in some way or another that cannot
be avoided by simply neglecting them. One way to overcome this problem is to
implement algorithms that approximate the divisions.

5.1.1 Thesis algorithm

In this thesis such an algorithm was produced. The algorithm can be split up in
two separate steps, where in some cases the first step is sufficient enough. For
demonstrative purposes, following division will be considered in the first step:

Eq. (5.1)

The dividend and divisor, and respectively, are known. The pressure is
determined via the ideal gas law which includes a rather simple division. The
first step is to rewrite Eq. (5.1):

Eq. (5.2)

The is the engine cylinder volume, where . Now, by using the fact that a
division on the form 1/(1 ) is:

Eq. (5.3) 1

This is however only true when 1 and 1.

The equation Eq. (5.2) can easily be rewritten on the form that is given by
Eq. (5.3):

Eq. (5.4)

This approximately equals:

Eq. (5.5) 1 1 1 1 1

The approximation given by Eq. (5.5) is true due to the fact that 1 / 0 and
1 / 1 . The parameter is known beforehand and thus the division
1/ can be set a priori. The introduction of , on a first sight, might be thought
of as a pure mathematical trick to ensure that Eq. (5.3) holds. However, this can

Page 44 (80)
also be seen as an introduction of a damper which will have an effect on the
convergence rate of the algorithm.

First and foremost, the algorithm can be rewritten on an arbitrary form as:

Eq. (5.6)

The divisor is known to be constrained by a . This could for example be the

maximal cylinder volume, in the case of the ideal gas law, as previously stated.
ensures the fulfillment of Eq. (5.3). Then the division, as stated in Eq. (5.6),
is approximately:

Eq. (5.7) 1 1 1 1 1

The term is a rest term and is set a priori as:

Eq. (5.8)

However, as previously stated, the is a damper. This will have a huge impact
on the convergence rate. If is sufficiently close to the true , the algorithm
will be close enough by only using terms up to second degree in Eq. (5.7).

The good thing is that may very well be updated to be close to the real solution,
due to the intrinsic mechanical work of a combustion process (compression and
expansion). The simplest way to see this is to study the compression phase of
combustion and what this implies for Eq. (5.1). A good initial damper would be
the volume , but this will only be close to at initial states, when
compression starts. However, if the division 1/ is calculated where is
sufficiently close to , the damper at next time step can be set equal to 1/ .
This is due to the fact that the volume will decrease at compression and thus
. The total volume, , is always known a priori as:

Unfortunately, the divisor will not decrease in all cases during the combustion
and some divisors will even lack basic physical interpretation. Therefore a second
step can be added that utilize a Newton Raphson method. The idea is to use the
previous algorithm, as given in Eq. (5.7), to get an initial guess. Assuming that
the guess is , the true value is and the divisor is , the equation to be solved is:

Eq. (5.9) 1

Page 45 (80)
The function that will be used in the Newton method is written as:

Eq. (5.10) 1

Eq. (5.11)

The corresponding Newton method will be:

Eq. (5.12)

One problem in Eq. (5.12) is the division 1/q, which is in fact the division to be
solved. Approximating 1/q with gives:

Eq. (5.13) 2

At the first step, the guess is used, . The Newton Raphson method will
unfortunately be sensitive to poor initial guesses, therefore the first step is
essential in the implementation to get a good guess. It is not a good idea to use an
arbitrary guess. However, this method implies a large time consumption and is
thus not very time efficient. It is a far better idea to sometimes sacrifice accuracy
in favor of better calculation speed.

5.1.2 Modified thesis algorithm

There is also a possible alternative approach to a division calculation algorithm

which could be suitable for a C code implementation. Assuming that all float
point numbers are converted to fixed point numbers, using the representation 2 ,
the most significant bit becomes crucial in this modified algorithm. The most
significant bit of a number ( ), is the bit position having the greatest value. If
is calculated fast enough, the algorithm could be stated as:

Eq. (5.14) 1 2 2 2

Eq. (5.15) 1

The is the of number . This algorithm can handle all possible

divisions if the of is known. Eq. (5.14) can be further simplified by
algebraic calculations:

Eq. ( 5.16) 1 3 2

The bit shift will limit 2 , 2 which corresponds to [0.5, 1] in a float point
representation. Hence, the stability criteria is satisfied and the convergence rate
will be sufficient. A possible approach to find the is to use the parallel
structure of the FPGA. The can then be found fairly quickly by
implementing the algorithm suggested by Warren in Hackers Delight.20

Page 46 (80)
In a case when the divisor is 1 in a division operation A/B, a slight
modification can be done to the thesis algorithm to make it far more precise in
terms of convergence. Remembering that the fundamental part of the algorithm
is to rewrite it on a form 1/(1-x), following trick can be done in this case which
proves to be especially accurate for small divisors:

Eq. (5.17)

This can in fact be simplified as follows:

Eq. (5.18) 1 1 1 1 1 1

As can be seen in Eq. (5.18), the algorithm will be accurate when:

1
1
1

This of course implies that 0, which will be especially true for divisors of
small magnitudes.

Page 47 (80)
5.2 Tabulation of data

When running the model in MATLAB offline there are no time constraints that
have to be met. However, the case is quite different online on an FPGA. For
example, the non linear equations cannot be solved sufficiently fast. Also the
NOx formation equations, as given by the Zeldovich mechanism, are nonlinear
and gruesome to calculate. The easiest way to approach this particular problem is
to tabulate the solutions for different pressures and temperatures and to utilize
interpolations. In this sub section this is explained in short.

A multi dimensional least squares method is an appropriate approach in the

case of interpolating tabulated data.21

5.2.1 The least squares method

If there is a set of output data , , ,, with a corresponding input set

, , ,, and , , ,, , a least squares interpolation is
possible. For example, the output data may be the specific heat capacity at
constant pressure, and the input data could be EGR percentage and temperature.
To find a least square fit, the error function has to be minimized:

Eq. (5.19)

Where , , , , , and , , , , , 1 . Minimizing the

error is equivalent to solving the equation:

Eq. (5.20) 2 0

The indicates that the polynomial fit to Eq. (5.20) is parabolic. This will be a
sufficient degree of complexity for the purposes of this thesis. In this case the
equation can then be rewritten in a matrix form as:

Eq. (5.21)

For example, it is now easy to find a fit for the specific heat as a function of EGR
and temperature. The resulting polynomial will only require multiplications
when implemented in C. Multiplication is a lot faster operation than actually
calculating the specific heat from the NASA interpolations. These calculations
would require both divisions and tabulations of NASA thermodynamical data,
which are both time consuming and memory inefficient.

Page 48 (80)
5.2.2 Functions of a single variable

In other cases, where there is only one dependent variable, it is possible to use a
simple linear fit that can be retrieved from suitable MATLAB functions. This is
the case for the isentropic compression term given by:

Eq. (5.22)

If the sampling rate is known, discretizations of geometrical properties such as

surface area, volume and volume derivative, are fairly easy to tabulate and store
a priori.

5.2.3 Data storage in memory maps

One of the most common ways to solve the issues involving nonlinearities is to
use the possibility to retrieve and store data. This data can then be accessed by
an FPGA very fast.13

Page 49 (80)
6 Results

6.1 Influence of the EGR

As previously stated, EGR is used to lower temperature conditions in the burned

zone during combustion. This will have an effect on the NOx formation in such a
manner that the conditions would be poorer for NOx to form. The temperature
calculations reveal why this is the case. The increase in temperature can be
calculated with basic physics as dT/mcp. The effect of EGR on oxygen
concentration in the combustion zone is however omitted in the validation.

The increase of EGR during the intake stroke will have an effect on the specific
heat value. If only , , and are assumed to be present in the exhaust
gas, the NASA specific heat interpolation for each species will be:

Eq. (6.1)

The specific heat for carbon dioxide and water will be somewhat higher. This
implies a smaller increase in temperature if EGR is present during the
combustion stroke. Eq. (6.1), the NASA polynomial used in the tabulation, of
course holds for other molecular compounds as well.

at 500 K at 800 K at 1500 K at 2500 K

N2 29.6 31.4 34.8 36.6
O2 31.1 33.8 36.6 39.2
CO2 44.6 51.5 58.4 61.4
H 2O 35.3 38.7 47.7 57.9
Pure air 29.8 31.9 35.2 37.1
Air+10% EGR 30.1 32.2 35.6 37.6
Air+20% EGR 30.4 32.5 36.0 38.1
Table 6.1 Values of specific heats at different temperatures for N2, O2, CO2 and
H2O. The values are given in J/Kmole. The Air + EGR values are calculated for a
global mean = 1.3.

The interesting thing however is not the specific heat values of the species
individually, but rather the averaged . The averaged will be:

Eq. (6.2)

1
1 % 3.773 %
4.773

Eq. (6.3)

1 3.773
.

Page 50 (80)
The dependency of the mean specific heat on temperature and EGR percentage is
shown in Fig. 6.1 and Fig. 6.2. The main thing that can be concluded is that the
EGR will have a non linear influence on the mean specific heat for higher
temperatures. On the other hand, the level curves shown in Fig. 6.2 imply an
almost linear relation between EGR and temperature for lower temperatures.

Fig. 6.1 Mean specific heat capacity as function of temperature and EGR.

0.9

0.8
Exhaust gas recirculation / %

0.7

0.6

0.5

0.4

0.3

0.2

0.1

500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600
Temperature / Kelvin

Fig. 6.2 The corresponding level curves. Notice the nonlinear curves at higher
temperatures.

Page 51 (80)
The EGR will also have an impact on the specific heat ratio that will be relevant
in heat release calculations. However, if is estimated during the compression
phase the EGR will be taken into account. Because of calculation time issues,
was set to a constant value. This value will approximately be 1.3 at all times. The
deviation from this value will purely be due to temperature and EGR inflow
changes. The calculations described by Eq. (3.41) assume only one zone. This is
due o the fact that the heat release calculations also assume only one zone.

6.2 Calculation of heat release

In the calculation of heat release, pressure data and volume calculations must be
used. The volume calculations are purely theoretical (see section 3.5.6) and quite
accurate, although approximative. However, the same cannot be said about the
pressure trace used in the calculations. The combustion process itself will initiate
acoustic phenomena in the cylinder, with a certain frequency, due to abrupt
pressure changes. This frequency will be relatively high (3.7 kHz) and will cause
problems when the differential is calculated (see Eq. (3.7)). Differentiation of
noisy measurements is not a good idea, as can be seen in the calculation of net
heat release in Fig. 6.3.

Unfiltered HR data

Estimated HR
Measured HR
400

300

200
J/CAD

100

-100

-100 -50 0 50 100 150

CAD / degrees

Fig. 6.3 Calculations of heat release when data is unfiltered.

Several low pass filters were tested in MATLAB. The filter that proved to be the
most satisfying was a SavitzkyGolay smoothing filter22. This particular filter
performs a polynomial regression on a data collection in local sections. Another
interesting fact is that it does not matter when data is filtered. The filtering can

Page 52 (80)
be done directly on pressure trace data or on the calculated heat release. The
filtered result is shown in Fig. 6.4.

Filtered HR data

Estimated HR
350
Measured HR

300

250

200
J/CAD

150

100

50

-100 -50 0 50 100

CAD / degrees

Fig. 6.4 SavitzkyGolay filtered heat release data.

As can be seen, the filter is very effective in eliminating high frequency noise
without deteriorating the true heat release. The filter could also be applied to the
pressure trace directly, before any calculations of the heat release.
6
x 10
11
105
10.5
100
10
95
9.5
90
9
Pressure / bar

Pressure / Pa

85
8.5
80
8
75
7.5
70
7
65
6.5
60 6

5 10 15 20 25 30 -60 -40 -20 0 20 40 60 80

CAD CAD

Fig. 6.5 Unfiltered pressure data at CAD interval [0,35]. The measurements
contain very peaky resonances that will cause problems in the differentiation of
pressure. The graph to the right is the same pressure data when a Savitzky
Golay filter is applied.

Page 53 (80)
The first peak in pressure, as can be seen in Fig. 6.5, at 0 CAD is caused by the
isentropic compression performed during the compression stroke. It is in fact this
isentropic pressure trace data that is used when is estimated in Eq. (3.41).
-3
x 10
2.5

1.5
Volume / m3

0.5

0
-400 -300 -200 -100 0 100 200 300 400
CAD

Fig. 6.6 Cylinder volume as function of CAD. Notice that the low
volume at CAD = 0 has a corresponding pressure peak due to
isentropic compression.

Page 54 (80)
6.2.1 Estimation results of and the impact on heat release

The estimation of gamma, as proposed in Eq. (3.41), resulted in a mean value of

1.32. This is somewhat higher than the constant value used in most of the
simulations ( 1.3). The corresponding heat release is presented in Fig. 6.9.
Estimation of gamma
1.4
mean value: 1.32
gammavalue (-)

1.35

1.3

1.25
-55 -50 -45 -40 -35 -30 -25 -20 -15
CAD
7 Pressure trace
x 10
2
Pressure/ Pa

1.5

0.5

0
-400 -300 -200 -100 0 100 200 300 400
CAD

Fig. 6.8 Pressure trace and corresponding gamma estimate. The mean value
was determined to be 1.32.

Normalized accumulated Heat Release Normalized accumulated Heat release

1.2
Measured HR Calculated HR
Calculated HR 1 Measured HR
1

0.8
0.8
HR (normalized)

0.6
HR (normalized)

0.6

0.4
0.4

0.2
0.2

0
0

-0.2
-0.2
-100 -50 0 50 100 150
-40 -20 0 20 40 60 80 100 120
CAD
CAD

Normalized accumulated Heat release

Calculated HR
1 Measured HR

0.8
HR (normalized)

0.6

0.4

0.2

-40 -20 0 20 40 60 80 100 120

CAD

Fig. 6.9 Normalized accumulated heat release at the estimated gamma,

1.32 (top left), =1.3 (top right) and =1.35 (bottom left). CA50 = 9 in
top cases and CA50 = 10 when =1.35.

Page 55 (80)
Greater accuracy could be achieved by utilizing the fact that can be estimated
before and after combustion fairly accurately, assuming adiabatic compression.
Tunestl proposes a linear approximation inbetween to account for changing
during the combustion process itself.23 However, there is an upside by the simple
algorithm as given by Eq. (3.41). It does in fact account for any EGR that may
contribute to change in the polytropic exponent, remembering that .

In this thesis, , the specific heat value at constant pressure, is determined

using NASA Glenn interpolation tables. Remembering that the polytropic
exponent is in fact a ratio between the specific heat value at constant pressure
and at constant volume it is possible to have a temperature dependent estimate
of . The specific heat at constant volume can be determined from ,
and thus is also known. Results from heat release calculations when this
approach is used, is shown in Fig. 6.10 6.12.

Estimated heat release

350
Heat release
Heat release at gamma = 1.3
300

250
Heat Release J/CAD

200

150

100

50

-100 -50 0 50 100

Crank angle (CAD)

Fig. 6.10 Results of heat release calculations when the polytropic exponent is
made adaptive by using NASA Glenn tabulations of specific heats. Calculated
with EGR% 10%. Interpolated specific heats are used i.e.
1.3. . 6.25 .

Page 56 (80)
 Estimated heat release
350 Heat release
Heat release at gamma = 1.3
300

250
Heat Release J/CAD

200

150

100

50

-150 -100 -50 0 50 100

Crank angle (CAD)

Fig. 6.11 Results as described in Fig. 6.10 but for another set of data.
Calculated with EGR% 15%. 1.3. . 6.25 .

400
Estimated heat release
Heat release
350 Heat release at gamma = 1.3

300

250
Heat Release J/CAD

200

150

100

50

-50
-50 0 50 100
Crank angle (CAD)

Fig. 6.12 Results as described in Fig. 6.10 but for another set of data.
Calculated with EGR% 0% . 1.3. . 6.25 .

As can be seen in Figs. 6.106.12, the heat release estimated by using an

updating will have a somewhat smaller peak in all three cases. The case when
a constant 1.3 gives almost the exact heat release results as the adaptive .

Page 57 (80)
6.3 Temperature estimations

The model assumes two different zones during the combustion phase. One zone
will contain all burned products due to combustion of fuel, EGR included. The
other zone contains unburned fuel, air and EGR. The concentration of EGR is
assumed to be uniform in the whole cylinder. That means that the EGR
concentrations in both zones are in fact equal at all times. The temperature will
be essential in the Zeldovich mechanism and thus also essential to NOx
formation.

The unburned zone was assumed to change temperature purely due to isentropic
compression. This means that the heat transfers between the zones where
neglected. Although somewhat physically unjustified, this assumption is common
in zero dimensional modeling.

As can be seen in Fig. 6.13, the temperature in the burned zone peaks at around
2500 K. The temperature will in fact, as stated previously, vary with the EGR in
the gas content. The peak temperature in the unburned zone is much lower, at
around 1000 Kelvin. This emphasizes the unrealistic aspects of the assumption
that the temperature changes in the unburned zone is solely due to isentropic
compression. This would imply a temperature gradient, between the zones, that
would be too large in magnitude for it to hold in reality.

Additionally a temperature decrease was added due to radiative heat losses of

the combustion products. This was stated as:

The variable is used to calibrate the temperature loss.

Aside from the radiative temperature loss, there is also a decrease in

temperature due to cylinder expansion which is taken into account with a factor
for the isentropic compression. However, the temperature decrease due to
dissociation of carbon dioxide into carbon monoxide and water to hydrogen was
neglected. The relationship between dissociation and temperature decrease is
complex and hard to calibrate, especially with EGR present in the gas. A previous
study by Andersson1 indicates that the temperature decrease due to dissociation
is at average circa 5 % of the total flame temperature at the pressure range of
interest and at a close to 1.

Page 58 (80)
 burned zone
cylinder
2500 unburned zone

2000
Temperature (K)

1500

1000

500
0 20 40 60 80 100 120
CAD

Fig. 6.13 Temperature in burned and unburned zone with the mean gas
temperature included.

A closer observation of Fig. 6.13 implies that the generated temperature results
seems feasible. It would be expected that the burned zone is the hottest one,
followed by the mean gas temperature and finally the unburned zone
temperature which is the colder zone in the cylinder.

Page 59 (80)
6.4 Equilibrium concentrations

When the temperatures are known it is possible to calculate the equilibrium

concentration in the burned zone using the equations described in the theory
section. The variables that are dependent on temperature are in fact only the
equilibrium constants that are needed in two of the equations. The resulting
equilibrium constants at a specific temperature are shown in the figures below.

4
Kp value

0
1000 1500 2000 2500
Temperature / Kelvin

Fig. 6.14 Values of at different temperature (see R3.7). This constant shows
a near linear dependency on temperature in the given temperature range.
5
x 10
2.5

1.5

0.5

0
1400 1600 1800 2000 2200 2400 2600

Fig. 6.15 Values of at different temperatures (see R3.8). This equilibrium

constant is far more dependent on temperature than the values of .

Page 60 (80)
The equilibrium constants were interpolated with polynomials for easier and
faster calculations. was approximated with only one fourth degree
polynomial, but had to be split up in three separate sections due to its heavy
dependence on temperature. The resulting polynomials are given as:

Eq. (6.4)

2.6223 10 2.4793 10 7.8685 10 0.0064678 1.5537

Eq. (6.5)

2.9827 10 0.19986 447.42 3.3477 10 , 1500, 2250

2.9375 10 0.022093 55.554 46730, 2250, 2500

8.8 10 1.13 10 0.0054 11.4 9156 , 2500, 3500

The large magnitudes involved in the Jacobi matrix rendered the

Newton Raphson method(see section 3.6) unstable if a relaxation matrix was not
introduced. However, the method still needed a good initial guess. Fortunately,
the initial guess is simple to set if the fuel is known. At combustion without the
gas water shift reaction, it is possible to calculate the equilibrium of , ,
and analytically with some help from the stoichiometric relation given in R3.5.
A good initial guess for iso octane is:

Eq. (6.6) 8 0.1 13 1 0.1

There are a total of 6 concentrations that need to be determined. However, the

initial guess contains only 5 concentrations. This is due to the fact that is
trivial to determine. This can be seen in Eq. (3.95) by the fact that the
concentration of nitrogen is independent of all the other equilibrium
concentrations. The nitrogen is in the first step non reactive, because the model
assumes two subzones in the burned zone. The nitrogen will only participate in
chemical reactions in the post combustion zone where NOx is formed. This
simplifies the system of equations to a total of five equations.

The convergence of the Newton Raphson method will not depend on deviations
of lambda or normalized pressure. The largest contribution to stability issues are
large magnitudes of the equilibrium constant . In some cases the
Newton Raphson method can solve the equations in only 20 iterations with
fairly good accuracy. However, there is much more room for improvement when it
comes to the relaxation matrix as described in section 3.6. One suitable solution
would be to make it adaptive.

Page 61 (80)
The most suitable approach to the issue of solving these linear equations is
simply to solve them offline for variations in normalized pressure and
equilibrium constants. The algorithm used above and the MATLAB function
fsolve are too slow to motivate an implementation in C for real time use in an
embedded system. The solutions retrieved offline can then be tabulated and used
online in a table format.

The Newton Raphson method derived in this thesis was also validated for some
specific values of parameters. The Newton method results matched the fsolve
results fairly well in most cases, with a far lesser computational time. However,
even the Newton method, as it is today, is probably not time efficient enough for
real time executions and for large magnitudes of the equilibrium coefficients, the
method was unstable at times. The average time to solve the equations with
~1500 iterations took 35 milliseconds in MATLAB. This is of course highly
unsuitable results if a model ought to be implemented and run in real time
applications. In cases where the convergence was sufficient after only 20
iterations, it took 6 ms in MATLAB. However, the method is considerably faster
than fsolve, which utilize a GaussNewton method on the form min
. This criterion is minimized over , which will give the
solution .

Page 62 (80)
6.5 Model NOx output

The output of the model is NOx concentration, given in parts per million. As
stated in Egnell4, NOx is not formed during the premixed combustion phase.
Therefore this part was ignored during the calculations. Calculations of NOx
concentrations are initiated a number of samples before the temperature peak in
the burned zone. This approach gave the most accurate result at given calibration
values.

Measurement Measured NOx (ppm) Calculated NOx (ppm)

Point 1 (EGR% = 7.4%) 86 86
Point 2(EGR% = 0.0%) 2261 2879
Point 3(EGR% = 5.9%) 2201 1771
Point 4(EGR% = 9.4%) 1419 1747
Point 5(EGR% = 11.3%) 1342 1538
Point 6(EGR% = 13.3%) 1126 1394
Table 6.2 Measured NOx vs. calculated NOx from model at 1000 rates per
minute and 1 at a set calibration.

As can be seen by table 6.2 the model estimate of NOx is much larger than the
measured NOx in Point 2. This could be motivated by the fact that the calculated
temperature in the burned zone is somewhat higher than the real temperature.
This will have a large impact on the calculated NOx concentrations that cannot
be neglected. The assumption that the temperature drop due to dissociation can
be neglected will hold in cases where EGR is present. However, when EGR is not
present in the gas content, the assumption will introduce an error. The same
holds for the temperature drop due to radiative heat losses in the post
combustion zone.

3000 25
Measured NOx
Estimated NOx
2500
20

2000
C onc entration (ppm )

15
Error (% )

1500

10
1000

5
500

0 0
1 2 3 4 5 6 1 2 3 4 5 6
Measurement point
Fig. 6.16 Graphical comparison of estimated and measured NOx from Table
6.2. EGR% = 7.4%, 0.0%, 5.9%, 9.4%, 11.3% and 13.3%, respectively.
Page 63 (80)
 1200
Measured NOx
Estimated NOx
1000

Concentration (ppm) 800

600

400

200

0
1 2 3 4

Fig. 6.17 Model comparison to measurement values of NOx for a data set at rpm
= 1200, 1.2 and EGR 14.9%, 21.3%, 20.7% and 22.6%, respectively.

Fig. 6.18 An example of a CADresolved simulation results of NOx formation.

As can be seen by Fig. 6.18, the NOx concentration will, at some point in time,
have an overshoot and then reach a final stationary value. This simulation result
is consistent with other results that utilize the Zeldovich mechanism to describe
NOx formation.8

Page 64 (80)
Additionally, a revised version of the Zeldovich mechanism was used. The NOx
concentration as given in Eq. (3.26) assumes a constant volume. However, this is
hardly the case during combustion and therefore the equation was slightly
modified to account for volume changes12:

Eq. (6.7)

The calibration parameter of the NOx model was the temperature drop due to
radiative heat losses generated. This specific approach of model calibration is
supported by previous research conducted by Ericson.24

Page 65 (80)
6.6 The division algorithm

The division algorithm derived in this thesis was evaluated in MATLAB with
fixed point representation implemented with the MATLAB builtin function bit
shift. The method seems promising as long as it is possible to retrieve the most
significant bit in a reasonable amount of time. The error is displayed in Fig. 6.19
where the function 1/x was plotted with utilization of the thesis division
algorithm and compared to the real division. One thing should be underlined: the
error in Fig. 6.19 also takes quantization into account. This means that the graph
shows what the error would be if the algorithm was run real time in C.

14
2nd degree term
3rd degree term
12
4th degree term

10

8
Error (%)

0
0 200 400 600 800 1000 1200 1400 1600 1800 2000

Fig. 6.19 The error relative to the true value of division .The maximum error is
around 12% and most often below 8% for the 2nd degree case.. Errors caused by
quantization are also included in the plot. The xaxis shows the divisor. As can
be seen, the error is cut in half for each step higher in degree.

As can be seen, the error is reasonably low considering the degree of the
approximation (see Eq. 5.16). Accuracy, unfortunately, must be sacrificed in order
to get a good run time online. Another positive aspect of the algorithm is the fact
that its stable at all times. This is in contrast to the NewtonRaphson approach
to the problem, where the initial guess has to be very close to the true value.

Initial guess 0.1 0.5 0.75 0.9

Result 0.2857 675
0.2828 0.67 million
Table 6.3 Results with NewtonRaphson applied to the division 1/3.5 after 4
iterations. The true value is in this case 0.2857. The Newton can use a relaxation
value << 1, if more stability is desired. Higher yield greater stability issues. A
conclusion from this is that the Newton method is too poor with respect to
convergence and speed for a reasonable online utilization.

Page 66 (80)
For the case when the divisor is close to 0, the alternative algorithm renders an
error plot as shown in Fig. 6.20. The approximation used was:

Eq. (6.8) ~ 1 1

Error plot for approximation of 1/(1-x)

7

4
Error (%)

0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
x (dimensionless)

Fig. 6.20 The error plot when the alternative algorithm is used in a range x
[0.0,0.5].

Higher accuracy can of course be reached by using terms of higher degree. The
next, more accurate approximation would be (see Eq. 5.18):

Eq. (6.9) ~ 1 1 1

However, this is rather unnecessary in, terms of accuracy, if following

mathematical trick is done:

Eq. (6.10) ~ 1 1 1 1 1

Neglecting the term will give an approximation of the division as follows:

Eq. (6.11) ~ 1 1 1 1 1

The resulting error plot can be seen in Fig. 6.21. As can be seen in Fig. 6.21 the
maximum error is far smaller than the maximum error shown in Fig. 6.20. This
is of course not for free as the approximation in Eq. (6.11) needs 2 additional
operations (1 multiplication and 1 addition). However, the same accuracy can be
achieved by using the approximation in Eq. (6.8) and simply bit shifting into the
range where 0.0,0.25 . This seems like a promising approach to approximate
divisional operations.

Page 67 (80)
 Error plot for approximation of 1/(1-x)
1.6

1.4

1.2

Error (%) 0.8

0.6

0.4

0.2

0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
x (dimensionless)

Fig. 6.21 The corresponding error plot when 1 1 1 is used as

an approximation instead.

Page 68 (80)
6.7 Interpolation results

This section is supposed to illustrate how interpolation can simplify rather

gruesome and nonlinear calculations. To avoid unnecessary nonlinear
calculations some interpolations were made to enable fast calculations. For
example, the isentropic compression gives a relation on the form:

Eq. (6.12)

This is however neither desirable nor efficient to calculate on an FPGA. The

approximation of the relation was done by linear interpolation in specific
segments. As can be seen in Fig. 6.22, three segments give a decent fit even when
a simple linear interpolation is used.

3.5

2.5
Result

1.5

1
0 20 40 60 80 100 120 140 160 180 200
Ratio p/p0

Fig. 6.22 Isentropic compression calculations with a fixed gamma, 1.3. The
three segments match the analytical calculation quite decently, in the interval of
interest.

In the case when multivariable functions are interpolated, the approach was to
produce a least squares fit in several variables, as described in Section 5.2. For
example, the results for the specific heat is given in Fig. 6.25. As can be seen the
least square fit gives results that can be used in practice. Unfortunately, there
are times when the least squares fit gives inappropriate results. However, these
cases are extreme and do not occur during engine combustion. As an example, the
approximation of the Woschni heattransfer coefficient would give a negative
result when the pressure is at 1 bar and a temperature at [1500, 2000] K. This is
not a realistic condition during combustion in a diesel engine (see Fig. 6.26).

Page 69 (80)
 1000
analytical
interpolated
900

800
Temperature (K)

700

600

500

400
-100 -50 0 50 100 150
CAD

Fig. 6.23 Comparison between using exact isentropic equation and discretized.
The spikes in the interpolated version are due to discontinuities in the
interpolation, as can be seen in Fig. 6.17, but this can be smoothed out.

2600 interpolated
analytical
2400

2200

2000
Temperature (K)

1800

1600

1400

1200

1000

800
-20 0 20 40 60 80 100 120
CAD

Fig. 6.24 Interpolated temperature (red) of burned zone versus the analytical
calculation of the same. As can be seen the interpolation gives a result that
follows the analytical one very accurately.

Page 70 (80)
 38 0.5
S pec ific heat (J / K m ole)

0.4
36

0.3

Error (%)
34
0.2
32
0.1

30 0
1400 1400
1200 0.8
1200 0.8
1000 0.6
1000 0.6
0.4
800 0.4
0.2 800 0.2
Temperature (K) 600 0
EGR (%) Temperature (K) 600 0
EGR (%)

Fig. 6.25 Specific heat interpolation results together with an error plot. The
interpolation was done at a fixed global lambda, global = 1.3. A value which is
fairly usual in practice.

4
x 10

-1
4000
3000 3
2000 2
7
1000 1 x 10
Temperature (Kelvin) 0 0
Pressure (Pascal)

Fig. 6.26 Approximation of the factor . .

which is used in Woschni (see
Eq. (3.30)). The interpolation, approximated with least squares calculations, is
the one that is below the other. Alternatively, . and .
could be interpolated
separately due to the trivial multiplicative nature of the factor.

Page 71 (80)
7 Calculations of time efficiency

In this section a study of the time efficiency of the different algorithms and
calculations is conducted. Because of the fact that the final model will be
implemented on an FPGA, only calculation operations such as multiplication,
addition, divisions and subtractions are accounted in the analysis, and not actual
time per se. Another aspect to keep in mind is that the FPGA runs in a truly
parallel manner. For example, this implies that the calculation can be
done in only two clock cycles. How so? The parallel behavior of the FPGA gives
the possibility to calculate and at the same time before being added
together.

7.1 Division algorithm

The main purpose of the division algorithm was to avoid conventional division
on an FPGA, which is proven to be rather time consuming. A conventional
division consumes at least 20 clock cycles and also occupies a substantial amount
of memory25. The time efficiency (number of operations) for the different proposed
ways to approximate division is illustrated in Table 7.1. Calculations that can be
done in parallel are only counted once in the table, i.e considerations are taken to
the FPGAimplementation and its parallel execution.

Algorithm Additions Subtractions Multiplications Bitshifts

Eq. (5.16) First 2 1 2 2 7
Eq. (6.8) Improved 2 0 2 2 6
Eq. (6.11) Best fit 3 0 3 2 8
Table 7.1 Number of operations for each division approximation calculation. As
can be seen, the approximation described by Eq. (6.8) only needs a total of 6
mathematical operations.

The most promising one seems to be Eq. (6.8) in terms of speed and accuracy.
After all it is possible to bit shift into a range 0.0,0.25 , as has been stated
previously.

7.2 Specific heat interpolation

An approximative approach to calculate the specific heat at constant pressure is

presented in this thesis. The method was to utilize a multivariable least square
approach to fit a function , to a set of output data (i.e. specific heats).
The other approach would of course be to use NASAinterpolations and
calculations as suggested by Eq. (6.2) and Eq. (6.3). The results in section 6.7
indicate that the least square interpolation gives highly accurate fits. However,
another interesting aspect is of course the time saving of doing so. This is in fact
the main idea of the interpolation in the first place: saving calculation time.

Page 72 (80)
The interpolation function of the specific heats of interest, in the NASA Glenn
table, is on a form:

Eq. (7.1)

This interpolation function holds for all the molecular compounds used in the
calculations, which are , , and . There is also the calculations in Eq.
(6.2) and Eq. (6.3) that have to be counted. The time comparison between the
interpolation and the analytical calculation of the specific heat is given in
Table 7.2. Calculations that can be done a priori is of course neglected.

Method Subtractions Additions Multiplications Divisions

Analytical 2 30 55 9 96
Interpolation 0 5 8 0 13
Table 7.2 The number of operations needed in the two different approaches to
calculate the specific heat at constant pressure needed in temperature and
calculations.

Yet again, considerations to the fact that FPGA runs in a true parallel manner
should be taken into account. Thus the results in Table 7.2 do not correctly
mirror the possibility to, for example, calculate Eq. (7.2) in parallel for each
molecule: , , and . In an FPGA implementation, the number of
operations would instead a lot less, which can be seen in Table 7.3. However, it
seems reasonable to prefer the interpolation function for specific heat
calculations instead.

Method Subtractions Additions Multiplications Divisions

Analytical 2 12 22 3 39
Interpolation 0 5 2 0 7
Table 7.3 Number of operations needed in an FPGA case when true parallelism
can be utilized to save calculation time.

The interpolation still saves some considerable computation time, even in an

FPGA implementation. Additionally, less memory is needed to store constants
and coefficients.

Page 73 (80)
8 Discussion

The heat release calculations are coherent with the heat release calculation (in
this thesis the HR is incorrectly called measured HR) done at a Scania
predevelopment group, considering the fact that the automatic heat release
calculation at Scania is not 100% accurate with respect to the actual heat
release, as well. The SavitzkyGolay filter23 successfully removes the high
frequency calculation noise. In a real time implementation it would be wise to
add some kind of electrical filter to filter the pressure measurements. Another
alternative would be to smooth out the derivative part of the pressure in the
heat release equation. Also, introducing an adaptive did not contribute much
change to the heat release. At least not on the particular set of data used in the
validation process. However, there are reasons to believe that combustion at
small global with high EGR rates will make the adaptive more important in
heat release calculations. It is also quite easy to calculate the using the least
square fit polynomial of the specific heat at constant pressure, knowing that
and / .

The approach to calculate the temperature in this thesis seems promising and
easy enough for a future implementation on an FPGA. It also accounts the
impact of EGR on temperature in a trivial manner, through changing mean
specific heat. Of course, the approach assumes a homogeneous content in the
engine cylinder. The leastsquares interpolation of the mean specific heat heavily
reduces calculation time at a very low accuracy cost of below 0.5%. Such a small
percentage is in fact negligible compared to model errors, and quantization
effects in a fixed point implemented code. The temperature calculation results
show that the burned zone will be the warmest zone, followed by the cylinder
mean temperature and lastly, the unburned zone. In fact, this result is expected.
However, large temperature gradients should arise between the two zones which
would render heat transfers between the burned and unburned zone. One can
easily question the plausibility of assuming zero heat transfers between the two
zones in question. The assumption is however, as stated previously, quite
standard in combustion modeling.

The NOx results indicate that the model of NOx formation is at least acceptable,
even if it lacks a realistic physical interpretation. The errors are quite large at
over 20% and it would indicate that more work can be done on parameter tuning
and model modification.

When it comes to calculating equilibrium concentrations (see theory section 3.6),

no time efficient method was found for solving the equations.
The Newton Raphson method developed in this thesis needs a relaxation matrix

Page 74 (80)
which will slow down the algorithm as a whole, because more iterations will be
needed.

The suggested algorithms for approximation of the division operation seem to be

promising. Especially the Eq. (6.8) indicates two essential properties: it can be
highly accurate for small values of and it does not require a lot of operations. In
a NOx model, a large number of the quantities relevant are small magnitude. For
example, the total amount of substance in the engine cylinder is in the region 0.2
mole, also using SI units will give small values of volume, equilibrium
concentrations etc. This implies that the requirement 1 in Eq. (6.8) will be
fulfilled in many cases.

Page 75 (80)
9 Conclusions

In this thesis a NOx model is derived for the purposes of calculation of NOx
output from heavy duty diesel engines. The model is zero dimensional and it
also utilizes two zones during the combustion phase.

The reactions occurring in the engine cylinder are hard to describe in full
extension, especially if the purpose is to simplify things for future real time use.
However, considering the results, the three reactions in the Zeldovich mechanism
together with the assumption that the system reaches chemical equilibrium
seems to be a good approach to a simple NOx model.

The heat release is calculated by using the diagnostic ROHR equation. The
results are coherent with data received from Scania CV AB. Although, the gross
heat release is harder to calculate considering the fact that calibration of the
Woschni heat transfer equation is needed. One way to do this is of course to
compare the accumulated heat release with the energy content of the fuel
injected.

The division algorithm derived in this thesis can reduce the number of operations
needed for a division operation. However, an implementation and validation
procedure on an FPGA is needed for a final confirmation of this conclusion.

The leastsquares fitting of specific heats give a very accurate fit with errors
below 0.5%. This is negligible compared to other errors, such as tuning of the
parameters used in the Woschni heattransfer equation. The leastsquares
fitting of specific heats reduce the amount of operations and also less memory
usage is needed.

Page 76 (80)
10 Future work

The natural continuation in the process of creating a fast and sufficiently

accurate NOx model would be to:

(1) Modify the model to include radiative heat losses in the heat release
equation.
(2) Implement and utilize a predictive heat release model instead of the
current diagnostic one. This would of course need an additional model of
the combustion process and fuel injection. A suitable approach is described
extensively in Anderssons Licentiate dissertation.
(3) Produce memory maps of NOx concentrations as functions of temperature,
pressure and , instead of using the equations derived from the Zeldovich
mechanism.
(4) Calibrate the model using larger data sets to produce a more accurate
model of NOx formation.
(5) Implement the model on an FPGA.

Another interesting investigation would be to study the possibility to derive a

black box model of NOx formation. The input data that could be used in such a
model are: pressure trace, EGR rate and cylinder temperature. Alternatively, one
could implement a gray box model. After all, the heat release and temperature
can easily be calculated using pressure data and the diagnostic heat release
equation. Knowing these two, it would be possible to use system identification to
produce a NOx model with burned zone temperature, reaction rates and heat
release as input data.

Page 77 (80)
Nomenclature

Mass of element XX in kilograms ( )

Molar mass of element XX in kilograms ( / )

Pressure in Pascal ( )

Universal gas constant ( / )

Volume in cubic meters ( )

Displacement volume ( )

Compression volume ( )

Temperature in Kelvin (K)

Crank angle in CAD ( )

Stoichiometric air to fuel ratio ( )

Polytropic exponent ( )

Specific heat at constant pressure ( / )

Specific heat at constant volume ( / )

Amount of substance ( )

Cylinder bore ( )

Conrod length ( )

% Percentage of EGR in gas content ( )

Page 78 (80)
11 References

1: M. Andersson, Fast NOx Prediction in Diesel Engines, Licentiate Thesis,

Division of Combustion Engines, Department of Energy Sciences, Lund

2: Euro VI Regulations for emissions from Lorries and Buses, EEB position
paper, European Environmental Bureau, Contact person: Dragomira Raeva,
dragomira.raeva@eeb.org, 2008.

3: C.G Zander, SingleCylinder Diesel Engine Experiments, Modeling and In

Cycle Control with Heatrealease Emphasis, Division of Combustion Engines,
Department of Energy Sciences, Lund, 2010.

4: J.E Dec, A conceptual model of DI diesel combustion based on lasersheet

imaging, SAE technical paper, 1997.

5: J. B. Heywood, Internal Combustion Engine Fundamentals, 1988, ISBN 007

1004998.

6: http://www.tc.gc.ca/, accessed 21 June, 2011.

7: http://www.cormetech.com/catalystoverview.htm, accessed 23 march 2011.

8: www.afcd.energy.gov, accessed 20 June, 2011.

9: R. Egnell, On Zerodimensional Modelling of Combustion and NOx formation

in Diesel Engines, 2001, ISSN: 0282-1990.

10: Long Liang, Efficient Simulation of Diesel Engine Combustion Using

Realistic Chemical Kinetics in CFD, SAE International, 2010-01-0178, 2010.

11: http://www.ngk.co.jp/english/news/2008/image/0611_06.jpg, accessed March

24, 2011.

12: R. Wain., An overview of FPGAs and FPGA programming; Initial experience

at Daresbury, CCLRC Daresbury Laboratory, November 2006.

13: C. Wilhelmsson, Embedded Systems and FPGAs for implementation of

control oriented Models, November 2009, LUTMDN/TMHP09/1068--SE.

14: http://en.wikipedia.org/wiki/Field-programmable_gate_array, accessed 24th

March, 2011.

15: http://www.intertechnology.com/Kistler/pdfs/Pressure_Model_6013C.pdf,
accessed 6th September, 2011.

16:http://energihandbok.se/x/a/i/10198/NOx-bildningen-vid-forbranning.html,
accessed, 30 march, 2011.

Page 79 (80)
17: G. Woschni, A universally applicable equation for instantaneous heat
transfer coefficient in the internal combustion engine, SAE technical Papers
670931, 1967.

18: B. Johansson, Frbrnningsmotorer, Division of Combustion Engines,

Department of Energy Sciences, Lund, 2006.

19: B.J. McBridge et al., Nasa Glenn Coefficients for calculating Thermodynamic
properties of Individual Species, TP-2002-211 556, 2002.

20: H. S. Warren, Hackers delight, 2002, ISBN: 978-0-201-91465-8.

21: D. Eberly, Leastsquares Fitting of Data, 2008,

http://www.geometrictools.com/, accessed, June, 2011.

22: A. Savitzky & M.J Golay, Smoothing and differentiation of data by simplified
least squares procedures, Analytical Chemistry, 36, p. 16271639.

23: P. Tunestl et al., Validation of a Self Tuning Gross Heat Release Algorithm,
SAE paper, 2008-01-1672, 2008.

24: C. Ericson, Model Based Optimization of a Complete Diesel Engine/SCR

System, 2009. ISBN: 9789162877248.

25: C.G Zander, Internal document (document name/label code is missing), 2011.

Page 80 (80)