Combustion-response mapping procedure for internal-combustion engine emissions

Applied Energy 156 (2015) 149–158

Contents lists available at ScienceDirect

Applied Energy journal homepage: www.elsevier.com/locate/apenergy

Combustion-response mapping procedure for internal-combustion engine emissions T. Korakianitis a,⇑, S. Imran b,c, N. Chung d, Hassan Ali c, D.R. Emberson c, R.J. Crookes c a

Parks College of Engineering, Aviation and Technology, Saint Louis University, St. Louis, MO 63103, USA Faculty of Mechanical Engineering, University of Engineering and Technology, Lahore, Pakistan c School of Engineering and Materials Science, Queen Mary University of London, Mile End Road, E1 4NS, UK d New Green Way Pte. Ltd., 51 Goldhill Plaza, #7-10/11, Singapore 308900, Singapore b

h i g h l i g h t s  Prediction of emissions by using CFD and detailed chemical kinetics reactions.  Decoupling of CFD calculations during the combustion process.  Pre-calculations of chemical kinetics for different thermodynamic conditions.  Experimental validation of the technique for a CI and an SI engine.

a r t i c l e

i n f o

Article history: Received 20 November 2014 Received in revised form 12 May 2015 Accepted 17 June 2015

Keywords: Combustion response mapping Decoupling chemical kinetics Combustion CFD

a b s t r a c t This paper describes a new method to predict emissions in internal combustion (IC) engines. The method couples a multi-dimensional engine modeling program with pre-integrated non-equilibrium chemical kinetics reaction results. Prior to engine simulation, detailed chemical kinetics reactions of air/fuel mixture at different temperatures, pressures, and compositions, are calculated using SENKIN, a subprogram in the CHEMKIN-II computer package. The reaction results are decoupled from their chemical eigenvalue (order of about 1010 s), then integrated and saved in physical time scale (order of about 105 s) in a database file. In the database reaction results of different initial conditions (temperature, pressure, and species composition) are stored in different zones; the zones are indexed using their respective reaction conditions. Fluid dynamics and thermal dynamics of the movement of piston and valves, and spray droplets interaction are simulated by KIVA-3 V. Instead of calculating directly the non-equilibrium chemical reactions of the air/fuel mixture, reaction results are obtained from the database file via an interpolating subroutine, which returns temperature, heat release, and species concentrations after reaction to the main program. The approach avoids direct time consuming calculation of detailed chemical reactions as well as the errors introduced by coupling the physical and chemical processes. Emissions are predicted accurately since reaction of air/fuel mixture is calculated using the detailed chemical kinetics mechanism. The approach is applied to model a Caterpillar 3401 direct injection compression ignition (CI) diesel engine. In addition we carried out experimental tests on a Toledo 1500 SI gasoline engine, and those results are compared with the proposed computational approach. In all cases the predicted results agree well with the experimental data. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The continued use of the internal combustion (IC) engine is highly dependent on new designs and new fuels: meeting ever stricter emissions regulation; and meeting customer demands for ⇑ Corresponding author. E-mail addresses: [email protected] (T. Korakianitis), chung.nguyen@ greenway.sg, [email protected] (N. Chung). http://dx.doi.org/10.1016/j.apenergy.2015.06.039 0306-2619/Ó 2015 Elsevier Ltd. All rights reserved.

greater efficiencies reducing fueling costs. A greater understanding of the physics and chemistry taking place in the IC engine is key to meeting these demands. The field of engine research is very active with a number of strategies being employed in order to achieve this insight into engine operation. Research mainly consists of the development of an accurate and efficient method to simulate and model the internal combustion engine with a variety of configurations and fuels. This field can further be broken down into those conducting experimental work with real engines and those

150

T. Korakianitis et al. / Applied Energy 156 (2015) 149–158

Nomenclature Latin G Nc Nr t

arbitrary variable number of species in reaction number of reaction time

Greek

r

difference operator

Superscripts 0 standard state, reference point c chemical n step index Subscripts 1 species 1 (fuel) b backward

constructing accurate mathematical models to simulate the engine. Both of these fields feed into each other with the purpose to increase accuracy and efficiency at the design stage. Alternative fueling of IC engines is predominated with experimental work in laboratory based engines with fuels such as RME and DME [1–3] used extensively along with natural gas, often in a compression ignition design. These works such as those by McTaggart Cowan [4] employ direct measurement of emissions and any performance analysis is firmly routed in the first law of thermodynamics frame work, with engine pressures and the heat release rate ubiquitous in the field. There is a need to model these alternative fuels on a mathematical basis without the need to resort to empirical data all the time and to address emissions along with performance. The coupling of first law analysis to empirically collected data is used extensively to provide quick and easy to interpret data concerning an engines performance. The majority of these first law models, such as that developed by Payri et al. [5] are often described as zero-dimensional, with the combustion chamber treated as a single or number of zones. Within each zone the first law of thermodynamics is applied with respect to time (dt). Whilst this treatment allows for a computationally quick determination of heat releases, temperatures or pressures, emissions are often neglected from such a model or may be inferred as is demonstrated by Kegl [6]. Emissions predictions models have been previously employed based around these zero-dimensional models [7,8], here the combustion chamber is divided into 3 zones: unburned gas; burned gas; and burned gas zone adjacent to the combustion chamber walls. The burned gas and unburned gas zones are separated by the flame front. It is assumed that each zone, burned or unburned is completely and instantaneously mixed (uniform in temperature, pressure, and composition). The combustion is calculated based on the propagation speed of the flame front. The nitric oxide emissions are calculated using the extended Zeldovich kinetic scheme [9] with the steady-state assumptions made for N concentration and equilibrium values used for H, O, O2, N2, and OH concentrations. Because of the simplified assumptions made, this method is normally used only for parametric studies. Multidimensional modeling of an IC engine really tries to address the complex fluid flows taking place in the engine throughout a cycle or over a number of cycles. Here the temporal and spatial variations of the flow field, temperature, composition, pressure and turbulence within the combustion chamber are dealt

