Evaluation of a combustion model for the simulation of hydrogen spark-ignition engines using a CFD code

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

Available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/he

Evaluation of a combustion model for the simulation of hydrogen spark-ignition engines using a CFD code C.D. Rakopoulos a,*, G.M. Kosmadakis a, E.G. Pariotis b a

Internal Combustion Engines Laboratory, Thermal Engineering Department, School of Mechanical Engineering, National Technical University of Athens, 9 Heroon Polytechniou St., Zografou Campus, 15780 Athens, Greece b Laboratory of Naval Propulsion Systems, Section of Naval Architecture and Marine Engineering, Department of Naval Sciences, Hellenic Naval Academy, End of Hatzikyriakou Ave., Hatzikyriakio, 18539 Piraeus, Greece

article info

abstract

Article history:

The present work deals with the evaluation of a combustion model that has been devel-

Received 12 July 2010

oped, in order to simulate the power cycle of hydrogen spark-ignition engines. The moti-

Received in revised form

vation for the development of such a model is to obtain a simple combustion model with

1 September 2010

few calibration constants, applicable to a wide range of engine configurations, incorporated

Accepted 1 September 2010

in an in-house CFD code using the RNG ke3 turbulence model. The calculated cylinder

Available online 22 September 2010

pressure traces, gross heat release rate diagrams and exhaust nitric oxide (NO) emissions are compared with the corresponding measured ones at various engine loads. The engine

Keywords:

used is a Cooperative Fuel Research (CFR) engine fueled with hydrogen, operating at

Hydrogen

a constant engine speed of 600 rpm. This model is composed of various sub-models used

Combustion

for the simulation of combustion of conventional fuels in SI engines; it has been adjusted in

Spark-ignition

the current study specifically for hydrogen combustion. The basic sub-model incorporated

CFD model

for the calculation of the reaction rates is the characteristic conversion time-scale method,

Laminar

meaning that a time-scale is used depending on the laminar conversion time and the

Turbulent

turbulent mixing time, which dictates to what extent the combustible gas has reached its chemical equilibrium during a predefined time step. Also, the laminar and turbulent combustion velocity is used to track the flame development within the combustion chamber, using two correlations for the laminar flame speed and the Zimont/Lipatnikov approach for the modeling of the turbulent flame speed, whereas the (NO) emissions are calculated according to the Zeldovich mechanism. From the evaluation conducted, it is revealed that by using the developed hydrogen combustion model and after adjustment of the unique model calibration constant, there is an adequate agreement with measured data (regarding performance and emissions) for the investigated conditions. However, there are a few more issues to be resolved dealing mainly with the ignition process and the applicability of a reliable set of constants for the emission calculations. ª 2010 Professor T. Nejat Veziroglu. Published by Elsevier Ltd. All rights reserved.

Abbreviations: ABDC, after bottom dead center; ATDC, after top dead center; BBDC, before bottom dead center; BTDC, before top dead center;
CA, degrees of crank angle; CFD, computational fluid dynamics; CR, compression ratio; EGR, exhaust gas recirculation; EVO, exhaust valve opening; IT, ignition timing; IVC, inlet valve closure; NO, nitric oxide; NOx, nitrogen oxides; rpm, revolutions per minute; SI, spark-ignition; TDC, top dead center. * Corresponding author. Tel.: þ30 210 7723529; fax: þ30 210 7723531. E-mail address: [email protected] (C.D. Rakopoulos). 0360-3199/$ e see front matter ª 2010 Professor T. Nejat Veziroglu. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2010.09.002

12546

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

Nomenclature cp DT,u f h k Lt P Pr qW rk S4 Scr Sh SRNG Sct t T TW Tu u ul ut uT u0 ! u v

1.

specific heat capacity under constant pressure, J kg1 K1 thermal diffusivity of the unburned mixture, m2/s residual gas mass fraction, kgi/kg enthalpy, J/kg turbulent kinetic energy (per unit mass), m2/s2 turbulent integral length scale, mm pressure, N/m2 Prandtl number, e wall heat flux, W/m2 flame kernel radius, m source term source term due to crevice flows source term of the enthalpy equation source term due to RNG turbulence model turbulent Schmidt number, e time, s temperature, K wall temperature, K unburned gas temperature, K velocity along x-axis, m/s laminar flame speed, m/s turbulent flame speed, m/s friction velocity, m/s rms turbulent velocity, m/s velocity vector, m/s velocity along y-axis, m/s

Introduction

The strict emission legislation has forced the automotive industry to take into consideration the use of alternative fuels in internal combustion engines, in order to reduce exhaust emissions with the least impact on engine’s performance. Towards this direction, bio-fuels gain an increasing attention lately, and research is conducted in order to evaluate the performance and emissions from internal combustion engines using such promising fuels [1e4]. Another alternative fuel that seems to be promising for spark-ignition engines is hydrogen, either used in dedicated hydrogen engines [5] or blended with another fuel [6]. Its use in such engines has been experimentally tested by many researchers [7e9], and proper computational models have been developed so far [10e13]. To gain a fundamental understanding of the hydrogen combustion mechanism in internal combustion engines, detailed CFD simulation models should be used. The present study presents a hydrogen combustion model, which is incorporated in an in-house CFD code for the simulation of a hydrogen spark-ignition engine. The combustion model uses recently developed correlations for the hydrogen laminar flame speed [11], while a more detailed description of the spatial flame propagation is accomplished by using local gas properties. An advantage of this model is that only one constant is needed to be tuned within a quite small range, contrary to what it is done in most existing multi-dimensional combustion models [14].

V w Wgi y yþ

cylinder volume, m3 velocity along z-axis, m/s gross indicated work, J distance of the computational node from the wall, m non-dimensional distance from the wall, e

Greek symbols diffusion coefficient, kg m1 s1 G4 3 turbulent dissipation rate (per unit mass), m2/s3 meff effective viscosity, kg m1 s1 n kinematic viscosity, m2/s r density, kg/m3 rb burned gas density, kg/m3 ru unburned gas density, kg/m3 sh coefficient for the enthalpy equation, e coefficient for the dissipation of turbulent kinetic s3 energy equation, e coefficient for the turbulent kinetic energy sk equation, e local wall shear stress, kg m1 s2 sW sc characteristic conversion time, s laminar kinetics time, s sl turbulent mixing time, s st species mass fraction, kgi/kg Yi species chemical equilibrium mass fraction, kgi/kg Yi* f generalized variable 4 equivalence ratio, e

For the evaluation of this model, the calculated results regarding engine performance and nitric oxide emissions are compared with the corresponding measured ones at various engine loads. From the analysis it is concluded that the presented combustion model manages to predict adequately the engine performance within the range of examined equivalence ratios. As far as the nitric oxide exhaust emissions are concerned, the model manages to capture qualitatively correct the trend as equivalence ratio varies, while the absolute values present a small difference from the measured ones.

2.

Experimental data

The experimental data used in the present study, for the evaluation of the developed combustion model, are provided by another research group [8]. These data concern measurements on a Cooperative Fuel Research (CFR) engine, operating at a constant engine speed equal to 600 rpm. Fuel is admitted in the engine through a gas-carburetor in the air-inlet manifold and the fuel rate is controlled with a valve, in order to adjust the equivalence ratio. This means that the load of the engine is altered with the variation of the equivalence ratio (quantified load control). All measurements are conducted with a wide-open throttle (WOT). The main specifications of this engine can be seen in Table 1, where the valve timing events are the ones referenced in [8].

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

Table 1 e CFR engine specifications. Engine model and type

CFR engine, single cylinder, naturally aspirated, four-stroke, water-cooled

Bore Stroke Swept volume Connecting rod length Compression ratio

82.55 mm 114.2 mm 0.6117 l 254 mm Variable (data provided for CR: 7e9.5) 2 (unshrouded) 600 rpm (constant) Variable 5 18
CA ATDC 25
CA ABDC 34
CA BBDC

Valves Engine speed Ignition timing Number of piston rings Valve Inlet valve opening Inlet valve closure timing Exhaust valve opening Exhaust valve closure

Table 2 e Variables that represent the generalized variable f in Eq. (1).

1 u v w h k 3 Yi

3.1.

CFD model

7
CA ATDC

