Internal Combustion Engine
Simulation
[Link] Babu
Head, Centre for Energy Studies
I.I.T Delhi
Internal Combustion Engine
Internal Combustion (IC) engines have
completely revolutionized transportation,
power generation and have perhaps altered
the way the society operates forever.
I.C engines are classified as
Spark Ignition Engine
Compression Ignition Engine
Simplest Model for an [Link]
The system is closed.
Air is the working fluid and is modeled as an ideal gas
throughout the cycle.
Compression and expansion processes are isentropic.
A reversible heat transfer process characterizes the combustion
of fuel and air.
Heat rejection takes place reversibly and at constant volume.
Otto cycle is used to model a basic Spark Ignition engine
Diesel Cycle is the basic model for the Compression Ignition
engine.
Spark Ignition Engines (Otto
Cycle)
Spark-ignition engines are the most common type used in cars.
Larger engines operate using a four-stroke cycle, while smaller
engines operate on a two-stroke cycle.
In a simple four-stroke cycle, a combustible mixture of air and
fuel is drawn into a cylinder during the intake stroke, and the
temperature and pressure of the mixture is raised during the
compression stroke.
At near the maximum compression, a spark initiates combustion
of the mixture, raising its temperature and forcing expansion.
The expanding gases do work on the piston during the power
stroke and then the burnt gases are purged during the exhaust
stroke.
Typically 3000 or more such cycles are repeated in a minute.
Otto Cycle
The Otto cycle is an air-standard model of the actual cycle.
In addition to the air-standard assumptions listed above, the
combustion process is modelled as a reversible constant volume
heat addition process.
The four steps of the air-standard Otto cycle are outlined below:
(1-2) Isentropic compression (Compression Stroke)
(2-3) Constant-volume, reversible heat addition (Ignition)
(3-4) Isentropic expansion (Power Stroke)
(4-1) Reversible, constant-volume heat rejection (Exhaust)
Otto Cycle
Otto Cycle
Ideal Otto Cycle
Compression Stroke. Adiabatic compression of gas
fuel mixture in the cylinder.
Ignition of gas fuel mixture. Take place rapidly at top
of the compression stroke while the volume is
essentially constant.
Expansion Stroke. Adiabatic, isentropic expansion of
gases in the cylinder after fuel mixture is ignited. This
is the part of the cycle that does positive work.
Exhaust of the spent gases and the intake of a new
fuel mixture into the cylinder.
The volume is the same at beginning and ending of
the exhaust and intake stroke.
Diesel Cycle
Diesel Cycle
Ideal Diesel Cycle
Compression Stroke. Adiabatic compression of gas
fuel mixture in the cylinder.
Ignition of gas fuel mixture. Fuel is ignited by high
temperature due a large compression. Burning takes
places while the pressure is essentially constant.
Expansion Stroke. Adiabatic, isentropic expansion of
gases. This is the part of the cycle that does positive
work.
Exhaust of the spent gases and the intake of a new
fuel mixture into the cylinder.
The volume is the same at beginning and ending of
the exhaust and intake stroke.
Two Zone Formulation
Two Zone Formulation
Two Zone Formulation Basic
Relations
Derivatives for the two temperatures,
pressure and the two volumes may be
derived
The result of the derivation is the
following five ordinary differential
equations:
Three-Zone Formulation
The next step is to divide the burned zone into the boundary layer and
adiabatic core zones.
Figure 2 is a schematic of the boundary layer and adiabatic zones.
Initially, the boundary layer has zero mass, and the adiabatic core is
equal to the burned mass.
As shown, the boundary layer increases in mass by entraining mass
from the adiabatic core.
The rate of mass entrainment of the boundary layer is dictated by the
thermodynamics and temperature relations (described below).
All the burned gas heat transfer is assigned to the boundary layer.
For these additional two zones, the formulation of the energy
equation is similar to the above for the overall system.
The first law, then, for the adiabatic and boundary layer zones,
respectively, are
Three Zone Formulation
Boundary layer
The growth of the boundary layer is dictated
by the mass and energy conservation
relations for the burned zone.
A choice at this point is how the average
boundary layer temperature is defined.
One classical definition for the average
boundary layer temperature is based on a
logarithmic relation of the wall temperature
and the adiabatic zone temperature
Boundary Layer Equations
Three Zone Formulation
To summarize, the new unknowns for the boundary layerand
adiabatic core zones are Ta, Tbl, Va, Vbl, ma and mbl.
A number of different sets of differential equations may be
formulated for these variables.
Because these differential equations are highly nonlinear and
exhibit complex interrelations, not all combinations of the above
equations were successful.
After extensive trials, the following six (6) ordinary differential
equations for these variables were found to be the most robust
in terms of obtaining stable solution
Boundary Layer Temperature
Three Zone Equations
Items for the Gas Temperature
Derivative Equation
To solve the above governing
differential equations (eq. 2,10 and 18),
several items are needed.
These items include the thermodynamic
properties, the mass fraction burned,
the convective heat transfer, and the
mass flow rates.
These items are described next
Thermodynamic
Properties:
The thermodynamic properties needed for
solving the first law include the instantaneous
specific internal energy, specific enthalpy,
specific gas constant, molecular mass, and six
(6) property derivatives
All the properties are needed as a function of
time (crank angle) for the cylinder contents
for the three zones, and for any matter
entering the cylinder
Thermodynamic Properties
Depending on the specific processes during a cycle,
the thermodynamic properties may be for different
mixtures of air, fuel vapor, and combustion products.
The concentrations of the combustion products may
be “frozen” for the lower temperatures, or these
concentrations may be based on an instantaneous
determination of chemical equilibrium (see e.g.,
Olikara and Borman, 1975) for higher temperatures.
Complete descriptions of the algorithms used for
determining the compositions are presented in a
number of references
Thermodynamic Properties
Once the composition is known, the individual thermodynamic
properties may be determined.
The properties of each species in the mixture are first
determined for the given temperature and pressure by the use
of polynomial curve fits (Heywood, 1988) to the thermodynamic
data.
The overall mixture properties are then determined by suitable
expressions.
Mass Fraction Burned: The mass fraction of fuel (or fuel
mixture) burned (xb) may be determined from a first-law
analysis of the measured cylinder pressure
The result of this analysis is almost always an “S” shape type
curve. The mass fraction burned is defined as
Mass Burnt Fraction
Wiebe Constants
and where “a” and “m” are parameters that
are selected to provide a match with
experimental information.
The values that provide a good match with
the data depend on complex functions of the
turbulence, chemistry, fuel, equivalence ratio,
chamber geometry, and a number of other
features.
The values used here are based on the work
of Heywood et al. a = 5.0; m = 2.2 (22)
Heat Transfer
Heat Transfer
During the combustion process, the total heat transfer is
allocated to the various zones (but the total is still given by eq.
23).
This allocation may be completed in a number of ways. For this
work, the heat transfer is assumed to be proportional to the
surface area of each zone.
The appropriate surface area is assumed to be proportional to
the volume raised to the 2/3 power.
In addition, the heat transfer is proportional to the temperature
difference.
For example, since the unburned zone temperature is lower, the
proportion of the heat transfer allocated to the unburned zone
should be lower as well.
Heat Transfer
Heat Transfer
Since the burned zone is divided into an
adiabatic core and a boundary layer, the
burned zone heat transfer is assigned in total
to the boundary layer.
A number of correlations for the
instantaneous convective heat transfer
coefficient ( hc ) have been proposed and
used throughout the years.
The current work has used the Woschni’s
correlation
Mass Flow Rate
The mass flow rates are based on
assumptions and approximations which
include that the flow is quasi-steady, one-
dimensional, adiabatic and reversible.
These assumptions and approximations are
corrected by use of an empirical discharge
coefficient.
For purposes of this study, the discharge
coefficient has been assumed constant with a
value of 0.7 (average value)
Flow Area
Finally, the open flow area is needed which
depends on the valve diameter and
instantaneous valve lift.
Since the valve lift profiles are not necessarily
known, this work has used a standard
assumption for the valve lift based solely on
the maximum lift and the valve open
duration, and using a sinusoidal shape as
recommended by Sherman and Blumberg
(1977).
Mass Conservation
In addition to the energy equation described
above, mass conservation must be satisfied
throughout the computation.
The instantaneous change in the cylinder
mass is equal to the net result of the mass
flow rates:
Numerical Solution
The instantaneous cylinder conditions (temperatures, pressure,
volumes, masses, and thermodynamic properties) as a function of
crank angle are obtained by the simultaneous numerical integration of
the various differential equations.
Since a number of the differential equations are themselves dependent
on the other variables, an iterative approach was chosen.
That is, the set of differential equations were evaluated, and then
these values were used to determine the “final” values of the
derivatives.
For the two zone portion (eq. 10), two iterations were sufficient for
convergence.
For the three zone portion (eq. 18), four iterations were necessary.
Numerical Solution
A number of numerical methods are available
to solve a set of differential equations like the
ones outlined above.
Both the Runge-Kutta technique and a simple
Euler technique were evaluated.
The two methods provided surprisingly similar
results.
Since the application of the Euler technique
was more straight ward, it was selected for
this work.
Numerical Solution
To complete the required input information, the boundary
conditions for the inlet (temperature and pressure) and for the
exhaust (pressure) are specified.
To begin a particular engine cycle calculation, several
parameters are not known.
The initial amount of exhaust gases (residual), and the initial
cylinder gas temperature and pressure must be assumed.
The complete calculation is repeated until the final values agree
with the initial values.
Depending on the specified tolerance, this procedure usually
finds convergence within about three (3) complete cycles.
Ref A and B are Heywood’s publications
Input Engine Data
Table 1 lists the engine specifications as provided by
Heywood et al. (1979), and as used in the present
work (current).
Using the listed valve timings resulted in excessive
back flow and unrealistic flow characteristics during
the valve overlap period.
To avoid excess back-flow, the amount of valve
overlap was reduced.
As shown, the intake valve open and close crank
angles, and the exhaust valve close crank angle were
modified.
The overall valve overlap was decreased from the
apparent 40°to 14°.
This provided much more realistic flow rates which
agreed with those published (see following).
The reasons for the need for this valve timing
modification included the fact that the original
modeling work (Heywood, et al., 1979) used actual
cam profiles whereas the current work used a
“sinusoidal” cam profile (since the actual cam profiles
were not available).
Further, the actual valve timings used in the original
modeling work (Heywood,et al., 1979) may have
differed from those which were listed as the engine
specifications
Heat Transfer
Finally, the table lists an empirical heat
transfer multiplier which was set to 1.33 to
provide good agreement with the published
results.
The need for this “multiplier” may be due to a
difference in the implementation of the
original heat transfer algorithm, or a similar
constant used by the original investigators,
but not reported.
Using this “multiplier,” the final results for the
heat transfer are in good agreement with the
published results.
Engine Input Data
Table 2 lists other parameters (and their source)
which were needed in this work.
Although most of the significant engine parameters
were duplicated in an exact fashion by the current
study, a few parameters were estimated, inferred or
interpreted from figures, or modified.
The parameters affected included items such as the
valve timings, the valve lift profiles, the exhaust and
inlet pressures, the unused fuel, and the cylinder wall
temperature.
In general, modest changes in these parameters will
not significantly change the results of the current
study.
RESULTS AND
DISCUSSION
Results from this simulation are described in
the following four subsections:
(1) comparison of basic results,
(2) numerical accuracy results,
(3) basic engine performance results, and
(4) detailed thermodynamic results
Basic Results
With the above input parameters, the simulation was executed
for the part load base-case condition: 1400 rpm and 52 kPa
(inlet pressure).
Table 3 compares the major output parameters from the current
calculations with those listed in Heywood et al. (1979).
In general, the agreement is good.
In particular, the global engine performance parameters such as
torque, power, mean effective pressure (MEP), specific fuel
consumption, thermal efficiency, mass averaged exhaust
temperature, and the overall air and fuel flow rates were in
good agreement for both brake and indicated values.
In addition, the distribution of the energy between brake work,
friction work, heat loss, exhaust energy, and unused fuel energy
was adequately duplicated
Results
Figure 3 shows four of the computed gas temperatures as a
function of the crank angle for the base line conditions
Also shown are results (symbols) from Heywood et al. (1979)
for these temperatures.
As shown, the computed values of the current study are in good
agreement with the values obtained from the figure from
Heywood et al. (1979).
The major difference between the current computations and
those of Heywood et al. (1979) are a slight over-prediction of
the adiabatic zone temperatures between 20° and 40°aTDC.
These differences may be due to interpretation of the values
from the figure, due to different valve timings (see above), due
to different valve lift profiles, and possibly due to slight
differences in the formulation.
In light of these possible discrepancies, the level of agreement
provides a high level of confidence that the current simulation is
correctly constructed and implemented
Gas temperature
The overall burned gas temperature (Tb) is
the average of the adiabatic core temperature
(Ta) and the boundary layer temperature
(Tbl).
This is an energy-based average; that is, the
energy represented by total burned zone is
the same as the sum of the energy in the
adiabatic core and in the boundary layer.
Results
Figure 4 shows the burned gas temperature, the unburned gas
temperature, the average cylinder gas temperature, and cylinder
pressure.
Also, the figure includes values from Heywood et al. (1979) for
the cylinder pressures.
Again, the agreement is good.
The average gas temperature (Tavg) is the equivalent of a one-
zone temperature and results in the same total energy that
exists in the burned and unburned zones.
The average gas temperature transitions from the unburned gas
temperature to the burned gas temperature during the course
of the combustion process.
Results
A heat transfer rate parameter may be defined as the
instantaneous heat transfer rate divided by the total fuel energy
per combustion event.
For the base case, the fuel energy per combustion event was
0.936 kJ.
Figure 5 shows this heat transfer rate parameter as a function
of crank-angle for the base case.
Results for this parameter from both the current work and from
Heywood et al. (1979) are shown.
The highest heat transfer rate occurs during the latter stages of
the combustion process at about 20°aTDC.
Results
A physical picture of the developing boundary
layer is useful.
Figure 6 shows a schematic illustration of the
rough relation of the boundary layer thickness
and the cylinder dimensions at 90° aTDC.
As shown, the boundary layer thickness is a
small faction of the bore dimension, but due
to the lower temperature of the boundary
layer, the mass in the boundary layer will be
greater than that indicated by the simple
volume fraction.
Numerical Accuracy Results
One of the specific goals of this work was to implement a high level of
accuracy in the overall simulation.
A “computational error” was determined for each case2.
This “error” was determined by summing the various energy transfers
throughout the calculation.
These energy transfers included the work, heat loss,inlet energy,
exhaust energy, and any unused fuel.
Since total (chemical and sensible) energy was used in the simulation,
the sum of the various energy transfers should equal zero.
Any value different than zero was deemed a “computational error.”
During the development of the simulation, the computational error
helped signal formulation or implementation problems.
Once the simulation was performing correctly, the computational error
was monitored to insure appropriate choices of computational
variables.
The characteristics of the computational variables are described next.
Numerical Accuracy Results
Two of the major computational variables are the
computational step size (number of calculations per
crank angle) and the minimum temperature at which
the equilibrium composition was used.
As mentioned above, this minimum temperature was
called the “switch” temperature, since the
computations would switch to using or not using the
equilibrium composition at this temperature.
In addition to a computational error, a relative
computational time was also determined.
This relative computational time was the ratio of the
actual computational time to the computational time
for a reference case.
Numerical Accuracy Results
The computational error as a function of the “switch” temperature was
determined for the base case conditions for sixteen
(16) calculations per crank angle (a computational step size of 0.0625 crank
angles).
As the “switch” temperature decreased from 2000 to 1000 K, the computational
error decreased. For a “switch” temperature of 1000 K, the computational error
was about 0.007%.
The complex behavior of the computational error is largely due to the
importance of the partial derivatives which have finite values even down to 1000
K.
On the other hand, as shown elsewhere the thermodynamic properties for the
most part are essentially unchanged for temperatures below about 1700 K.
Also, note that the relative calculation time increases as the “switch”
temperature decreases.
The calculation time increases by a factor of about 2 for a decrease in the
“switch” temperature from 1400 to 900 K.
Numerical Accuracy Results
Figure 7 shows the absolute computational error as a function of the
number of calculations completed per crank angle for a “switch”
temperature of 1000 K for the base case condition.
As shown in this figure, the computational error decreases as the
number of calculations increase.
For the fewest number of calculations per crank angle (1), the error is
about 0.6%, while for two calculations per crank angle, the error is
about 0.4%.
Sixteen (16) calculations per crank angle resulted in the lowest error
(0.007%) for the conditions studied.
Also shown is the increase in the relative computational time.
As expected, the relative computation time increases linearly with
increases in the number of calculations per crank angle.
For example, increasing from 1 to 32 calculations per crank angle
increases the computational time by a factor of 2
Basic Engine Performance
Results
This subsection will present basic engine
performance results In particular the brake
specific fuel consumption (bsfc), the thermal
efficiency, and the mass average exhaust
temperature will be presented as functions of
spark timing, equivalence ratio, compression
ratio, and burn duration.
Most of the definitions of
Basic Engine Performance
Results
This subsection will present basic engine performance results In
particular the brake specific fuel consumption (bsfc), the
thermal efficiency, and the mass average exhaust temperature
will be presented as functions of spark timing, equivalence ratio,
compression ratio, and burn duration.
Most of the definitions of these parameters are well known, but
the use of the mass average exhaust temperature may need to
be defined.
The exhaust gas temperature that is the most representative of
the exhaust gas energy is an average temperature based on the
energy of the gases flowing out the exhaust port.
Exhaust Gas Temperature
If the energy is
characterized by the
enthalpy, then the
average exhaust
temperature may be
determined from
Exhaust Gas Temperature
The temperatures that result from this calculation will
typically be higher than measured exhaust gas
temperatures of the exhaust.
This is because measurements typically are more
closely related to time-averaged values, and are
generally obtained downstream of the exhaust valve
exit plane.
A more complete description of the differences
between the mass-averaged, time averaged, and
thermocouple exhaust temperatures may be found
elsewhere.
Results
Before examining basic engine performance results,
the cylinder pressure is considered.
Figure 8 shows the logarithm of the cylinder pressure
as a function of the logarithm of the cylinder volume
for the base case.
Also shown on this figure are the locations of the four
valve openings and closings, and the start and end of
the combustion process.
This type of figure clearly shows the pressure
changes during the exhaust and intake stokes.
Further, the nearly constant slopes of the curves
during compression and expansion are shown.
Results
The next set of results are based on
examining the effect of one major variable
while the other variables are held constant.
In particular, the engine load is constant.
To maintain constant load, the intake
manifold pressure (throttle position) is
adjusted.
Furthermore, the following results are for
MBT timing (except where spark timing is the
independent variable).
Effect of [Link]
First, the effect of compression ratio on several of the
engine performance parameters will be examined.
Note that the following results for compression ratio
changes are based on the assumption that the
engine does not knock.
Figure 9 shows the brake specific fuel consumption
(bsfc) as a function of compression ratio for four burn
durations for the base case.
As shown, the bsfc decreases as compression ratio
increases, but the differential gains are decreasing as
compression ratio continues to increase.
Also, the bsfc increases for shorter burn durations.
Effect of Comp. Ratio
Figure 10 shows the brake thermal efficiency as a function of
compression ratio for four combustion durations.
These results are proportional to the inverse of the results in figure 9.
Figure 11 shows the exhaust gas temperature as a function of
compression ratio for three combustion durations for the base case
Although the agreement is not as good as before, the maximum
difference is less than 3%.
This appears to be a result of some differences in the computed
instantaneous exhaust flow rates.
Since the overall engine performance parameters are in good
agreement, this slight discrepancy is tolerable.
The exhaust temperatures decrease as compression ratio increases
since more of the thermal energy is used for work.
Similarly the EGT decreases as the comb. duration is decreased as
more energy is used for work
Results
The next figure will examine some of the details of
the combustion process, and in particular, the
differences of the three temperatures as one of the
major parameters is varied.
Figure 12 shows the maximum temperatures for the
adiabatic core (Ta), for the burned gas (Tb), and for
all the cylinder gases (Tavg) as a function of
compression ratio for a combustion duration of 60°
for the base case.
As shown, the maximum adiabatic zone gas
temperature is higher than the maximum burned
zone gas temperature by about 40 K, and this
difference is roughly constant for the different
compression ratios.
Results
As mentioned above, one implication of these
higher temperatures is the effect on nitric
oxide predictions.
Although nitric oxide calculations were
outside the scope of the current work,an
estimate may be completed to illustrate the
impact of these higher temperatures.
From Heywood (1988), the rate of nitric oxide
production roughly scales as
The approximate increase of the nitric oxide production rate
for
temperature of 2600 K compared to 2560 K is about 52%.
This is an instantaneous value, and would not apply to the
whole
combustion event. Nevertheless, the computation of the
adiabatic
core temperatures is, important for accurate determinations
of nitric oxide concentrations
Effect of Spark Timing
The relative spark timing is the timing relative to MBT timing for
each case.
Figure 13 shows the brake specific fuel consumption (bsfc) as a
function of the relative spark timing for three combustion
durations for the base case.
Also shown are the results for these cases from Heywood et al.,
1979.
The agreement is quite satisfactory.
Bsfc increases for both retard and advance timings relative to
the MBT timing.
Further, the bsfc increases as the combustion duration
increases.
Note that the increase from a combustion duration of 60° to
100° is greater that the increase from 20° to 60°.
Results
Figure 14 shows the exhaust gas temperature as a function of
the relative spark timing for three values of the combustion
duration for the base case.
The agreement is acceptable.
As the relative spark timing is retarded from MBT timing the
exhaust temperature increases since less of the energy was
used in producing work.
As the relative spark timing is advanced from MBT timing, the
exhaust temperature decreases in spite of the fact that less
energy is used to produce work (see figure 12).
This is because, for advanced timings, the higher combustion
temperatures result in higher heat losses and hence lower
energy content in the exhaust gases
Summary of simple combustion
model
In simple model, different heat transfer
surface areas are not considered
Simple combustion model has been employed
(Wiebe)
Emission prediction is not complete
Gas exchange process is not undertaken
Number of assumptions made in the model
are quite large
This type of model is quite useful for
preliminary investigations
Improved Combustion Models
The combustion of fuel in an engine is a chemical
process influenced by many parameters.
One of these is the ratio between air and fuel.
If more air is available than required to burn the fuel
completely the combustion is called lean.
The opposite is called rich combustion.
The ratio between air and fuel at which neither
unburned fuel nor air remains after combustion is the
stoichiometric air fuel ratio.
The following equation for the stoichiometric air
requirement specifies how much air is required for a
complete combustion of 1 kg fuel:
Variables
Hu lower heating value
c mass fraction of carbon in the fuel
h mass fraction of hydrogen in the fuel
o mass fraction of oxygen in the fuel
s mass fraction of sulfur in the fuel
n mass fraction of nitrogen in the fuel
w mass fraction of water in the fuel
In rich combustion, the total heat supplied during the cycle is
limited by the amount of air in the cylinder.
The fuel is totally converted to combustion products even if the
amount of air available is less than the amount of stoichiometric
air.
Gas Composition
However, the composition of the combustion products is
different if fuel is burned under rich or lean conditions.
The composition itself depends on the type of fuel used, the air
fuel ratio, pressure and temperature.
It is always the same if sufficient time is available to reach
chemical equilibrium.
It is well known that under real engine conditions, complete
combustion as assumed above can never be achieved.
This is very important for excess air ratios close to 1.0 (the
excess air ratio is defined as the ratio between the amount of
air in the cylinder and the amount required for stoichiometric
combustion).
Quasi-dimensional
Combustion Models
The quasi-dimensional combustion model for SI
engines implements the rate of heat release in a
homogeneous charge engine.
Thereby the influence of the following parameters is
considered
. The combustion chamber shape
• The spark plug location and spark timing
• The composition of the cylinder charge (residuals,
recirculated exhaust gas, air and fuel vapor)
• The macroscopic charge motion and turbulence
level
The thermodynamics of the two zone combustion
model is outlined earlier
Flame Propagation
A spherical propagation of the flame from the
spark plug through the combustion chamber
is assumed.
With this assumption the instantaneous flame
radius, the flame area and the wetted piston,
head and liner surface of the burned and
unburned zone can be determined from
purely geometrical considerations.
Entrainment and depletion
However, experimental evidence shows that
both entrainment and depletion of fuel do not
strictly follow geometrical correlations.
This is corrected by an approach by Bargende
, which improves the heat release versus
crankangle.
Therefore, the following correction is
introduced into Equation 2.1.22.