f m

forward species m

Abbreviations AC adaptive chemistry CFD computational fluid dynamics CFR fuel co-operative research CI compression ignition CPU central processing unit DAC dynamic adaptive chemistry DME dimethyl ether HCCI homogeneous charge compression ignition IC internal combustion ILDM intrinsic low dimensional manifold ODE ordinary differential equation RME rape methyl ester SI spark ignition QSSA quasi-steady-state assumptions

with [10]. These models may range from a less complex one dimensional model [11] to full three dimensional models examining the complex in-cylinder air motion characteristics of squish, swirl and turbulence [12]. This has become a major part of the engine design and development stage as both CI and SI engines have moved towards direct injection. This paper deals with how complex multidimensional modeling is used to model an event involving combustion, as we do with IC engines. Researchers and designers have employed a number of techniques and with the rapid progress in computer technology, ever more detailed modeling techniques are available including: full CFD to model the flow field characteristics; and detailed chemical kinetics to detail the species formation during combustion. There have been numerous methods to integrate the fluid dynamic to the combustion simulation. In the single zone approach [13,14], during the combustion process the combustion chamber is assumed to be a homogeneous volume (temperature, pressure, and composition are uniformly distributed over the considered volume). Modeling starts from the beginning of ignition to the end of combustion. The combustion of the air/fuel mixture is handled by detailed chemical kinetics, with a number of chemical kinetics modeling packages available, such as HCT [15] and SENKIN, a subprogram in CHEMKIN computer package [16] The state of fluid in the cylinder chamber at the beginning of combustion is specified using different methods. A CFD analysis of the cylinder describes the cylinder air/fuel mixture from the beginning of compression to the point of ignition. The results of this step are then used as input for the detailed chemical kinetics. An example of this two step scheme can be seen in the work of Ogink and Golovitchev [14]. Heat transfer between the mixture and the cylinder wall in the single-zone models is computed using the Woschni method [17]. Since combustion of the air/fuel mixture is calculated for the entire combustion chamber volume, this approach has several disadvantages: (1) fluid dynamics is not taken into account when calculating combustion and (2) the simplicity of the heat transfer model during the compression and expansion strokes does not describe the complexity of the thermal processes. The approach predicted well the trends (qualitative results) of emissions, but failed to predict emissions quantitatively. A multi-zone approach [18] is really an extension of the single zone approach. Here the combustion chamber is divided into three typical zones: (1) constant volume zone representing crevices; (2) zone representing the thermal boundary layer; and (3) core zones.

T. Korakianitis et al. / Applied Energy 156 (2015) 149–158

Each zone is assumed to be homogeneous. The zones have different temperatures, pressures, and composition. The advantage of this approach over the single-zone approach is that temperature and pressure in the combustion chamber are specified for each zone. The disadvantages are similar to that of the single-zone approach. The predicted CO2 is 60% and predicted CO is 80% of the experimental data, with worse accuracy for NOx emissions [18]. It is clear that with both the single and multi zone approach, the multi dimensional fluid dynamic simulations are decoupled from the detailed chemical kinetics once combustion has commenced. Direct coupling of the detailed chemical kinetics and the detailed multi dimensional flow field modeling is desirable to meet modern IC engines requirements. There are a number of obstacles to overcome to fully realize this technique. The detailed chemical reaction mechanisms required in the modeling of the combustion taking place in an IC engine are formally described as systems of ordinary differential equations (ODE). These coupled equations frequently have widely disparate characteristic time scales. The chemical timescales cover several orders of magnitude (1010 –100 s, [19]) referred to as ‘‘stiffness’’ [20]. The shear number of species and reactions involved in the chemical mechanism make the solution of these stiff ODE in the CFD framework extremely CPU intensive. For example, the oxidation of n-heptane contains 560 species and 2539 reactions [21]. Previous studies have demonstrated that even with mechanisms with species less than 50, over 90% of the CPU time was consumed in solving the ODE systems in the reactive flow calculation [22]. A solution to this problem is to reduce the reaction mechanism. A great deal of work has been applied to this reduction strategy, with a number of techniques introduced in the following section.

2. Mechanism reduction Mechanism reduction aims to identify unimportant species and reactions in the detailed chemical kinetic mechanism and remove these from any calculation. This must be carried out with no loss of validity and maintain the important features of the calculation. A number of methods have been proposed to do this. A detailed review of some of these techniques are described by Law [23] and Androulakis [24]. Here we see that the majority of the methods developed fall into two categories: skeletal reduction that eliminates unimportant species and the reactions; and time-scale reduction that moderates the stiffness. Global reduction approaches include: globally lumping techniques; sensitivity analysis; optimization-based approaches; and element flux analysis approaches [25]. These global reduction techniques lead to a skeletal kinetic mechanism, which is used for the entire simulation. There are accuracy compromises with these techniques due to the range of reaction space where the skeletal mechanism are considered valid. The single skeletal model may represent chemical activity well in one region of the flow field, but it’s accuracy reduces over the entire composition/temperature space of interest. An approach developed to overcome this problem has been the so called adaptive chemistry approach [26]. Most of the adaptive chemistry techniques pre-generate a library of mechanisms. A range of the overall thermochemical conditions expected to be encountered in a given reacting flow calculation decides about validity of any sub-mechanism included in the library. At each step of the CFD calculation, the most suitable sub-mechanism is selected to match the current local conditions. Dynamic adaptive chemistry (DAC) scheme has been proposed to overcome some of the challenges faced by adaptive chemistry scheme. The DAC scheme is based on directed relation graph with error propagation (DRGEP) method. Liang et al. [27,28] have recently applied the DAC technique with CHEMKIN and KIVA for

