- Email: [email protected]
Towards large eddy simulation of combustion in spark ignition engines
Proceedings of the
Proceedings of the Combustion Institute 31 (2007) 3059–3066
Combustion Institute www.elsevier.com/locate/proci
Towards large eddy simulation of combustion in spark ignition engines S. Richard a, O. Colin a,*, O. Vermorel a, A. Benkenida a, C. Angelberger a, D. Veynante b b
a IFP, 1 et 4 avenue de Bois-Pre´au, 92852 Rueil-Malmaison Cedex, France Laboratoire EM2C, CNRS Ecole Centrale Paris, 92295 Chaˆtenay-Malabry Cedex, France
Abstract Internal combustion engine simulations are commonly performed using the RANS (Reynolds averaged Navier–Stokes) approach. It gives a correct estimates of global quantities but is by nature not adapted to describe phenomena strongly linked to cyclic variations. On the other hand, large eddy simulation (LES) is a promising technique to determine successive engine cycles. This work demonstrates the feasibility of LES engine cycles simulation by using a flame surface density (FSD) approach. This approach, presented in a first section, combines an Eulerian spark ignition model derived from the RANS AKTIM model [J.M. Duclos, O. Colin, Arc and Kernel Tracking Ignition Model for 3D SI Engines Calculations, Comodia, Nagoya, Japan, 2001, pp. 343–350] and a Coherent Flame Model (CFM) [S. Candel, T. Poinsot, Combust. Sci. Tech. 70 (1990) 1–15; O. Colin, A. Benkenida, C. Angelberger, Oil & Gas Sci. Techn.—Rev. IFP 58 (1) (2003) 47–32] describing the flame propagation. The CFM model, commonly used in RANS simulations, is here formulated in a LES context. In a second part, the whole ignition-combustion model is validated against an experiment relative to the turbulent ignition and flame propagation of a stoichiometric propane-air mixture [B. Renou, A. Boukhalfa, Combust. Sci. Tech. 162 (2001) 347–371]. Finally, LES engine cycles simulations are performed on a real engine configuration. First, the sensitivity of the model to the ^ is examined, showing a weak dependence of the modelling approach to D. ^ LES combustion filter size D Then results are compared to those obtained with the algebraic model for the FSD proposed by Boger et al. [M. Boger, D. Veynante, H. Boughanem, A. Trouve´, Proc. Combust. Inst. 27 (1998) 917–925] and the need for non-equilibrium combustion models is demonstrated. 2006 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: LES; Combustion; CFM; Spark ignition; Engine cycle simulation
1. Introduction
*
Corresponding author. Fax: +33 0 147527068. E-mail addresses: [email protected] (S. Richard), [email protected] (O. Colin), [email protected] (D. Veynante).
The RANS (Reynolds averaged Navier– Stokes) approach provides today correct predictions of mean quantities such as pressure or pollutant emissions in internal combustion engines but does not describe unsteady phenomena such as cyclic variations or engine transients. On the other
1540-7489/$ - see front matter 2006 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.proci.2006.07.086
3060
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
hand, large eddy simulation (LES), where resolved flow motions are explicitly computed, is well adapted to unsteady problems. However, the effects of sub-grid-scale (sgs) motions should be modelled. Previous premixed combustion modelling approaches in LES have been mainly based on the artificial thickening of a flame [6] or on the so-called G-equation [7]. The laminar flamelet approach appears as an interesting alternative for piston engines, according to the combustion diagram from Borghi [8] and considering typical Damko¨hler and Karlovitz numbers of about Da 10–100 and Ka 0.01–1. Algebraic models have already been proposed for the flame surface density (FSD) in LES applications [5]. Balance equations for the wrinkling factor [9] or the FSD [10] have also been investigated. In piston engines, where the flame is initially laminar at the spark plug, then grows under the action of thermal expansion, curvature and turbulence and is finally quenched when interacting with walls, the correct description of the wrinkling evolution is crucial, giving an advantage to the balance equation formalism. The present work is devoted to the presentation and the validation of a FSD equation and a spark ignition model dedicated to piston engine applications. These developments have been implemented in the AVBP code [11] specifically dedicated to LES combustion calculations and have already been successfully tested on simple test cases [12]. The first section summarises the modelling approach. In a second part, the whole ignition–combustion model is validated against an experiment relative to the ignition and flame propagation [4]. Finally, simulations are performed on a real engine configuration. The sensitivity of the model to the LES filter size is first examined. Results are then compared to those obtained with the algebraic closure of FSD proposed by Boger et al. [5]. 2. The CFM-LES combustion model Assuming a flamelet regime, a single step chemistry and a unit Lewis number, the thermochemistry of the reacting flow may be described by the progress variable c of the reaction (c = 0 in fresh gases and c = 1 in fully burnt gases) according to P
oq~c þ r ðq~u~cÞ þ r ðquc q~u~cÞ ¼ qS d jrcj; ot ð1Þ
where Sd is the local propagation speed. Q denotes ~ a mass weighted filtered a filtered quantity and Q ~ ¼ qQ). The filtered FSD is defined as quantity (qQ R ¼ jrcj [5] and introducing the weighted flame
surface average as hQis ¼ QR=R, the source term in Eq. (1) is written as qS d jrcj ¼ hqS d is R. In practice, ÆqSdæs is approximated by quSl [13] where qu is the unburned gases density and Sl the laminar flame speed. Equation (1) contains two unclosed terms: an unresolved transport term and the filtered FSD. In the CFM-LES approach, the filtered FSD is computed through a transport equation. An exact unclosed equation for the local FSD was first proposed by Pope [14] or Candel and Poinsot [2] and filtered, leading to the following equation [10–15]: oR þ r huis R ¼ hr u nn : ruis R ot r ðhS d nis RÞ þ hS d r nis R; ð2Þ where n is the local normal to the flame front pointing towards the fresh gases. The first term on the right-hand side corresponds to the straining of the flame by the flow. The second and third terms represent, respectively the planar propagation of the flame at speed Sd and the propagation of a curved front. The unresolved transport term of the ~c equation is commonly closed under a gradient assumption and is written as a turbulent diffusion. Following this type of formulation, this term vanishes when turbulence decreases. However, Rymer [16] showed that in a laminar case the unresolved transport is non-zero and corresponds to nongradient transport (NGT) which can be written as r ðquc q~ u~cÞ ¼ qu S l r ððc ~cÞNÞ, where N ¼ rc=jrcj is the normal to the iso-surface of the filtered progress variable. Rymer then proposed the following general formulation: r ðquc q~u~cÞ ¼ qu S l r ððc ~cÞNÞ r
qtt r~c ; Sct ð3Þ
where tt is the turbulent viscosity and Sct a turbulent Schmidt number. Unfortunately, NGT terms may lead to numerical instabilities because they are anti-diffusive and are difficult to handle in practical applications. Consequently, this term is usually omitted and the equation becomes oq~c qtt þ r ðq~ u~cÞ ¼ r r~c þ qu S l R~c ot Sct
ð4Þ
with R~c defined as: R~c ¼ R þ r ððc ~cÞNÞ. A phenomenological transport equation for R~c is proposed here, assuming that this equation satisfies the same main properties as Eq. (2), i.e., it should ensure: (i) the laminar propagation of the flame in multidimensional cases, (ii) a reasonable behaviour—in particular correct turbulent flame speeds—in case of high turbulence levels and (iii)
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
the control of the flame brush thickness. The CFM-LES, unlike RANS formulations, is a spatially filtered equation and it takes into account the effects of the resolved flow on the flame front. For this purpose the FSD transport equation terms are split into resolved and unresolved parts oR~c þ T res þ T sgs ¼ S res þ S sgs þ C res þ C sgs þ P ; ot ð5Þ where Tres, Sres, Cres and P are, respectively, the transport, strain, curvature and propagation terms due to resolved flow motions, and Tsgs, Ssgs and Csgs are, respectively, the unresolved transport, strain and curvature terms. 2.1. Modelling of the flame propagation and curvature The propagation and resolved curvature terms are physically linked and ensure the laminar flame propagation, when the sgs turbulence is low. The proposed formulation of these terms is based on the normal to the iso-surface of the filtered progress variable and reads: P ¼ r ðS d NR~c Þ and C res ¼ S d ðr NÞR~c . Sd is defined by S d ¼ ð1 þ s~cÞS l , where s = qu/qb 1 is the thermal expansion rate across the flame front and qb is the burned gases density. DNS simulations [13] have shown that the sgs curvature acts as a sink term towards burned gases and as a source term towards fresh gases. Hence, Rymer proposed to express sgs curvature effects as: C sgs ¼ bS l ðc cÞ=ðcð1 cÞÞ ðR jrcjÞR, where b = 4/3 and c* = 0.5 are modelling constants. This formulation is retained here, corrected for the R~c equation: C sgs ¼ bS l ðc cÞ=ðcð1 cÞÞðR~c R~lam c is c ÞR~c , where estimated from the usual BML expression as: c ¼ ~cq=qb and R~lam ¼ jr~cj þ ðc ~cÞr N. c 2.2. Modelling of the flame strain rate The resolved strain is S res ¼ ðr ~u NN : r~uÞR~c and corresponds to the flame strain rate due to the resolved structures. The turbulent strain rate in RANS is classically modeled as at = C(u 0 /Sl, lt/dl)u 0 /lt, where C is an efficiency function taking into account the ability of all vortices to wrinkle the flame, often defined by the ITNFS function of Meneveau and Poinsot [17]. In LES, Colin et al. [6] showed that resolved eddies smaller than the flame thickness defined by ^ ¼ nres Dx (nres 5–10 for finite volume solvers) D are not able to wrinkle the flame front. Consequent^ (and not Dx like ly the combustion filter size is D in the momentum equations) and we consider now that the FSD Eq. (5) is filtered at this scale. The sgs strain is accordingly defined by ^ l Þ^u0 =D ^ R~c where ^u0 is the S sgs ¼ at R~c ¼ Cð^u0 =S l ; D=d
3061
^ As a first turbulent velocity fluctuation at scale D. approach ^ u0 is estimated from the sub-grid scale turbulence intensity u 0 assuming a Kolmogorov cas^ x Þ1=3 where u 0 is obtained from cade: ^ u0 ¼ u0 ðD=D the Smagorinsky model by a simple dimensional relation: u 0 = tt/(CDx) with C = 0.12 in practice. The ITNFS expression requires as input the integral length scale lt, even in the context of LES [6]. As in LES lt is unknown, the efficiency function proposed by Charlette et al. [18] is retained here as ^ and ^ based only on SGS quantities D u0 . 2.3. Modelling of the FSD transport In the present paper, the unresolved transport is simply closed under a gradient assumption and is written as a turbulent diffusion. As R~c is fil^ the turbulent viscosity ^tt ¼ C^ ^ tered at scale D, u0 D is used instead of tt. Besides, this viscosity is also retained in the transport equation of the progress variable as jr~cj is linked to R~c . The FSD transport terms then read: T res ¼ r ð~ uR~c Þ and T sgs ¼ r ðð^tt =Sct ÞrR~c Þ. 2.4. Control of the flame brush thickness Like in the TFLES approach [6], the flame brush thickness is controlled by a balance between the turbulent transport and source terms of the FSD equation. Considering a 1D steady flame and integrating the FSD equation on the whole computational domain, the natural flame brush thickness dcn is dcn ¼
2 b 4S l ðNeq 1Þ 1þs ð1 þ sÞ 12 þ 3s c 1 þ s þ s3 at
;
ð6Þ
where Neq is the equilibrium wrinkling factor given by a KPP analysis [19] sffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^tt at 2 : ð7Þ Neq ¼ S T KPP =S l ¼ 1 þ bc S l 1 1þs Figure 1 indicates that dcn evolves with sgs turbulence: under engine operating conditions, it can vary from more than 50Dx to less than 5Dx, which corresponds respectively to over-resolution and under-resolution of the flame front. In order to ensure a constant resolution on nres grid points, a correction factor rc is introduced in the expressions of the diffusivity and sgs strain: ^mt is replaced by rc^tt and at by at/rc, where rc is defined as rc ¼ nres Dx =dcn . It should be noted that: • rc does not modify the turbulent flame speed (Eq. (7)). • rc is of the order of one for high levels of turbulence which means the sgs turbulent stretch is weakly modified by rc.
3062
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
3. The spark ignition model
Fig. 1. Evolution of dcn (Eq. (6)) and rc with turbulence ^ x ¼ 1, — rc for rc for D=D intensity. dcn, ^ x ¼ 2, - - - rc for D=D ^ x ¼ 5, – – rc for D=D ^ x ¼ 10. D=D Operating conditions are: P = 20 · 105 Pa, Tfresh gasfuel/air ratio = 1, Sl = 1 m s1, s = 8, es = 700 K, dl = 0.1 mm, Dx = 1 mm.
• When turbulence decreases towards zero, rc^tt and at/rc tend towards zero which means that the flame thickness is no more controlled. This is clearly a drawback of the model, but in practical LES calculations the sgs turbulence level is expected to reach only locally and temporarily low values, so that the flame resolution may not be altered. Moreover, when the sgs strain is low, the FSD is simply convected at the laminar flame speed and the resolved thickness may be conserved. • In practice, to avoid modifying the mixing processes involved in the species and energy transport equations outside the flame front, a dynamic procedure is adopted: the diffusivity rc^mt is only used inside the resolved flame front. Out of the reaction zone, tt = Cu 0 Dx, corresponding to the standard Smagorinsky diffusivity, is retained. Finally, the proposed CFM-LES model reads oq~c q^tt r~c þ qu S l R~c þ r ðq~u~cÞ ¼r rc Sct ot oR~c þ r ð~uR~c Þ ¼ðr ~u NN : r~uÞR~c ot r ðS d NR~c Þ þ S d r ðNÞR~c ^tt þ r rc rR~c Sct ! ^ ^u0 ^ u0 D ; R =r þC ^ ~c c S l dl D c c R~c R~lam R~c c cð1 cÞ
þ bS l
ð8Þ
In SI engines, an electrical spark generated with a spark plug initiates the combustion. This phenomenon may be simply modelled by imposing, at the ignition time tign, an initial profile of R~c close to the spark plug. However, the initial amount of burned gases is very small: typically, the diameter of the initial sphere of burned gases is of the order of the inter-electrode distance, i.e. 1 mm, close to the mesh size. Consequently, such a small flame front cannot be resolved on the mesh and the total reaction rate should be correctly estimated as long as the volume occupied by the burned gases is too small to be resolved. For this purpose, a spark ignition model based on the RANS model AKTIM (Arc and Kernel Tracking Ignition Model) [1] is retained here. The main difference with this model is that no lagrangian tracers are used for describing the flame kernels. Instead, the Eulerian progress variable field is chosen to track the initial flame kernel. The proposed SI model, combines three phases. First, an initial spherical profile of the progress variable is imposed at spark timing tign: ~cðx; tign Þ ¼ c0 =2ð1 tan hðjx xspk j=dign ÞÞ, where ^ in the vicinity of the dign is the mean value of D spark plug, xspk is the position of the spark plug electrodes and c0 is an arbitrary small parameter. The initiation hardly depends on the value of c0 as long as c0 remains small, that is, as long as the initial burned gases mass is negligible compared to the total mass in the chamber. After time tign, the reaction rate (or FSD) has to be defined for the progress variable, Eq. (4). This FSD cannot be given by the proposed model Eq. (5) as the flame front is not fully established. Instead, an ignition FSD called Rign is reconstructed. As in the AKTIM model, the mean flame kernel surface Smean is assumed to keep a roughly spherical shape during the first stage of combustion described by the ignition model. Smean may then be related to the local progress variable by Z 3 4 cðx; tÞ dV ¼ p rb ðtÞ ; V b ðtÞ ¼ 3 X 2 S mean ðtÞ ¼ 4p rb ðtÞ ; where Vb is the total volume occupied by burned gases on the computational domain and rb is the radius of the sphere of volume Vb. Assuming that the wrinkling of the reacting front by turbulent eddies remains small during ignition, the total surface of the flame Stot is given by S tot ðtÞ ¼ Nspk ðtÞS mean ðtÞ;
ð10Þ
where Nspk is the mean wrinkling factor over this surface. As the flame is initially laminar, Nspk should satisfy Nspk(tign) = 1. Then, Nspk increases due to the turbulent stretching of the flame.
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
Consequently, Nspk cannot be given by an equilibrium algebraic expression such as (7). For this reason a model evolution equation for Nspk is retained dNspk Neq Nspk ¼ hat is Nspk ; ð11Þ dt Neq 1 s where ÆQæs = XQRign(x,t) dV/XRign(x,t) dV represents the average of Q over the flame surface. Equations (9)–(11) allow to compute the total flame surface Stot(t) at each instant. This surface now needs to be distributed spatially over the computational grid, that is, the local FSD Rign(x, t) used in Eq. (4) needs to be defined. For this purpose the FSD is expressed as in Boger [5] Rign ðx; tÞ ¼ a~cðx; tÞð1 ~cðx; tÞÞ;
ð12Þ
where a is a global coefficient which is determined by imposing the integral of Rign(x, t) over the domain to be equal to the total R surface Stot(t) defined previously : a ¼ S tot ðtÞ= X ~cðx; tÞð1 ~cðx; tÞÞ dV . The FSD Eq. (8) is theoretically valid only when a fully established flame front is considered. Therefore the reaction rate is defined by Rign until the progress variable reaches the value 1 somewhere in the domain. At this instant ttransition, the flame surface density Rign is transferred to the FSD transport equation: R~c ðttransition ; xÞ ¼ Rign ðttransition ; xÞ. For t > ttransition, the FSD is given by Eq. (8) and the spark ignition model is not used anymore.
3063
defined as: rs = (3Vb/4p)1/3 and rp = (Stot/4p)1/2 where Vb is the total volume occupied by burned gases R and Stot the total flame surface S tot ¼ X R~c dV . The mean resolved Rwrinkling Nres = Sres/Smean is defined by S res ¼ X jrcj dV and S mean ¼ 4pr2s . Surfaces are finally related by: Stot = NsgsNresSmean. Simulations were performed ^ ¼ 3 mmÞ correfor two filter sizes. Case A (D sponds to a typical LES resolution. Case B ^ ¼ 6:5 mm lt Þ is equivalent to a RANS calcu(D lation, and therefore no resolved turbulent velocity field was used and ^ u0 ¼ 0:9S l was imposed. For both cases, the mean and total radii presented in Fig. 2 correctly match the experimental evolutions. Figure 3 shows that the initial laminar flame is progressively wrinkled by turbulence. Case A is close to a DNS calculation since the flame wrinkling is entirely resolved (Nsgs 1). In case B the filtered flame remains perfectly spherical (Nres = 1), and the wrinkling is totally at the
4. Validation on an experimental configuration The ignition-combustion models presented above were validated on the experiment of Renou and Boukhalfa [4]. A propane–air mixture is injected through a turbulence grid into a channel where a thin spark plug allows to ignite the mixture. The turbulence is found to be nearly isotropic (PIV and laser tomography imaging) and hardly decays during the flame kernel growth. Laminar spherical propagation cases have also been analysed. The computations have been carried out on a 20 · 20 · 20 mm cubic mesh composed of regular hexahedra with a characteristic size of 0.5 mm. This size corresponds to the present target resolution for engine cycle calculations. For simplicity, the computational domain is assumed to follow the convected homogeneous and isotropic turbulence (HIT), so that the flame centre remains roughly at the centre of the domain. In practice, a HIT is simply initialized according to a Passot–Pouquet spectrum and outlet boundary conditions [20] are adopted. A stoichiometric propane-air flame propagates in a turbulent field characterized by u 0 /Sl = 0.9 and lt = 6.5 mm. For the need of comparison with experiments, mean and total radii are respectively
Fig. 2. Evolution of mean and total flame radii, s rs experiment, rs Case A, — rs Case B, h rp experiment, rp Case A, - - - rp Case B.
Fig. 3. Flame wrinkling evolution, s experiment, Ntot, — Nres, - - -Nsgs.
3064
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
sgs level. The CFM-LES model therefore shows a very low dependency to the combustion filter size, even for such extreme cases. 5. Towards engine applications The proposed ignition and combustion models have been tested on a real engine configuration. For this purpose, intake, compression and combustion strokes have been computed on the four-valve pent-roof PSA XU10 engine. This configuration has been investigated experimentally and is characterised by strong cycle-to-cycle variations [21]. The main engine specifications are summarised in Table 1. The Arbitrary Lagrangian Eulerian formalism was previously introduced in the AVBP solver to account for moving meshes [11]. Turbulent transports are described with the classical Smagorinsky model. NSCBC [20] inflow/outflow boundary conditions (BC) were implemented to impose the experimental pressure in inlet and exhaust ports. Isothermal BC are adopted to compute wall heat fluxes and wall/flow interactions are treated with a ‘‘law of the wall’’ based on the linear/log laws. Finally wall/flame interactions are taken into account extending the model of Bruneaux et al. [22] to LES. Computations have been performed on unstructured meshes composed of hexahedra
Table 1 Engine specifications Bore · Stroke
86 · 86 mm
Compression ratio Engine speed (RPM) Fuel/air ratio Volumetric efficiency Spark timing
10 2000 0.7 (propane) 0.35 40 CA deg
elements (Fig. 4). The number of elements varies with the movement of piston and valves in order to keep a roughly constant mesh resolution of the order of Dx = 0.5 mm (330 000 elements at top dead centre and 660 000 at bottom dead centre). The main parameters of the calculations are summarised in Table 2. A sensitivity analysis of ^ is the models to the LES combustion filter size D first performed by modifying the resolution parameter nres and by using a more refined mesh with a typical resolution of Dx = 0.25 mm. Results are also compared pffiffiffiffiffiffiffito ffi the Boger’s algebraic model ^ where Neq is cho[5]: R~c ¼ Neq 4 6=p~cð1 ~cÞ=D, sen equal to the equilibrium wrinkling factor of the CFM-LES model, Eq. (7). The cylinder pressure computed with the CFM-LES model is correctly situated in the experimental envelope (Fig. 5). On the contrary, the Boger’s model overestimates the pressure increase because the SGS wrinkling is assumed to be in equilibrium with turbulence (Fig. 6), which is not confirmed by CFM-LES calculations where Nsgs reaches its equilibrium value after about 20 CAD. Figure 7 shows that the total flame surface using CFM-LES is nearly independent of the ^ when D ^ increases (calculations flame filter size D: C to A to B), the resolved flame surface decreases
Table 2 Parameters of the different simulations ^ Case Model nres =DðmmÞ
Mesh
A B C D E
Reference Reference Refined Reference Reference
CFM-LES CFM-LES CFM-LES Boger CFM-LES
5/2.5 10/5 5/1.25 5/2.5 –/2.5
Crank angles are relative to the combustion top dead centre.
Fig. 4. Computational mesh during the intake stroke.
Fig. 5. Cylinder pressure vs crank angle, – – experimenCase B, — Case C, tal envelope, s Case A, Case D, h Case E. Transitions between SI and combustion models are represented by arrows.
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
Fig. 6. Sub-grid-scale and equilibrium wrinkling evoluÆNsgsæs Case A, — ÆNeqæs Case A, – – . ÆNsgsæs tions, Case D.
3065
to unity during the calculation and the growing rate of the flame surface is not strongly modified. This hypothesis is confirmed by case E performed setting rc = 1 which behaves like case A (Fig. 7). The spatial distribution of flame surface is also similar to case A (Fig. 8) which indicates that starting from the initial profile of R~c with a thick^ the flame propagation time in the ness close to D, engine is too short to allow the flame front to evolve towards its natural thickness dcn. This figure also shows higher peak values of R~c for case E compared to case A. The reason is that rc increases with turbulence intensity u 0 (Fig. 1), then, for case A, rc tends to decrease the turbulent stretch at/rc in regions of high turbulence and vice versa. On the contrary, this smoothing effect is absent for case E. 6. Conclusions
Fig. 7. Resolved and total flame surfaces vs crank angle, Total Case A, s Resolved Case A, Total Case B, h Resolved Case B, — Total Case C, + Resolved Case C, - - - Total Case E, D Resolved Case E.
(Figs. 7 and 8) and the SGS flame surface increases accordingly. This feature is essential for a LES model to ensure that the solution is nearly independent of the mesh resolution. The thickness correction parameter rc (not shown) remains close
This work validates the ability of the CFMLES and spark-ignition models to simulate SI engine cycles. First, these models correctly reproduce the flame wrinkling observed in turbulent ignition-propagation experiments of a propaneair flame [4]. Then, simulations on a real engine configuration have been presented. Results obtained with CFM-LES correctly match the experimental data, while those obtained with the Boger et al.’s algebraic closure [5] overestimate the pressure growth during the first stages of the combustion. This finding shows the ability of the FSD balance equation to handle non-equilibrium situations between flame surface production and destruction when an algebraic model assumes such an equilibrium. The global results obtained with the CFM-LES model are also independent of the flame front resolution: a decrease of the resolved flame front wrinkling, due to a larger combustion filter size, is correctly balanced by an increase of the sub-grid scale flame front wrinkling. This point is very important to ensure that the overall reaction rate does not depend on the computational mesh resolution. Multi-cycle simulations are currently under investigation on the XU10 engine in order to check the ability of the
Fig. 8. FSD fields and isocontours of progress variable. Pictures are taken at 0 deg on a cut plane situated 3 mm under the spark plug (represented by +). (a) Case B; (b) Case A; (c) Case C; (d) Case E.
3066
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
proposed LES approach to recover experimental cyclic variations. Acknowledgments The authors are grateful to Drs. B. Renou and B. Lecordier from COmplexe de Recherche Interprofessionnel en Ae´rothermochimie (CORIA Rouen - France) for providing experimental data of the configuration studied in this paper. This work was partially funded by the European Commission within the LESSCO2 RTD project (project number NNE5-2001-00495). References [1] J.M. Duclos, O. Colin, Arc and Kernel Tracking Ignition Model for 3D SI Engines Calculations, Comodia, Nagoya, Japan, 2001, pp. 343–350. [2] S. Candel, T. Poinsot, Combust. Sci. Tech. 70 (1990) 1–15. [3] O. Colin, A. Benkenida, C. Angelberger, Oil & Gas Sci. Tech.—Rev. IFP 58 (1) (2003) 32–47. [4] B. Renou, A. Boukhalfa, Combust. Sci. Tech. 162 (2001) 347–371. [5] M. Boger, D. Veynante, H. Boughanem, A. Trouve´, Proc. Combust. Inst. 27 (1998) 917–925. [6] O. Colin, F. Ducros, D. Veynante, T. Poinsot, Phys. Fluids 12 (7) (2002) 1843–1863. [7] N. Peters, Turbulent Combustion, Cambridge University Press, 2000. [8] R. Borghi, On the Structure and Morphology of Turbulent Premixed Flames, Rec. Adv. Aerosp. Sci. (1985) 117–138.
[9] H. Weller, G. Tabor, A. Gosman, C. Fureby, Proc. Combust. Inst. 27 (1998) 899–907. [10] E.R. Hawkes, R. Cant, Combust. Flame 126 (2001) 1617–1629. [11] V. Moureau, G. Lartigue, Y. Sommerer, C. Angelberger, O. Colin, T. Poinsot, J. Comp. Phys. 202 (2005) 710–736. [12] S. Richard, O. Colin, O. Vermorel, A. Benkenida, D. Veynante, in: ECCOMAS Thematic Conference on Computational Combustion, Lisbon, 2005. [13] A. Trouve´, T. Poinsot, J. Fluid Mech. 278 (1994) 1– 31. [14] S. Pope, Int. J. Eng. Sci. 26 (5) (1988) 445– 469. [15] D. Veynante, T. Poinsot, RANS and LES modeling for turbulent combustion, New tools in turbulent modelling, Les e´ditions de physique, Springer-Verlag, 1997. [16] G. Rymer, Analyse et mode´lisation du taux de re´action moyen et des me´canismes de transport en combustion turbulente pre´me´lange´e, PhD thesis, Ecole Centrale de Paris-FR, 2001. [17] C. Meneveau, T. Poinsot, Combust. Flame 86 (1991) 311–332. [18] F. Charlette, C. Meneveau, D. Veynante, Combust. Flame 131 (2002) 150–180. [19] J.M. Duclos, D. Veynante, T. Poinsot, Combust. Flame 95 (1993) 101–117. [20] T. Poinsot, S. Lele, J. Comp. Phys. 101 (1992) 104– 129. [21] A. Torres, S. Henriot, 3D Modeling of Combustion in Lean Burn Four-valve Engines: Influence of Intake Configuration, Comodia, Yokohama, Japan, 1994, pp.139–156. [22] G. Bruneaux, T. Poinsot, J. Ferziger, J. Fluid Mech. 349 (1997) 199–219.
Proceedings of the Combustion Institute 31 (2007) 3059–3066
Combustion Institute www.elsevier.com/locate/proci
Towards large eddy simulation of combustion in spark ignition engines S. Richard a, O. Colin a,*, O. Vermorel a, A. Benkenida a, C. Angelberger a, D. Veynante b b
a IFP, 1 et 4 avenue de Bois-Pre´au, 92852 Rueil-Malmaison Cedex, France Laboratoire EM2C, CNRS Ecole Centrale Paris, 92295 Chaˆtenay-Malabry Cedex, France
Abstract Internal combustion engine simulations are commonly performed using the RANS (Reynolds averaged Navier–Stokes) approach. It gives a correct estimates of global quantities but is by nature not adapted to describe phenomena strongly linked to cyclic variations. On the other hand, large eddy simulation (LES) is a promising technique to determine successive engine cycles. This work demonstrates the feasibility of LES engine cycles simulation by using a flame surface density (FSD) approach. This approach, presented in a first section, combines an Eulerian spark ignition model derived from the RANS AKTIM model [J.M. Duclos, O. Colin, Arc and Kernel Tracking Ignition Model for 3D SI Engines Calculations, Comodia, Nagoya, Japan, 2001, pp. 343–350] and a Coherent Flame Model (CFM) [S. Candel, T. Poinsot, Combust. Sci. Tech. 70 (1990) 1–15; O. Colin, A. Benkenida, C. Angelberger, Oil & Gas Sci. Techn.—Rev. IFP 58 (1) (2003) 47–32] describing the flame propagation. The CFM model, commonly used in RANS simulations, is here formulated in a LES context. In a second part, the whole ignition-combustion model is validated against an experiment relative to the turbulent ignition and flame propagation of a stoichiometric propane-air mixture [B. Renou, A. Boukhalfa, Combust. Sci. Tech. 162 (2001) 347–371]. Finally, LES engine cycles simulations are performed on a real engine configuration. First, the sensitivity of the model to the ^ is examined, showing a weak dependence of the modelling approach to D. ^ LES combustion filter size D Then results are compared to those obtained with the algebraic model for the FSD proposed by Boger et al. [M. Boger, D. Veynante, H. Boughanem, A. Trouve´, Proc. Combust. Inst. 27 (1998) 917–925] and the need for non-equilibrium combustion models is demonstrated. 2006 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: LES; Combustion; CFM; Spark ignition; Engine cycle simulation
1. Introduction
*
Corresponding author. Fax: +33 0 147527068. E-mail addresses: [email protected] (S. Richard), [email protected] (O. Colin), [email protected] (D. Veynante).
The RANS (Reynolds averaged Navier– Stokes) approach provides today correct predictions of mean quantities such as pressure or pollutant emissions in internal combustion engines but does not describe unsteady phenomena such as cyclic variations or engine transients. On the other
1540-7489/$ - see front matter 2006 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.proci.2006.07.086
3060
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
hand, large eddy simulation (LES), where resolved flow motions are explicitly computed, is well adapted to unsteady problems. However, the effects of sub-grid-scale (sgs) motions should be modelled. Previous premixed combustion modelling approaches in LES have been mainly based on the artificial thickening of a flame [6] or on the so-called G-equation [7]. The laminar flamelet approach appears as an interesting alternative for piston engines, according to the combustion diagram from Borghi [8] and considering typical Damko¨hler and Karlovitz numbers of about Da 10–100 and Ka 0.01–1. Algebraic models have already been proposed for the flame surface density (FSD) in LES applications [5]. Balance equations for the wrinkling factor [9] or the FSD [10] have also been investigated. In piston engines, where the flame is initially laminar at the spark plug, then grows under the action of thermal expansion, curvature and turbulence and is finally quenched when interacting with walls, the correct description of the wrinkling evolution is crucial, giving an advantage to the balance equation formalism. The present work is devoted to the presentation and the validation of a FSD equation and a spark ignition model dedicated to piston engine applications. These developments have been implemented in the AVBP code [11] specifically dedicated to LES combustion calculations and have already been successfully tested on simple test cases [12]. The first section summarises the modelling approach. In a second part, the whole ignition–combustion model is validated against an experiment relative to the ignition and flame propagation [4]. Finally, simulations are performed on a real engine configuration. The sensitivity of the model to the LES filter size is first examined. Results are then compared to those obtained with the algebraic closure of FSD proposed by Boger et al. [5]. 2. The CFM-LES combustion model Assuming a flamelet regime, a single step chemistry and a unit Lewis number, the thermochemistry of the reacting flow may be described by the progress variable c of the reaction (c = 0 in fresh gases and c = 1 in fully burnt gases) according to P
oq~c þ r ðq~u~cÞ þ r ðquc q~u~cÞ ¼ qS d jrcj; ot ð1Þ
where Sd is the local propagation speed. Q denotes ~ a mass weighted filtered a filtered quantity and Q ~ ¼ qQ). The filtered FSD is defined as quantity (qQ R ¼ jrcj [5] and introducing the weighted flame
surface average as hQis ¼ QR=R, the source term in Eq. (1) is written as qS d jrcj ¼ hqS d is R. In practice, ÆqSdæs is approximated by quSl [13] where qu is the unburned gases density and Sl the laminar flame speed. Equation (1) contains two unclosed terms: an unresolved transport term and the filtered FSD. In the CFM-LES approach, the filtered FSD is computed through a transport equation. An exact unclosed equation for the local FSD was first proposed by Pope [14] or Candel and Poinsot [2] and filtered, leading to the following equation [10–15]: oR þ r huis R ¼ hr u nn : ruis R ot r ðhS d nis RÞ þ hS d r nis R; ð2Þ where n is the local normal to the flame front pointing towards the fresh gases. The first term on the right-hand side corresponds to the straining of the flame by the flow. The second and third terms represent, respectively the planar propagation of the flame at speed Sd and the propagation of a curved front. The unresolved transport term of the ~c equation is commonly closed under a gradient assumption and is written as a turbulent diffusion. Following this type of formulation, this term vanishes when turbulence decreases. However, Rymer [16] showed that in a laminar case the unresolved transport is non-zero and corresponds to nongradient transport (NGT) which can be written as r ðquc q~ u~cÞ ¼ qu S l r ððc ~cÞNÞ, where N ¼ rc=jrcj is the normal to the iso-surface of the filtered progress variable. Rymer then proposed the following general formulation: r ðquc q~u~cÞ ¼ qu S l r ððc ~cÞNÞ r
qtt r~c ; Sct ð3Þ
where tt is the turbulent viscosity and Sct a turbulent Schmidt number. Unfortunately, NGT terms may lead to numerical instabilities because they are anti-diffusive and are difficult to handle in practical applications. Consequently, this term is usually omitted and the equation becomes oq~c qtt þ r ðq~ u~cÞ ¼ r r~c þ qu S l R~c ot Sct
ð4Þ
with R~c defined as: R~c ¼ R þ r ððc ~cÞNÞ. A phenomenological transport equation for R~c is proposed here, assuming that this equation satisfies the same main properties as Eq. (2), i.e., it should ensure: (i) the laminar propagation of the flame in multidimensional cases, (ii) a reasonable behaviour—in particular correct turbulent flame speeds—in case of high turbulence levels and (iii)
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
the control of the flame brush thickness. The CFM-LES, unlike RANS formulations, is a spatially filtered equation and it takes into account the effects of the resolved flow on the flame front. For this purpose the FSD transport equation terms are split into resolved and unresolved parts oR~c þ T res þ T sgs ¼ S res þ S sgs þ C res þ C sgs þ P ; ot ð5Þ where Tres, Sres, Cres and P are, respectively, the transport, strain, curvature and propagation terms due to resolved flow motions, and Tsgs, Ssgs and Csgs are, respectively, the unresolved transport, strain and curvature terms. 2.1. Modelling of the flame propagation and curvature The propagation and resolved curvature terms are physically linked and ensure the laminar flame propagation, when the sgs turbulence is low. The proposed formulation of these terms is based on the normal to the iso-surface of the filtered progress variable and reads: P ¼ r ðS d NR~c Þ and C res ¼ S d ðr NÞR~c . Sd is defined by S d ¼ ð1 þ s~cÞS l , where s = qu/qb 1 is the thermal expansion rate across the flame front and qb is the burned gases density. DNS simulations [13] have shown that the sgs curvature acts as a sink term towards burned gases and as a source term towards fresh gases. Hence, Rymer proposed to express sgs curvature effects as: C sgs ¼ bS l ðc cÞ=ðcð1 cÞÞ ðR jrcjÞR, where b = 4/3 and c* = 0.5 are modelling constants. This formulation is retained here, corrected for the R~c equation: C sgs ¼ bS l ðc cÞ=ðcð1 cÞÞðR~c R~lam c is c ÞR~c , where estimated from the usual BML expression as: c ¼ ~cq=qb and R~lam ¼ jr~cj þ ðc ~cÞr N. c 2.2. Modelling of the flame strain rate The resolved strain is S res ¼ ðr ~u NN : r~uÞR~c and corresponds to the flame strain rate due to the resolved structures. The turbulent strain rate in RANS is classically modeled as at = C(u 0 /Sl, lt/dl)u 0 /lt, where C is an efficiency function taking into account the ability of all vortices to wrinkle the flame, often defined by the ITNFS function of Meneveau and Poinsot [17]. In LES, Colin et al. [6] showed that resolved eddies smaller than the flame thickness defined by ^ ¼ nres Dx (nres 5–10 for finite volume solvers) D are not able to wrinkle the flame front. Consequent^ (and not Dx like ly the combustion filter size is D in the momentum equations) and we consider now that the FSD Eq. (5) is filtered at this scale. The sgs strain is accordingly defined by ^ l Þ^u0 =D ^ R~c where ^u0 is the S sgs ¼ at R~c ¼ Cð^u0 =S l ; D=d
3061
^ As a first turbulent velocity fluctuation at scale D. approach ^ u0 is estimated from the sub-grid scale turbulence intensity u 0 assuming a Kolmogorov cas^ x Þ1=3 where u 0 is obtained from cade: ^ u0 ¼ u0 ðD=D the Smagorinsky model by a simple dimensional relation: u 0 = tt/(CDx) with C = 0.12 in practice. The ITNFS expression requires as input the integral length scale lt, even in the context of LES [6]. As in LES lt is unknown, the efficiency function proposed by Charlette et al. [18] is retained here as ^ and ^ based only on SGS quantities D u0 . 2.3. Modelling of the FSD transport In the present paper, the unresolved transport is simply closed under a gradient assumption and is written as a turbulent diffusion. As R~c is fil^ the turbulent viscosity ^tt ¼ C^ ^ tered at scale D, u0 D is used instead of tt. Besides, this viscosity is also retained in the transport equation of the progress variable as jr~cj is linked to R~c . The FSD transport terms then read: T res ¼ r ð~ uR~c Þ and T sgs ¼ r ðð^tt =Sct ÞrR~c Þ. 2.4. Control of the flame brush thickness Like in the TFLES approach [6], the flame brush thickness is controlled by a balance between the turbulent transport and source terms of the FSD equation. Considering a 1D steady flame and integrating the FSD equation on the whole computational domain, the natural flame brush thickness dcn is dcn ¼
2 b 4S l ðNeq 1Þ 1þs ð1 þ sÞ 12 þ 3s c 1 þ s þ s3 at
;
ð6Þ
where Neq is the equilibrium wrinkling factor given by a KPP analysis [19] sffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^tt at 2 : ð7Þ Neq ¼ S T KPP =S l ¼ 1 þ bc S l 1 1þs Figure 1 indicates that dcn evolves with sgs turbulence: under engine operating conditions, it can vary from more than 50Dx to less than 5Dx, which corresponds respectively to over-resolution and under-resolution of the flame front. In order to ensure a constant resolution on nres grid points, a correction factor rc is introduced in the expressions of the diffusivity and sgs strain: ^mt is replaced by rc^tt and at by at/rc, where rc is defined as rc ¼ nres Dx =dcn . It should be noted that: • rc does not modify the turbulent flame speed (Eq. (7)). • rc is of the order of one for high levels of turbulence which means the sgs turbulent stretch is weakly modified by rc.
3062
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
3. The spark ignition model
Fig. 1. Evolution of dcn (Eq. (6)) and rc with turbulence ^ x ¼ 1, — rc for rc for D=D intensity. dcn, ^ x ¼ 2, - - - rc for D=D ^ x ¼ 5, – – rc for D=D ^ x ¼ 10. D=D Operating conditions are: P = 20 · 105 Pa, Tfresh gasfuel/air ratio = 1, Sl = 1 m s1, s = 8, es = 700 K, dl = 0.1 mm, Dx = 1 mm.
• When turbulence decreases towards zero, rc^tt and at/rc tend towards zero which means that the flame thickness is no more controlled. This is clearly a drawback of the model, but in practical LES calculations the sgs turbulence level is expected to reach only locally and temporarily low values, so that the flame resolution may not be altered. Moreover, when the sgs strain is low, the FSD is simply convected at the laminar flame speed and the resolved thickness may be conserved. • In practice, to avoid modifying the mixing processes involved in the species and energy transport equations outside the flame front, a dynamic procedure is adopted: the diffusivity rc^mt is only used inside the resolved flame front. Out of the reaction zone, tt = Cu 0 Dx, corresponding to the standard Smagorinsky diffusivity, is retained. Finally, the proposed CFM-LES model reads oq~c q^tt r~c þ qu S l R~c þ r ðq~u~cÞ ¼r rc Sct ot oR~c þ r ð~uR~c Þ ¼ðr ~u NN : r~uÞR~c ot r ðS d NR~c Þ þ S d r ðNÞR~c ^tt þ r rc rR~c Sct ! ^ ^u0 ^ u0 D ; R =r þC ^ ~c c S l dl D c c R~c R~lam R~c c cð1 cÞ
þ bS l
ð8Þ
In SI engines, an electrical spark generated with a spark plug initiates the combustion. This phenomenon may be simply modelled by imposing, at the ignition time tign, an initial profile of R~c close to the spark plug. However, the initial amount of burned gases is very small: typically, the diameter of the initial sphere of burned gases is of the order of the inter-electrode distance, i.e. 1 mm, close to the mesh size. Consequently, such a small flame front cannot be resolved on the mesh and the total reaction rate should be correctly estimated as long as the volume occupied by the burned gases is too small to be resolved. For this purpose, a spark ignition model based on the RANS model AKTIM (Arc and Kernel Tracking Ignition Model) [1] is retained here. The main difference with this model is that no lagrangian tracers are used for describing the flame kernels. Instead, the Eulerian progress variable field is chosen to track the initial flame kernel. The proposed SI model, combines three phases. First, an initial spherical profile of the progress variable is imposed at spark timing tign: ~cðx; tign Þ ¼ c0 =2ð1 tan hðjx xspk j=dign ÞÞ, where ^ in the vicinity of the dign is the mean value of D spark plug, xspk is the position of the spark plug electrodes and c0 is an arbitrary small parameter. The initiation hardly depends on the value of c0 as long as c0 remains small, that is, as long as the initial burned gases mass is negligible compared to the total mass in the chamber. After time tign, the reaction rate (or FSD) has to be defined for the progress variable, Eq. (4). This FSD cannot be given by the proposed model Eq. (5) as the flame front is not fully established. Instead, an ignition FSD called Rign is reconstructed. As in the AKTIM model, the mean flame kernel surface Smean is assumed to keep a roughly spherical shape during the first stage of combustion described by the ignition model. Smean may then be related to the local progress variable by Z 3 4 cðx; tÞ dV ¼ p rb ðtÞ ; V b ðtÞ ¼ 3 X 2 S mean ðtÞ ¼ 4p rb ðtÞ ; where Vb is the total volume occupied by burned gases on the computational domain and rb is the radius of the sphere of volume Vb. Assuming that the wrinkling of the reacting front by turbulent eddies remains small during ignition, the total surface of the flame Stot is given by S tot ðtÞ ¼ Nspk ðtÞS mean ðtÞ;
ð10Þ
where Nspk is the mean wrinkling factor over this surface. As the flame is initially laminar, Nspk should satisfy Nspk(tign) = 1. Then, Nspk increases due to the turbulent stretching of the flame.
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
Consequently, Nspk cannot be given by an equilibrium algebraic expression such as (7). For this reason a model evolution equation for Nspk is retained dNspk Neq Nspk ¼ hat is Nspk ; ð11Þ dt Neq 1 s where ÆQæs = XQRign(x,t) dV/XRign(x,t) dV represents the average of Q over the flame surface. Equations (9)–(11) allow to compute the total flame surface Stot(t) at each instant. This surface now needs to be distributed spatially over the computational grid, that is, the local FSD Rign(x, t) used in Eq. (4) needs to be defined. For this purpose the FSD is expressed as in Boger [5] Rign ðx; tÞ ¼ a~cðx; tÞð1 ~cðx; tÞÞ;
ð12Þ
where a is a global coefficient which is determined by imposing the integral of Rign(x, t) over the domain to be equal to the total R surface Stot(t) defined previously : a ¼ S tot ðtÞ= X ~cðx; tÞð1 ~cðx; tÞÞ dV . The FSD Eq. (8) is theoretically valid only when a fully established flame front is considered. Therefore the reaction rate is defined by Rign until the progress variable reaches the value 1 somewhere in the domain. At this instant ttransition, the flame surface density Rign is transferred to the FSD transport equation: R~c ðttransition ; xÞ ¼ Rign ðttransition ; xÞ. For t > ttransition, the FSD is given by Eq. (8) and the spark ignition model is not used anymore.
3063
defined as: rs = (3Vb/4p)1/3 and rp = (Stot/4p)1/2 where Vb is the total volume occupied by burned gases R and Stot the total flame surface S tot ¼ X R~c dV . The mean resolved Rwrinkling Nres = Sres/Smean is defined by S res ¼ X jrcj dV and S mean ¼ 4pr2s . Surfaces are finally related by: Stot = NsgsNresSmean. Simulations were performed ^ ¼ 3 mmÞ correfor two filter sizes. Case A (D sponds to a typical LES resolution. Case B ^ ¼ 6:5 mm lt Þ is equivalent to a RANS calcu(D lation, and therefore no resolved turbulent velocity field was used and ^ u0 ¼ 0:9S l was imposed. For both cases, the mean and total radii presented in Fig. 2 correctly match the experimental evolutions. Figure 3 shows that the initial laminar flame is progressively wrinkled by turbulence. Case A is close to a DNS calculation since the flame wrinkling is entirely resolved (Nsgs 1). In case B the filtered flame remains perfectly spherical (Nres = 1), and the wrinkling is totally at the
4. Validation on an experimental configuration The ignition-combustion models presented above were validated on the experiment of Renou and Boukhalfa [4]. A propane–air mixture is injected through a turbulence grid into a channel where a thin spark plug allows to ignite the mixture. The turbulence is found to be nearly isotropic (PIV and laser tomography imaging) and hardly decays during the flame kernel growth. Laminar spherical propagation cases have also been analysed. The computations have been carried out on a 20 · 20 · 20 mm cubic mesh composed of regular hexahedra with a characteristic size of 0.5 mm. This size corresponds to the present target resolution for engine cycle calculations. For simplicity, the computational domain is assumed to follow the convected homogeneous and isotropic turbulence (HIT), so that the flame centre remains roughly at the centre of the domain. In practice, a HIT is simply initialized according to a Passot–Pouquet spectrum and outlet boundary conditions [20] are adopted. A stoichiometric propane-air flame propagates in a turbulent field characterized by u 0 /Sl = 0.9 and lt = 6.5 mm. For the need of comparison with experiments, mean and total radii are respectively
Fig. 2. Evolution of mean and total flame radii, s rs experiment, rs Case A, — rs Case B, h rp experiment, rp Case A, - - - rp Case B.
Fig. 3. Flame wrinkling evolution, s experiment, Ntot, — Nres, - - -Nsgs.
3064
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
sgs level. The CFM-LES model therefore shows a very low dependency to the combustion filter size, even for such extreme cases. 5. Towards engine applications The proposed ignition and combustion models have been tested on a real engine configuration. For this purpose, intake, compression and combustion strokes have been computed on the four-valve pent-roof PSA XU10 engine. This configuration has been investigated experimentally and is characterised by strong cycle-to-cycle variations [21]. The main engine specifications are summarised in Table 1. The Arbitrary Lagrangian Eulerian formalism was previously introduced in the AVBP solver to account for moving meshes [11]. Turbulent transports are described with the classical Smagorinsky model. NSCBC [20] inflow/outflow boundary conditions (BC) were implemented to impose the experimental pressure in inlet and exhaust ports. Isothermal BC are adopted to compute wall heat fluxes and wall/flow interactions are treated with a ‘‘law of the wall’’ based on the linear/log laws. Finally wall/flame interactions are taken into account extending the model of Bruneaux et al. [22] to LES. Computations have been performed on unstructured meshes composed of hexahedra
Table 1 Engine specifications Bore · Stroke
86 · 86 mm
Compression ratio Engine speed (RPM) Fuel/air ratio Volumetric efficiency Spark timing
10 2000 0.7 (propane) 0.35 40 CA deg
elements (Fig. 4). The number of elements varies with the movement of piston and valves in order to keep a roughly constant mesh resolution of the order of Dx = 0.5 mm (330 000 elements at top dead centre and 660 000 at bottom dead centre). The main parameters of the calculations are summarised in Table 2. A sensitivity analysis of ^ is the models to the LES combustion filter size D first performed by modifying the resolution parameter nres and by using a more refined mesh with a typical resolution of Dx = 0.25 mm. Results are also compared pffiffiffiffiffiffiffito ffi the Boger’s algebraic model ^ where Neq is cho[5]: R~c ¼ Neq 4 6=p~cð1 ~cÞ=D, sen equal to the equilibrium wrinkling factor of the CFM-LES model, Eq. (7). The cylinder pressure computed with the CFM-LES model is correctly situated in the experimental envelope (Fig. 5). On the contrary, the Boger’s model overestimates the pressure increase because the SGS wrinkling is assumed to be in equilibrium with turbulence (Fig. 6), which is not confirmed by CFM-LES calculations where Nsgs reaches its equilibrium value after about 20 CAD. Figure 7 shows that the total flame surface using CFM-LES is nearly independent of the ^ when D ^ increases (calculations flame filter size D: C to A to B), the resolved flame surface decreases
Table 2 Parameters of the different simulations ^ Case Model nres =DðmmÞ
Mesh
A B C D E
Reference Reference Refined Reference Reference
CFM-LES CFM-LES CFM-LES Boger CFM-LES
5/2.5 10/5 5/1.25 5/2.5 –/2.5
Crank angles are relative to the combustion top dead centre.
Fig. 4. Computational mesh during the intake stroke.
Fig. 5. Cylinder pressure vs crank angle, – – experimenCase B, — Case C, tal envelope, s Case A, Case D, h Case E. Transitions between SI and combustion models are represented by arrows.
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
Fig. 6. Sub-grid-scale and equilibrium wrinkling evoluÆNsgsæs Case A, — ÆNeqæs Case A, – – . ÆNsgsæs tions, Case D.
3065
to unity during the calculation and the growing rate of the flame surface is not strongly modified. This hypothesis is confirmed by case E performed setting rc = 1 which behaves like case A (Fig. 7). The spatial distribution of flame surface is also similar to case A (Fig. 8) which indicates that starting from the initial profile of R~c with a thick^ the flame propagation time in the ness close to D, engine is too short to allow the flame front to evolve towards its natural thickness dcn. This figure also shows higher peak values of R~c for case E compared to case A. The reason is that rc increases with turbulence intensity u 0 (Fig. 1), then, for case A, rc tends to decrease the turbulent stretch at/rc in regions of high turbulence and vice versa. On the contrary, this smoothing effect is absent for case E. 6. Conclusions
Fig. 7. Resolved and total flame surfaces vs crank angle, Total Case A, s Resolved Case A, Total Case B, h Resolved Case B, — Total Case C, + Resolved Case C, - - - Total Case E, D Resolved Case E.
(Figs. 7 and 8) and the SGS flame surface increases accordingly. This feature is essential for a LES model to ensure that the solution is nearly independent of the mesh resolution. The thickness correction parameter rc (not shown) remains close
This work validates the ability of the CFMLES and spark-ignition models to simulate SI engine cycles. First, these models correctly reproduce the flame wrinkling observed in turbulent ignition-propagation experiments of a propaneair flame [4]. Then, simulations on a real engine configuration have been presented. Results obtained with CFM-LES correctly match the experimental data, while those obtained with the Boger et al.’s algebraic closure [5] overestimate the pressure growth during the first stages of the combustion. This finding shows the ability of the FSD balance equation to handle non-equilibrium situations between flame surface production and destruction when an algebraic model assumes such an equilibrium. The global results obtained with the CFM-LES model are also independent of the flame front resolution: a decrease of the resolved flame front wrinkling, due to a larger combustion filter size, is correctly balanced by an increase of the sub-grid scale flame front wrinkling. This point is very important to ensure that the overall reaction rate does not depend on the computational mesh resolution. Multi-cycle simulations are currently under investigation on the XU10 engine in order to check the ability of the
Fig. 8. FSD fields and isocontours of progress variable. Pictures are taken at 0 deg on a cut plane situated 3 mm under the spark plug (represented by +). (a) Case B; (b) Case A; (c) Case C; (d) Case E.
3066
S. Richard et al. / Proceedings of the Combustion Institute 31 (2007) 3059–3066
proposed LES approach to recover experimental cyclic variations. Acknowledgments The authors are grateful to Drs. B. Renou and B. Lecordier from COmplexe de Recherche Interprofessionnel en Ae´rothermochimie (CORIA Rouen - France) for providing experimental data of the configuration studied in this paper. This work was partially funded by the European Commission within the LESSCO2 RTD project (project number NNE5-2001-00495). References [1] J.M. Duclos, O. Colin, Arc and Kernel Tracking Ignition Model for 3D SI Engines Calculations, Comodia, Nagoya, Japan, 2001, pp. 343–350. [2] S. Candel, T. Poinsot, Combust. Sci. Tech. 70 (1990) 1–15. [3] O. Colin, A. Benkenida, C. Angelberger, Oil & Gas Sci. Tech.—Rev. IFP 58 (1) (2003) 32–47. [4] B. Renou, A. Boukhalfa, Combust. Sci. Tech. 162 (2001) 347–371. [5] M. Boger, D. Veynante, H. Boughanem, A. Trouve´, Proc. Combust. Inst. 27 (1998) 917–925. [6] O. Colin, F. Ducros, D. Veynante, T. Poinsot, Phys. Fluids 12 (7) (2002) 1843–1863. [7] N. Peters, Turbulent Combustion, Cambridge University Press, 2000. [8] R. Borghi, On the Structure and Morphology of Turbulent Premixed Flames, Rec. Adv. Aerosp. Sci. (1985) 117–138.
[9] H. Weller, G. Tabor, A. Gosman, C. Fureby, Proc. Combust. Inst. 27 (1998) 899–907. [10] E.R. Hawkes, R. Cant, Combust. Flame 126 (2001) 1617–1629. [11] V. Moureau, G. Lartigue, Y. Sommerer, C. Angelberger, O. Colin, T. Poinsot, J. Comp. Phys. 202 (2005) 710–736. [12] S. Richard, O. Colin, O. Vermorel, A. Benkenida, D. Veynante, in: ECCOMAS Thematic Conference on Computational Combustion, Lisbon, 2005. [13] A. Trouve´, T. Poinsot, J. Fluid Mech. 278 (1994) 1– 31. [14] S. Pope, Int. J. Eng. Sci. 26 (5) (1988) 445– 469. [15] D. Veynante, T. Poinsot, RANS and LES modeling for turbulent combustion, New tools in turbulent modelling, Les e´ditions de physique, Springer-Verlag, 1997. [16] G. Rymer, Analyse et mode´lisation du taux de re´action moyen et des me´canismes de transport en combustion turbulente pre´me´lange´e, PhD thesis, Ecole Centrale de Paris-FR, 2001. [17] C. Meneveau, T. Poinsot, Combust. Flame 86 (1991) 311–332. [18] F. Charlette, C. Meneveau, D. Veynante, Combust. Flame 131 (2002) 150–180. [19] J.M. Duclos, D. Veynante, T. Poinsot, Combust. Flame 95 (1993) 101–117. [20] T. Poinsot, S. Lele, J. Comp. Phys. 101 (1992) 104– 129. [21] A. Torres, S. Henriot, 3D Modeling of Combustion in Lean Burn Four-valve Engines: Influence of Intake Configuration, Comodia, Yokohama, Japan, 1994, pp.139–156. [22] G. Bruneaux, T. Poinsot, J. Ferziger, J. Fluid Mech. 349 (1997) 199–219.