c
TU Delft, The Netherlands, 2006
APPLYING ILDM TECHNIQUE WITH TIME SCALING
FOR CFD SIMULATION OF DIESEL ENGINE
Ravindra H. Aglave, J¨urgen Warnatz†
Universit¨at Heidelberg, Interdisciplinary Center for Scientific Computing INF 368, 69120, Germany
e-mail: aglave@iwr.uni-heidelberg.de † e-mail: warnatz@iwr.uni-heidelberg.de
web page: http://www.iwr.uni-heidelberg.de/groups/reaflow
Key words: ILDM,Computational Fluid Dynamics, Combustion, Reduced Mechanisms Abstract. Detailed mechanisms describing ignition, flame propagation and pollutant formation typically involve several hundred species and elementary reactions, prohibiting their use in practical three-dimensional engine simulations. Conventionally reduced mech-anisms often fail to predict minor radicals and pollutant precursors. The ILDM-method is an automatic reduction of a detailed mechanism, which assumes local equilibrium with respect to the fastest time scales identified by a local eigenvector analysis. In the reactive flow calculation, the species compositions are constrained to these manifolds. Conserva-tion equaConserva-tions are solved for only a small number of reacConserva-tion progress variables, thereby retaining the accuracy of detailed chemical mechanisms. This gives an effective way of coupling the details of complex chemistry with the time variations due to turbulence.
A standard and RNG k-e model is used to model the turbulent flow field. Turbulence-chemistry interactions are taken into account by integrating the chemical source terms over a presumed probability density function (pdf ). The standard KIVA III code along with the above models for chemistry and turbulence were used to simulate the flow in a single cylinder of a DI diesel engine. Since no equilibrium or partial-equilibrium approach is involved, accurate concentrations of O radicals and soot-precursors (C2H2 & C3H3) can
be obtained from the ILDM.
1 INTRODUCTION
The simulation of diesel engines involves modeling several processes, viz. spray dynam-ics, auto-ignition, chemistry, turbulence and heat transfer. One also needs to consider the interactions between chemistry and turbulence. Accurate pollutant prediction requires accurate models for all these processes. In a previous paper1, basic development and
chemistry were carried out. The main advantage of these models is that they are based on detailed chemical mechanisms and hence, do not contain conditional parameters, which need to be adjusted as the engine operating conditions are changed. To recapitulate the study carried out in the previous paper, following can be said. ILDM method is used as the automatic reduction method for the rest of the chemistry modeling. A model fuel (n-heptane) with 1-d ILDM was used with CO2 as the progress variable. The chemistry
turbulence interactions were accounted for, by using a pdf of progress variables and mix-ture fraction, both during chemistry as well as ignition NOx was calculated using the Zeldovich mechanism. Soot calculation were performed using a two equation phenomeno-logical model2. These models were implemented in a modified version of the standard
three-dimensional CFD code, KIVA III3 that is capable of simulating two-phase engine
flows. KIVA solves the three-dimensional Favre-averaged Navier-Stokes equations with a k- turbulence model. The spray dynamics was described by means of a discrete-droplet-model (DDM) with standard sub-discrete-droplet-models for collision and evaporation. The breakup discrete-droplet-model used was a Kelvin-Helmholz type wave breakup model4. The resulting code was used to
model a Caterpillar direct-injection diesel engine for which experimental results were avail-able5. Good agreement was observed between simulated and experimental quantities for
pressure. In this paper, this work is advanced in two ways. The ILDM table used for calculation of detail chemistry is recreated using the improved numerical schemes6.
Sec-ond, the correlation between pdfs of temperature, mixture fraction and reaction progress variables is altered by calculating a Damkohler number (a ratio of chemistry to flow time scales). In the following sections, some of the complex phenomenons modeled by this code are described briefly. In the following sections, the modeling and results due to the new ILDM table and two progress variables are presented.
2 IGNITION
Ignition delay and location has great impact on the performance of the engine. Low temperature reactions are responsible for formation of a radical pool7 from spray droplets
and vapor in auto-ignition. During the ignition delay period, the concentration of this pool increases, mainly due to chain reactions and finally leads to ignition. Although it’s possible to model the detail chemistry of the ignition process8, it’s impractical to use such
a calculations in a CFD code. Instead, ignition is tracked by a representative species viz. CO. However the reaction rates are obtained from a detail model9 in a stirred reactor
environment. A single representative species, viz. CO, is used to track the concentration of radical pool. This species has negligible concentration in the ignition delay period, while it rises rapidly and monotonically during ignition. The transport equation solved for the CO concentration is given by;
Where ρ is the Reynolds-averaged density, gYCO is the Favre-averaged CO mass fraction, ~u
is the velocity vector, DT is the turbulent diffusion coefficient and gY˙COC is the mean source
term due to chemistry. The laminar chemical source term was obtained from the reaction rates for CO, which is calculated for the case of a perfectly stirred reactor. This pre-calculation is performed using a detailed mechanism containing about 200 species involved in about 1200 elementary reactions and is proven to accurately reproduce ignition delay times10 and flame velocities11 in a heptane-air mixture for a range of equivalence ratios.
These rates are tabulated as a function of the mixture fraction (ξ ), the temperature (T), the pressure (p) and the CO concentration itself (YCO ). They are integrated over a
probability density function (pdf) to account for turbulent fluctuations. A presumed pdf method is used, where the shape of the pdf is assumed a priori and the parameters of the pdf are calculated from the moments of the variable (its mean and variance). The mean source term is therefore given by:
g˙ CC CO = Z ˙ YCOC (ξ, T, p, YCO) ˜P (ξ, T, p, YCO) dξ · dT · dp · d (YCO) (2)
Assuming negligible pressure fluctuations, and assuming statistical independence of the variables, the pdf was split into the product of one-dimensional pdf, i.e.,
˜
P (ξ, T, YCO) = ˜P (ξ) · ˜P (T ) · ˜P (YCO) (3)
and beta functions12 were used for the three pdfs. The ranges of the pdf variables (CO
mass fraction and temperature) are those that are tabulated. The CO mass fraction is pre-tabulated up to its critical value of 0.1 and the temperature is tabulated up to 1500 K from a stirred reactor simulation. These values are recalled when the pdf integration is performed during the engine simulation.
3 TURBULENCE MODELING
has an added term ’cǫ3’ and a constant. ∂ρeǫ ∂t + ∇ · ρ˜~u˜ǫ= − 2 3cǫ1 − cǫ3− extra term z }| { 2 3cµcη ek eǫ∇ · ˜~u ρeǫ∇ · ˜~u + ˜ǫ ˜ k " cǫ1− extra constant z}|{c η ! = σ: ∇˜~u − cǫ2ρ˜ǫ + csW˙ S # + ∇ · µef f P rǫ ∇eǫ (4) where cη = η (1 − η/η0 ) 1 + βη3 (5) η = S˜k ǫ (6) S = (2SijSij)1/2 (7) Sij = 1 2 ∂ ˜u~i ∂xj +∂ ˜u~j ∂xi ! (8) η is the ratio of the turbulent to mean strain rate and S is the magnitude of the mean strain Sij. The values of all constants are give in table 1. For rapid distortion in a
compressible, low Mach number flow, following equations give the remaining constants13.
cǫ3 = −1 + 2Cǫ1− 3m (n − 1) + (−1) δ√ 6cµcηη 3 (9) cηη = −1 βη0 (10) m = 0.514 (11) δ = 1 if∇ · ~u < 0 δ = 0 if∇ · ~u > 0
The sign of cǫ3 changes with the sign of the velocity dilatation and the polytropic process,
thus improving the model.Here, ’n’ is the index for a polytropic process, generally taken as 1.4 for engine studies15. The RNG k-ǫ model especially addresses the appreciable ratios
Constant cǫ1 cǫ2 cǫ3 cs P rk P rǫ η0 β cµ
Standard k-ǫ 1.44 1.92 -1.0 1.5 1.0 1.3 - - 0.09
RNG k-ǫ 1.42 1.68 Eq.12 1.5 0.719 0.719 4.38 0.012 0.0845
Table 1: Values of the constants in the turbulence models3,15.
4 TURBULENCE - CHEMISTRY INTERACTIONS
Effects of turbulence on reaction zone structure in non-premixed as well as partially premixed flames have significant repercussions on combustion models. A complete un-derstanding is sought to elucidate the response of mass fractions of species, especially intermediates and minor but critical species, to turbulence16. There can be sever errors in
estimating the mean reaction rates if they can calculated from the mean species concentra-tions. One needs to use the instantaneous species concentrations instead. The approach using multidimensional joint pdf involves solving the velocity-composition joint pdf or at least the composition joint pdf approach (when turbulence is modeled)17,18,19. Alternately,
it is possible to assume the shape of the pdf as an analytical function. The parameters of the pdf are determined from the moments calculated from the transport equations12.
The commonly used assumed shapes are the Gaussian functions20, beta functions21 or
the Dirac’s delta function22. A dynamical approach, in contrast to the conventional
ap-proaches(equilibrium, steady-state assumption or partial equilibrium assumption), which tries to find out the directions in which the source term vector will rapidly reach a steady-state is the ILDM approach8,23,24. The directions are not associated with any particular
reaction or species and are not fixed throughout the combustion process.
5 CHEMISTRY MODELING
The intrinsic low-dimensional manifold (ILDM) method8,23 is an automatic reduction
do not need to be identified a priori. An eigenvalue analysis of the detailed chemical mechanism is carried out which identifies the fast processes in dynamic equilibrium with the slow processes. The computation of ILDM points can be expensive, and hence an in-situ tabulation procedure6 was used, which enables the calculation of only those points
that are needed during the CFD calculation.
For the diesel engine calculation, heptane was used as a model diesel fuel, and a detailed chemical mechanism consisted of 43 species and 393 elementary reactions. Although it is possible to keep the pressure variable in the ILDM method, in the range of pressures important for combustion in diesel engines (50 - 100 bar), the chemistry does not show a significant variation. Hence, the pressure in the mechanism was kept constant, with 80 bars chosen as the representative pressure in the engine. This mechanism was proven to accurately reproduce flame velocities in laminar premixed flames using heptane as the fuel.11 A one-dimensional ILDM with the CO
2 concentration as the progress variable was
used. The performance of the 1-d ILDM was seen to be excellent in comparison with the detailed mechanism in homogeneous reactors, as shown in Figure 1. & Figure 2
This chemistry model was used in KIVA III for engine simulation. The reaction rates and the species concentrations obtained from ILDM are integrated over a presumed pdf in order to use it in turbulent flows. The mean reaction rate of the progress variable is given by: ]˙ YC CO2 = Z ˙ YCOC 2(ξ, T, YCO2) ˜P (ξ, T, YCO2) dξ · dT · d (YCO2) (12)
In accordance with the assumptions made in the ignition model, ξ, T and YCO2 were
assumed to be statistically independent. Therefore, the three variable join pdf is split into three single variable pdfs. Beta pdfs were used for ξ and T . Since YCO2 is highly
dependent on ξ and T , the pdf P (YCO2) was normalised into a pdf with a new variable b,
where b is defined as:
b = YCO2
YCO2,eq
(13)
˜
P (ξ, T, b) = ˜P (ξ) · ˜P (T ) · ˜P (b) (14) Where YCO2,eq is the CO2 concentration at equilibrium. Since YCO2,eq is itself dependent
5.1 RECOMPILING ILDM
The ILDM code has undergone several modifications in terms of its numerical schemes. As a result of this, the species concentrations can be calculated over a wider range of Temperature and mixture fraction ranges. This has removed zones of composition spaces where model points were required to be used instead of ILDM points6,14. This is especially
important in the low temperature regions as well as fuel rich and fuel lean areas, as these conditions are common during ignition phase. Figure1 and Figure 2 shows the comparison for the case of a stoichiometric mixture (ξ =0.062, equivalence ratio = 1.0). The blue regions are where ILDM points are found. The red regions indicate where only model points could be found.
As can be seen from Figure1 and Figure 2, in the case of stoichiometric mixture, the area of red region has decreased. (Areas showing no blue or red square mean that points are yet to be calculated in this region). Almost of the range required for a diesel engine simulation now shows existence of ILDM points. Similar enhancements are seen for lean and rich regions as well.
6 SOOT FORMATION
A semi-empirical model is used, where the formation and oxidation of soot is described via global steps. These are: 1. Particle formation via soot precursors viz. polycyclic aromatic hydrocarbons (PAH). PAH’s are formed via C3H3 , which can form the first
benzene ring. Concentration of C3H3 is obtained from ILDM table. 2. Particle growth,
which includes surface growth due to C2H2 and coagulation of existing particles. 3. Soot
oxidation, where oxidizing species are O, OH or O2. Concentration of O, OH and O2 are
also obtained from ILDM table.
The temporal rate of change of the soot volume fraction fv and soot concentration Cs is described using certain parameters that are calculated from the variables obtain from ILDM table.
7 COMPARISON OF PHYSICAL AND CHEMICAL TIME SCALES
Mean scalar dissipation and mean rate of strain was calculated in each computational cell to evaluate the relative ratio of the flow time scale to that of the chemical time scale. The chemical time scale is obtained from the ILDM computation. It is important to note that the maximum reaction progress variable concentration (CO2,eq and H2Oeq) used in
the unphysical states even if a non-equilibrium chemistry model is used. Effectively, the value of the CO2,eq and H2Oeq used for normalisation is scaled by the ratio of flow to
chemistry time scales. When the ratio of flow time to chemical time is smaller than 1 for a given cell, the pdf integration calculation using chemistry source terms from ILDM is skipped and a simpler one step reaction is used to calculate the chemistry and heat source terms. The Damkohler number therefore introduces a correlation between the pdf’s which are assumed to be independent. The results differ slightly when the inverse of the mean rate of strain is used to determine the flow time scale as against the scalar dissipation rate.
7.1 Scalar Dissipation
Local instantaneous scalar dissipation rate, which describes the rate of molecular mix-ing of fuel and oxidizer is considered as the most important parameter in non-premixed combustion15. The structure of the reaction zone is strongly coupled to the underlying
strain rate field through the influence of fluctuating strain on the scalar dissipation rate16
. Scalar Dissipation differs from energy dissipation because its less uniform in both space and time - a phenomenon called intermittancy; described as intense localized fluctuations of any quantity in a turbulent flow. This can result in localized extinction and reignition of combustion processes17.
8 RESULTS
The above ignition and chemistry models were implemented in KIVA III and were used to model a Caterpillar engine, for which experimental results were available5. The
engine specifications are listed in Table 8 Since the combustion chamber geometry and the six-hole injector configuration are symmetrical, the entire computational domain was divided into 6 equal sectors, and the computational domain actually simulated was one-sixth (60o) of the total chamber. The numerical mesh used contained about 35,000 cells
at BDC and about 12000 cells at TDC. Experimental data were available for 5 different injection timings. Figure 3 shows the behavior of the representative species (CO) during the ignition delay period for the case with SOI -7o ATDC. Its concentration remains near
zero during this period, and shows a sharp increase, indicating ignition.
Bore 137.16mm Stroke 165.1mm ConRod Length 263mm
Engine Speed 1600 rpm Number of Nozzle Orifices 6
Injection Timing -7, -4, -1, 2, 5o ATDC
Duration of Injection 19.75o
Fuel Injected 0.168 g/cycle Cylinder Wall Temperature 433K
Piston Wall Temperature 553K Head Temperature 523K Initial Gas Temperature 361±15K
Spray Temperature 341K
Initial Gas Composition (g/cm3)
O2 N2 CO2 H2O
4.6012× 10−4 1.5337× 10−3 2.8579× 10−6 1.2579× 10−6
Table 2: Caterpillar 3401 Engine Specifications
Literature reports that when more than one reaction contributes to heat release, the most reactive mixture fraction (defined loosely as the mixture fraction where ignition occurs), varies with time29. The mixture fraction of the cells where ignition is reported is
seen to vary with time and is in a range of few points on the leaner and richer side of the stoichiometric mixture fraction. It predominantly lies in the leaner side for late injection cases and rich side for early injection cases.
OH is one of the most important radical formed by the reaction of H atoms with molecular O2. 80% of the O2 is consumed via this route30 . Figure 5, a through f show
a representative flame propagation in terms of an iso-surface of the mass fraction of OH radical. One can see that in the early stages, the flame is on the outer edges of (see a) where the spray particles are (not shwon in figure). As the crank progresses, the flame spreads rapidly to engulf the entire space and licks the cylinder (see f ) head at crank angle of 8+. The iso-surface is colored by the temperature, where one can notice the highest temperatures on the leading edge.
9 CONCLUIONS
An ignition model and a chemistry model based on detailed chemistry were incorpo-rated in a modified version of KIVA III and used to simulate a direct injection Caterpillar engine. An ILDM table with improved numerical schemes was created. A time-scaled approach was used to calculated the chemistry in the ILDM cells. A good agreement is seen between experimental and simulated pressure, temperature as well as that of the heat release rate. The end of run soot concentrations also compare well.
REFERENCES
[1] C. Correa, H. Niemann, B. Schramm and J. Warnatz, Use of ILDM Reduced Chem-istry in Direct Injection Diesel Engines, Thermo- and Fluid-Dynamic Processes in Diesel Engines, Selected Papers from the THIESEL 2000 Conference, Valen-cia, Spain, Sept. 13-15, p353-362, (2002). Thermo- and Fluid-Dynamic Processes in Diesel Engines, Selected Papers from the THIESEL 2000 Conference, Valencia, Spain, Sept. 13-15, pp 353-362, (2002).
[2] K.J. Syed, J.B. Moss and C.D. Stewart, Modeling soot formation and thermal radi-ation in buoyant turbulent diffusion flames, Twenty-fourth (Intl.) Symposium of the Combustion Institute 24, p 1533, (1990).
[3] A.A. Amsden, P.J O’Rourke and T.D. Butler, KIVA-II: A Computer Program for Chemically Reactive Flows with Sprays, Los Alamos National Laboratory, LA-11560-MS, (1989).
[4] J. Xin, L. Ricart and R.D. Reitz, Computer Modeling of Diesel Spray Atomization and Combustion, Combust. Sci. Tech. 137, 171-194 (1998).
[5] R.D. Reitz, Personal Communications of Chrys Correa, Engine Research Center, University of Wisconsin, Personal Communications of Chrys Correa (1999).
[6] B.S. Schramm, Automatische Reduktion chemischer Reaktionsmechanismen am Beispiel der Oxidation von h¨oheren Kohlenwasserstoffen und deren Verwendung in reaktiven Str¨omungen, Ph.D Thesis, IWR, University of Heidelberg, (2003).
[7] J. Warnatz, U. Maas and R.W. Dibble, Combustion, Springer Verlag, Berlin Heidel-berg, 2nd edition, (1999).
[8] U. Maas and S.B. Pope, Simplifying chemical kinetics: Intrinsic low-dimensional manifolds in composition space, Comb. Flame 88, 239, (1992).
[10] M. Nehse, J. Warnatz and C. Chevalier, In Proc. Combust. Inst., Vol. 26, 173, The Combustion Institue, Pittsburgh, PA. (1996).
[11] C. Correa, H. Niemann, B. Schramm and J. Warnatz, Reaction mechanism reduction for higher hydrocarbons by the ILDM method, In Proc. of Comb. Inst., Vol. 28, pp 607–1614, Pittsburg, PA. (2000).
[12] S.S. Girimaji. Assumed β-pdf model for turbulent mixing: Validation and extension to multiple scalar mixing. Combust. Sci. Technol. 78, p177, (1991).
[13] Z. Han and R.D. Reitz, Turbulence modeling of internal combustion engines using RNG k-ǫ models, Combust. Sci. and Tech. 106, pp 267–295 (1995).
[14] G.N Coleman and N.N. Mansour, Modeling the rapid spherical compression of isotropic turbulence, Phys. Fluids A 3, p 2225, (1991).
[15] R.B. Bird, W.E. Stewart and E.N. Lightfoot (Ed.), Transport Phenomenon, John Wiley & Sons, New York, (1960).
[16] R.S. Barlow and J.H. Frank, Effects of Turbulence on Species Mass Fractions in Methan/Air Jet Flames, In Twenty-Seventh Symposium (Intl.) on Combustion, Vol. 28, pp 1087-1095, The Combustion Institute,Pittsburg, PA, (1998).
[17] S.B. Pope, Monte Carlo calculations of premixed turbulent flames, Eighteenth Sym-posium (Intl.) on Combustion 18, pp 1001-1010, (1981).
[18] S.B. Pope, PDF methods for turbulent reactive flows, Prog. Energy Combust. Sci. 11, p119, (1985).
[19] S.B. Pope. Turbulent Flows, Cambridge University Press, Cambridge, UK, (2000). [20] F.C. Lockwood and A.S. Naguib, Aspects of combustion modeling in engineering
turbulent diffusion flames, J. of the Institute of Fuel 49, pp 218-223, (1976).
[21] J. Janicka and W. Kollmann, A two-variables formalism for the treatment of chemical reactions in turbulent hydrogen-air diffusion flames. Seventeenth Symposium (Intl.) on Combustion 17, pp421–430 (1981).
[22] E.E. Khalil, D.B. Spalding and J.H. Whitelaw, Calculation of local flow properties in two-dimensional furnaces, Int. J. Heat Mass Transfer 18, pp 775–791, (1975). [23] U. Maas, Automatische Reduktion von Reaktionsmechanismen zur Simulation
[24] U. Maas and S.B. Pope, Implementation of simplified chemical kinetics based on in-trinsic low-dimensional manifolds, In Proc. Combust. Inst., Vol. 24, 103, Pittsburgh, PA, USA, (1993).
[25] H. Niemann, Automatische Reduktion chemischer Reaktionsmechanismen, Ph.D The-sis, IWR, University of Heidelberg, (2000).
[26] H. Pitsch, C. M. Cha and S. Fedotov, Interacting flamelet model for non-premixed turbulent combustion with local extinction and re-ignition, Center for Turbulence Research, Stanford University, Annual Research Briefs, pp 65–77, (2001).
[27] N. Peters, Laminar Diffusion Flamelet Models in Non-Premixed Turbulent Combus-tion, Prog. Energy Combust. Sci., 10, p 319, (1984).
[28] Author Unknown, Creation of the whorled, National Energy Research Scientific Com-puting Center, Annual Report, http://www.nersc.gov/news/annual reports/ annrep04/ annrep04.pdf, pp 14-18, (2004).
[29] T. Ishiyama, K. Miwa and O. Horikoshi, A study on ingition process of diesel en-gines,COMODIA 94, Tokyo, Japan, pp 337-342, (1994).
Temperature (K)
C
O
2
M
a
s
s
F
ra
c
ti
o
n
1500
0
2000
2500
3000
0.02
0.04
0.06
0.08
0.1
Figure 1: ILDM points for stoichiometric fuel conditions - New Case with 1 rpv (Blue area indicates good points). 1500 2000 2500 3000 Temperature (K) 0.00 0.04 0.08 0.12 0.16 C O2 m a s s fr a c ti o n
Figure 3: Rise in CO concentration leading to Ignition for the case of -4 degree SOI in a Caterpillar Engine. Tp 450 430 410 390 370 350 Y X Z
Y X Z T 2400 2120 1840 1560 1280 1000 (a) Crank=-4.2 Y X Z T 2400 2120 1840 1560 1280 1000 (b) Crank=-3.9 Y X Z T 2400 2120 1840 1560 1280 1000 (c) Crank=-2.9 Y X Z T 2400 2120 1840 1560 1280 1000 (d) Crank=-1.9 Y X Z T 2400 2120 1840 1560 1280 1000 (e) Crank=+2.1 Y X Z T 2400 2120 1840 1560 1280 1000 (f) Crank=+8.1
Crank Angle (Degrees) P re ssu re (b a rs) -100 -50 0 50 100 20 40 60 80 Experiment Simulation
Figure 6: Temperature Profile for SOI=-1 ATDC with CO2 as Progress Variable.
Crank Angle (degrees)
T
e
m
p
e
ra
tu
re
(K
)
-100 -50 0 50 100 400 600 800 1000 1200 1400 1600 1800 Experiment SimulationCrank Angle (Degrees) H e a t R e le a se R a te (J/ C ra n k A n g le ) -50 0 50 100 0 50 100 150 200 250 300 350 400 Experiment Simulation
Figure 8: Heat Release Rate for SOI=-1 ATDC with CO2as Progress Variable.
SOI (Degree ATDC)
S o o t C o n ce n tr a ti o n (m o ls/ cc) -5 0 5 0 5E-11 1E-10 1.5E-10 2E-10 Experiment Simultaion Old Simulation