151

a simulated HCCI engine and surrogate fuels. The method showed excellent agreement with the pressure and species mass fraction results obtained using the full mechanisms for n-heptane mechanism (578 species) and a reduced set (178 species) as parent mechanisms. The DAC method gave a 30- and 4-fold reduction in CPU times. These methods have become known as an ‘‘on the fly’’ method and have been implemented in various ways by a number of researchers incorporating element flux analysis and identification of flux magnitudes [25], all aiming to reduce CPU intensiveness. Stiffness of the ODE is tackled by a number techniques such as quasi-steady-state assumptions (QSSA) and partial equilibrium assumptions. Combination of the quasi-steady state approximation (QSSA) and the partial equilibrium assumption (PA) attempt to replace the stiff ODEs of the steady-state species by their corresponding algebraic equations [24]. These techniques have been built upon using methods that automatically identify the fast and slow components of the chemical reacting mechanism by analysis of the timescales of chemical systems and a decoupling of the mechanisms. The intrinsic low dimensional manifolds method (ILDM) introduced by Maas and Pope in 1992 [29] is one such technique. The method has been used by a number of researchers [30]3119 and shown considerable improvements in CPU usage. In another study conducted by Maas and Pop [32], another technique named as ‘‘Trajectory Generated Low Dimensional Manifolds (TGLDM)’’ was proposed to simplify the chemical kinetics. The technique shared many features of ILDM. Reaction trajectory was presented as a main theme and was defined by the following equation.

d /ðtÞ ¼ Sð/ðtÞÞ dt

ð1Þ

Spatial homogeneity was assumed and it was governed by Eq. (1). At t = 0, /ð0Þ ¼ / is any point in a defined realizable region CR . The solution of the Eq. (1) after time t provides a new value for / that is denoted by /R ðt; / Þ. The curve developed by /R : t ½0; 1Þ is called reaction trajectory generated by / and the coordinates of this trajectory are specified by /R ðt; / Þ. The current study uses the concept of reaction indexing and polynomial fitting for the production rates as a function of the reaction conditions (as explained in the text below) instead of mapping as proposed by Maas and Pop [32]. In the database, reaction results of different initial conditions (temperature, pressure, and composition at the beginning of an arbitrary time interval) are stored in different zones; the zones are indexed using their reaction conditions. In another study In Situ Adaptive Tabulation (ISAT), Pope [33] presented a technique to achieve computationally efficient chemistry. By using a mechanism containing 16 species and 41 reactions, he modeled a non-premixed methane – air mixture in statistically homogeneous reactor. The author reported very good control of the tabulation errors with respect to a specified error tolerance as well as the calculations were speeded up by a factor of 1000 when compared to the time taken by the direct integration of the reaction equations. It was indeed a breakthrough so far as the application of the probability density function (PDF) methods were concerned in the calculation of turbulent combustion. Matorakos and his co-workers have made large contributions in the development and application of different techniques that proved helpful in integrating chemistry into turbulent flows generally encountered in the environments similar to IC engines. This also includes large eddy simulation/conditional moment closure (LES–CME) [34] and a recent article examining the influence of turbulence–chemistry interaction for n-heptane spray combustion under diesel engine conditions with emphasis on soot formation and oxidation [35]. Kraft and his co-investigators have investigated

152

T. Korakianitis et al. / Applied Energy 156 (2015) 149–158

different solver strategies and have reported quality work on the turbulence–chemistry interaction [36,37]. The chemistry–spray interaction as implemented in KIVA codes can be found in Ref. [38]. 3. Numerical method In the present paper we propose a pre-integrated non-equilibrium combustion response mapping approach for prediction of emissions of IC engines. The method has previously been employed by Korakianitis et al. [39,40] to predict the emissions in a gas turbine engine. In this previous work an algebraic response model has been formulated to take the place of the ODE. The method is to couple a detailed chemical kinetics model with a multi-dimensional fluid dynamics model while at the same time keeping all the chemical information from most reactive species in a manner that this information is not directly coupled to the high eigenvalue (fast components). 3.1. The pre-integrated response mapping approach Prior to the simulation of the engine SENKIN is used to calculate detailed kinetics reaction of the air/fuel mixture in a variety of initial reaction conditions (temperatures, pressures, and compositions). The results are saved in a database file in the physical timescale. During the simulation process, KIVA-3 V is used to calculate the motion of piston and valves, injection, collision, breakup and evaporation of spray droplets, heat transfer, fluid movement, and diffusion of species. During the combustion process of air/fuel mixture, instead of directly calculating detailed chemical kinetics reaction of the mixture, the reaction results are retrieved from the database file. The procedure to couple SENKIN and KIVA-3 V

in the pre-integrated approach is illustrated in Fig. 1. There are two main requirements for the pre-integrated approach: (1) results of the needed reaction must always be found in the database file and (2) those results must be ready to use in the fluid dynamic calculation process. To satisfy the first requirement, the database file must contain reaction results of a very wide range of reaction conditions. To fulfill this SENKIN is executed in batch mode to calculate reaction of air/fuel mixture at different temperatures, pressures, and compositions. The reaction results (temperature, pressure, species concentrations, etc.) are saved with their absolute values in discrete finite moments determined in the input file. For convenience in the retrieval process, the reaction results are saved in the form of production rates so that they can be obtained as production rate times for each computational time step. All the reaction results related to species concentration are also normalized to the total density of the mixture. This process is illustrated in Fig. 2. For an arbitrary reaction result G (it could be production rates of species, or rate of increase of temperature, etc.), which has a value of G0 at the beginning of a filtering cycle t 0c , where subscript c denote chemical and the superscript is the step index. After one step of calculation, at the time t 1c , SENKIN returns the result with value of G1 , and the duration of reaction is Dt1c ¼ t1c  t0c . The filter compares Dt1c to a physical time scale t p (order of about 105 —104 depending on the particular application), if Dt1c P t p the production rate of G is calculated as:

