- Email: [email protected]
Modelling cycle to cycle variations in an SI engine with detailed chemical kinetics
Combustion and Flame 158 (2011) 179–188
Contents lists available at ScienceDirect
Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e
Modelling cycle to cycle variations in an SI engine with detailed chemical kinetics Jonathan Etheridge a, Sebastian Mosbach a, Markus Kraft a,⇑, Hao Wu b, Nick Collings b a b
Department of Chemical Engineering and Biotechnology, University of Cambridge, New Museums Site Pembroke Street, Cambridge, CB2 3RA, UK Department of Engineering, University of Cambridge, Trumpington Street, Cambridge, CB2 1PZ, UK
a r t i c l e
i n f o
Article history: Received 8 January 2010 Received in revised form 26 May 2010 Accepted 16 August 2010 Available online 24 September 2010 Keywords: SI engine modelling Detailed chemistry Cycle to cycle variation
a b s t r a c t This paper presents experimental results and a new computational model that investigate cycle to cycle variations (CCV) in a spark ignition (SI) engine. An established stochastic reactor model (SRM) previously used to examine homogeneous charge compression ignition (HCCI) combustion has been extended by spark initiation, flame propagation and flame termination sub-models in order to simulate combustion in SI engines. The model contains a detailed chemical mechanism but relatively short computation times are achieved. The flame front is assumed to be spherical and centred at the spark location, and a pent roof and piston bowl geometry are accounted for. The model is validated by simulating the pressure profile and emissions from an iso-octane fuelled single cylinder research engine that showed low CCV. The effects of key parameters are investigated. Experimental results that show cycle to cycle fluctuations in a four-cylinder naturally aspirated gasoline fuelled SI engine are presented. The model is then coupled with GT-Power, a one-dimensional engine simulation tool, which is used to simulate the breathing events during a multi-cycle simulation. This allows an investigation of the cyclic fluctuations in peak pressure. The source and magnitude of nitric oxide (NO) emissions produced by different cycles are then investigated. It was found that faster burning cycles result in increased NO emissions compared with cycles that have a slower rate of combustion and that more is produced in the early stages of combustion compared with later in the cycle. The majority of NO was produced via the thermal mechanism just after combustion begins. Ó 2010 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
1. Introduction Nitric oxides (NOx), unburned hydrocarbons (UHC) and carbon monoxide (CO) are major components of air pollution caused by Internal Combustion Engines (ICE). Increasingly stringent restrictions on vehicle emissions have led to much research into reducing these emissions. Spark ignition (SI) engine NOx emissions have been studied both experimentally [1] and with models of varying complexity [2–5]. The high temperatures reached in SI engines cause the conversion of atmospheric nitrogen to nitric oxide (NO). The three accepted mechanisms responsible for NO formation are the thermal, prompt and nitrous oxide mechanisms. The majority of NO is formed in the burned gases behind the flame front, where the temperatures remain high for enough time to produce NO by the thermal mechanism that has relatively slow rate constants. The prompt mechanism produces NO in the flame, although the amount is thought to be much lower than in the end gas.
⇑ Corresponding author. E-mail address: [email protected] (M. Kraft). URL: http://como.cheng.cam.ac.uk (M. Kraft).
Empirical equations exist to predict NO emissions from SI engines, however these are extremely simplified. SI engine NO emissions have been modelled previously by assuming that equilibrium is reached instantaneously behind the flame front for all species, with the exception of NO due to its slow rate of formation. Chemical kinetics are then used in the hot end gas to calculate NO production [6]. Improved predictions have been made by splitting the end gas into multiple zones that have burnt at different times, taking into account the difference in temperature of the charge that burns earlier and later in the engine cycle [4,7]. As well as different points in a cycle being responsible for varying amounts of NO, fluctuations between cycles will result in cycle to cycle variations in NO emissions [8]. Cycle to cycle variations (CCV) in SI engine combustion are a common problem which reduce performance due to variations in combustion phasing, and, in some circumstances, can lead to knocking or misfire. Several causes of CCV are known and the magnitude of their effects varies between engine design and operating conditions. In-cylinder charge motion can vary between cycles, altering the rate at which unburned mixture is entrained and the shape of the turbulent flame front. The total mass of fuel, air and exhaust gas recirculation (EGR) can vary between cycles, as well as local compositions in the cylinder, which are especially
0010-2180/$ - see front matter Ó 2010 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.combustflame.2010.08.006
180
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
List of symbols and abbreviations Symbols C CADspark lT NB NE NU NTot Pmax Rf;n SL t tspark tstop Dt i u uT VB VE Vf V Unext
qi qu sb sE
smax E entrainment rate constant (–) spark timing angle (CAD ATDC) characteristic length scale (m) number of burned particles (–) number of entrained particles (–) number of unburned particles (–) total number of particles at IVC (–) peak pressure (bar) flame radius at nth CAD step (m) laminar flame speed (m s1) time after IVC (s) spark time (s) time at EVO (s) time step (s) mean inlet gas speed (m s1) characteristic flame speed (m s1) burned zone volume (m3) entrained zone volume (m3) flame volume calculated from the radius (m3) volume of the next particle to be entrained (m3) inlet air density (kg m3) unburned zone density (kg m3) characteristic burn time (s) entrained zone mixing time (s)
important near the spark location at the time of ignition. Variations in local turbulence at the time of spark can lead to a movement of the flame kernel altering the interactions between the flame and the piston and cylinder walls. Many studies of these effects have been completed. In [9], variations in the peak pressures of two different engines were simulated by random fluctuations, with a normal distribution, in the rms turbulent velocity. Fluctuations in injected fuel mass were also given as a possible cause of the measured CCV. Experimental results have also suggested that CCV is not very dependent on spray fluctuations, and the in-cylinder flow just after the spark timing has a greater effect [10]. Convection of the flame kernel has also been investigated and proposed to be a key factor [11,12]. Movement of the flame kernel and the time delay before the flame becomes fully turbulent have also been discussed in [13] as effects of variations in the rms turbulent velocity. Combustion in SI engines takes place in a turbulent flow field. The structure and speed of the flame depend on charge motion, composition, temperature and chamber geometry. Detailed SI models exist that predict the structure and propagation of the turbulent flame [14,15]. To run many cycles of a computational fluid dynamics (CFD) model with detailed chemistry however would be extremely computationally expensive. Less complex models exist that assume the flame is spherical and propagates at the turbulent flame speed from the spark [16–18]. In these models the cylinder is divided into unburned and burned zones and the burned zone may be subdivided into several zones that have burned at different times. The burned zone composition is frequently calculated by assuming the mixture behind the flame front reaches equilibrium. In this paper we present a new probability density function (PDF) based SI model that uses a detailed chemical mechanism throughout the cycle. The model has been created from a previously developed stochastic reactor model (SRM) successfully used to model Homogeneous Charge Compression Ignition (HCCI) combustion. The model has been validated by simulating a single-cylinder research engine which demonstrated low cycle to cycle variation. The model was then coupled with GT-Power, a
maximum entrained zone mixing time (s)
Abbreviations BTDC before top dead centre CA50 crank angle at 50% heat release CAD crank angle degree CCV cycle to cycle variation CFD computational fluid dynamics CR compression ratio EGR exhaust gas recirculation EMST Euclidean minimum spanning tree EVO exhaust valve open FID flame ionisation detector GDI gasoline direct injection HCCI homogeneous charge compression ignition IVC inlet valve close NDIR non-dispersive infrared PDF probability density function PRF primary reference fuel RON research octane number RPM revolutions per minute SI spark ignition SRM stochastic reactor model UHC unburned hydrocarbons
one-dimensional engine simulation tool, which was used to simulate the breathing events during multi-cycle simulation. Experimental results from a four-cylinder naturally aspirated gasoline fuelled SI engine were modelled and the cyclic fluctuations in pressure and NO emissions investigated. The paper is structured as follows. In the next section the model details and implementation is reported. A description of both engine set-ups is given at the start of Section 3. Calibration of the model using experimental data from a single-cylinder research engine is then presented. Subsequently the model is calibrated using experimental data from fast, slow and average cycles from a fourcylinder SI engine. Results from a multi-cycle simulation used to model cyclic fluctuations are contained in the same section along with experimental measurements of emissions. Finally conclusions are drawn and future work outlined. 2. Model details The model was built into the framework of the SRM used to simulate HCCI combustion [19–23]. The SRM is based on the PDF transport equation, and models convective heat transfer [24], turbulent mixing, piston movement and detailed chemistry. The PDF is approximated by dividing the charge into a notional ensemble of NTot stochastic particles. For SI simulation the particle ensemble is divided into three zones: unburned, entrained and burned, which contain NU, NE and NB particles respectively. Mixing occurs within each zone, but not between zones and all particles are initially in the unburned zone. The Euclidean minimum spanning tree (EMST) mixing model was used in the unburned and burned zones and is described in [25]. A Curl based mixing model was used in the entrained zone and will be described in more detail later [26]. The entrained zone is subdivided into unburned and burning particle categories to keep track of the particles’ states, as will be discussed later. Chemistry occurs at every CAD step in every particle, as in the HCCI model. The primary reference fuel (PRF) mechanism used
181
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
contains 157 species and 1552 reactions [19–23]. The model could be used to investigate knocking as in [27], however, similar work is not presented here. The main SRM operator splitting loop, with the addition of spark initiation, flame propagation and flame termination, was implemented as follows: 1. Initialise t = 0, Dt, CAD = IVC, CADspark, spark energy, Rf = 0, NU = NTot, NE = NB = 0 all particles are in the unburned zone. Determine temperature, composition, mass, volume and pressure of particle ensemble. 2. Progress in time t ´ t + Dt. If CAD P EVO or t P tstop then save the detailed exhaust composition as input EGR and stop. 3. Perform volume change due to piston movement. 4. Perform the mixing step. 5. Perform stochastic heat transfer. 6. Perform the pressure equilibration step [27]. 7. If CAD P CADspark and NU = NTot perform spark initiation (as described below) else if CAD P CADspark and NU > 0 perform flame propagation (as described below) else go to (8). 8. Perform the chemistry step. If NE > 0 perform flame termination (as described below). 9. Perform the pressure equilibration step [27]. 10. Perform the mixing step. 11. Go to step (2). 2.1. Spark initiation During flame propagation, particles are moved from the unburned zone to the entrained zone, according to the required volume of the flame. This is a discrete process, and the volume of the entrained and burned particles is unlikely to match the calculated flame volume exactly. During spark initiation, the initial flame volume is very small. If 100 particles are used, the initial entrained volume is much lower than the volume represented by a single particle. If the energy in the spark were to be distributed uniformly throughout a volume represented by a particle, then the temperature would not be sufficient to initiate combustion. For this reason, the first particle randomly chosen to enter the entrained zone is not moved directly into the zone, but the required entrained volume is removed from the particle and added to a newly created particle in the entrained zone. Once the calculated entrained volume has increased above the particle’s volume, it is reunited into a single particle in the entrained zone. The diameter of the entrained zone at the time of spark was set to 1.0 mm, the rough size of the spark gap. The spark is simulated by setting the temperature of the mass added to the entrained zone, at each time step, to 2000 K, until all of the spark energy has been added. This low temperature was chosen to prevent any problems with the chemical mechanism, which was not developed to represent the high temperatures experienced at the spark. The spark initiation was modelled as follows: 1. If its the first call to the spark initiation subroutine, continue to (2) else go to (7). 2. Set the flame radius to Rf = 0.5 mm. 3. Select a particle at random with a uniform distribution. 4. Remove the spherical volume with 0.5 mm radius from the particle and add to a newly created particle in the entrained zone. NE = 1, V E ¼ 4=3pR3f . 5. Increase the temperature of the entrained particle to 2000 K by adding energy at constant pressure. Remove the energy added from the initial spark energy. 6. Go to the chemistry step. 7. Calculate the new flame volume, Vf, as described in the flame propagation section.
8. If V f > V UNext þ V E go to (13) else go to (9). 9. Remove Vf VE from the unburned particle being entrained. 10. If there is spark energy remaining, increase the temperature of the freshly entrained mass to a maximum of 2000 K by adding energy at constant pressure. 11. Reunite the freshly entrained mass with the entrained particle. 12. Go to the chemistry step. 13. Remove the unburned particle being entrained from the unburned zone, NU = NTot 1. 14. If there is spark energy remaining, increase the temperature of the freshly entrained mass to a maximum of 2000 K by adding energy at constant pressure. 15. Reunite the freshly entrained mass with the entrained particle. 16. Select the next particle to entrain, at random with a uniform distribution. Go to step (2) in flame propagation section described below. 2.2. Flame propagation A randomly chosen particle is moved from the unburned zone to the entrained zone, when the calculated flame volume exceeds the burned and entrained zone volumes by more than half the volume of the next particle to be entrained. The volume of the flame is calculated from the geometry of the piston and cylinder, and the flame radius, assuming it is spherical and centred at the spark location. The flame radius at the nth time step, Rf;n, is obtained from
ðt t spark Þ Rf ;n ¼ Rf ;n1 þ uT 1 exp þ SL Dt:
sb
ð1Þ
The increase in flame radius is calculated as the time step multiplied by the entrainment velocity. The flame radius is also updated after the pressure equilibration step, to take into account the expansion of the flame. The radius initially increases at the laminar flame speed, SL, taken from [28], until a wrinkled turbulent flame front develops after a characteristic burn time, sb [29,28]. The characteristic burn time is given by
sb ¼
lT ; SL
ð2Þ
where lT is the characteristic length scale [28]. One of the most desirable features of the EMST mixing model is that it preferentially mixes particles which are close in composition space. This is based on the reasoning that fluid parcels which are close in physical space should be similar in terms of composition. The flame brush is relatively thin with high temperature and concentration gradients across it, so EMST is not a good mixing model to use in the entrained zone. A mixing model based on Curl’s has been used in the entrained zone [26], whereby two particles are selected to mix. The difference being that only unburned entrained particles are chosen to mix with entrained and burning particles. The resulting particle compositions are not identical, as the particles are only moved to their mean composition by a factor of 10%, and the mixing time is accordingly adjusted to reflect this. This method was chosen to prevent the entrained zone from becoming unphysically homogeneous. The mixing time in the entrained zone, sE, is increased at the same exponential rate as the transition from a laminar to a turbulent flame speed. The mixing time increases to the value smax once E the turbulent flame has fully developed. Before this, the entrained zone mixing time is smaller causing faster mixing and heat release. This was included as the rate of burning is initially the rate of entrainment before a turbulent flame front has developed containing unburned pockets [13].
182
sE
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
ðt t spark Þ : ¼ smax 1 exp E
sb
ð3Þ
The rate of burning is controlled by the mass entrained, and the rate it mixes with the hot, reacting entrained material, and is calculated with the detailed chemical mechanism. The maximum entrained zone mixing time was smaller than the other zones, to prevent the volume of the entrained zone growing unrealistically large, and to propagate heat through the entrained zone at the correct rate. Flame propagation was implemented in the model as follows: 1. Increase the flame radius using Eq. (1) and calculate the new flame volume, Vf. 2. Update the entrained zone mixing time according to Eq. (3). P E E PNB B 3. If V f > Ni¼1 V i þ i¼1 V i þ 1=2V UNext go to (4) else go to (6). 4. Add the particle to the entrained zone and select the next particle to be added randomly with a uniform distribution. NU = NU 1, NE = NE + 1. Flag the newly entrained particle as unburned and entrained. 5. Go to step (3). 6. Perform the chemistry step. The correlation used to obtain the characteristic flame speed, uT, was taken from [18]:
i uT ¼ 0:08C u
qu qi
1=2 ;
ð4Þ
i is the mean inlet gas speed, qu is the unburned gas density where u and qi is the inlet air density. The correlation used is relatively old, however this work was focused on producing a detailed chemistry SI model incorporated into the SRM, to enable multi-cycle simulation at low computational cost. There were no measurements of in-cylinder turbulence, hence the constant C was used to fit the model to the experimental data by varying the rate of entrainment. The model was intended to predict not only the pressure profile, but also temperature and composition. The cylinder geometry model contains a pent roof and the piston contains a bowl. The volume of the spherical flame is calculated exactly in the bowl and cylinder areas. In the pent roof region, numerical integration of the flame cross sectional area at different heights is used to obtain the flame volume. 2.3. Flame termination In this section we describe how particles are moved from the entrained zone to the burned zone. Particles move from the entrained and unburned zone to the entrained and burning zone once their temperature has risen above 1500 K. During the chemistry step the rate of heat release in each entrained particle is calculated. When a particle’s heat release rate has dropped below 2% of its maximum value, it is moved to the burned zone. When no particles remain in the unburned zone, the entrained zone mixing time is set to smax as Eq. (2) can no longer be evaluated. Flame termination E was modelled as follows: 1. If any unburned entrained particle’s temperature has risen above 1500 K, flag the particle as entrained and burning. 2. If any entrained and burning particle’s rate of heat release has dropped below 2% of its maximum, move the particle to the burned zone. NE = NE 1, NB = NB + 1. 3. If NU = 0 then sE ¼ smax E . In summary, an SI model has been integrated into the SRM. The particle ensemble is divided into three zones: unburned, entrained and burned. Mixing occurs within each zone, but not between zones. Chemistry occurs at every CAD step in every particle, as in
the HCCI model. An initial temperature and flame radius of 2000 K and 0.5 mm are set at the spark timing respectively. The flame radius is updated using an entrainment velocity. The flame radius and cylinder and piston geometry are used to obtain the flame volume. Particles are moved from the unburned zone to the entrained zone when the calculated flame volume exceeds the burned and entrained zone volumes by more than half the volume of the next particle to be entrained. Entrained particles, initially in the entrained and unburned zone, are moved to the entrained and burning zone once their temperature has risen above 1500 K. Entrained and burning particles are moved to the burned zone when their heat release rate has dropped below 2% of its maximum value. 3. Results and discussion In this section we present and compare experimental and simulation results from two different engines: a single-cylinder naturally aspirated direct injection iso-octane fuelled research engine and a four-cylinder naturally aspirated GDI engine. The engine specifications and operating conditions are given in Table 1. Exhaust emissions of CO and UHC, as well as the in-cylinder pressure, were measured in the single-cylinder research engine. The CO volume fraction was measured with a non-dispersive infrared (NDIR) sensor. UHC were measured 200 mm downstream of the exhaust valve with a fast flame ionisation detector (FID). The pressure in each cylinder of the four-cylinder GDI SI engine was monitored with Kistler 6123 piezoelectric transducers and recorded. Emissions of NOx were measured using a Horiba EXSA1500 analyser. More details of the engine can be found in [30]. 3.1. Model calibration The model was initially calibrated using results from the single cylinder research engine. As fuel injection occurred early during intake, the cylinder charge was assumed to be initially homogeneous in terms of composition and temperature, hence all particles were set to be identical at Inlet Valve Closing (IVC). The model was calibrated to find values for C and smax E . The values of these parameters are given in Table 2. The pressure profile was matched to the experimental data shown in Fig. 1. The emissions matched experimental data reasonably well, as shown in Table 3. No NOx measurements were taken. Fig. 2 shows the particles remained in the entrained zone for between 0 and 1.75 ms. The PDF was obtained by taking each value from the ensemble as the mean of a normal distribution, with a r of 0.05 ms, summing the distributions and dividing by the total number. From experimental observations, the burn time in the turbulent flame front has been estimated at values ranging from 1.0 to 1.8 ms [13]. The PDF becomes smoother when the particle number
Table 1 SI mode engine specifications and operating conditions. Engine
Research
Production
Cylinders Fuel Bore (mm) Stroke (mm) Con. rod length (mm) Disp. volume (cm3/cyl) CR Speed (rpm) Air/fuel equiv. ratio EGR (%) Spark timing (CAD BTDC) Start of injection (CAD BTDC)
1 iso-octane 89.0 90.3 148.9 562 11 1500 1.0 15.0 35 280
4 95 RON gasoline 87.5 83.0 146.3 499 12 1500 1.0 28.8 40 308
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188 Table 2 Parameters used in simulation of singlecylinder research engine. No. of particles
100
CAD step (CAD) C smax (ms) E
0.2 9.0 1.7
Fig. 1. Comparison of experimental and simulated pressure profiles for the singlecylinder research engine.
Table 3 Emission results for single-cylinder research engine.
UHC (ppm) CO (ppm) NOx (ppm)
Exp.
Sim.
600 2700 –
746 3185 671
183
is increased. The PDF of entrained time is skewed. The shape is caused by the time particles are entrained for, before they begin to mix. This is affected by the entrained zone mixing time, which increases during the simulation, so particles entrained later in the cycle are likely to spend longer in the entrained zone. Fig. 3 shows the average rate of heat release against CAD, calculated by centering each particles’ peak heat release rate at 0CAD and averaging over the 100 particles. The rate of heat release increased rapidly, and the majority of the chemical energy was released within 2CAD. The number of particles entrained and burning remains roughly constant, however, the number of particles entrained and unburned grows and then falls off (see Fig. 4a). This indicates a thickening flame brush, followed by consumption of the remaining fuel as the flame reaches the walls. The temperatures remain roughly constant throughout (see Fig. 4b). The unburned entrained region is slightly warmer than the unburned zone, due to some particles having mixed with the burning entrained particles. The burning entrained particles are slightly cooler than the burned particles, as all the chemical energy has not been released in this zone, and particles can mix with the cooler unburned entrained particles. A spark temperature of 2000 K was chosen, as the mechanism was not created for the extremely high combustion temperatures that are reached within a spark plug gap during the passage of a spark. Fig. 5a suggests this is an acceptable approach, as there does not appear to be a trend between the spark temperature and the pressure profile at the conditions investigated. Fig. 5b shows that increasing the initial flame radius caused a slight advance in combustion phasing. The effect is not significant, given the relatively large range investigated. As variations in this parameter caused minor changes in the pressure profile, the initial flame radius assumption was taken as acceptable and a detailed spark initiation sub-model, e.g. [31], was not implemented. The simulation was repeated with n-heptane as the fuel. This was done without changing the laminar flame speed constants, in order to investigate the sensitivity of the combustion to the entrained zone composition. The initial burning rate was faster, and also led to knocking, as shown in Fig. 6. 3.2. Multi-cycle simulation with CCV A four-cylinder SI engine displayed high CCV at the speed and load investigated. Measurements suggested an EGR ratio of 30%
Fig. 2. Probability density of particles’ time in the entrained zone.
Fig. 3. Average rate of heat release from each particle against CAD, normalised so that all peak rates were at 0 CAD.
184
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
(a)
(b)
Fig. 4. Number of particles and the average temperature of each zone against CAD.
(a)
(b)
Fig. 5. Comparison of pressure profiles using different spark temperatures and initial flame radii.
Fig. 6. Comparison of in-cylinder pressure profiles with iso-octane and n-heptane fuel.
at this operating point. The experimental method uses the ratio of CO2 in the exhaust and compression stroke, and is described in [30]. Large amounts of internal EGR have been reported to increase CCV [28]. Details of the engine are given in Table 1. A mixture of 95% iso-octane and 5% n-heptane was used to simulate the 95 Research Octane Number (RON) gasoline. The model was calibrated to fit an experimental pressure profile from a cycle with a peak pressure close to the average as shown in Fig. 7. The model parameters used are given in Table 4. Fig. 7 shows experimental pressure against crank angle, and the simulation results when C was varied and smax was constant. The E pressure traces match the experimental data well, with C values of 11 and 22 used to match the slow and fast pressure traces respectively. The cycles in which the flame initially entrains mass at a faster rate, continue to entrain mass at a fast rate throughout the cycle. This suggests variations in the flame surface area occur early in the cycle resulting, in different rates of unburned gas entrainment and combustion. This could be caused by the flame kernel being convected towards a wall, reducing the available area, or due to turbulence altering the shape of the kernel early in the cycle, resulting in a higher surface area. Fig. 8 shows the PDF of peak pressure obtained from the experimental data, and a normal distribution fitted to the data. The peak
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
185
Table 4 Parameters used in simulation of four-cylinder engine. C
smax (ms) E
16.0 1.8
Fig. 9. C against simulated peak pressure.
Fig. 7. Fast, medium and slow cycles obtained by varying C compared with experiment.
using Eq. (5). The detailed exhaust gas composition was stored at the end of each cycle, and used in the following cycle. The EGR mass was assumed to remain constant. Fifty cycles were simulated and the values of peak pressure against the crank angle they occurred are plotted in Fig. 10, along with the experimental results. The simulation matches the experimental data well and the peak pressures do not fall on a straight line, but a range of peak pressures were obtained at the same crank angle. This was because the temperature, pressure and composition at IVC are all affected by the previous cycle, and go on to affect the current cycle. This is a good result, due to the detail required, along with short computation times to make a multi-cycle simulation feasible. The average NOx emissions were 528 ppm and 794 ppm in the experiment and simulation respectively. The simulated NOx emissions reproduced the experimental result well, suggesting the temperature paths taken by the particle ensemble were a good representation of the physical process’ in the experiment. Fig. 11 shows the variation in NOx emissions from cycles with different peak pressures. The faster burning cycles with higher peak pressures produced higher NOx emissions. The emissions increased
Fig. 8. Normal distribution fit to experimental distribution of peak pressures.
pressures measured in the experiment fit a normal distribution well. Using the values from the C sweep, a relationship between peak pressure, in bar, and C was obtained for the four-cylinder engine. Fig. 9 shows C against peak pressure in the cylinder and the curve with Eq. (5).
C ¼ 0:0015P3max 0:054P2max þ 1:83Pmax 10:49:
ð5Þ
The engine model was coupled with GT-Power, a one-dimensional engine simulation tool. This is a two-way coupling as GT-Power provides the SRM with IVC conditions, and the SRM returns the temperatures and pressures after IVC, up to EVO. The normal distribution fitted to the experimental peak pressures shown in Fig. 8, with a mean of 19.7 bar and standard deviation 2.5 bar, was used to generate a random peak pressure at the beginning of each simulated cycle. The peak pressure was used to obtain C for the cycle
Fig. 10. Peak pressures against crank angle at which they occur in 50 simulated and 96 experimental cycles.
186
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
Fig. 11. NOx emissions and peak temperature against the peak pressure during the cycle they occurred in, for 50 simulated cycles.
almost linearly with peak pressure. This was caused by increased peak combustion temperatures, as was expected.
(a)
In Fig. 12a the rates of the thermal NO reactions, in a particle that burns close to CA50, are plotted against crank angle. The zone containing the particle is also shown. The reaction rates are highest when the particle is in the burning entrained zone. Fig. 12b shows the particle’s peak temperature is not reached until it has moved into the burned zone. The high reaction rates in the burning entrained zone can be explained by the higher concentrations of O, N and OH. Formation of NO continues, at a slower rate, when the particle is in the burned zone. The other NO formation reactions occurred at much lower rates compared with the thermal reactions. Fig. 13a shows the evolution of NO concentration in each particle during combustion in the average peak pressure simulation case. The particles were ordered by the CAD they were entrained, with the first to be entrained labelled as particle 1. It can be seen that more NO is produced in particles that are entrained earlier in the cycle, compared with those entrained later in the cycle. This is due to the higher temperatures reached in these particles. The temperature evolution of each particle in the order they were entrained is shown in Fig. 13b. The majority of NO was produced via the thermal mechanism. The reaction rates for this pathway peaked as the particles moved from the entrained and unburned to the entrained and burning sub-division of the entrained zone. It can be seen in Fig. 13a that the level of NO in the particles entrained early in the cycle peaks and then decreases. This effect is mainly due to mixing.
(b)
Fig. 12. Rates of thermal NO reactions and temperature and NO mass fraction in a particle that burns close to CA50 and the zone the particle is in against CAD.
(a)
(b)
Fig. 13. NO mass fraction and temperature evolution during the average SI case in each particle, ordered from 1st particle to be entrained to the last.
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
(a)
187
(b)
Fig. 14. Heat released in the fast and slow cycles and the peak heat release rate in each particle against the CAD at which they occurred.
Fig. 14a compares the heat release in the fast and slow simulated cycles shown in Fig. 7. Although heat release occurred later in the slow cycle, the total heat released was almost the same as in the fast cycle. Fig. 14b compares the peak heat release rate of each particle and the CAD it occurred in the fast and slow cycles. The peak heat release rates reached by particles in the faster burning cycle were higher than the slow cycle. This was because the rate of heat release was lower in the individual particles in the slow cycle, due to the lower temperature of the entrained zone.
4. Conclusions A new PDF based stochastic reactor model for simulating spark ignition engines has been developed. The particle ensemble was split into unburned, entrained and burned zones. Mixing occurred within each zone but not between zones. A detailed mechanism was used to perform chemistry in each of the zones, resulting in an exhaust gas composing of 157 species. Experimental results from a single-cylinder research engine, that showed low CCV, were compared with the simulation. Key parameters in the calibrated model were then investigated. The temperature used at the spark did not affect the simulation results at the conditions studied. The size of the initial flame kernel also had low impact on the simulation over a wide range of values, suggesting the spark initiation sub-model was sufficient. Experiments on a four-cylinder SI engine resulted in cycle to cycle fluctuations. The new SI model was coupled with GT-Power, to enable multi-cycle simulation of the four-cylinder engine. The variation of peak pressures and crank angles they occurred matched experimental data well. The average NOx emissions also matched the experiment, suggesting correct combustion temperatures in the model. The levels of NOx produced in faster burning cycles was much larger than in slower cycles. If CCV could be controlled, then NOx emissions could be reduced at the same speed and load. The amount of NO, and reactions contributing to its production, in cycles that burned at different rates were analysed. It was found the highest rate of NO production occurred in the entrained and burning zone, via the thermal mechanism. The same reactions continued at a slower rate when the particle moved to the burned zone, until the temperature decrease caused the reactions to freeze. Investigation into an individual cycles NO production revealed the majority of NO is formed at the beginning of combustion. After
combustion occurred in the early burning particles, a slight decrease in NO mass fraction was observed. The decrease was mainly due to mixing with particles with a lower NO mass fraction. Acknowledgment Funding from the EPSRC (UK), Grant number EP/D068703/1, is gratefully acknowledged. Experimental data has been provided by the Internal Combustion Engine Group, University of Oxford. References [1] R. Stevens, P. Ewart, H. Ma, C.R. Stone, Combust. Flame 148 (2007) 223–233. [2] L. Dodge, J. Kubesh, D. Naegeli, P. Campbell, Modelling NOx Emissions from Lean-Burn Natural Gas Engines, SAE Paper No. 981389, 1998. [3] H. Safari, S.A. Jazayeri, R. Ebrahimi, Int. J. Hydrogen Energy 34 (2009) 1015– 1025. [4] R. Raine, L. Wyszynski, R. Stone, Proc. Inst. Mech. Eng. 216 (5) (2002) 403–412. [5] S. Narayan, S. Rajan, Int. J. Engine Res. 2 (1) (2001) 34–45. [6] A. Ibrahim, S. Bari, Fuel 87 (2008) 1824–1834. [7] C.D. Rakopoulos, C.N. Michos, Energy Convers. Manage. 49 (2008) 2924–2938. [8] J.K. Ball, C.R. Stone, N. Collings, Proc. Inst. Mech. Eng. 213 (2) (1999) 175–189. [9] E. Abdi Aghdam, A. Burluka, T. Hattrell, K. Liu, C. Sheppard, J. Neumeister, N. Crundwell, Study of Cyclic Variation in an SI Engine Using Quasi-Dimensional Combustion Model, SAE Paper No. 2007-01-0939, 2007. [10] J. Serras-Pereira, P. Aleiferis, D. Richardson, S. Wallace, Mixture Preparation and Combustion Variability in a Spray-Guided DISI Engine, SAE Paper No. 2007-01-4033, 2007. [11] A. Brown, C. Stone, P. Beckwith, Cycle-by-Cycle Variations in Spark Ignition Engine Combustion – Part I Flame Speed and Combustion Measurements and a Simplified Turbulent Combustion Model, SAE Paper No. 960612, 1996. [12] C. Stone, A. Brown, P. Beckwith, Cycle-by-Cycle Variations in Spark Ignition Engine Combustion – Part II Modelling of Flame Kernel Displacements as a Cause of Cycle-by-Cycle Variations, SAE Paper No. 960613, 1996. [13] H. Daneshyar, P.G. Hill, Prog. Energy Combust. Sci. 13 (1987) 47–73. [14] J. Ewald, N. Peters, Proc. Combust. Inst. 31 (2007) 3051–3058. [15] Z. Tan, R.D. Reitz, Combust. Flame 145 (2006) 1–15. [16] H. Bayraktar, Renew. Energy 32 (2007) 758–771. [17] M. Baratta, A. Catania, E. Spessa, A. Vassallo, Development and Assessment of a Multizone Combustion Simulation Code for SI Engines Based on a Novel Fractal Model, SAE Paper No. 2006-01-0048, 2006. [18] J. Keck, Turbulent flame structure and speed in spark-ignition engines, in: Nineteenth Symposium (International) on Combustion/The Combustion Institute, 1982, pp. 1451–1466. [19] M. Kraft, P. Maigaard, F. Mauss, M. Christensen, B. Johansson, Proc. Combust. Inst. 28 (2000) 1195–1201. [20] H. Su, A. Vikhansky, S. Mosbach, M. Kraft, A. Bhave, K. Kim, T. Kobayashi, F. Mauss, Combust. Flame 147 (2006) 118–132. [21] A. Bhave, M. Kraft, L. Montorsi, F. Mauss, Combust. Flame 144 (2006) 634–637. [22] S. Mosbach, A. Aldawood, M. Kraft, Combust. Sci. Technol. 180 (2008) 1263– 1277. [23] J. Etheridge, S. Mosbach, M. Kraft, H. Wu, N. Collings, SAE Int. J. Fuels Lubr. 2 (1) (2009) 13–27. [24] A. Bhave, M. Balthasar, M. Kraft, F. Mauss, Int. J. Engine Res. 5 (1) (2004) 93– 104.
188
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
[25] S. Subramaniam, S.B. Pope, Combust. Flame 115 (1998) 487–514. [26] R. Curl, AIChE J 9 (2) (1963) 175181. [27] A. Gogan, B. Sundén, A. Lehtiniemi, F. Mauss, Stochastic Model for the Investigation of the Influence of Turbulent Mixing on Engine Knock, SAE Paper No. 2004-01-2999, 2004. [28] J.B. Heywood, Internal Combustion Engine Fundamentals, McGraw-Hill, New York, 1988.
[29] J. Tagalian, J. Heywood, Combust. Flame 64 (1986) 243–246. [30] S. Karagiorgis, N. Collings, K. Glover, N. Coghlan, A. Petridis, Residual Gas Fraction and Estimation on a Homogeneous Charge Compression Ignition Engine Utilizing the Negative Valve Overlap Strategy, SAE Paper No. 2006-013276, 2006. [31] T. Kravchik, E. Sher, Combust. Flame 99 (1994) 635–643.
Contents lists available at ScienceDirect
Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e
Modelling cycle to cycle variations in an SI engine with detailed chemical kinetics Jonathan Etheridge a, Sebastian Mosbach a, Markus Kraft a,⇑, Hao Wu b, Nick Collings b a b
Department of Chemical Engineering and Biotechnology, University of Cambridge, New Museums Site Pembroke Street, Cambridge, CB2 3RA, UK Department of Engineering, University of Cambridge, Trumpington Street, Cambridge, CB2 1PZ, UK
a r t i c l e
i n f o
Article history: Received 8 January 2010 Received in revised form 26 May 2010 Accepted 16 August 2010 Available online 24 September 2010 Keywords: SI engine modelling Detailed chemistry Cycle to cycle variation
a b s t r a c t This paper presents experimental results and a new computational model that investigate cycle to cycle variations (CCV) in a spark ignition (SI) engine. An established stochastic reactor model (SRM) previously used to examine homogeneous charge compression ignition (HCCI) combustion has been extended by spark initiation, flame propagation and flame termination sub-models in order to simulate combustion in SI engines. The model contains a detailed chemical mechanism but relatively short computation times are achieved. The flame front is assumed to be spherical and centred at the spark location, and a pent roof and piston bowl geometry are accounted for. The model is validated by simulating the pressure profile and emissions from an iso-octane fuelled single cylinder research engine that showed low CCV. The effects of key parameters are investigated. Experimental results that show cycle to cycle fluctuations in a four-cylinder naturally aspirated gasoline fuelled SI engine are presented. The model is then coupled with GT-Power, a one-dimensional engine simulation tool, which is used to simulate the breathing events during a multi-cycle simulation. This allows an investigation of the cyclic fluctuations in peak pressure. The source and magnitude of nitric oxide (NO) emissions produced by different cycles are then investigated. It was found that faster burning cycles result in increased NO emissions compared with cycles that have a slower rate of combustion and that more is produced in the early stages of combustion compared with later in the cycle. The majority of NO was produced via the thermal mechanism just after combustion begins. Ó 2010 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
1. Introduction Nitric oxides (NOx), unburned hydrocarbons (UHC) and carbon monoxide (CO) are major components of air pollution caused by Internal Combustion Engines (ICE). Increasingly stringent restrictions on vehicle emissions have led to much research into reducing these emissions. Spark ignition (SI) engine NOx emissions have been studied both experimentally [1] and with models of varying complexity [2–5]. The high temperatures reached in SI engines cause the conversion of atmospheric nitrogen to nitric oxide (NO). The three accepted mechanisms responsible for NO formation are the thermal, prompt and nitrous oxide mechanisms. The majority of NO is formed in the burned gases behind the flame front, where the temperatures remain high for enough time to produce NO by the thermal mechanism that has relatively slow rate constants. The prompt mechanism produces NO in the flame, although the amount is thought to be much lower than in the end gas.
⇑ Corresponding author. E-mail address: [email protected] (M. Kraft). URL: http://como.cheng.cam.ac.uk (M. Kraft).
Empirical equations exist to predict NO emissions from SI engines, however these are extremely simplified. SI engine NO emissions have been modelled previously by assuming that equilibrium is reached instantaneously behind the flame front for all species, with the exception of NO due to its slow rate of formation. Chemical kinetics are then used in the hot end gas to calculate NO production [6]. Improved predictions have been made by splitting the end gas into multiple zones that have burnt at different times, taking into account the difference in temperature of the charge that burns earlier and later in the engine cycle [4,7]. As well as different points in a cycle being responsible for varying amounts of NO, fluctuations between cycles will result in cycle to cycle variations in NO emissions [8]. Cycle to cycle variations (CCV) in SI engine combustion are a common problem which reduce performance due to variations in combustion phasing, and, in some circumstances, can lead to knocking or misfire. Several causes of CCV are known and the magnitude of their effects varies between engine design and operating conditions. In-cylinder charge motion can vary between cycles, altering the rate at which unburned mixture is entrained and the shape of the turbulent flame front. The total mass of fuel, air and exhaust gas recirculation (EGR) can vary between cycles, as well as local compositions in the cylinder, which are especially
0010-2180/$ - see front matter Ó 2010 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.combustflame.2010.08.006
180
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
List of symbols and abbreviations Symbols C CADspark lT NB NE NU NTot Pmax Rf;n SL t tspark tstop Dt i u uT VB VE Vf V Unext
qi qu sb sE
smax E entrainment rate constant (–) spark timing angle (CAD ATDC) characteristic length scale (m) number of burned particles (–) number of entrained particles (–) number of unburned particles (–) total number of particles at IVC (–) peak pressure (bar) flame radius at nth CAD step (m) laminar flame speed (m s1) time after IVC (s) spark time (s) time at EVO (s) time step (s) mean inlet gas speed (m s1) characteristic flame speed (m s1) burned zone volume (m3) entrained zone volume (m3) flame volume calculated from the radius (m3) volume of the next particle to be entrained (m3) inlet air density (kg m3) unburned zone density (kg m3) characteristic burn time (s) entrained zone mixing time (s)
important near the spark location at the time of ignition. Variations in local turbulence at the time of spark can lead to a movement of the flame kernel altering the interactions between the flame and the piston and cylinder walls. Many studies of these effects have been completed. In [9], variations in the peak pressures of two different engines were simulated by random fluctuations, with a normal distribution, in the rms turbulent velocity. Fluctuations in injected fuel mass were also given as a possible cause of the measured CCV. Experimental results have also suggested that CCV is not very dependent on spray fluctuations, and the in-cylinder flow just after the spark timing has a greater effect [10]. Convection of the flame kernel has also been investigated and proposed to be a key factor [11,12]. Movement of the flame kernel and the time delay before the flame becomes fully turbulent have also been discussed in [13] as effects of variations in the rms turbulent velocity. Combustion in SI engines takes place in a turbulent flow field. The structure and speed of the flame depend on charge motion, composition, temperature and chamber geometry. Detailed SI models exist that predict the structure and propagation of the turbulent flame [14,15]. To run many cycles of a computational fluid dynamics (CFD) model with detailed chemistry however would be extremely computationally expensive. Less complex models exist that assume the flame is spherical and propagates at the turbulent flame speed from the spark [16–18]. In these models the cylinder is divided into unburned and burned zones and the burned zone may be subdivided into several zones that have burned at different times. The burned zone composition is frequently calculated by assuming the mixture behind the flame front reaches equilibrium. In this paper we present a new probability density function (PDF) based SI model that uses a detailed chemical mechanism throughout the cycle. The model has been created from a previously developed stochastic reactor model (SRM) successfully used to model Homogeneous Charge Compression Ignition (HCCI) combustion. The model has been validated by simulating a single-cylinder research engine which demonstrated low cycle to cycle variation. The model was then coupled with GT-Power, a
maximum entrained zone mixing time (s)
Abbreviations BTDC before top dead centre CA50 crank angle at 50% heat release CAD crank angle degree CCV cycle to cycle variation CFD computational fluid dynamics CR compression ratio EGR exhaust gas recirculation EMST Euclidean minimum spanning tree EVO exhaust valve open FID flame ionisation detector GDI gasoline direct injection HCCI homogeneous charge compression ignition IVC inlet valve close NDIR non-dispersive infrared PDF probability density function PRF primary reference fuel RON research octane number RPM revolutions per minute SI spark ignition SRM stochastic reactor model UHC unburned hydrocarbons
one-dimensional engine simulation tool, which was used to simulate the breathing events during multi-cycle simulation. Experimental results from a four-cylinder naturally aspirated gasoline fuelled SI engine were modelled and the cyclic fluctuations in pressure and NO emissions investigated. The paper is structured as follows. In the next section the model details and implementation is reported. A description of both engine set-ups is given at the start of Section 3. Calibration of the model using experimental data from a single-cylinder research engine is then presented. Subsequently the model is calibrated using experimental data from fast, slow and average cycles from a fourcylinder SI engine. Results from a multi-cycle simulation used to model cyclic fluctuations are contained in the same section along with experimental measurements of emissions. Finally conclusions are drawn and future work outlined. 2. Model details The model was built into the framework of the SRM used to simulate HCCI combustion [19–23]. The SRM is based on the PDF transport equation, and models convective heat transfer [24], turbulent mixing, piston movement and detailed chemistry. The PDF is approximated by dividing the charge into a notional ensemble of NTot stochastic particles. For SI simulation the particle ensemble is divided into three zones: unburned, entrained and burned, which contain NU, NE and NB particles respectively. Mixing occurs within each zone, but not between zones and all particles are initially in the unburned zone. The Euclidean minimum spanning tree (EMST) mixing model was used in the unburned and burned zones and is described in [25]. A Curl based mixing model was used in the entrained zone and will be described in more detail later [26]. The entrained zone is subdivided into unburned and burning particle categories to keep track of the particles’ states, as will be discussed later. Chemistry occurs at every CAD step in every particle, as in the HCCI model. The primary reference fuel (PRF) mechanism used
181
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
contains 157 species and 1552 reactions [19–23]. The model could be used to investigate knocking as in [27], however, similar work is not presented here. The main SRM operator splitting loop, with the addition of spark initiation, flame propagation and flame termination, was implemented as follows: 1. Initialise t = 0, Dt, CAD = IVC, CADspark, spark energy, Rf = 0, NU = NTot, NE = NB = 0 all particles are in the unburned zone. Determine temperature, composition, mass, volume and pressure of particle ensemble. 2. Progress in time t ´ t + Dt. If CAD P EVO or t P tstop then save the detailed exhaust composition as input EGR and stop. 3. Perform volume change due to piston movement. 4. Perform the mixing step. 5. Perform stochastic heat transfer. 6. Perform the pressure equilibration step [27]. 7. If CAD P CADspark and NU = NTot perform spark initiation (as described below) else if CAD P CADspark and NU > 0 perform flame propagation (as described below) else go to (8). 8. Perform the chemistry step. If NE > 0 perform flame termination (as described below). 9. Perform the pressure equilibration step [27]. 10. Perform the mixing step. 11. Go to step (2). 2.1. Spark initiation During flame propagation, particles are moved from the unburned zone to the entrained zone, according to the required volume of the flame. This is a discrete process, and the volume of the entrained and burned particles is unlikely to match the calculated flame volume exactly. During spark initiation, the initial flame volume is very small. If 100 particles are used, the initial entrained volume is much lower than the volume represented by a single particle. If the energy in the spark were to be distributed uniformly throughout a volume represented by a particle, then the temperature would not be sufficient to initiate combustion. For this reason, the first particle randomly chosen to enter the entrained zone is not moved directly into the zone, but the required entrained volume is removed from the particle and added to a newly created particle in the entrained zone. Once the calculated entrained volume has increased above the particle’s volume, it is reunited into a single particle in the entrained zone. The diameter of the entrained zone at the time of spark was set to 1.0 mm, the rough size of the spark gap. The spark is simulated by setting the temperature of the mass added to the entrained zone, at each time step, to 2000 K, until all of the spark energy has been added. This low temperature was chosen to prevent any problems with the chemical mechanism, which was not developed to represent the high temperatures experienced at the spark. The spark initiation was modelled as follows: 1. If its the first call to the spark initiation subroutine, continue to (2) else go to (7). 2. Set the flame radius to Rf = 0.5 mm. 3. Select a particle at random with a uniform distribution. 4. Remove the spherical volume with 0.5 mm radius from the particle and add to a newly created particle in the entrained zone. NE = 1, V E ¼ 4=3pR3f . 5. Increase the temperature of the entrained particle to 2000 K by adding energy at constant pressure. Remove the energy added from the initial spark energy. 6. Go to the chemistry step. 7. Calculate the new flame volume, Vf, as described in the flame propagation section.
8. If V f > V UNext þ V E go to (13) else go to (9). 9. Remove Vf VE from the unburned particle being entrained. 10. If there is spark energy remaining, increase the temperature of the freshly entrained mass to a maximum of 2000 K by adding energy at constant pressure. 11. Reunite the freshly entrained mass with the entrained particle. 12. Go to the chemistry step. 13. Remove the unburned particle being entrained from the unburned zone, NU = NTot 1. 14. If there is spark energy remaining, increase the temperature of the freshly entrained mass to a maximum of 2000 K by adding energy at constant pressure. 15. Reunite the freshly entrained mass with the entrained particle. 16. Select the next particle to entrain, at random with a uniform distribution. Go to step (2) in flame propagation section described below. 2.2. Flame propagation A randomly chosen particle is moved from the unburned zone to the entrained zone, when the calculated flame volume exceeds the burned and entrained zone volumes by more than half the volume of the next particle to be entrained. The volume of the flame is calculated from the geometry of the piston and cylinder, and the flame radius, assuming it is spherical and centred at the spark location. The flame radius at the nth time step, Rf;n, is obtained from
ðt t spark Þ Rf ;n ¼ Rf ;n1 þ uT 1 exp þ SL Dt:
sb
ð1Þ
The increase in flame radius is calculated as the time step multiplied by the entrainment velocity. The flame radius is also updated after the pressure equilibration step, to take into account the expansion of the flame. The radius initially increases at the laminar flame speed, SL, taken from [28], until a wrinkled turbulent flame front develops after a characteristic burn time, sb [29,28]. The characteristic burn time is given by
sb ¼
lT ; SL
ð2Þ
where lT is the characteristic length scale [28]. One of the most desirable features of the EMST mixing model is that it preferentially mixes particles which are close in composition space. This is based on the reasoning that fluid parcels which are close in physical space should be similar in terms of composition. The flame brush is relatively thin with high temperature and concentration gradients across it, so EMST is not a good mixing model to use in the entrained zone. A mixing model based on Curl’s has been used in the entrained zone [26], whereby two particles are selected to mix. The difference being that only unburned entrained particles are chosen to mix with entrained and burning particles. The resulting particle compositions are not identical, as the particles are only moved to their mean composition by a factor of 10%, and the mixing time is accordingly adjusted to reflect this. This method was chosen to prevent the entrained zone from becoming unphysically homogeneous. The mixing time in the entrained zone, sE, is increased at the same exponential rate as the transition from a laminar to a turbulent flame speed. The mixing time increases to the value smax once E the turbulent flame has fully developed. Before this, the entrained zone mixing time is smaller causing faster mixing and heat release. This was included as the rate of burning is initially the rate of entrainment before a turbulent flame front has developed containing unburned pockets [13].
182
sE
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
ðt t spark Þ : ¼ smax 1 exp E
sb
ð3Þ
The rate of burning is controlled by the mass entrained, and the rate it mixes with the hot, reacting entrained material, and is calculated with the detailed chemical mechanism. The maximum entrained zone mixing time was smaller than the other zones, to prevent the volume of the entrained zone growing unrealistically large, and to propagate heat through the entrained zone at the correct rate. Flame propagation was implemented in the model as follows: 1. Increase the flame radius using Eq. (1) and calculate the new flame volume, Vf. 2. Update the entrained zone mixing time according to Eq. (3). P E E PNB B 3. If V f > Ni¼1 V i þ i¼1 V i þ 1=2V UNext go to (4) else go to (6). 4. Add the particle to the entrained zone and select the next particle to be added randomly with a uniform distribution. NU = NU 1, NE = NE + 1. Flag the newly entrained particle as unburned and entrained. 5. Go to step (3). 6. Perform the chemistry step. The correlation used to obtain the characteristic flame speed, uT, was taken from [18]:
i uT ¼ 0:08C u
qu qi
1=2 ;
ð4Þ
i is the mean inlet gas speed, qu is the unburned gas density where u and qi is the inlet air density. The correlation used is relatively old, however this work was focused on producing a detailed chemistry SI model incorporated into the SRM, to enable multi-cycle simulation at low computational cost. There were no measurements of in-cylinder turbulence, hence the constant C was used to fit the model to the experimental data by varying the rate of entrainment. The model was intended to predict not only the pressure profile, but also temperature and composition. The cylinder geometry model contains a pent roof and the piston contains a bowl. The volume of the spherical flame is calculated exactly in the bowl and cylinder areas. In the pent roof region, numerical integration of the flame cross sectional area at different heights is used to obtain the flame volume. 2.3. Flame termination In this section we describe how particles are moved from the entrained zone to the burned zone. Particles move from the entrained and unburned zone to the entrained and burning zone once their temperature has risen above 1500 K. During the chemistry step the rate of heat release in each entrained particle is calculated. When a particle’s heat release rate has dropped below 2% of its maximum value, it is moved to the burned zone. When no particles remain in the unburned zone, the entrained zone mixing time is set to smax as Eq. (2) can no longer be evaluated. Flame termination E was modelled as follows: 1. If any unburned entrained particle’s temperature has risen above 1500 K, flag the particle as entrained and burning. 2. If any entrained and burning particle’s rate of heat release has dropped below 2% of its maximum, move the particle to the burned zone. NE = NE 1, NB = NB + 1. 3. If NU = 0 then sE ¼ smax E . In summary, an SI model has been integrated into the SRM. The particle ensemble is divided into three zones: unburned, entrained and burned. Mixing occurs within each zone, but not between zones. Chemistry occurs at every CAD step in every particle, as in
the HCCI model. An initial temperature and flame radius of 2000 K and 0.5 mm are set at the spark timing respectively. The flame radius is updated using an entrainment velocity. The flame radius and cylinder and piston geometry are used to obtain the flame volume. Particles are moved from the unburned zone to the entrained zone when the calculated flame volume exceeds the burned and entrained zone volumes by more than half the volume of the next particle to be entrained. Entrained particles, initially in the entrained and unburned zone, are moved to the entrained and burning zone once their temperature has risen above 1500 K. Entrained and burning particles are moved to the burned zone when their heat release rate has dropped below 2% of its maximum value. 3. Results and discussion In this section we present and compare experimental and simulation results from two different engines: a single-cylinder naturally aspirated direct injection iso-octane fuelled research engine and a four-cylinder naturally aspirated GDI engine. The engine specifications and operating conditions are given in Table 1. Exhaust emissions of CO and UHC, as well as the in-cylinder pressure, were measured in the single-cylinder research engine. The CO volume fraction was measured with a non-dispersive infrared (NDIR) sensor. UHC were measured 200 mm downstream of the exhaust valve with a fast flame ionisation detector (FID). The pressure in each cylinder of the four-cylinder GDI SI engine was monitored with Kistler 6123 piezoelectric transducers and recorded. Emissions of NOx were measured using a Horiba EXSA1500 analyser. More details of the engine can be found in [30]. 3.1. Model calibration The model was initially calibrated using results from the single cylinder research engine. As fuel injection occurred early during intake, the cylinder charge was assumed to be initially homogeneous in terms of composition and temperature, hence all particles were set to be identical at Inlet Valve Closing (IVC). The model was calibrated to find values for C and smax E . The values of these parameters are given in Table 2. The pressure profile was matched to the experimental data shown in Fig. 1. The emissions matched experimental data reasonably well, as shown in Table 3. No NOx measurements were taken. Fig. 2 shows the particles remained in the entrained zone for between 0 and 1.75 ms. The PDF was obtained by taking each value from the ensemble as the mean of a normal distribution, with a r of 0.05 ms, summing the distributions and dividing by the total number. From experimental observations, the burn time in the turbulent flame front has been estimated at values ranging from 1.0 to 1.8 ms [13]. The PDF becomes smoother when the particle number
Table 1 SI mode engine specifications and operating conditions. Engine
Research
Production
Cylinders Fuel Bore (mm) Stroke (mm) Con. rod length (mm) Disp. volume (cm3/cyl) CR Speed (rpm) Air/fuel equiv. ratio EGR (%) Spark timing (CAD BTDC) Start of injection (CAD BTDC)
1 iso-octane 89.0 90.3 148.9 562 11 1500 1.0 15.0 35 280
4 95 RON gasoline 87.5 83.0 146.3 499 12 1500 1.0 28.8 40 308
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188 Table 2 Parameters used in simulation of singlecylinder research engine. No. of particles
100
CAD step (CAD) C smax (ms) E
0.2 9.0 1.7
Fig. 1. Comparison of experimental and simulated pressure profiles for the singlecylinder research engine.
Table 3 Emission results for single-cylinder research engine.
UHC (ppm) CO (ppm) NOx (ppm)
Exp.
Sim.
600 2700 –
746 3185 671
183
is increased. The PDF of entrained time is skewed. The shape is caused by the time particles are entrained for, before they begin to mix. This is affected by the entrained zone mixing time, which increases during the simulation, so particles entrained later in the cycle are likely to spend longer in the entrained zone. Fig. 3 shows the average rate of heat release against CAD, calculated by centering each particles’ peak heat release rate at 0CAD and averaging over the 100 particles. The rate of heat release increased rapidly, and the majority of the chemical energy was released within 2CAD. The number of particles entrained and burning remains roughly constant, however, the number of particles entrained and unburned grows and then falls off (see Fig. 4a). This indicates a thickening flame brush, followed by consumption of the remaining fuel as the flame reaches the walls. The temperatures remain roughly constant throughout (see Fig. 4b). The unburned entrained region is slightly warmer than the unburned zone, due to some particles having mixed with the burning entrained particles. The burning entrained particles are slightly cooler than the burned particles, as all the chemical energy has not been released in this zone, and particles can mix with the cooler unburned entrained particles. A spark temperature of 2000 K was chosen, as the mechanism was not created for the extremely high combustion temperatures that are reached within a spark plug gap during the passage of a spark. Fig. 5a suggests this is an acceptable approach, as there does not appear to be a trend between the spark temperature and the pressure profile at the conditions investigated. Fig. 5b shows that increasing the initial flame radius caused a slight advance in combustion phasing. The effect is not significant, given the relatively large range investigated. As variations in this parameter caused minor changes in the pressure profile, the initial flame radius assumption was taken as acceptable and a detailed spark initiation sub-model, e.g. [31], was not implemented. The simulation was repeated with n-heptane as the fuel. This was done without changing the laminar flame speed constants, in order to investigate the sensitivity of the combustion to the entrained zone composition. The initial burning rate was faster, and also led to knocking, as shown in Fig. 6. 3.2. Multi-cycle simulation with CCV A four-cylinder SI engine displayed high CCV at the speed and load investigated. Measurements suggested an EGR ratio of 30%
Fig. 2. Probability density of particles’ time in the entrained zone.
Fig. 3. Average rate of heat release from each particle against CAD, normalised so that all peak rates were at 0 CAD.
184
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
(a)
(b)
Fig. 4. Number of particles and the average temperature of each zone against CAD.
(a)
(b)
Fig. 5. Comparison of pressure profiles using different spark temperatures and initial flame radii.
Fig. 6. Comparison of in-cylinder pressure profiles with iso-octane and n-heptane fuel.
at this operating point. The experimental method uses the ratio of CO2 in the exhaust and compression stroke, and is described in [30]. Large amounts of internal EGR have been reported to increase CCV [28]. Details of the engine are given in Table 1. A mixture of 95% iso-octane and 5% n-heptane was used to simulate the 95 Research Octane Number (RON) gasoline. The model was calibrated to fit an experimental pressure profile from a cycle with a peak pressure close to the average as shown in Fig. 7. The model parameters used are given in Table 4. Fig. 7 shows experimental pressure against crank angle, and the simulation results when C was varied and smax was constant. The E pressure traces match the experimental data well, with C values of 11 and 22 used to match the slow and fast pressure traces respectively. The cycles in which the flame initially entrains mass at a faster rate, continue to entrain mass at a fast rate throughout the cycle. This suggests variations in the flame surface area occur early in the cycle resulting, in different rates of unburned gas entrainment and combustion. This could be caused by the flame kernel being convected towards a wall, reducing the available area, or due to turbulence altering the shape of the kernel early in the cycle, resulting in a higher surface area. Fig. 8 shows the PDF of peak pressure obtained from the experimental data, and a normal distribution fitted to the data. The peak
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
185
Table 4 Parameters used in simulation of four-cylinder engine. C
smax (ms) E
16.0 1.8
Fig. 9. C against simulated peak pressure.
Fig. 7. Fast, medium and slow cycles obtained by varying C compared with experiment.
using Eq. (5). The detailed exhaust gas composition was stored at the end of each cycle, and used in the following cycle. The EGR mass was assumed to remain constant. Fifty cycles were simulated and the values of peak pressure against the crank angle they occurred are plotted in Fig. 10, along with the experimental results. The simulation matches the experimental data well and the peak pressures do not fall on a straight line, but a range of peak pressures were obtained at the same crank angle. This was because the temperature, pressure and composition at IVC are all affected by the previous cycle, and go on to affect the current cycle. This is a good result, due to the detail required, along with short computation times to make a multi-cycle simulation feasible. The average NOx emissions were 528 ppm and 794 ppm in the experiment and simulation respectively. The simulated NOx emissions reproduced the experimental result well, suggesting the temperature paths taken by the particle ensemble were a good representation of the physical process’ in the experiment. Fig. 11 shows the variation in NOx emissions from cycles with different peak pressures. The faster burning cycles with higher peak pressures produced higher NOx emissions. The emissions increased
Fig. 8. Normal distribution fit to experimental distribution of peak pressures.
pressures measured in the experiment fit a normal distribution well. Using the values from the C sweep, a relationship between peak pressure, in bar, and C was obtained for the four-cylinder engine. Fig. 9 shows C against peak pressure in the cylinder and the curve with Eq. (5).
C ¼ 0:0015P3max 0:054P2max þ 1:83Pmax 10:49:
ð5Þ
The engine model was coupled with GT-Power, a one-dimensional engine simulation tool. This is a two-way coupling as GT-Power provides the SRM with IVC conditions, and the SRM returns the temperatures and pressures after IVC, up to EVO. The normal distribution fitted to the experimental peak pressures shown in Fig. 8, with a mean of 19.7 bar and standard deviation 2.5 bar, was used to generate a random peak pressure at the beginning of each simulated cycle. The peak pressure was used to obtain C for the cycle
Fig. 10. Peak pressures against crank angle at which they occur in 50 simulated and 96 experimental cycles.
186
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
Fig. 11. NOx emissions and peak temperature against the peak pressure during the cycle they occurred in, for 50 simulated cycles.
almost linearly with peak pressure. This was caused by increased peak combustion temperatures, as was expected.
(a)
In Fig. 12a the rates of the thermal NO reactions, in a particle that burns close to CA50, are plotted against crank angle. The zone containing the particle is also shown. The reaction rates are highest when the particle is in the burning entrained zone. Fig. 12b shows the particle’s peak temperature is not reached until it has moved into the burned zone. The high reaction rates in the burning entrained zone can be explained by the higher concentrations of O, N and OH. Formation of NO continues, at a slower rate, when the particle is in the burned zone. The other NO formation reactions occurred at much lower rates compared with the thermal reactions. Fig. 13a shows the evolution of NO concentration in each particle during combustion in the average peak pressure simulation case. The particles were ordered by the CAD they were entrained, with the first to be entrained labelled as particle 1. It can be seen that more NO is produced in particles that are entrained earlier in the cycle, compared with those entrained later in the cycle. This is due to the higher temperatures reached in these particles. The temperature evolution of each particle in the order they were entrained is shown in Fig. 13b. The majority of NO was produced via the thermal mechanism. The reaction rates for this pathway peaked as the particles moved from the entrained and unburned to the entrained and burning sub-division of the entrained zone. It can be seen in Fig. 13a that the level of NO in the particles entrained early in the cycle peaks and then decreases. This effect is mainly due to mixing.
(b)
Fig. 12. Rates of thermal NO reactions and temperature and NO mass fraction in a particle that burns close to CA50 and the zone the particle is in against CAD.
(a)
(b)
Fig. 13. NO mass fraction and temperature evolution during the average SI case in each particle, ordered from 1st particle to be entrained to the last.
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
(a)
187
(b)
Fig. 14. Heat released in the fast and slow cycles and the peak heat release rate in each particle against the CAD at which they occurred.
Fig. 14a compares the heat release in the fast and slow simulated cycles shown in Fig. 7. Although heat release occurred later in the slow cycle, the total heat released was almost the same as in the fast cycle. Fig. 14b compares the peak heat release rate of each particle and the CAD it occurred in the fast and slow cycles. The peak heat release rates reached by particles in the faster burning cycle were higher than the slow cycle. This was because the rate of heat release was lower in the individual particles in the slow cycle, due to the lower temperature of the entrained zone.
4. Conclusions A new PDF based stochastic reactor model for simulating spark ignition engines has been developed. The particle ensemble was split into unburned, entrained and burned zones. Mixing occurred within each zone but not between zones. A detailed mechanism was used to perform chemistry in each of the zones, resulting in an exhaust gas composing of 157 species. Experimental results from a single-cylinder research engine, that showed low CCV, were compared with the simulation. Key parameters in the calibrated model were then investigated. The temperature used at the spark did not affect the simulation results at the conditions studied. The size of the initial flame kernel also had low impact on the simulation over a wide range of values, suggesting the spark initiation sub-model was sufficient. Experiments on a four-cylinder SI engine resulted in cycle to cycle fluctuations. The new SI model was coupled with GT-Power, to enable multi-cycle simulation of the four-cylinder engine. The variation of peak pressures and crank angles they occurred matched experimental data well. The average NOx emissions also matched the experiment, suggesting correct combustion temperatures in the model. The levels of NOx produced in faster burning cycles was much larger than in slower cycles. If CCV could be controlled, then NOx emissions could be reduced at the same speed and load. The amount of NO, and reactions contributing to its production, in cycles that burned at different rates were analysed. It was found the highest rate of NO production occurred in the entrained and burning zone, via the thermal mechanism. The same reactions continued at a slower rate when the particle moved to the burned zone, until the temperature decrease caused the reactions to freeze. Investigation into an individual cycles NO production revealed the majority of NO is formed at the beginning of combustion. After
combustion occurred in the early burning particles, a slight decrease in NO mass fraction was observed. The decrease was mainly due to mixing with particles with a lower NO mass fraction. Acknowledgment Funding from the EPSRC (UK), Grant number EP/D068703/1, is gratefully acknowledged. Experimental data has been provided by the Internal Combustion Engine Group, University of Oxford. References [1] R. Stevens, P. Ewart, H. Ma, C.R. Stone, Combust. Flame 148 (2007) 223–233. [2] L. Dodge, J. Kubesh, D. Naegeli, P. Campbell, Modelling NOx Emissions from Lean-Burn Natural Gas Engines, SAE Paper No. 981389, 1998. [3] H. Safari, S.A. Jazayeri, R. Ebrahimi, Int. J. Hydrogen Energy 34 (2009) 1015– 1025. [4] R. Raine, L. Wyszynski, R. Stone, Proc. Inst. Mech. Eng. 216 (5) (2002) 403–412. [5] S. Narayan, S. Rajan, Int. J. Engine Res. 2 (1) (2001) 34–45. [6] A. Ibrahim, S. Bari, Fuel 87 (2008) 1824–1834. [7] C.D. Rakopoulos, C.N. Michos, Energy Convers. Manage. 49 (2008) 2924–2938. [8] J.K. Ball, C.R. Stone, N. Collings, Proc. Inst. Mech. Eng. 213 (2) (1999) 175–189. [9] E. Abdi Aghdam, A. Burluka, T. Hattrell, K. Liu, C. Sheppard, J. Neumeister, N. Crundwell, Study of Cyclic Variation in an SI Engine Using Quasi-Dimensional Combustion Model, SAE Paper No. 2007-01-0939, 2007. [10] J. Serras-Pereira, P. Aleiferis, D. Richardson, S. Wallace, Mixture Preparation and Combustion Variability in a Spray-Guided DISI Engine, SAE Paper No. 2007-01-4033, 2007. [11] A. Brown, C. Stone, P. Beckwith, Cycle-by-Cycle Variations in Spark Ignition Engine Combustion – Part I Flame Speed and Combustion Measurements and a Simplified Turbulent Combustion Model, SAE Paper No. 960612, 1996. [12] C. Stone, A. Brown, P. Beckwith, Cycle-by-Cycle Variations in Spark Ignition Engine Combustion – Part II Modelling of Flame Kernel Displacements as a Cause of Cycle-by-Cycle Variations, SAE Paper No. 960613, 1996. [13] H. Daneshyar, P.G. Hill, Prog. Energy Combust. Sci. 13 (1987) 47–73. [14] J. Ewald, N. Peters, Proc. Combust. Inst. 31 (2007) 3051–3058. [15] Z. Tan, R.D. Reitz, Combust. Flame 145 (2006) 1–15. [16] H. Bayraktar, Renew. Energy 32 (2007) 758–771. [17] M. Baratta, A. Catania, E. Spessa, A. Vassallo, Development and Assessment of a Multizone Combustion Simulation Code for SI Engines Based on a Novel Fractal Model, SAE Paper No. 2006-01-0048, 2006. [18] J. Keck, Turbulent flame structure and speed in spark-ignition engines, in: Nineteenth Symposium (International) on Combustion/The Combustion Institute, 1982, pp. 1451–1466. [19] M. Kraft, P. Maigaard, F. Mauss, M. Christensen, B. Johansson, Proc. Combust. Inst. 28 (2000) 1195–1201. [20] H. Su, A. Vikhansky, S. Mosbach, M. Kraft, A. Bhave, K. Kim, T. Kobayashi, F. Mauss, Combust. Flame 147 (2006) 118–132. [21] A. Bhave, M. Kraft, L. Montorsi, F. Mauss, Combust. Flame 144 (2006) 634–637. [22] S. Mosbach, A. Aldawood, M. Kraft, Combust. Sci. Technol. 180 (2008) 1263– 1277. [23] J. Etheridge, S. Mosbach, M. Kraft, H. Wu, N. Collings, SAE Int. J. Fuels Lubr. 2 (1) (2009) 13–27. [24] A. Bhave, M. Balthasar, M. Kraft, F. Mauss, Int. J. Engine Res. 5 (1) (2004) 93– 104.
188
J. Etheridge et al. / Combustion and Flame 158 (2011) 179–188
[25] S. Subramaniam, S.B. Pope, Combust. Flame 115 (1998) 487–514. [26] R. Curl, AIChE J 9 (2) (1963) 175181. [27] A. Gogan, B. Sundén, A. Lehtiniemi, F. Mauss, Stochastic Model for the Investigation of the Influence of Turbulent Mixing on Engine Knock, SAE Paper No. 2004-01-2999, 2004. [28] J.B. Heywood, Internal Combustion Engine Fundamentals, McGraw-Hill, New York, 1988.
[29] J. Tagalian, J. Heywood, Combust. Flame 64 (1986) 243–246. [30] S. Karagiorgis, N. Collings, K. Glover, N. Coghlan, A. Petridis, Residual Gas Fraction and Estimation on a Homogeneous Charge Compression Ignition Engine Utilizing the Negative Valve Overlap Strategy, SAE Paper No. 2006-013276, 2006. [31] T. Kravchik, E. Sher, Combust. Flame 99 (1994) 635–643.