The CFD code developed by the authors can simulate threedimensional curvilinear domains using the finite volume method in a collocated grid. It incorporates the RNG ke3 turbulence model with some slight modifications to introduce the compressibility of a fluid in generalized coordinates, as described in [15]. It solves the following transport equations for the conservation of mass, momentum and energy (Eq. (1)), as suggested in [16]. vðrfÞ þ V,ðr! u fÞ ¼ V,ðGf VðfÞÞ þ Sf þ Scr vt

G4

Continuity u-momentum v-momentum w-momentum Enthalpy Turbulent kinetic energy Dissipation of turbulent kinetic energy Species (i) mass fraction

0 meff meff meff meff/sh meff/sk meff/s3 meff/Sct

the mesh generation model, as well as the detailed evaluation of the developed CFD model under motoring conditions, can be found in previous published papers by the authors [19,20].

The available experimental data concern the inlet mass flow rates of air and hydrogen, the cylinder pressure traces, together with the NOx emissions, when varying the compression ratio, equivalence ratio and ignition timing. Further details concerning the implementation of the experimental facilities can be found in [8].

Numerical model

Equation

Generalized variable (f)

3.2.

3.

12547

(1)

In Table 2, all the variables corresponding to the general transport property f are shown. The source term Scr represents the effect of the crevice model used [17]. Furthermore, for the turbulence model the source terms of the turbulent kinetic energy and its dissipation are the same ones used in [17], where an extra term SRNG has been added to this version of turbulence model [18]. All the constants used in the implementation of the RNG ke3 turbulence model, the enthalpy equation and the calculation of the effective viscosity are the same as in [18]. Moreover, the volumetric heat of combustion, which is inserted in the source term of the enthalpy equation Sh, takes also into account the dissociation of the combustion products. Further details of the in-house CFD model used, concerning the discretization process, the pressure-correction algorithm used, the boundary conditions, the turbulence model,

Heat transfer model

Hydrogen-fueled engines require special attention in the description of the wall heat transfer, due to their thinner thermal boundary layer (shorter quenching distance) and their much higher laminar flame speed, compared to the ones of the conventional-fueled engines [8]. The heat transfer mechanism in this kind of engines has been enlightened by some recent experimental studies, such as the ones in [21,22]; they have measured the wall heat fluxes in the cylinder of hydrogen spark-ignition engines, and concluded that the peak heat flux is increased in this kind of engines as well as the cumulative heat loss to the cylinder walls especially for stoichiometric operation. In the present study the wall heat flux correlation developed by the authors [23] is used for the simulation of the heat transfer mechanism. Its formulation is based on a compressible version of the standard law-of-the-wall [24] and includes the unsteady pressure term, which has been also shown in [25] to give more accurate results. The wall heat flux (in W/m2), which is inserted into the source term of the energy equation Sh for the boundary cells, is shown in Eq. (2). This expression is used for every non-dimensional distance from the wall pffiffiffi pffiffiffiffiffiffiffiffiffiffiffi ky=n the friction velocity. ( yþ ¼ yuT/n) with uT ¼ sW =r ¼ C1=4 m !  .  dP n yþ  40 ruT cp T ln TW T  þ 117:31 dt uT 0:4767 þ 1=Pr      qW ¼ 1 1 1 þ ln y þ  ln 40 þ þ 10:2384 0:4767 0:4767Pr 0:4767Pr (2) Further details concerning the development, implementation and validation of this heat transfer model used in the in-house CFD code can be found in previous works of the authors [23].

3.3.

Crevice model

For the simulation of the crevices a quite simpler phenomenological model than the one found in [26] has been developed and incorporated into the in-house CFD code. The methodology followed is in accordance to the one originally presented

12548

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

in [26], but for the calculation of the mass flow rates between the crevice regions simpler expressions are used. The ring motion is neglected (axial and twist), as well as the flow through the ring grooves. For the calculation of the mass flow rates in every inter-ring region, isentropic compressible flow is assumed and an equivalent surface is introduced [27]. Therefore, the only input required is the number of the piston rings and the estimated/approximated dimensions of the crevice regions, which are actually the equivalent discharge surface and the volume of every inter-ring region. Fewer dimensions are needed than in the original model presented in [26], which makes the implementation of the present model easier. Further details concerning the description, implementation and evaluation of this crevice model can be found in Ref. [17].

3.4.

Combustion model

In the next subsections the combustion model developed by the authors will be presented, as applied for the simulation of hydrogen-fueled, spark-ignition engines. There are various sub-models used in this combustion model to describe the laminar flame speed, the ignition process etc., and for each one of them an analysis follows in order to justify its use.

3.4.1.

Hydrogen laminar burning velocity

There are numerous experimental investigations concerning the hydrogen laminar burning velocity at atmospheric conditions for flame speeds that are stretch-corrected [28,29] or not [30,31]. On the other hand, measurements of the laminar burning velocities for elevated pressures and temperatures (relevant to engine conditions) are quite scarce, due to issues concerning the flame stability in such conditions [8]. Table 3 shows a representative set of published experimental studies at high pressures and temperatures regarding the laminar flame speed [8,30,32e34]. From the various experimental studies shown in Table 3, algebraic correlations have been derived, which can be easily implemented in an engine simulation code [35]. Another approach for obtaining a correlation describing the hydrogen laminar flame speed is through the use of one-dimensional chemical kinetics codes, incorporating one of the many reaction mechanisms of hydrogen-air chemistry [36,37] that can be found in literature [38]. In that case, a stretch-free laminar flame speed is calculated. The authors’ intention is to use in the expression of the turbulent flame speed (see Section 3.4.2) a correlation for the laminar flame speed, which includes stretch effects. In this way it will be avoided to include such effects externally, with the use of appropriate

formulas, as for example the one found in Ref. [39]. For this reason, stretch-free expressions of the hydrogen laminar flame speed are not taken into consideration here. In the developed combustion model of the in-house CFD code a correlation for the laminar flame speed has been incorporated. The authors investigated the features of the available ones in the literature, and selected for further investigation the correlations, which have been derived from measurements under engine-like conditions (wide range of equivalence ratio, high pressure and temperature) [40]. A brief discussion of the correlations shown in Table 3 follows. Beginning with the correlation of Milton and Keck [32], this is rarely used in engine simulations since it corresponds only to stoichiometric operation. The expression of Liu and MacFarlane [30] is also not suitable to be used in CFD simulations, since it corresponds only to ambient pressure and, as it has been shown in Ref. [8], the magnitude of the hydrogen burning velocity depends to some extent on the mixture pressure. Therefore, their use will not be opted and the rest three studies shown in Table 3 are further examined herein as being more suitable, although the correlation of Gerke et al. [34] is derived under conditions matching better the ones observed in real engines. From these experimental studies, correlations regarding the laminar flame speed are obtained, which take into account the stretched flame. The numerical formulas of these three correlations can be found in the Appendix. As discussed in Ref. [12], a correlation should adequately calculate the varying effect of all parameters involved in engines, namely the equivalence ratio, unburned gas temperature, pressure, and residual gas mass fraction. Since the correlation proposed by Iijima and Takeno [33] does not take into account the residual gas mass fraction, it has been modified by adding the relevant term suggested by Verhelst [8]. It should be noted that this formula is the only reliable one, concerning the effect of dilution on the hydrogen laminar burning velocity that can be found in the literature [12]. The effect of the equivalence ratio, the unburned gas temperature and pressure on hydrogen laminar flame speed is shown in Fig. 1aec. These parameters are varied in a range similar to the one encountered in spark-ignition engines. The residual gas mass fraction is kept constant and equal to 10%, which is a realistic value for hydrogen-fueled engines in the absence of external EGR. The correlation of Iijima and Takeno [33] has been derived when the flame is in the cellular region, and includes flame stretch effects, as described in Ref. [5]. Therefore, increased values for the flame speed are expected [41], as shown in Fig. 1aec. On the other hand, the correlation of Verhelst [8]

Table 3 e Representative experimental studies concerning the hydrogen laminar flame speed. Study Milton and Keck [32] Liu and MacFarlane [30] Iijima and Takeno [33] Verhelst [8] Gerke et al. [34]

Equivalence ratio range

Temperature range (K)

Pressure range

f¼1 0.23  f  1.9 0.5  f  4 0.3  f  1 0.4  f  2.5