G_ ¼ ðG1  G0 Þ=Dt 1c and this value is stored in the result file along with its reaction condition. If Dt 1c < tp no result is stored in the result file. The computation process with SENKIN continues until step tnc is reached, at which Dt nc P tp corresponding to the values Gn with Dt nc ¼ t nc  t 0c . At this moment production rate of G is stored as:

ðGn  G0 Þ=Dt nc and a new filtering cycle starts. This procedure ensures all chemical reaction results are stored in the physical timescale t p , so that they are ready to use in the fluid dynamics process. Computations are expedited by reducing the number of species accounted for in the calculation. Only N c species (the most active species and species of interest, in our application 12 species are used in KIVA-3 V) are chosen to represent the behavior of reaction involving N c species. The conservation of mass is satisfied by an

Fig. 1. Pre-integrated non-equilibrium combustion response mapping calculation algorithm.

Fig. 2. Physical timescale filtering algorithm.

153

T. Korakianitis et al. / Applied Energy 156 (2015) 149–158

equivalent N2 technique, in which all the remaining species are added in the nitrogen molecular concentration. This equivalent N2 adding technique causes a change in the mixture structure and affects the reaction results. However, the overall dominance of N2 molecule concentration in the mixture is large enough to minimize this effect. After all possible reactions are computed, a set of relations between reaction conditions and their results are created. Although these conditions and results are discretely stored in the result file, the result file is extremely large. It is computationally inefficient to find a particular reaction result in a very large database. To solve this problem the reaction results are stored as functions of reaction condition (i.e. Fourier or polynomial form). Because of its simplicity, the polynomial form is chosen to represent the relation between conditions and results of reactions. The coefficients of polynomials are stored in the database file in matrix form, along with the their reaction conditions. Once the database is created, engine simulation and emissions prediction are implemented as follows. An interpolating subroutine, which retrieves reaction results from the database file is added in the KIVA-3 V code. At the beginning of the simulation KIVA-3 V reads in geometry, operating condition, and valve shift profiles of the engine. The state of fluid in the cylinder chamber, heat transfer, fuel injection, collision, coalescence, and breakup of droplet are calculated. Transport of the fluid is also computed. At the moment the spark plug ignites (or temperature and pressure in the combustion chamber are high enough for the mixture to auto-ignite), KIVA-3 V sends a request, which contains reaction condition (composition, temperature, pressure of the mixture, and reaction time) to the interpolating subroutine. Based on these indices, the interpolating subroutine locates the position of reaction result in the database file (specifies where to find the coefficients for the fitted functions of the requested reaction). The normalized production rate is converted back to absolute value, and the reaction result is:

G2 ¼ G1 þ G_  Dt;

ð2Þ

where Dt is requested reaction time, G1 is the value sent to the interpolating subroutine to request reaction results (initial condition for this reaction), and G2 is the new value after reaction. Then these reaction results are sent back to KIVA-3 V to compute the new state of the fluid in the combustion chamber. The motion of piston and valves, and the new volume of the combustion chamber are calculated for a new computational cycle. The reaction result retrieval is illustrated in Fig. 3. 3.2. SENKIN detailed kinetics reaction modeling program SENKIN is a subprogram in CHEMKIN-II computer package written by Lee et al. [16] in Sandia National Laboratories. SENKIN can be used to predict the time-dependent chemical kinetics behavior of homogeneous gas mixtures in a closed adiabatic system. The results of SENKIN are later used to provide chemical kinetics information for each individual computation cell in KIVA-3 V. KIVA-3 V accounts for changes in cell volume, and for fluid dynamics, heat transfer, and species diffusion processes between its cells. SENKIN requires two input files to calculate the chemical reactions of a mixture: (1) an input file, which contains the initial condition at the beginning of reaction (temperature, pressure, species, concentrations of species, reaction time, critical temperature of the mixture, etc.) and (2) a binary linking file, which consists of all thermal properties and chemical reaction mechanisms of involving species. The linking file is created by an interpreter program based on the thermodynamic data and reaction mechanism of the considered species, which are stored in the thermodynamic database and

Fig. 3. Reaction results interpolation algorithm.

reaction mechanism files, respectively. For reactions of a mixture involving m species vm and r reaction with stoichiometric coefficients m0 and m00 : Nc X

m0m;r vm ()

m¼1

Nc X

m00m;r vm ðr ¼ 1; . . . ; Nr Þ

ð3Þ

m¼1

SENKIN solves the equation system of conservation of species and energy:

dY m _ m W m m ¼ 1; . . . ; Nc ¼ vx dt Nc X dT _ mWm ¼ 0 cv þv em x dt m¼1

ð4Þ

where Y m is the mass fraction of species m; W m is the molecular weight of species m; N c is total number of species in the system, v is the specific volume, e is the specific energy, p is the pressure, cv is specific heat at constant volume, and T is temperature. The pro_ m is: duction rate of species m; x

x_ m ¼

Nr Nc Nc X Y Y 0 00 ðm00m;r  m0m;r Þ kf ;r ½vm mm;r  kb;r ½vm mmr r¼1

m¼1

! ð5Þ

m¼1

where ½vm  is the molar concentration of species m. Forward and backward rate coefficients of reaction r; kf ;r ; kb;r are calculated by Arrhenius law.

kf ;r ¼ Ar T br exp kb;r ¼

  Er RT

kf ;r K c;r

where K c;r is obtained from:

ð6Þ

154

