Table Of ContentSixteenth International Multidimensional Engine Modeling User’s Group Meeting
at the SAE Congress, April 2, 2006, Detroit, Michigan
A G-equation Combustion Model Incorporating Detailed Chemical
Kinetics for PFI/DI SI Engine Simulations
Long Liang∗, Rolf D. Reitz
Engine Research Center, University of Wisconsin-Madison, Madison, WI 53706,USA
Jianwen Yi, Claudia O. Iyer
Ford Research and Advanced Engineering, Dearborn, MI 48124,USA
A G-equation-based combustion model incorporating detailed chemical kinetics has been developed and
implementedinKIVA-3VforSpark-Ignition(SI)enginesimulationsforbetter predictionsofflamepropa-
gation and pollutant formation. A progressvariable concept is introduced into the turbulent flame speed
correlation to account for the laminar to turbulent evolution of the spark kernel flame. The flame front
in the spark kernel stage is tracked using the Discrete Particle Ignition Kernel (DPIK) model. In the
G-equation model, it is assumed that after the flame front has passed, the mixture within the mean
flame brushtends to localequilibrium. The subgrid-scaleburnt/unburntvolumes of the flame-containing
cells are tracked for the primary heat release calculation. An iso-octane kinetic mechanism coupled with
a reduced NOX mechanism is used to describe the chemical processes in the post-flame region and the
potential heat release from the end gas. The integrated model was used to simulate the combustion pro-
cess in a Fordfour-valve single-cylinder SI engine, which is equipped with both Port-Fuel-Injection(PFI)
and Direct-Injection (DI) fuel systems. For both PFI and DI operational modes, good agreement with
experimentalin-cylinder pressure,heatreleaseratesandengine-outNOX wasobtainedfor differentspark
timings and internal residual levels.
1 Introduction mental engine combustion process and to further im-
prove the versatility of multidimensional models, at-
The in-cylinder turbulent combustion in SI engines tention is being given to models incorporating com-
is a complicated aero-thermo-chemical process es- prehensive elementary chemical kinetic mechanisms.
pecially due to the turbulence and chemistry in- The objective of the current work is to incorporate
teractions on tremendously different time-scale and detailed chemical kinetics into the G-equation-based
length-scale levels. In this paper, we present a turbulent combustion model which was implemented
G-equation-based flamelet combustion model incor- into the KIVA-3V code by Tan et al. [5] [7]. Specif-
porating detailed chemical kinetics for both Port- ically, detailed fuel oxidation mechanisms coupled
Fuel-Injection(PFI)andDirect-Injection(DI)Spark- with a reduced NOX mechanism are applied behind
Ignition (SI) engine combustion simulations. The the mean flame front for modeling post flame com-
level set method is a powerful tool for describing in- bustion and NOX formation. The chemical kinetic
terface evolution. With its application to combus- mechanismsarealsoappliedinfrontoftheflamefront
tion,Williams [1]firstsuggestedatransportequation for potential capability of predicting the compression
of a non-reactive scalar, G, for laminar flame propa- autoignition of the end-gas. In the course of coupling
gation. Peters [2] [3] subsequently extended this ap- detailed chemistry with the G-equation combustion
proach to the turbulent flame regime. The turbulent model for the primary heat release calculation within
G-equation concept has been successfully applied to the flame front, it was required to revisit and im-
SIenginecombustionsimulationsbyDekenaetal.[4], prove laminar and turbulent flame speed correlations
Tan [5] and Ewald et al. [6]. for better description of the turbulent flame propa-
In recent years, to better understand the funda- gation process.
∗Corresponding author. Email address: [email protected]
1
2 Model Formulation the present model to account for the laminar to tur-
bulent evolution of the spark kernel flame. This term
2.1 G-equation description of turbulent
can be obtained by solving Eq. (2) assuming a uni-
flame propagation form turbulenceprofile [3]. Physically, the additional
exponential term can be interpreted as a progress
In the flamelet modeling theory of premixed turbu-
variable which accounts for the increasing disturbing
lentcombustionbyPeters[3],tworegimesofpractical
effect of the surrounding eddies on the flame front
interest, i.e., the corrugated flamelet regime and the
surface as the ignition kernel grows from the laminar
thinreactionzoneregime,weredescribedbythesame
flame stage into the fully developed turbulent stage.
group of level set equations, including the transport
equationsfortheFavreaveraged,G(cid:1),anditsvariance, In the present study, Cm2 in Eq. (3) is selected as a
(cid:2) tunable model constant for different engines consid-
G(cid:1)(cid:1)2, and a model equation for the turbulent/laminar
ering uncertainties due to other sub-models and/or
flame surface area ratio, which in turn gives an alge-
mesh resolution. However, for the same Ford engine
braic solution for steady-state planar turbulentflame
speed, S0. The set of equations used in the current studied in this paper, Cm2 is fixed over all operating
T
conditions.
implementation is [5] [7]:
Laminar flame speed S0 is one of the key scaling
L
∂G(cid:1) ρ factors in Eq. (3). Metghalchi et al. [9] suggested a
∂t +(v(cid:3)f−vvertex)·(cid:1)G(cid:1) = ρuST0|(cid:1)G(cid:1)|−DT(cid:1)κ|(cid:1)G(cid:1)| (1) correlation for SL0 as a function of equivalence ratio
φ, temperature and pressure of the unburnt mixture
given by:
(cid:2)
∂G(cid:1)(cid:1)2 (cid:2) T p
∂t +(v(cid:3)f −vvertex)·(cid:1)G(cid:1)(cid:1)2 SL0 = SL0,ref(T u )α(p )β(1−2.1·Ydil) (4)
u,ref ref
ρ (cid:1)(cid:4)
= (cid:1)(cid:2)·( ρuDT(cid:1)(cid:2)G(cid:2)(cid:1)(cid:1)2)+2DT((cid:1)G(cid:1))2−cs(cid:1)kG(cid:2)(cid:1)(cid:1)2 (2) where the subscript ref means the reference condi-
tion of 298 K and 1 atm. Ydil is the mass fraction of
diluent. The S0 and fuel-type independent expo-
L,ref
S0 (cid:4) nents α and β are correlated as functions of φ as:
ST0 = 1+ 1−exp(−Cm2·t/τ)
L (cid:5) (cid:6) (cid:7) SL0,ref = BM +B2(φ−φM)2 (5)
a b2 l a b2 l u(cid:1)l
· − 24b13lF + ( 24b13lF)2+a4b23SL0lF (3) βα==−20..1186−+00..82(2φ(φ−−1)1) ((76))
where (cid:1)(cid:2) is the tangential gradient operator, (cid:3)vf is Since Eq. (5) is invalid for very lean and very rich
the fluid velocity vector, vvertex is the vertex moving mixtures (by predicting negative flame speed val-
velocity. DT is the turbulent diffusivity, and a4,b1,b3 ues), which is unacceptable in modelingDI operating
and cs are constants from the turbulence model or conditions, a formula proposed by Gu¨lder [10] was
experiment (cf. [3]). (cid:1)(cid:4) and (cid:1)k are the Favre mean tur- adopted in this study:
bulence kinetic energy and its dissipation rate from S0 = ωφη ·exp(−ξ(φ−σ)2) (8)
theRNGk-(cid:4) model[8]. u(cid:1) is theturbulenceintensity. L,ref
SL0 is the unstretched laminar flame speed. l and lF Due to the relatively coarse mesh resolution in
are the turbulence integral length scale and laminar engine simulations, the growth of the ignition kernel
flame thickness, respectively. κ(cid:1) is the mean flame is tracked by using the DPIK model [5], where the
front curvature. These equations together with the flame front position is marked by Lagrangian parti-
Reynolds averaged Navier-Stokes equations and the cles. The kernel growth rate is:
turbulencemodelingequationsformacompletesetto dr ρ
describepremixedturbulentflamefrontpropagation. dtk = ρu(Splasma+ST) (9)
k
One significant advantage of the G-equation for-
where rk is the kernel radius, ρu and ρk are the local
mulation of turbulent premixed flames is the absence
unburnt gas density and the gas density inside the
of chemistry source terms in the transportequations.
Asaconsequence,theturbulentflamespeedS0 plays kernel, respectively. The plasma velocity Splasma is
T
given as [5]:
a crucial role as a predetermined input. Compared
wexipthontehneticaolrlyreilnatciroenasfionrgStTe0rmderisivaedddbeydPinetEerqs.[2(3],)ainn Splasma = 4πrk2[ρu(uQk˙s−pkh·uη)ef+fp·ρu/ρk] (10)
2
where Q˙spk is the electrical energy discharge rate,
ηeff is the electrical energy transfer efficiency due
to heat loss to the spark plug. hu is the enthalpy
of unburnt mixture, uk is the internal energy of the
mixture inside the kernel. To account for turbulent
strain and curvature effects on the kernel flame, the
unstretched laminar flame speed S0 in Eq. (4) was
L
multiplied by a stretch factor I0, which takes the fol-
Figure 1: Numerical descriptions of the turbulent flame
lowing form according to Herweg et al. [11]:
structure and the flame containing cells.
I0 = 1−(lF/15l)1/2(u(cid:1)/SL0)3/2−2·(lF/rk)(ρu/ρk) As shown in Fig. 1, in this method, in order to
(11) predict Yi,u for all the species in the flame-containing
cells, the sub-grid scale unburnt/burnt volumes par-
where the second and third terms on the right hand
titioned by the mean flame front are tracked based
side represent the contributions due to turbulent
on geometrical information. As the mean flame front
strain and due to the geometrical curvature of the
sweepsforward,themixturebehindthesweepingvol-
kernel, respectively. Note that the mean curvature
ume tends to local equilibrium following a constant
effects are also considered in theG-equation combus-
tionmodelbythelasttermofEq. (1). Thetransition pressure, constant enthalpy process. Yi,u can be cal-
culated as follows:
from the kernel model to the turbulent G-equation
(1) Determine the equilibrium species mass frac-
combustion model follows the same criterion as the
one used in the previous work by Tan and Reitz [5], tions Yi,b and the equilibrium flame temperature. In
this study, an element potential method-based code
i.e., the transition is controlled by a comparison of
by Pope [13] was used for the chemical equilibrium
the kernel radius with a critical size which is pro-
calculation.
portional to the locally averaged turbulence integral
(2) Calculate the burntgas density and the burnt
length scale,viz.,
species densities based on the equation of state and
rk ≥ Cm1·l = Cm1·0.16k3/2/(cid:4) (12) the Yi,b from step (1):
where Cm1 is a model constant. ρi,b = Yi,b·(pMWmix,b)/((cid:2)Tb) (14)
Although elements of the above G-equation de- where MWmix,b is the average molecular mass of the
scription were originally developed for premixed burnt mixture, and (cid:2) is the universal gas constant.
flames, it is also successfully applied to partially pre- (3) Calculate the unburnt species densities ρi,u
mixed flames in the DI operating mode in this study.
based on species mass conservation:
2.2 Primary heat release within the tur- ρi,u = (ρiVi4−ρi,bVb)/Vu (15)
bulent flame brush
(4) Finally, determine the unburnt species mass
In the present implementation of the G-equation fractions: (cid:8)
model, it is assumed that after the flame front has Yi,u = ρi,u/ ρi,u (16)
i
passed, the mixture within the turbulentflame brush
In KIVA, the heat release rate due to the chem-
tends to the local and instantaneous thermodynamic
istrysourcetermisdirectlyrelatedtothespeciescon-
equilibrium. Based on this assumption, an updated
version rate [12]. It needs to be noted that the four
method is suggested to calculate the species conver-
sion rate and the associated primary heat release at species associated with the NOX formation mecha-
the flame front, viz., nism, i.e., NO, NO2, N, and N2O, are excluded from
the equilibrium calculation due to their relatively
dρi/dt = ρu(Yi,u−Yi,b)ST0Af,i4/Vi4 (13) short residence time within the flame front, and the
relative slow rate of the NOX chemical reactions.
wherei4isthecomputationalcellindexinKIVA[12].
Yi,u and Yi,b are the mass fractions of species i with 2.3 Post-flame heat release and pollutant
respect to the unburnt and burnt mixtures, respec-
formation
tively. Af is the mean flame front area and V is the
cell volume. In this study, the computational cells ahead and be-
hindthepropagating flamefrontaremodeledas Well
3
StirredReactors(WSR).AdetailedPRFchemicalki- constants Cm1 = 2.0 and Cm2 = 1.0 were held fixed
netic mechanism was applied to account for the fur- (see Eqs. (12) and (3)).
ther oxidation of CO and other intermediate species, The predicted in-cylinder pressure traces match
such as small hydrocarbon molecules and the species the measured data well in terms of peak pressure
intheH2-O2 system. Toconsidertheeffectsofturbu- and combustion phasing in the PFI operating mode
lent mixing, the reaction rates could be adjusted by with spark timing sweeps, as shown in Fig. 2. Fig-
considering the eddy turnover time as an additional ure 3 shows the evolution of the mean turbulent
timescale, and by combining this timescale with the flame front (G(cid:1) = 0 iso-surface) in the PFI mode with
kinetic timescale. However, this was not done in the spark timing = -40 ◦ATDC and engine speed 1500
present study, i.e., only kinetic rates were used. A rev/min). As seen, in the PFI mode, the flame prop-
nine-reaction reduced NOX mechanism was coupled agates throughout the whole cylinder reaching the
with the hydrocarbon oxidation mechanism for pre- wall at about 20 ◦ATDC.
dicting the formation of NO and NO2 [14].
3 Results and Discussion
TheenginestudiedisaFordfour-valvesingle-cylinder
SI gasoline engine which features a pentroof combus-
tion chamber and a converging-shaped piston bowl.
The engine is equipped with both PFI and DI fuel
systems. The DI system employs the wide-spacing
arrangement, with a centrally mounted spark plug
and an intake-side-mounted swirl-type injector. The
test data used in this work was reported by Mun˜oz
et al. [15]. The specifications of the engine and the
modeled operating conditions are listed in Table 1.
Figure 2: Measured (EXPT) and predicted (SIMU) in-
cylinder pressure in PFI mode.
Table 1: Ford engine specifications and operating condi-
tions [15].
Bore / Stroke 89 mm / 79.5 mm
Compression Ratio 12
Engine Speed 1500 rev/min
PFI Mode
(a) (b)
◦
Spark Timings ( ATDC) -44, -40, -36, -32
Internal Residual 28%
MAP 65kPa
DI Mode
◦
Spark Timings ( ATDC) -32, -28, -24, -20
Internal Residual 6%
MAP 75kPa
(c) (d)
◦
End-of-Injection ( ATDC) -72
Figure 3: Evolution of the mean turbulent flame front
Thecomputationalmeshcontainsaround170,000 (shaded G(cid:1) = 0 iso-surface) in PFI mode (Spark timing =
cells, including the intake and exhaust manifolds -40 ◦ATDC, Engine speed=1500 rev/min).
and the cylinder. A reduced 25-species, 51-reaction Compared to PFI cases, it is more challenging to
iso-octane mechanism [16] including the NOX reac- accurately predictthepressureevolution andheatre-
tions was used to model the post-flame chemistry lease in the DI mode. The equivalence ratio of the
and the low temperature chemistry in the end gas. stratified charged mixture varies from very rich to
CHEMKINIIwas usedto solvethedetailed chemical very lean. Therefore the laminar flame speed corre-
kinetic equations. In all simulated cases, the model lation needs to bereliable over a widerange of equiv-
4
alence ratios, and also as a function of temperature the dark lines denote stoichiometric conditions, the
and pressure of the unburnt mixture. As seen in Fig. bright lines represent the flame front surfaces. Ac-
4, the pressureand heat release rate predicted by the cording to the simulation, most of the NO is formed
present model agree with the measured data reason- around the stoichiometric lines, while CO is mainly
ably well for all spark timings. generated within the fuel rich region, as expected.
Figure 6 shows the comparison of normalized
engine-out NOXin both the PFI and DI modes. Al-
though therearediscrepancies inabsolute values, the
general trends as functions of spark timing are well
captured.
(a) (b)
(a) (b)
Figure 6: Measured and predicted engine-out NOX. (a)
PFI mode; (b) DI mode.
(c) (d)
Figure 4: Comparisons of in-cylinder pressure and heat
release rate in the DI mode.
(a) (b)
Figure 7: Calculated evolution of kernel radius. (a) PFI
mode; (b) DI mode.
During the calibration of the present ignition and
combustion models, it was found that the stretch ef-
fects on the kernel flame, as described by Eq. (11),
Figure 5: In-cylinderNOandCOmassfractioncontours played animportantroleinthepredictionofcombus-
at 20 ◦ATDC (DI, spark timing = -32 ◦ATDC). tion phasing. ReferringtoEq. (11), thestretch effect
Figures5(a)and(b)showthespatialdistribution due to flame curvature (the last term on the right
ofNOandCOspeciesmassfractionsat20◦ATDCfor hand side) decays very quickly with the growth of rk,
◦
the DI case with spark timing = -32 ATDC, where while the strain effect dueto strong turbulence and a
5
thicklaminarflamestructure(thesecondtermonthe [6] J. Ewald and N. Peters. A level set based
rightside)mightsignificantlydelay thekernelflame’s flamelet model for the prediction of combus-
fastpropagation,asinthePFIcasewithsparktiming tion in spark ignition engines. 15th Interna-
◦
= -44 ATDC shown in Fig. 7 (a). (Here for the first tional Multidimensional Engine Modeling User’s
◦
10 CA from spark timing, the kernel flame was sus- Group Meeting, Detroit, MI, 2005.
tained mainly by the plasma velocity). Comparing
[7] Z. Tan, S.-C. Kong, and R. D. Reitz. Modeling
Figs. 7 (a) and (b), the stretch effects in the DI cases
premixed and direct injection si engine combus-
aregenerally less significantcompared tothose inthe
tion using the G-equation model. JSAE Paper
PFI cases, mainly because of thelower internal resid-
2003-01-1843, 2003.
ual fractions and correspondingly higher unstretched
laminar flame speeds in the DI cases.
[8] Z. Han and R. D. Reitz. Turbulencemodelingof
internal combustion engines using rng k-e mod-
els. Comb. Sci. Tech., 106:267–295, 1995.
4 Conclusions
[9] M. Metghalchi and J. C. Keck. Burning veloc-
AG-equation-basedturbulentcombustionmodelwas
ities of mixtures of air with methanol, isooc-
coupled with detailed chemistry, so that the post-
tane, and indolene at high pressures and tem-
flame heat release and the low temperature heat
peratures. Combustion and Flame, 48:191–210,
release in the end gas could be modeled based
1982.
on detailed chemical kinetics. A sub-grid scale
burnt/unburnt volume-based method was developed [10] O. L. Gu¨lder. Correlations of laminar combus-
for calculating the primary heat release and species tion data for alternative s.i. engine fuels. SAE
conversion within the turbulent flame brush. Paper 841000, 1984.
Both the laminar and turbulent flame speed cor-
[11] R. Herweg and R. R. Maly. A fundamental
relationswereupdatedformoreaccuratedescriptions
model for flame kernel formation in s.i. engines.
of kernel flame evolution and turbulent flame propa-
SAE Paper 922243, 1992.
gation in stratified mixtures.
A Ford SI Gasoline engine was modeled in both [12] A. A. Amsden. A block-structured kiva pro-
PFI and DI operating modes with spark timing gram for engines with vertical or canted valves.
sweeps. Good agreement with the measured pressure TechnicalReportLA-13313-MS,LosAlamosNa-
traces and engine-out NOXwas achieved in all cases tional Lab, 1997.
with fixed values of the model constants.
[13] S. B. Pope. Ceq: A fortran li-
brary to compute equilibrium composi-
References
tions using gibbs function continuation.
http://eccentric.mae.cornell.edu/∼pope/CEQ,
[1] F. A. Williams. Turbulent Combustion. SIAM,
2003.
1985.
[14] S.-C. Kong, Y. Sun, and R. D. Reitz. Modeling
[2] N. Peters. The turbulent burning velocity for diesel spray flame lift-off, sooting tendency and
large scale and small scale turbulence. J. Fluid noxemissionsusingdetailedchemistrywithphe-
Mech., 384:107–132, 1999. nomenologicalsootmodel. Proceedings of ASME
ICES2005, 2005.
[3] N. Peters. Turbulent Combustion. Cambridge
University Press, 2000. [15] R. H. Mun˜oz, Z. Han, B. A. VanDerWege,
and J. Yi. Effect of compression ratio on
[4] M. Dekena and N. Peters. Combustion model- stratified-charge direct-injection gasoline com-
ingwith theG-equation. Oil & Gas Science and bustion. SAE Paper 2005-01-0100, 2005.
Technology-Rev., IFP, 54(2):265–270, 1999.
[16] S. Tanaka, F. Ayala, and J. C. Keck. A reduced
[5] Z. Tan and R. D. Reitz. Modeling ignition chemical kinetic model for hcci combustion of
and combustion in spark-ignition engines using primary reference fuels in a rapid compression
a level set method. SAE Paper 2003-01-0722, machine. Combustion and Flame, 133:467–481,
2003. 2003.
6