300  Tu  550 296  Tu  523 291  Tu  500 300  Tu  430 350  Tu  700

0.5  P  7 atm P ¼ 1 atm 0.5  P  25 atm 1  P  10 bar 5  P  45 bar

12549

Hydrogen laminar flame speed (m/s)

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

Correlations Iijima and Takeno [33] (0.5<φ<4 , 291
16

a

14 12

Tu = 800 K P = 20 bar f = 10%

10 8

b

P = 20 bar φ = 0.6 f = 10%

c

Tu = 800 K φ = 0.6 f = 10%

6 4 2 0 0.2

0.4

0.6

0.8

Equivalence ratio

1

400

600

800

1000

1200

5

10

Unburned gas temperature (K)

15

20

25

30

35

40

Pressure (bar)

Fig. 1 e Effect of (a): equivalence ratio, (b): unburned gas temperature, and (c): pressure, on the hydrogen laminar flame speed calculated from various correlations under engine-relevant conditions.

predicts lower laminar flame speed as the equivalence ratio (Fig. 1a) and pressure (Fig. 1c) varies, while this is also true for temperatures higher than approximately 500 K (Fig. 1b). In Fig. 1c, it is observed that beyond the lower pressure validity limit of Gerke et al. [34] expression, a quite steep increase of the laminar flame speed occurs. The authors decided not to use the correlation of Iijima and Takeno [33], in the combustion model, since it does not take into consideration some recent findings concerning the hydrogen laminar flame speed, such as the onset of cellularity [42]. Concerning the last two correlations, in Section 4.3.1, a parametric investigation will take place concerning the simulation of the hydrogen-fueled engine, in order to better identify which correlation of the hydrogen laminar flame speed is more suitable to be used within the particular CFD code developed.

3.4.2.

Hydrogen turbulent burning velocity

There are many expressions found in the literature concerning the turbulent burning velocity of a homogeneous mixture [43e48]. In their majority, they are not developed specifically for hydrogen flame propagation studies. Nevertheless, some general findings of the aforementioned apply to the hydrogen turbulent burning velocity as well. Additionally, in order to decide which expression better matches the specific needs of simulating hydrogen combustion in spark-ignition engines, one should have in mind the numerous relevant review works found in literature. Some of those are of general purpose, such as [43,46,49e51], where various turbulent flame speed correlations are compared with experimental test cases, in order to identify which ones can qualitatively capture the effect of some crucial parameters, such as pressure, Lewis number, turbulence level u0 , etc. Specifically for hydrogen, more details can be found in Refs. [8,11,52].

In Ref. [8], many turbulent flame speed formulas have been implemented in a quasi-dimensional code for the calculation of the closed part of the engine’s cycle. Among the correlations examined is the one proposed by Zimont/Lipatnikov and Chomiak [47,48]. The same correlation is also examined in Ref. [11], where it has been incorporated in a CFD code. In both these works the correlation of Zimont/Lipatnikov and Chomiak is considered adequate to be used, since it is quite simple with few uncertainties in its formulation. More specifically, this turbulent flame speed model has only one calibration factor in its expression, making its implementation and use in a CFD code easier and less time consuming for its tuning. Thus, the authors opted to use this correlation in the hydrogen combustion model developed. The expression describing the turbulent flame speed is given by Eq. (3). ut ¼ Aðu0 Þ

3=4

1=4 1=4 ul1=2 DT;u Lt þ ul

(3)

where A is the only calibration constant with reported values from 0.40 to almost 1 [11,48], u’ the rms turbulent velocity pffiffiffiffiffiffiffiffiffiffiffi given by u0 ¼ 2k=3, ul the stretched laminar flame speed, DT,u the thermal diffusivity of the unburned mixture, and Lt the turbulent integral length scale given by Lt ¼ 0:37ðu0 Þ3 =3. Fig. 2a shows the dependence of the ratio of the turbulent to laminar flame speed (ut/ul) for different equivalence ratios and turbulent intensities u0 under engine-relevant conditions, including also the effect of the residual gas mass fraction f. According to Ref. [53], for near-to-stoichiometric mixtures this ratio takes quite low values as also depicted in Fig. 2a. On the other hand, for lean or very lean mixtures the ratio of turbulent to laminar flame speed increases significantly, and this trend is also captured by the formula of the turbulent flame speed (Eq. (3)) used in the present study. The same correlation concerning the turbulent flame speed has been applied for different unburned gas temperatures, when varying the turbulent intensity. The calculated results are shown in Fig. 2b. The effect of temperature is less

12550

Ratio of turbulent to laminar flame speed (-)

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

u' = 3 m/s u' = 2 m/s u' = 1 m/s u' = 0 m/s (ut=ul)

28 24 20

a

b

Tu = 800 K P = 20 bar f = 10% A = 0.7 (Eq. (3))

16 12

c

φ = 0.4 Tu = 800 K f = 10% A = 0.7 (Eq. (3))

φ = 0.6 P = 20 bar f = 10% A = 0.7 (Eq. (3))

8 4 0 0.4

0.6

0.8

Equivalence ratio

1

400

600

800

1000

1200 5

Unburned gas temperature (K)

10

15

20

25

30

35

40

Pressure (bar)

Fig. 2 e Effect of (a): equivalence ratio, (b): unburned gas temperature, (c): pressure and turbulent intensity, on the ratio of turbulent to laminar flame speed under engine-relevant conditions.

pronounced than the one of the equivalence ratio. This minor dependence is due to the variation of the laminar flame speed with temperature (see. Fig. 1b) and of the mixture thermal diffusivity included in the expression of the turbulent flame speed (see Eq. (3)). The next parameter investigated is the effect of pressure on the turbulent flame speed, which has been pointed out by many researchers to play an important role on the absolute value of the turbulent flame speed [49,53,54]. In Fig. 2c the effect of pressure on the ratio of turbulent to laminar flame speed is depicted corresponding to very lean mixtures (4 equal to 0.4), since for higher equivalence ratios this effect is negligible as also reported in Ref. [53]. It is observed that for low turbulent intensity and elevated pressure the variation of the ratio under investigation is negligible as well, whereas for higher turbulence levels it becomes quite important. It is also concluded that for low pressure (independent of the turbulence level), the effect of pressure is much more pronounced, according to the experimental findings in [53]. From the above it is concluded that Eq. (3) describes adequately the effect of equivalence ratio, unburned gas temperature and pressure, at least qualitatively.

3.4.3.

Ignition phase

A simple approach has been selected for the simulation of the ignition process. It has been shown that during the ignition period, the flame kernel propagates approximately with the laminar flame speed, since its length is considered to be smaller than the prevailing turbulence eddies [11,55,56]. The laminar development of the flame kernel during the ignition phase is simulated here using Eq. (4), as proposed in Ref. [11]. drk ru ¼ ul dt rb

(4)

where rk the radius of the flame kernel, ru the unburned gas density, rb the burned gas density, and ul the laminar flame speed. The temporal discretization of Eq. (4) leads to Eq. (5), which is used to describe the ignition phase of the developed CFD code.

rt ¼ dt ut utl þ rtk rtþ1 k rb

(5)

where (t þ 1) denotes values corresponding to the current time step, while (t) to the previous one, and dt is the time step used. The initial condition used for rk at the start of ignition (t ¼ tign) is equal to 0.5 mm (initial flame kernel diameter equal to d0k ¼ 2r0k ¼ 1 mm [14]), which is approximately equal to the spark-plug gaps more often used in spark-ignition engines. The flame kernel radius is used, in order to track in which computational cells the propagating flame has reached. More details about the process of cell tracking are presented in Section 3.4.4, where the turbulent flame development is given. The end of the ignition process and the transition to the turbulent combustion is determined by the time instant that the flame kernel radius exceeds a specific value (critical radius). According to the two most widely used methods found in the literature dealing with the transition of the combustion process from laminar to turbulent mode, this critical radius is either a constant value [11] or is related to the integral length scale. In the latter case, it is usually taken to be equal to two or three times the integral length scale [14,57], which means that another calibration factor is added in the simulation code. Nevertheless, the calculated values when using both methods are within a small range. In the present study, a steady value for the transition to the turbulent combustion has been applied, since it seems that this method is quite simpler and equally efficient. The value of the critical radius used is equal to rk,crit ¼ 4 mm [11]. According to a parametric study that has been conducted, using the inhouse CFD code, it was revealed that by varying the critical radius, from 3 to 5 mm, no significant influence on the overall combustion process, the performance and the heat release rate was observed. The ignition energy supplied by the spark-plug was simulated by adding energy to the computational cell located at the point where the spark-plug is placed, according to Refs. [58,59]. This energy deposition follows a specific pattern concerning its duration and profile.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

For the present simulation of the CFR engine, which has a constant low rotational speed (600 rpm), the spark duration is equal to 5
CA (almost 1.5 ms). This value has been obtained from similar spark-ignition systems widely used [39,56,60]. The supplied spark energy to the combustible mixture used is 10 mJ, being a quite low value, since it incorporates also the heat energy loss to the electrodes; no additional model has been implemented to simulate this process, due to the great uncertainty of the exact shape, temperature and transient effects near the spark electrodes. Although the value used is not based on any measured value for the specific engine under investigation, a sensitivity analysis that has been conducted revealed that the variation of this value does not have any significant influence on the engine’s performance, and the influence on the heat release rate was actually negligible [60]. Various spark power profiles have also been investigated in the present study, with the one finally used shown in Fig. 3, which is a typical one resembling to some measured ones [60,61].

3.4.4.

Turbulent flame development phase

As was discussed in Section 3.4.3, as the initial flame kernel expands and its radius reaches a critical value (in this case equal to 4 mm), the ignition phase terminates and the main phase of flame propagation initiates. Then, the flame speed is governed mainly by turbulent phenomena, since the turbulent flame speed is quite higher than the laminar one. From that point and then, Eqs. (4) and (5) are not applied. Instead, Eq. (6) is used [11] incorporating the turbulent flame speed ut, whose correlation was presented in a previous subsection (see Eq. (3)). drk ru ¼ ut dt rb

3.4.5.

(7)

The basic sub-model incorporated in the developed CFD code is the characteristic conversion time-scale method. This submodel uses a time-scale, which depends on the laminar conversion time and the turbulent mixing time, and dictates to what extent the combustible gas of each computational cell has reached its chemical equilibrium in a predefined time step.

Spark power (W)

10

Spark power Spark energy

8

8

6

6

4

4

2

2

Total spark energy: 10 mJ Spark duration: ~ 5 0CA

0

0 0

0.25

0.5

0.75

1

1.25

1.5

1.75

Reaction rate calculation

Cumulative spark energy (mJ)

When applying Eq. (7), the radius of the flame is calculated in order to track in which computational cells the flame exists. It should be stressed that for each cell, which is located at the outer part of the flame brush, a different turbulent flame speed is calculated according to the local properties (e.g.

10

pressure, temperature, species mass fractions etc.). By doing so, the flame does not keep its exact initial spherical shape as the flame radius increases [57]. Nevertheless, as it will be shown in a next subsection, the flame shape is close to a spherical one for each case examined, as also noticed experimentally in Ref. [62]. Having determined the local flame propagation speed, it is possible to calculate the burned portion of each cell, defined as the ratio of burned cell volume to the total cell volume. For a better understanding of the methodology used, Fig. 4 gives schematically the main idea of this sub-grid approach implemented in the developed CFD code. The burned ratio of the cells is identified in this figure. To simplify the calculation of the burned volume of the cell, the flame surface within the cell is supposed to be flat. This is not actually true, especially for cells near the spark-plug, because at those locations the flame kernel radius is low. But for cells far away from the spark-plug, this is a fair approximation. This burned part of the cell (corresponding to the calculated burned volume) is taken into account in the source term of the species conservation equations for the computational cells that are located at the flame brush, instead of the whole cell’s volume. This is done, in order to make weaker the effect of grid size in the calculations and, at the same time, to enable the combustion model to calculate numerical solutions quite independent of the total number of computational cells used. As mentioned earlier, the same methodology for flame tracking is also followed during the ignition phase. In that case, the calculated radius taking into consideration only the laminar flame speed is used to track in which cells the flame has reached.

(6)

The temporal discretization of Eq. (6) leads to Eq. (7). rt ¼ dt ut utt þ rtk rtþ1 k rb

12551

2

Time after spark (ms) Fig. 3 e Spark power and energy profile used during the ignition phase.

Fig. 4 e Computational cells with unburned and burned regions.

12552

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

The methodology followed originates from the work of Ref. [63], which has gone through some major improvements ever since [14] and has also been extended to other types of engines. This combustion modeling methodology has been successfully implemented in diesel engines [64,65], in HCCI engines [66,67], as well as in spark-ignition engines [14,58]. Next, the combustion model is presented as was implemented in the developed CFD code. The species considered for the combustion of hydrogen in a spark-ignition engine are: H2, O2, N2, H2O, H, O, N, OH, NO. At each time step, the reaction rates of each of these nine species is calculated from Eq. (8), similar to Ref. [63], except for the nitric oxide one, which is kinetically controlled in engine flows; the calculation of the latter is shown in Section 3.4.6.  t tþ1  ðYi Þt  Yi dYi ¼ dt sc

(8)

where (Yi)t is the species mass fraction at the previous time step, (Yi*)t is the species mass fraction corresponding to the chemical equilibrium value at the previous time step, sc is the characteristic conversion time, and (dYi/dt)tþ1 is the reaction rate of each species at the current time step, which is inserted into the source term of the species transport equations for the cells in which the flame has propagated. The explicit version of Eq. (8) has been chosen in order to calculate at the beginning of every time step the reaction rates, which are then held constant for the calculations during this time step. Therefore, the computational time step during the combustion process should be substantially decreased, in order for the aforementioned methodology to produce reliable results as is also the case followed here. The species concentrations corresponding to chemical equilibrium Yi* are calculated using a developed subroutine [68]. This subroutine has been greatly adjusted and partially reformulated for the specific needs of hydrogen combustion calculations, by omitting all species involved containing carbon atoms, such as CO2 and CO. Its input is the local values of pressure, temperature and species concentrations, while its output is the chemical equilibrium values of all species considered. The characteristic conversion time sc, which is actually the time needed for the burning mixture to reach its chemical equilibrium state, is the same for all nine species considered in this study; it is the sum of a laminar kinetics time sl and a turbulent mixing time st. This characteristic time is calculated from Eq. (9), according to Ref. [63]. sc ¼ sl þ st

(9)

(at medium load) is equal to 0.054, as also suggested in Refs. [14,58], while for leaner mixtures it is higher and for richer lower. The delay coefficient fd found in the expression of the turbulent mixing time (Eq. (11)) is given by Eq. (12), as reported in [14].

fd ¼ 1  eðtts Þ=sd

(12)

where (t  ts) is the time after the spark, and sd ¼ Cm1Lt/ul is the time required by the laminar flame to propagate by two to three times the turbulent integral length scale. The estimated value of the required times is expressed through the Cm1 coefficient, and for the current study it is equal to a constant value of 2.5 for all cases examined, as proposed in Ref. [63]. A parametric investigation has revealed that by increasing the Cm1 coefficient, the timing of the peak pressure is advanced, while it has a weaker effect on the magnitude of the peak pressure. Nevertheless, no tuning of this coefficient has been attempted here.

3.4.6.

NO formation modeling

In Section 3.4.5, it was discussed how the characteristic conversion time needed for a species to reach its chemical equilibrium state is calculated. The aforementioned procedure is followed for all species considered, except for the nitric oxide, for which it has been shown that its reaction rate is kinetically controlled, according to three kinetic reactions, known as the extended Zeldovich mechanism [69]. The methodology followed is the same as the one presented in detail in Ref. [70]. The most usual set of constants for the reaction rates of the three chemical equations considered is also used in the present study, which is the one provided by Ref. [69]. Nevertheless, these constants were derived for hydrocarbon fuels, and their uncertainty (which is already quite high) increases when they are applied for the calculation of the nitric oxide reaction rates in the case of a fuel being carbon-free [12]. The chemical reactions taken into consideration for the calculation of the nitric oxide can be seen in Table 4, together with the constants used.

4.

Results and discussion

4.1.

Test cases examined

The test cases examined in the current study are shown in Table 5. Throughout the Cases 1e3, the equivalence ratio is varied in order to evaluate the combustion model. On the

The laminar kinetics time is given by Eq. (10). sl ¼

DT;u u2l

(10)

Chemical reactions

while the turbulent one is given by Eq. (11). st ¼ cm;f fd

k 3

Table 4 e Chemical reactions and constants used for NO calculations [69].

(11)

The coefficient cm,f shown in Eq. (11) is a function of the local equivalence ratio, as is also explained in Ref. [14], in order to capture better the combustion effects of variable mixture compositions. In the present study, the reference value of cm,f

Forward reaction rates: f f kfi ¼ Afi Tbi eðEi =RTÞ ; i ¼ 1  3 Af (mol m3 s1)

1. 2. 3.

N þ NO 4 N2 þ O N þ O2 4 NO þ O N þ OH 4 NO þ H

3.3  106 6.4  103 4.1  107

bf (e) 0 1 0

Ef (J/mol) 0 26.2  103 0

12553

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

Table 5 e Test cases examined in the present study. Equivalence ratio (4)

Compression ratio

0.71 0.59 0.50 0.59

8 8 8 9.5

1 2 3 4

other hand, Case 4 (with higher compression ratio) is used to comparatively evaluate the two alternative laminar flame speed correlations [8,34], implemented into the combustion model. All cases correspond to constant engine speed (equal to 600 rpm) as well as ignition timing (equal to 15
CA BTDC).

Cylinder pressure (bar)

Case

Measured Calculated (Correlation: [8]) Calculated (Correlation: [34])

60 50

CR=9.5 φ=0.59 IT=15 0CA BTDC

40 30 20 10 0

120 Computational details

In the present study some initial and boundary conditions are provided from the experimental measurements, such as the pressure at IVC, and the inlet air and hydrogen flow rates. For the rest of the required conditions, reliable estimations are used. The wall temperature values are held constant throughout the calculations and their values are estimated according to a recent work in the same CFR engine [21]. The initial flow field is calculated by considering an initial swirl ratio equal to 1.5, since the inlet valve is unshrouded [66]. The initial turbulence properties are the same as the ones derived from a relevant work concerning the same engine [66], and are held constant for every case examined here, since the engine speed is kept constant as well. More specifically, the initial u0 at IVC is equal to 5.16 m/s, giving a rms turbulent velocity equal to around 1 m/s at TDC, although this value slightly changes with engine load (from 0.9 m/s for the lower load to 1.2 m/s for the higher load). The initial gas temperature at IVC is held constant for all cases examined and equal to 370 K [8]. The estimation of the residual gas mass fraction plays an important role in the value of the laminar flame speed, as described in the correlations found in the Appendix. At IVC the initial mixture is composed by H2, O2, N2, and H2O, and the mass content of each of these species is affected by the residual gas mass fraction. Estimations concerning the residual gas and its trend are obtained from the works reported in Refs. [71,72]. No overlapping period of the inlet and exhaust valve exists in the specific engine, making the estimation of the residual gas quite simpler. The exact value of the residual gas mass fraction, wall temperature and initial pressure, at each case examined, is shown in Table 6. Concerning the computational mesh used, previous relevant studies of the present authors [17,23], together with

Table 6 e Initial and boundary conditions for the four cases examined. Case 1 2 3 4

Wall Initial pressure Residual gas temperature (K) at IVC (bar) mass fraction (%) 430 425 415 430

1.079 1.054 1.050 1.059

12.04 9.25 6.7 7.5

150 180 210 240 Crank angle degrees (ABDC)

Fig. 5 e Comparison of the calculated with the measured pressure trace, using two laminar flame speed correlations [8,34].

a parametric study using the firing version of the developed CFD code, have shown that the most reliable mesh, providing adequate accuracy and keeping the required CPU time relative low, is the one having (40  40  40) grid lines along each axis respectively. The methodology for adding or removing mesh layers is the same to the one described in Ref. [23]. The computational time step used is kept low and equal to 0.5
CA, since it has been shown that it is small enough to provide time-step independent solutions [23], even for low rotational speeds. During the combustion period, initiating from the ignition timing and afterwards, the time step is further decreased and becomes equal to 0.1
CA, as in Ref. [11]. This occurs in order to increase the accuracy of the computational results, but also to increase the stability of the code, since during each time step the variation of the cylinder pressure and temperature is large.

90

Gross heat release rate (J/0CA)

4.2.

Experimental Calculated (Correlation: [8]) Calculated (Correlation: [34])

80 70 60

CR=9.5 φ=0.59 IT=15 oCA BTDC

50 40 30 20 10 0 -10 160

170

180

190

200

Crank angle degrees (ABDC) Fig. 6 e Comparison of the calculated gross heat release rates, using two laminar flame speed correlations, with the experimental one.

12554

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

Table 7 e Calculated NO emissions and gross indicated work using the two correlations of the laminar flame speed, and their comparison with the measured values.

Measured Calculated (correlation: [34]) Calculated (correlation: [8])

NO emissions (ppm)

Gross indicated work (J)

2273 2356 2430

308.8 334.2 339.7

4.3. Evaluation of the combustion model with comparison to experimental data Prior to the evaluation of the combustion model, a comparative analysis of the two alternative laminar flame speed correlations described previously is made. This is done in order to decide which of these two correlations is more suitable for our combustion model.

4.3.1.

Investigation of laminar flame speed correlations

The two alternative laminar flame speed correlations considered are the ones proposed in Refs. [8,34]. To comparatively evaluate these two expressions, the CFR engine is examined with an increased compression ratio (Case 4), since in that case the pressure and temperature reach higher values (compared to the rest of the available cases) and the demand for a reliable correlation is more intense. Fig. 5 shows the predicted pressure traces using the two correlations, and their comparison with the measured one. When using the Verhelst correlation, the calibration constant of the turbulent flame speed correlation (Eq. (3)) is set equal to 1.6, which is substantially higher than the corresponding value when the Gerke et al. correlation is used (equal to 0.8), in order to better match the measured in-cylinder pressure trace. The calibration constant value used in the Verhelst correlation is also much higher than the reported values from other researchers, which are held below unity. Nevertheless, the calculated pressure history is well captured using both correlations.

Fig. 6 shows a comparison of the predicted gross heat release rates using the two alternative laminar flame speed correlations with the corresponding experimental one. The calculation of the experimental and numerical gross heat release rate is based on the cylinder pressure traces, using a single-zone approach [21]. The CFD code using both correlations predicted similar gross heat release rates, although the calibration constant has very different values. In general, the predicted gross heat release rate is close to the experimental one (qualitatively), while the timing of the peak gross heat release is slightly advanced in comparison to the measured one. Also, the calculated combustion duration is well captured with both correlations used, although it is slightly overestimated (approximately by 3
CA). The disadvantage observed when using the Verhelst correlation, is that the combustion rate during the ignition process is greatly underestimated, increasing afterwards significantly during the main flame propagation period, until approximately 175
CA ABDC, due to the high value of the tuning constant used. The change of the slope of the calculated gross heat release rate during the expansion stroke will be explained in Section 4.3.2. Moreover, Table 7 shows the exhaust NO emissions and the gross indicated work when using the two correlations, and their comparison with the measured values. The gross indicated work is calculated during the closed part of the cycle (from the IVC until the EVO) and is given by Eq. (13). ZEVO Wgi ðJÞ ¼

pdV

(13)

IVC

where V the instantaneous cylinder volume. The calculated nitric oxide emissions using both correlations give a value very close to the measured one (discrepancy smaller than 7%). On the other hand, using both correlations, the calculated gross indicated work is higher than the measured one, since during the expansion stroke the cylinder pressure is slightly over-predicted in comparison to the measured one (see Fig. 5). This may be attributed to the lower heat loss calculated during this time period [23]. Nevertheless,

Cylinder pressure (bar)

50

Measured Calculated

40

30

φ=0.59

φ=0.71

φ=0.50

20

10

0 120

150

180

210

240 120

150

180

210

240 120

150

180

210

240

Crank angle degrees (ABDC) Fig. 7 e Comparison of calculated pressure traces with the measured ones at various engine loads (CR [ 8, IT [ 15
CA BTDC).

12555

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

Gross heat release rate (J/oCA)

100 Experimental Calculated

80 60

φ=0.71

φ=0.59

φ=0.50

40 20 0 160

170

180

190

200 160

170

180

190

200 160

170

180

190

200

Crank angle degrees (ABDC) Fig. 8 e Comparison of calculated gross heat release rates with the experimental ones for three equivalence ratios (CR [ 8, IT [ 15
CA BTDC).

4.3.2.

Simulation of engine performance and emissions

In the current subsection, Cases 1e3 of Table 5 will be elaborated using the Gerke et al. correlation for the laminar burning velocity [34]. Fig. 7 shows the comparison of the calculated with the measured cylinder pressure traces at various engine loads. As shown, the peak pressure is well captured in all cases examined, as well as its timing, with the only exception for the richer mixture case (4 ¼ 0.71) where this timing is estimated almost 3
CA later. It should be stressed that the model manages to predict adequately the absolute value of the cylinder pressure at the IT, in all cases examined, without any constant tuning. This is crucial for the proper estimation of the cylinder pressure during the combustion phase, and is partly attributed to the further development of the CFD code by including a crevice model and a detailed wall-function formulation [17,23]. An aspect that needs further investigation is the slight pressure over-prediction during the expansion stroke, which becomes more pronounced for the richer mixtures and, to some extent, can be attributed to the relatively low wall heat flux observed during the expansion stroke, as also described in Ref. [23]. In Fig. 8, the calculated gross heat release rate is compared with the experimental one for Cases 1e3. As the cylinder charge becomes leaner, the peak gross heat release is

decreased rapidly, while the combustion duration is increased. The latter is expected, owing to the decreased hydrogen turbulent combustion velocity for lean mixtures. The change of the slope of the calculated gross heat release rate, during the expansion stroke, may be attributed to the retarded combustion of the outer computational cells adjacent to the cylinder liner. It has been observed that the propagating flame reaches these regions just before the end of combustion, where the prevailing unburned gas temperature and turbulence level is low. Thus, combustion of these specific cells is delayed, leading to rapid ignition and combustion of all these cells later on when the local conditions are proper. This causes an increase of the heat release for a few degrees CA. Nevertheless, further investigation is needed to understand better the reason for the change of slope in the gross heat release rate. This comparison reveals the capability of the CFD code not only to capture the effect of equivalence ratio on the cylinder pressure history of a spark-ignition engine, but also to simulate reliably the combustion processes taking place.

6000

Measured Calculated

5000

NO emissions (ppm)

this discrepancy is lower than 10%. It is to be understood in this work that the calculated values concern nitric oxide (NO), while the available measured ones nitrogen oxides (NOx); however, as known, their difference at a spark-ignition engine exhaust is extremely small. From the former extended evaluation, it has been revealed that both correlations can be used for reliable engine simulations using the developed CFD code. Nevertheless, the use of the correlation of Gerke et al. [34] is preferred for the combustion model, since it is derived from conditions more relevant to real engine operation, the calibration factor of the turbulent flame speed formula is within the reported values [11], and the ignition process is simulated better as shown previously.

4000 CR=8 IT=15 0CA BTDC 3000 2000 1000 0 0.45

0.5

0.55

0.6

0.65

0.7

0.75

Equivalence ratio Fig. 9 e Calculated NO emissions and comparison with measurements, when varying the equivalence ratio.

12556

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

Fig. 10 e Flame images for the three equivalence ratios examined (CR [ 8, IT [ 15
CA BTDC).

Fig. 9 shows the comparison of the calculated nitric oxide emissions with the measured tailpipe ones for all cases examined. As shown, the effect of varying the equivalence ratio is well captured by the simulation model, although in absolute values there is a small discrepancy, approximately 26% at high load and 50% at low load. However, in the latter case, the absolute value of the measured nitric oxide emissions is very low. The calculated

results would be more accurate, if a specially derived set of constants for carbon-free fuels were used for the Zeldovich mechanism, according to Ref. [12]. The flame development images are shown in Fig. 10, depicting how the flame propagates inside the cylinder after the ignition timing for each case examined. It is interesting to notice the difference of the flame development for the various

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

12557

Fig. 11 e NO mass fraction spatial distribution during the combustion period for the three equivalence ratios (Cases 1e3).

equivalence ratios. It should be noted that the spark-plug is offset from the cylinder axis, as one can guess from the location of the ignition kernel. Owing to the high burning velocity of hydrogen, the flame has reached the liner for the rich mixture case (Case 1) at around TDC, while for the leaner cases (Cases 2 and 3) the flame at that moment has also developed significantly. Moreover, the flame has covered the whole combustion chamber at around TDC for the richer mixture case (Case 1), at 10
CA ATDC for the mixture of Case 2, and at 15
CA ATDC for the leaner mixture of Case 3. Fig. 11 shows the spatial distribution of the NO mass fraction inside the cylinder, at various time instants, for all equivalence ratios examined at three planes (the top-right

plane crosses the spark-plug). The contour values are given as mass fractions (kgNO/kgmixture) and their range is different for every case, in order to better visualize the emissions production process. As shown in Fig. 11, as the equivalence ratio is increased NO emissions are radically increased. In all cases examined, during the first stages of combustion, NO emissions are produced near the spark-plug where the local temperature is also higher [70]. However, as combustion continues, the core of the high NO concentration moves towards the center of the cylinder, primarily due to the cooling of the outer regions and secondarily to the prevailing velocity field. This is an interesting finding revealing the importance of the proposed CFD code, since in the simple quasi-dimensional or multi-zone

12558

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

models the total heat losses are equally distributed to the incylinder zones, irrespective of their location (close or far away from the cylinder walls), and also no convection phenomena are considered there.

5.

Summary and conclusions

In the present work a combustion model has been developed, which has been incorporated in a CFD code for the simulation of a hydrogen-fueled spark-ignition engine. This combustion model takes into consideration some recent findings concerning the hydrogen laminar flame speed, and it uses a robust turbulent flame speed correlation in order to track the flame-containing computational cells. In these cells the characteristic conversion time-scale method is applied, in order to calculate the combustion rate. By doing so, the flame propagation is obtained within the engine cylinder in order to calculate the performance and NO emissions of the engine. The study showed that both hydrogen laminar flame speed correlations investigated in detail, proposed by Verhelst [8] and Gerke et al. [34], could be incorporated in the CFD code for the simulation of the hydrogen engine, giving reasonably acceptable results. However, the use of the latter correlation has been preferred, since it was based on a dataset in a wider pressure and temperature range closer to real engine operating conditions. Also, using this correlation, the ignition process is better captured. The combustion model has been applied to simulate the power cycle of the same engine at various engine loads. The calculated results (regarding cylinder pressure traces, gross heat release rates, and exhaust nitric oxide emissions) have been compared with the measured ones showing a good match. It should be noted that throughout this investigation, the only calibration constant used in the model has preserved its constancy of value. However, some aspects need further investigation as follows:  The set of constants used for the calculation of the nitric oxide emissions.  The underestimated gross heat release rate at the beginning of combustion, and the change of the slope of the gross heat release rate diagram near the end of combustion.  The slightly overestimated cylinder pressure trace during the expansion stroke. In any case, the results obtained so far are encouraging, by showing that the simulation of hydrogen-fueled engines using a detailed simulation numerical tool, such as the one

presented in the current work, can yield reliable results, without the need for recourse to a large set of calibration constants. However, a more extended investigation should be conducted in the future covering a wider range of operating conditions, to justify the validity of the model and its ability to capture combustion phenomena without adjusting the value of the calibration constant.

Acknowledgements The authors wish to thank Prof. S. Verhelst (Univ. of Gent, Belgium) for the provision of the experimental data. His contribution is greatly acknowledged. Also, G.M. Kosmadakis wishes to thank the Greek State Scholarships Foundation for granting him a Ph.D. research scholarship.

Appendix. The correlations for the hydrogen laminar flame speed mentioned in the present study are given here. In every correlation the effect of the residual gas mass fraction has been added, according to Verhelst [8], where f is the residual gas mass fraction. 1) Iijima and Takeno correlation [33], valid for 0.5  f  4, 291  Tu  500 K, 0.5  P  25 atm: h  i a Tu ð1  gf Þ ul ¼ ul0 1 þ b log PP0 T0 ul0 ¼ 2:98  ðf  1:7Þ2 þ0:32ðf  1:7Þ3 a ¼ 1:54 þ 0:026ðf  1Þ b ¼ 0:43 þ 0:003ðf  1Þ g ¼ 2:715  0:5f P0 ¼ 1 atm; T0 ¼ 291 K 2) Verhelst correlation [8], valid for 0.3  f  1, 300  Tu  430 K, 1  P  10 bar:

ul ¼ ul0

 a  b Tu T0

P P0

ð1  gf Þ

ul0 ¼ 4:77f þ 8:65f2  0:394f  0:296 a¼ 1:232 2:9f3  6:69f2 þ 5:06f  1:16 f < 0:6 b¼ 0:0246f þ 0:0781 f  0:6 g ¼ 2:715  0:5f P0 ¼ 5 bar; T0 ¼ 365 K 3

3) Gerke et al. correlation [34], valid for 0.4  f  2.5, 350  Tu  700 K, 10  P  45 bar:

 a  b P ul ¼ ul0 TTu0 ð1  gf Þ P0

6 0:25f  3:4774f5 þ 18:498f4  46:525f3 þ 52:317f2  13:976f þ 1:2994 ul0 ¼ 10:41096  .  1:272384ðf  2:5Þ f > 2:5 a ¼ 0:0163 1 f þ 2:2937  . b ¼ 0:2037 1 f  0:575 g ¼ 2:715  0:5f P0 ¼ 20 bar; T0 ¼ 600 K

0:4  f  2:5

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

references

[1] Agarwal AK. Biofuels (alcohols and biodiesel) applications as fuels for internal combustion engines. Prog Energy Combust Sci 2007;33(3):233e71. [2] Rakopoulos CD, Michos CN. Generation of combustion irreversibilities in a spark ignition engine under biogasehydrogen mixtures fueling. Int J Hydrogen Energy 2009;34(10):4422e37. [3] Rakopoulos CD, Rakopoulos DC, Hountalas DT, Giakoumis EG, Andritsakis EC. Performance and emissions of bus engine using blends of diesel fuel with bio-diesel of sunflower or cottonseed oils derived from Greek feedstock. Fuel 2008;87(2):147e57. [4] Tsolakis A, Megaritis A, Yap D. Application of exhaust gas fuel reforming in diesel and homogeneous charge compression ignition (HCCI) engines fuelled with biofuels. Energy 2008;33(3):462e70. [5] Verhelst S, Wallner T. Hydrogen-fueled internal combustion engines. Prog Energy Combust Sci 2009;35(6):490e527. [6] Dimopoulos P, Bach C, Soltic P, Boulouchos K. Hydrogenenatural gas blends fuelling passenger car engines: combustion, emissions and well-to-wheels assessment. Int J Hydrogen Energy 2008;33(23):7224e36. [7] Kirchweger W, Haslacher R, Hallmannsegger M, Gerke U. Applications of the LIF method for the diagnostics of the combustion process of gas-IC-engines. Exp Fluids 2007;43 (2e3):329e40. [8] Verhelst S. A study of the combustion in hydrogen-fuelled internal combustion engines. Ph.D. thesis, Ghent University, Ghent, Belgium; 2005. [9] Rottengruber H, Berckmueller M, Elsaesser G, Brehm N, Schwarz C. Direct-injection hydrogen SI-engine operation strategy and power density potentials. Trans SAE J Fuels Lubricants 2004;113:1749e61 [SAE Paper no. 2004-01-2927]. [10] Dimopoulos P, Rechsteiner C, Soltic P, Laemmle C, Boulouchos K. Increase of passenger car engine efficiency with low engine-out emissions using hydrogenenatural gas mixtures: a thermodynamic analysis. Int J Hydrogen Energy 2007;32(14):3073e83. [11] Gerke U. Numerical analysis of mixture formation and combustion in a hydrogen direct-injection internal combustion engine. Ph.D. thesis, Diss. ETH no. 17477, Cuvillier Go¨ttingen, ETH Zurich, Switzerland; 2007. [12] Knop V, Benkenida A, Jay S, Colin O. Modelling of combustion and nitrogen oxide formation in hydrogen-fuelled internal combustion engines within a 3D CFD code. Int J Hydrogen Energy 2008;33(19):5083e97. [13] Verhelst S, Sierens R. A quasi-dimensional model for the power cycle of a hydrogen-fuelled ICE. Int J Hydrogen Energy 2007;32(15):3545e54. [14] Fan L, Reitz RD. Development of an ignition and combustion model for spark-ignition engines. Trans SAE J Engines 2000; 109:1977e89 [SAE Paper no. 2000-01-2809]. [15] Thakur S, Wright JR. A multiblock operator-splitting algorithm for unsteady flows at all speeds in complex geometries. Int J Numer Methods Fluids 2004;46(4):383e413. [16] Shyy W. Computational modeling for fluid flow and interfacial transport. Amsterdam: Elsevier Science; 1994. [17] Rakopoulos CD, Kosmadakis GM, Dimaratos AM, Pariotis EG. Investigating the effect of crevice flow on internal combustion engines using a new simple crevice model implemented in a CFD code. Appl Energy; 2010;. doi:10.1016/j. apenergy.2010.07.012. [18] Han Z, Reitz RD. Turbulence modeling of internal combustion engines using RNG ke3 models. Combust Sci Technol 1995;106(4):267e95.

12559

[19] Rakopoulos CD, Kosmadakis GM, Pariotis EG. Evaluation of a new computational fluid dynamics model for internal combustion engines using hydrogen under motoring conditions. Energy 2009;34(12):2158e66. [20] Rakopoulos CD, Kosmadakis GM, Pariotis EG. Investigation of piston bowl geometry and speed effects in a motored HSDI diesel engine using a CFD against a quasi-dimensional model. Energy Convers Manage 2010;51(3):470e84. [21] Demuynck J, Raes N, Zuliani M, De Paepe M, Sierens R, Verhelst S. Local heat flux measurements in a hydrogen and methane spark ignition engine with a thermopile sensor. Int J Hydrogen Energy 2009;34(24):9857e68. [22] Shudo T, Nabetani S. Analysis of degree of constant volume and cooling loss in a hydrogen fuelled SI engine. Trans SAE J Fuels Lubricants 2001;110:1911e7 [SAE Paper no. 2001-01-3561]. [23] Rakopoulos CD, Kosmadakis GM, Pariotis EG. Critical evaluation of current heat transfer models used in CFD incylinder engine simulations and establishment of a comprehensive wall-function formulation. Appl Energy 2010;87(5):1612e30. [24] Launder BE, Spalding DB. The numerical computation of turbulent flows. Comput Methods Appl Mech Eng 1974;3(2): 269e89. [25] Han Z, Reitz RD. A temperature wall function formulation for variable-density turbulent flows with application to engine convective heat transfer modeling. Int J Heat Mass Transfer 1997;40(3):613e25. [26] Namazian M, Heywood JB. Flow in the pistonecylinderering crevices of a spark-ignition engine: effect on hydrocarbon emissions, efficiency and power. Trans SAE (Section 1) 1982; 91:261e88 [SAE Paper no. 820088]. [27] Keribar R, Dursunkaya Z, Flemming MF. An integrated model of ring pack performance. Trans ASME J Eng Gas Turbines Power 1991;113:382e9. [28] Vagelopoulos CM, Egolfopoulos FN, Law CK. Further considerations on the determination of laminar flame speeds with the counterflow twin-flame technique. Symp (Int) Combust [Proc] 1994;25(1):1341e7. [29] Kwon OC, Faeth GM. Flame/stretch interactions of premixed hydrogen-fueled flames: measurements and predictions. Combust Flame 2001;124(4):590e610. [30] Liu DDS, MacFarlane R. Laminar burning velocities of hydrogeneair and hydrogeneairesteam flames. Combust Flame 1983;49(1e3):59e71. [31] Koroll GW, Kumar RK, Bowles EM. Burning velocities of hydrogeneair mixtures. Combust Flame 1993;94(3):330e40. [32] Milton B, Keck J. Laminar burning velocities in stoichiometric hydrogen and hydrogenehydrocarbon gas mixtures. Combust Flame 1984;58(1):13e22. [33] Iijima T, Takeno T. Effects of temperature and pressure on burning velocity. Combust Flame 1986;65(1):35e43. [34] Gerke U, Steurs K, Rebecchi P, Boulouchos K. Derivation of burning velocities of premixed hydrogen/air flames at engine-relevant conditions using a single-cylinder compression machine with optical access. Int J Hydrogen Energy 2010;35(6):2566e77. [35] Verhelst S, Woolley R, Lawes M, Sierens R. Laminar and unstable burning velocities and Markstein lengths of hydrogeneair mixtures at engine-like conditions. Proc Combust Inst 2005;30(1):209e16. [36] Marinov NM, Curran HJ, Pitz WJ, Westbrook CK. Chemical kinetic modeling of hydrogen under conditions found in internal combustion engines. Energy Fuels 1998;12:78e82. [37] O’Conaire M, Curran HJ, Simmie JM, Pitz WJ, Westbrook CK. A comprehensive modeling study of hydrogen oxidation. Int J Chem Kinetics 2004;36:603e22. [38] Konnov AA. Remaining uncertainties in the kinetic mechanism of hydrogen combustion. Combust Flame 2008;152(4):507e28.

12560

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 5 ( 2 0 1 0 ) 1 2 5 4 5 e1 2 5 6 0

[39] Herweg R, Maly RR. A fundamental model for flame kernel formation in S.I engines. SAE paper no. 922243; 1992. [40] Hu E, Huang Z, He J, Miao H. Experimental and numerical study on laminar burning velocities and flame instabilities of hydrogeneair mixtures at elevated pressures and temperatures. Int J Hydrogen Energy 2009;34(20):8741e55. [41] Chen JH, Im HG. Stretch effects on the burning velocity of turbulent premixed hydrogen/air flames. Proc Combust Inst 2000;28(1):211e8. [42] Bradley D, Lawes M, Liu K, Verhelst S, Woolley R. Laminar burning velocities of lean hydrogeneair mixtures at pressures up to 1.0 MPa. Combust Flame 2007;149(1e2): 162e72. [43] Duclos JM, Veynante D, Poinsot T. A comparison of flamelet models for premixed turbulent combustion. Combust Flame 1993;95(1e2):101e17. [44] Weller HG, Uslu S, Gosman AD, Maly RR, Herweg R, Heel B. Prediction of combustion in homogeneous-charge sparkignition engines. In: Proceedings of 3rd international symposium on “Diagnostics and modelling of combustion in internal combustion engines (COMODIA 1994)”, Yokohama, Japan; July 1994. p. 163e9. [45] Peters N. The turbulent burning velocity for large-scale and small-scale turbulence. J Fluid Mechanics 1999;384(1): 107e32. [46] Guelder OL. Turbulent premixed flame propagation models for different combustion regimes. Symp (Int) Combust [Proc] 1990;23(1):743e50. [47] Zimont VL. Gas premixed combustion at high turbulence. Turbulent flame closure combustion model. Exp Therm Fluid Sci 2000;21(1e3):179e86. [48] Lipatnikov AN, Chomiak J. A simple model of unsteady turbulent flame propagation. Trans SAE J Engines 1997;106: 2441e52 [SAE Paper no. 972993]. [49] Lipatnikov AN, Chomiak J. Modeling of pressure and nonstationary effects in spark ignition engine combustion: A comparison of different approaches. Trans SAE J Engines 2000;109:1833e50 [SAE Paper no. 2000-01-2034]. [50] Dinkelacker F, Hoelzler S. Investigation of a turbulent flame speed closure approach for premixed flame calculations. Combust Sci Technol 2000;158(1):321e40. [51] Lipatnikov AN. Testing premixed turbulent combustion models by studying flame dynamics. Int J Spray Combust Dynamics 2009;1(1):39e66. [52] Abdel-Gayed RG, Bradley D. Dependence of turbulent burning velocity on turbulent Reynolds number and ratio of laminar burning velocity to r.m.s. turbulent velocity. Symp (Int) Combust [Proc] 1976;16:1725e35. [53] Kitagawa T, Nakahara T, Maruyama K, Kado K, Hayakawa A, Kobayashi S. Turbulent burning velocity of hydrogen-air premixed propagating flames at elevated pressures. Int J Hydrogen Energy 2008;33(20):5842e9. [54] Brandl A, Pfitzner M, Mooney JD, Durst B, Kern W. Comparison of combustion models and assessment of their applicability to the simulation of premixed turbulent combustion in IC-engines. Flow Turbul Combust 2005;75 (1e4):335e50. [55] Willems H, Sierens R. Modeling the initial growth of the plasma and flame kernel in SI engines. Trans ASME J Eng Gas Turb Power 2003;125(2):479e84. [56] Pischinger S, Heywood JB. A model for flame kernel development in a spark-ignition engine. Symp (Int) Combust [Proc] 1991;23(1):1033e40.

[57] Tan Z, Reitz RD. An ignition and combustion model based on the level-set method for spark-ignition engine multidimensional modeling. Combust Flame 2006;145(1e2): 1e15. [58] Kuo T-W. Multidimensional port-and-cylinder gas flow, fuel spray, and combustion calculations for a port-fuel-injection engine. SAE paper no. 920515; 1992. [59] Yossefi D, Maskell SJ, Ashcroft SJ, Belmont MR. Ignition source characteristics for natural-gas-burning vehicle engines. Proc Instn Mech Engrs Part D J Automob Eng 2000; 214(2):171e80. [60] Arcoumanis C, Bae C- S. Correlation between spark ignition characteristics and flame development in a constant-volume combustion chamber. Trans SAE J Engines 1992;101:556e70 [SAE Paper no. 920413]. [61] Aleiferis PG, Taylor AMKP, Ishii K, Urata Y. The relative effects of fuel concentration, residual-gas fraction, gas motion, spark energy and heat losses to the electrodes on flame-kernel development in a lean-burn spark ignition engine. Proc Inst Mech Engrs Part D J Automob Eng 2004;218 (4):411e25. [62] Heywood JB, Vilchis FR. Comparison of flame development in a spark-ignition engine fueled with propane and hydrogen. Combust Sci Technol 1984;38(5):313e24. [63] Abraham J, Bracco FV, Reitz RD. Comparison of computed and measured premixed charge engine combustion. Combust Flame 1985;60(3):309e22. [64] Kong S-C, Han Z, Reitz RD. The development and application of a diesel ignition and combustion model for multidimensional engine simulations. Trans SAE J Engines 1995;104:502e18 [SAE Paper no. 950278]. [65] Masood M, Ishrat MM, Reddy AS. Computational combustion and emission analysis of hydrogen-diesel blends with experimental verification. Int J Hydrogen Energy 2007;32(13): 2539e47. [66] Kong S-C, Ayoub N, Reitz RD. Modeling combustion in compression ignition homogeneous charge engines. Trans SAE J Engines 1992;101:896e911 [SAE Paper no. 920512]. [67] Kim S, Ito K, Yoshihara D, Wakisaka T. Application of a genetic algorithm to the optimization of rate constants in chemical reaction submodels for engine combustion simulation. In: Proceedings of 6th international symposium on “Diagnostics and modelling of combustion in internal combustion engines (COMODIA 2004)”, Yokohama, Japan, August 2004, pp. 43e51. [68] Rakopoulos CD, Hountalas DT, Tzanos EI, Taklis GN. A fast algorithm for calculating the composition of diesel combustion products using 11 species chemical equilibrium scheme. Adv Eng Software 1994;19(2):109e19. [69] Lavoie GA, Heywood JB, Keck JC. Experimental and theoretical study of nitric oxide formation in internal combustion engines. Combust Sci Technol 1970;1(4):313e26. [70] Rakopoulos CD, Michos CN, Giakoumis EG. Thermodynamic analysis of SI engine operation on variable composition biogasehydrogen blends using a quasi-dimensional, multizone combustion model. SAE Int J Engines 2009;2(1):880e910 [SAE Paper no. 2009-01-0931]. [71] Fox JW, Cheng WK, Heywood JB. A model for predicting residual gas fraction in spark-ignition engines. Trans SAE J Engines 1993;102:1538e44 [SAE Paper no. 931025]. [72] Senecal PK, Xin J, Reitz RD. Predictions of residual gas fraction in IC engines. Trans SAE J Engines 1996;105:2243e54 [SAE Paper no. 962052].