T. Korakianitis et al. / Applied Energy 156 (2015) 149–158 Nc patm X ðm00  m0m;r Þ RT m¼1 m;r ! DS0r DH0r ¼ exp  R RT

K c;r ¼ K p;r K p;r

ð7Þ

Entropy of the mixture is:

S¼

  Nc  X p ½v m  S0m  Rln½vm   Rln patm m¼1

ð8Þ

Closure of the system is attained when specific heat, enthalpy, and entropy of the mixture are computed using polynomial coefficients from the thermodynamic database file. 4. Experimental setup and the simulation parameters This section details the specifications of the engines used for the validation of the results. It also include a range of initial and boundary conditioned chosen to run the simulations for the two engines selected. 4.1. Specification and modeling parameters for the diesel engine The pre-integrate approach was applied to model a Caterpillar 3401 direct injection CI diesel engine [41]. The computational meshes for all the modeling cases were created by ICEM-CFD, a mesh generation package with ANSYS. Database files were created for each of the two fuels on two Dell Optiplex GX20. SENKIN was used to calculate reaction of air/fuel mixture with equivalence ratio from 0.2 to 5.0 in increments of 0.05; the initial temperatures of reactions was from 700 K to 1800 K with steps of 10 K. Initial pressures of reactions varied from 10 bar to 140 bar with interval of 1 bar. The modeling parameters for the Caterpillar 3401 are listed in Table 1. The computational mesh consists of 92,000 cells, partly illustrated in Fig. 4(a). Because of the symmetry of the combustion chamber, the engine geometry is created for 60 degrees. In this model n-heptane was used to simulate reaction of diesel [42]. The simulation took 9.5 h.

created for half the engine cylinder due to its symmetry (Fig. 4(b)). A mixture of 90% iso-octane (C8H18) and 10% n-heptane (C7H16) was used to model the reaction of gasoline (this mixture is equivalent to gasoline, which has octane number of 90) [43,44]. The simulation took 7.3 h to complete. The reaction mechanism consists of 716 species and 2541 reactions. The equivalent N2 adding technique was used to reduce the number of mixture representing species to 12 (fuel, O2, N2, CO2, H2O, H, H2, O, N, CO, and NO). The chemical reactions were calculated for the mixture of gasoline and air with equivalence ratio varied from 0.2 to 5.0 in increments of 0.05; the initial temperatures of reactions was from 700 K to 1800 K with steps of 10 K. Initial pressures of reactions varied from 10 bar to 140 bar with interval of 1 bar. For experimental investigation, experiments were carried out on a Toledo 1500 four-stroke SI gasoline engine. The engine is fitted with an HD hydraulic dynamometer. The averaged air flow rate was measured using an air box with an orifice plate. The design of the air box follows BS 1042 (1964). The in-cylinder pressure was measured using a piezoelectric pressure transducer (Kistler 6061B) with measuring range from 0 to 250 bars. To reduce the effect of thermal shock on measured pressure, the transducer is water cooled. The signals (25pC  1 bar) from the transducer were amplified using a Kistler charge amplifier (5041E type) before being fed to a personal computer (PC). CO2, CO, and NO concentrations were measured using the sampling technique. The exhaust gas was sampled from a point 60 mm from the exhaust manifold outlet. The exhaust gas was pumped through the measuring system using a vacuum pump. Water vapor and particulates larger than 1 lm were filtered out. The exhaust gas was pumped to the sample cell of the nondispersive infrared gas analyzer, where CO, CO2, and unburned hydrocarbon (HC) were measured, then the sample gas was led through NOx and O2 electrochemical sensors, and driven out. The Toledo 1500 engine was tested at 2/3 WOT, ignition timing at h ¼ 30 BTDC, and the dynamometer was adjusted to vary engine speed from 1500 rpm to 3500 rpm.

5. Results 5.1. Comparison between calculated results and data in the literature

4.2. Specification and modeling parameters for the gasoline engine For further validation, numerical simulation and experimental investigation on the operation of a Toledo 1500 four-stroke gasoline SI engine were carried out, and the corresponding results were compared. The specifications of the engine are presented in Table 2. The engine operating conditions modelled are described in Table 3. In numerical simulation the pre-integrated approach is applied, on a mesh of 68,000 cells. The computational geometry was

Table 1 Simulation parameters of the Caterpillar 3401 diesel engine (from [25]). Rotational speed Cylinder wall temperature Cylinder head temperature Piston face temperature Volume displacement Inlet air pressure Inlet air temperature Exhaust gas pressure Exhaust gas temperature Injection pressure Fuel injected per cycle Injection crank angle Duration of injection pulse Number of spray particle

1600 rev/min 525 K 525 K 525 K 2.44 l 184 kPa 310 K 200 kPa 580 K 90 MPa 0.1622 g/cycle 349 21.5 9000

The calculated pressure profile of the Caterpillar 3401 is presented in Fig. 5. The calculated pressure profile agrees well with the measured data. The difference between the predicted and measured pressures lies within a range of 93–106%. From h ¼ 340 to h ¼ 355 the predicted pressure is lower than the measured one. There is a pressure shift at h ¼ 355 , where the calculated pressure increases faster than the measured one: the predicted pressure profile crosses, and from then on lies above the experimental profile. The pressure shift could be a result of adding the heavy intermediates into the equivalent nitrogen molecular concentration. The pressure then would increase as one molecule of a heavy intermediate species results to the addition of a few equivalent nitrogen molecules. The model predicts correctly the peak pressure with 5% higher than the measured value The model failed to capture the rate of rise in the in-cylinder pressure in the premixed combustion phase. The rate of heat released for the Ctaerpillar engine has also been presented. Both of the predicted heat release peaks stay higher when compared to the corresponding rate of the heat release curves calculated from the experimental data. The first heat release peaks for both the predicted as well as experimental value occur at the same crank angle position although there is slight delay in the predicted curves for the rate of heat release. Although the second peak for the predicted rate of heat release is higher than the corresponding rate of heat release obtained experimentally but the rate of rise in the experimental heat release is

155

T. Korakianitis et al. / Applied Energy 156 (2015) 149–158

Fig. 4. Modeling meshes used for Caterpillar 3401 diesel engine (a) and Toledo 1500 gasoline engine (b).

Table 2 Specifications of the Toledo 1500 SI gasoline engine.

Table 3 Simulation parameters of the Toledo 1500 gasoline engine.

Main specifications Bore Stroke Squish Connecting rod Compression ratio Cylinder head Piston head

73.66 mm 87.38 mm 10.00 mm 241.30 mm 9:1 Flat Flat

Cam timing IVO

60 BTDC

IVC

870 ABDC

Peak lift inlet

8.1 mm at 440 BBDC

EVO

350 BBDC

EVC

60 ATDC

Peak lift exhaust

8.1 mm at 760 ABDC

higher than the corresponding rate of heat release predicted computationally. The calculated in-cylinder NO history of the Caterpillar 3401 is illustrated in Fig. 6. The predicted results is shown per crank angle in the cylinder, and comparison is made with the available

Rotational speed Ignition timing Fuel Cylinder wall temperature Cylinder head temperature Piston face temperature Inlet air pressure Inlet air temperature Exhaust gas pressure Exhaust gas temperature

(1500–3500) rpm from 30° BTDC to 100 ATDC C8H18 and C7H16 580 K 580 K 580 K 103 kPa 320 K 135 kPa 550 K

experimental data, which was obtained from a surge tank at the engine exhaust. The higher values of the predicted NO can be explained as follows. The formation of No is a very strong function of temperature as well as the concentration of the oxidizer. As it is evident from the in-cylinder pressure as well as the rate of heat release rate curves, the predicted values of the in-cylinder pressure and the rate of heat release are higher when compared to the experimentally obtained values. The higher pressure in case of the predicted in-cylinder pressure results in higher temperature

156

T. Korakianitis et al. / Applied Energy 156 (2015) 149–158

Fig. 5. In-cylinder pressure history (bar) and rate of energy release (MJ) of the Caterpillar 3401 diesel engine (1600 rpm, injection from h ¼ 349 to h ¼ 354 ).

Fig. 6. In-cylinder NO history (grams) of the Caterpillar 3401 diesel engine (1600 rpm, injection from h ¼ 349 to h ¼ 354 ).

and hence higher NO. Although, the first peaks of the rate of heat release occurs at approximately same crank angle position but the predicted heat release peak is higher than the corresponding heat release peak obtained experimentally. 5.2. Comparison between calculated results and experimental data The predicted and measured in-cylinder pressures of the Toledo 1500 SI engine (operating at 2000 rpm, 2/3 WOT) are illustrated in Fig. 7. It can be seen that the maximum predicted pressure is about 8% higher than the maximum measured value. The model predicts the maximum pressure at h ¼ 375 , while the measured peak pressure occurs at h ¼ 378 . It is evident from

Fig. 7 that the model under predicts the in-cylinder pressure during the final stages of the compression stroke, but near the TDC, the rate of pressure rise for the computationally calculated in-cylinder pressure overtakes the rate of the experimentally calculated rate of pressure rise. Similarly, when the two heat release rate curves are compared, higher peak is predicted when compared to the heat release rate peak calculated form the experimental data. The variations of CO2, CO, and NO concentrations are illustrated in Figs. 8(a), (b) and 9, respectively. The increase of engine speed is associated with the lowering of load. From the speed of 1500 rpm, the residual gas in the cylinder reduces and engine revolution increases with lower load. The increase of engine speed relatively moves the ignition point towards TDC at constant ignition timing. The combustion process moves further into the expansion stroke. The combustion of the air/fuel mixture still completes close to TDC when the engine speed increases to 1900 rpm. On the other hand, the fuel consumption increases with increasing engine speed, while the cooling water flow rate is constant, leading to an increase of the in-cylinder temperature,which in turn accelerate the combustion rate. CO2 concentration is seen increasing from the speed of 1500 rpm to 1900 rpm. NO concentration is also observed increasing with engine speed within this range, and a maximum peak of NO concentration occurs at the speed of 1900 rpm. From the speed of 1900 rpm CO2 and NO concentrations decrease with increase of engine speed, and minimum peaks of CO2 and NO concentrations occur at the speed of 2200 rpm. This is explained as follows. The retarded ignition leads to reduction of the in-cylinder maximum temperature. On the other hand, the combined effect of reduction of the residual gas and the increase in fuel consumption leads to increase of the in-cylinder maximum temperature. However, in these cases the effect of the retarded ignition on temperature is greater than that of the residual gas and the fuel consumption. The overall effect is that the combustion of the unburned mixture moves further towards the expansion stroke, leading to decrease of the in-cylinder maximum temperature and pressure. Therefore the production rates of CO2 and NO reduce. For speeds above 2200 rpm the added energy from the fuel is larger than the energy taken away by the exhaust, cooling water, and radiation. More vigorous combustion occurs at high temperature, and the production rates of CO2 and NO monotonically increase with engine speed. The trends for different engine emissions can also be explained on the basis of equivalence ratio (/) or relative air fuel ratio (k). The change in engine changes the equivalence ratio. At higher speeds the engine is running hotter due to higher equivalence ratio resulting more complete combustion. These trends are consistent with the data reported in the literature. For a given compression ratio, the NO increases initially, reaches a maximum value and then drops off. This can be attributed to the fact that the NO formation is a strong function of combustion chamber temperature, mixture strength and available oxygen and the formation generally occurs in the post flame gases. With air fuel ratio getting leaner, the combustion temperature drops resulting in the weakening of the NO formation. For lower equivalence ratios, the thermal quenching during the NO formation process restricts the process and hence result in lower NO. In addition to the factor mentioned above there are other factors that may influence the formation of NO. The amount of heat being transferred to the cylinder walls also changes with changes in the equivalence ratio (being different at different speeds). Similarly the turbulence levels also vary with this variation in the engine speed and can affect the in-cylinder flow fields. Due to these moving flow fields, the center of the flame, that ideally lies at the spark plug gap, can move away.

T. Korakianitis et al. / Applied Energy 156 (2015) 149–158

157

Fig. 7. In-cylinder pressure and rate of energy release of the Toledo 1500 SI engine (2/3 WOT, ignition at 30° BTDC, 2000 rpm).

Fig. 8. Exhaust CO2 (a) and CO (b) variation with speed of the Toledo 1500 SI engine (2/3 WOT, ignition at h ¼ 30 BTDC).

dynamics processes of engine operation; SENKIN is used to pre-calculate reaction of air/fuel mixture. The approach couples the thermal and fluid dynamics processes with the chemical kinetics reaction process in the physical timescale. The chemical reaction is pre-calculated, the reaction results are decoupled from their chemical eigenvalue, and integrated in physical timescale, then stored in a database file. During the chemical reaction process instead of calculating the chemical reaction process directly, KIVA-3 V obtains reaction results from the pre-integrated database file. The pre-integrated approach predicts accurately emissions in IC engines using detailed chemical kinetics reaction, while the computer time is shorten. The approach is applied to predict emissions in a Caterpillar 3401 direct injection compression ignition (CI) diesel engine and d Toledo1500 SI gasoline engine. The predicted values of NOx are 5–15% higher than experimental values, and computer time is saved. Fig. 9. Exhaust NO variation with speed of the Toledo 1500 SI engine (2/3 WOT, ignition at h ¼ 30 BTDC).

6. Conclusions Many approaches have been proposed to predict emissions from IC engines. Each approach has its own advantages and disadvantages. To avoid direct coupling of the physical processes with the chemical process, the pre-integrated approach is developed. The approach couples the advantages of multi-dimensional engine modeling and of detailed chemical kinetics reaction modeling. KIVA-3 V is used to model the thermal dynamics and fluid

References [1] Korakianitis T, Namasivayam AM, Crookes RJ. Diesel and rapeseed methyl ester (RME) pilot fuels for hydrogen and natural gas dual-fuel combustion in compression–ignition engines. Fuel 2011;90(7):2384–95. [2] Korakianitis T, Namasivayam AM, Crookes RJ. Natural-gas fuelled sparkignition (SI) and compression–ignition (CI) engine performance and emissions. Prog Energy Combust Sci 2011;37(1):89–112. [3] Korakianitis T, Namasivayam AM, Crookes RJ. Hydrogen dual-fuelling of compression ignition engines with emulsified biodiesel as pilot fuel. Int J Hydrogen Energy 2010;35(24):13329–44. [4] McTaggart-Cowan GP, Jones HL, Rogak SN, Bushe WK, Hill PG. The effects of high-pressure injection on a compression–ignition, direct injection of natural gas engine. J Eng Gas Turb Power – Trans ASME 2007;129(2):579–88 [ISI

158

[5] [6] [7]

[8]

[9] [10] [11]

[12]

[13]

[14]

[15]

[16]

[17]

[18] [19] [20] [21]

[22] [23] [24]

T. Korakianitis et al. / Applied Energy 156 (2015) 149–158 Document Delivery No.: 159AA Times Cited: 2 Cited Reference Count: 31 McTaggart-Cowan, G.-P. Jones, H.L. Rogak, S.N. Bushe, W.K. Hill, P.G. ASMEAMER SOC MECHANICAL ENG NEW YORK]. Payri F, Olmeda P, Martn J, Garca A. A complete 0D thermodynamic predictive model for direct injection diesel engines. Appl Energy 2011;88(12):4632–41. Kegl B. Influence of biodiesel on engine combustion and emission characteristics. Appl Energy 2011;88(5):1803–12 [Times Cited: 1]. Hires SD, Ekchian A, Heywood JB, Tabaczynski RJ, Wall JC. Performance and NOx emissions modeling of a jet ignition pre-chamber stratified charge engine. SAE paper 760161, 1976. Heywood JB, Hinggins JM, Watts PA, Tabaczynski RJ. Development and use of a cycle simulation to predict SI engine efficiency and NOx emissions. SAE paper 790291, 1979. Warnatz J, Mass U, Dibble RW. Combustion. 3rd ed. New York: SpringerVerlag; 2000. ISBN 3-540-677518. Ramos JI. Internal combustion engine modelling. United States: Hemisphere Publishing; 1989. D’Errico G, Cerri T, Pertusi G. Multi-objective optimization of internal combustion engine by means of 1D fluid-dynamic models. Appl Energy 2011;88(3):767–77. Prasad B, Sharma CS, Anand TNC, Ravikrishna RV. High swirl-inducing piston bowls in small diesel engines for emission reduction. Appl Energy 2011;88(7):2355–67 [ISI Document Delivery No.: 749UW Times Cited: 0 Cited Reference Count: 46 Prasad, B.V.V.S.U. Sharma, C.S. Anand, T.N.C. Ravikrishna, R.V. ELSEVIER SCI LTD OXFORD]. Flowers D, Aceves S, Westbrook CK, Smith R, Dibble R. Detailed chemical kinetic simulation of natural gas HCCI combustion: gas composition effects and investigation of control strategies. ASME J Eng Gas Turb Power 2001;123:433–9. Ogink R, Golovitchev V. Gasoline HCCI modelling: computer program combining detailed chemistry and gas exchange process. SAE Paper 200101-3614, 2001. Lund CM. HCT-A general computer program for calculating time-dependent phenomena involving one-dimensional hydrodynamics, transport, and detailed chemical kinetics. Technical report UCRL-52504, Livermore: Lawrence Livermore National Laboratory; 1978. Kee RJ, Rupley FM, Miller JA. CHEMKIN-II: A Fortran chemical kinetics package for the analysis of gas-phase chemical kinetics. Technical report SAND898009, Livermore (CA, USA): Sandia National Laboratories; 1980. Woschni G. A universally applicable equation for the instantaneous heat transfer coefficient in the internal combustion engine. SAE paper 670931, 1967. Ogink R, Golovitchev V. Gasoline HCCI modeling: an engine cycle simulation code with a multi-zone combustion model. SAE paper 2002-01-1745, 2002. Nafe J, Maas U. Modeling of NO formation based on ILDM reduced chemistry. Proc Combust Inst 2002;29:1379–85. Westbrook CK, Mizobuchi Y, Poinsot TJ, Smith PJ, Warnatz E. Computational combustion. Proc Combust Inst 2005;30:125–57. Youngchul Ra, Reitz RD. A reduced chemical kinetic model for IC engine combustion simulations with primary reference fuels. Combust Flame 2008;155:713–38. Kaiyuan He, Ierapetritou MG, Androulakis IP. Integration of on-the-fly kinetic reduction with multidimensional CFD. AIChE J 2010;56:1305–14. Law Chung K. Combustion at a crossroads: status and prospects. Proc Combust Inst 2007;31(1):1–29. Androulakis I. store and retrieve representations of dynamic systems motivated by studies in gas phase chemical kinetics. Comput Chem Eng 2004;28:2141–55.

[25] Kaiyuan He, Androulakis IP, Ierapetritou MG. Numerical investigation of homogeneous charge compression ignition (HCCI) combustion with detailed chemical kinetics using on-the-fly reduction. Energy Fuels 2011;25:3369–76. [26] Bhattacharjee B, Schwer DA, Barton PI, Green WH. Optimally-reduced kinetic models: reaction elimination in large-scale kinetic mechanisms. Combust Flame 2003;135:191–208. [27] Long L, Stevens John G, Farrell John T. A dynamic adaptive chemistry scheme for reactive flow computations. Proc Combust Inst 2009;32:527–34 [Times Cited: 19]. [28] Liang L, Stevens John G, Raman S, Farrell John T. The use of dynamic adaptive chemistry in combustion simulation of gasoline surrogate fuels. Combust Flame 2009;156(7):1493–502 [Times Cited: 10]. [29] Maas U, Pope SB. Simplifying chemical-kinetics – intrinsic low-dimensional manifolds in composition space. Combust Flame 1992;88:239–64. [30] Lovas T, Amneus P, Mauss F, Mastorakos E. Comparison of automatic reduction procedures for ignition chemistry. Proc Combust Inst 2002;29:1387–93 [Times Cited: 21 Part 1 29th International Combustion Symposium JUL 21-26, 2002 HOKKAIDO UNIV, SAPPORO, JAPAN]. [31] Blasenbrey T, Maas U. ILDMs of higher hydrocarbons and the hierarchy of chemical kinetics. Proc Combust Inst 2000;28:1623–30. [32] Pop SB, Maas U. Simplifying chemical kinetics: trajectory generated low dimensional manifolds. FDA 1993:93–111. [33] Pope SB. Computationally efficient implementation of combustion chemistry using in situ adaptive tabulation. Combust Theory Model I 1997:41–63. [34] Zhang H, Garmory A, Cavaliere DE, Mastorakos E. Large Eddy Simulation/ Conditional Moment closure modeling of swirl-stabilized non-premixed flames with local extinction. Proc Combust Inst 2015;35:1167–74 [ISSN 1540-7489]. [35] Bolla M, Farrace D, Wright YM, Boulouchos K, Mastorakos E. Influence of turbulence–chemistry interaction for n-heptane spray combustion under diesel engine conditions with emphasis on soot formation and oxidation. Combust Theory Model 2014;18:330–60 [ISSN 1364-7830]. [36] Asproulias I, Revell A, Craft TJ. An investigation into solver strategies for the modelling of compressible turbulent flows. In: Proc. 28th int. symposium on shock waves, 2011 [eScholarID:189780]. [37] Klein T, Craft TJ, Iacovides H. Developments in two-time-scale turbulence models applied to non-equilibrium flows. In: Proc. 7th int. symposium on turbulence and shear flow phenomena, 2011 [eScholarID:189776]. [38] Golovitchev VI, Nordin N. Detailed chemistry spray combustion model for the KIVA code. In: International multidimensional engine modeling user’ group meeting at SAE congress, Detroit, January 2001. [39] Korakianitis T, Dyer RS, Subramanian N. Pre-integrated non-equilibrium combustion-response mapping for gas-turbine emissions. J Eng Gas Turb Power – Trans ASME 2004;126(2):300–5. [40] Dyer RS, Korakianitis T. Pre-integrated response map for inviscid propane–air detonation. Combust Sci Technol 2007;179(7):1327–47. [41] Uludogan A, Xin J, Reitz RD. Exploring the use of multiple injectors and split injection to reduce DI Diesel engine emissions. [42] Curran HJ, Gaffuri P, Pitz WJ, Westbrook CK. A comprehensive modelling study of n-heptane mechanism. Combust Flame 1998;114:149–77. [43] Curran HJ, Gaffuri P, Pitz WJ, Westbrook CK. A comprehensive modeling study of iso-octane oxidation. Combust Flame 2002;129:253–80. [44] Curran HJ, Pitz WJ, Westbrook CK, Callahan CV, Dryer m FL. Oxidation of automotive primary reference fuels at elevated pressures. Proc Combust Inst 1998;27:379–87.