**3.1 Introduction and overview**

274 Computational Simulations and Applications

The slight difference between the calculated and the measured brake torque at the transient end (Fig. 6f) can primarily be ascribed to an underestimation of the gas pressure

(Baratta & Spessa, 2009) modified the previous model with reference to the engine installation on a commercial vehicle for urban transportation. The model was extensively revised by modifying the pertinent pipe, bend and 'flowsplit' objects. In addition, the following differences were considered with respect to the dyno test-bed configuration: - The bus IC was of the air-to-air type, and displayed different volumes and temperatures than the water-to-air IC installed on the dyno test rig. This influenced the engine volumetric efficiency, the charge temperature and the wave propagation phenomena within the intake system. Therefore, it was necessary to revise the IC calibration, with specific reference to flow losses and heat transfer effects (surface roughness and


The new model was tuned on the basis of a specific test, in which the hydraulic torque converter of the considered urban bus was kept under stall conditions by means of the vehicle brakes, while a quick opening of the throttle valve was actuated. Throughout this transient process, the engine torque demand was thus proportional to N2. As an example, Fig. 7 compares the experimental (black solid line) and simulated (red dotted line) timehistories of the boost pressure, mass flow rate and engine speed for a load step from N / Nmax = 0.3 to 0.8. The same EM profile versus PR was used as in (Baratta et al., 2010).

The reliability of the 1-D approach can be improved if predictive 0-D combustion models are used to predict the heat-release rate within the engine combustion chamber. The turbulent combustion process is a complex phenomenon that involves many chemical, thermodynamic and fluid-dynamic aspects, which should be studied by adopting a three-dimensional approach. However, as discussed in great detail by (Lipatnikov & Chomiak, 2002), even in this case, the development of a fundamentally substantiated model, that is, a model which is based only on the application of 'first principles', is very difficult. A possible, practical solution is that of shifting from the first principles to phenomenology, i.e., in the use of well established experimental facts and approximate descriptions of selected combustion-process characteristics which are assumed to be the main controlling factors. Similarly, predictive 0-D combustion models, which are the topic of this section, are based on a phenomenological description of the turbulent combustion process of a premixed fuelair mixture. Although they generally need a preliminary tuning procedure, they can potentially predict the dependence of the heat-release rate on, among other factors, incylinder flow, combustion chamber geometry, mixture composition, thermodynamic state, and spark timing. Since the pioneering work of (Blizzard & Keck, 1974), a large number of papers have been published, which have focused on the development and/or the application of predictive combustion models to SI engines. A rather good review of the main aspects that have to be faced in a thermodynamic combustion model formulation can be

contribution to the friction mean effective pressure under full-load operations.

material, and the related multipliers).

estimation of the engine brake torque.

**3. Predictive 0-D combustion models** 

found in (Velherst & Sheppard, 2009).

The goal of a predictive combustion model is to predict the rate at which the unburned mixture is converted into burned gases. This allows the computation of the in-cylinder pressure through Eqs. (4-6). The different models are based on the definition of a 'turbulent burning velocity', Sb, and of a flame burning-front area Abf, whereas the flame-brush thickness is generally neglected. The flame area is often modelled assuming a spherical shape of the flame front, which gradually intersects the combustion chamber surfaces as it grows (see, among others, Baratta et al., 2008; Bozza et al., 2005; Wahiduzzaman et al., 1993). This assumption has been confirmed by experiments, at least for combustion chambers with sufficiently low swirl and tumble ratios. In a thermo-dynamic modelling approach, this is also the most reasonable a priori choice. In fact, a sub-model for the flame deformation by the in-cylinder flow would need detailed information on the flow motion characteristics, which is not compatible with the thermo-dynamic nature of the overall model.

The evolution equation for the burned-gas mass fraction xb can be derived adopting two different approaches (Baratta et al ., 2006).

Numerical Simulation Techniques for the

heat-transfer coefficient.

mean flow quantities:

combustion phase (Bozza et al., 2005).

Prediction of Fluid-Dynamics, Combustion and Performance in IC Engines Fuelled by CNG 277

used by the in-cylinder heat-transfer sub-model, in order to derive a value for the convective

A rather simple, but quite widely applied model is the K-k one (Poulos & Heywood, 1983), which is based on a zero-dimensional energy cascade from the mean flow to the viscous eddy dissipation. According to such a model, the rates of change in the mean flow kinetic energy ( <sup>2</sup> 1 2 *K mU cyl* ) and in the turbulent kinetic energy ( <sup>2</sup> 32 ' *cyl k mu* ) are, respectively:

*<sup>o</sup> i i*

*<sup>o</sup> cyl*

*cyl*

*cyl*

1 2

(18)

*i*

*L m* 

*dt <sup>m</sup>* (16)

(17)

'3

is the

*i u L* 

1 <sup>2</sup> 2

where: U is the average mean-flow velocity, u' is the rms turbulent velocity,

*dK <sup>m</sup> mv P K*

*dk <sup>m</sup> Pm k dt <sup>m</sup>* 

rate of turbulent kinetic energy dissipation per unit mass, Li is the characteristic size of the large-scale eddies, m is the mass in the cylinder, *mi* is the mass flow rate into the cylinder, *mo* is the mass flow rate out of the cylinder, *<sup>i</sup> v* is the jet velocity entering the chamber and P is the rate of turbulent kinetic energy production. Since the turbulence model does not spatially resolve the flow parameters, the production term P is empirically estimated from

0.3307

*K k P c*

where the turbulent dissipation constant c is adjusted to give the expected profiles of u' and U throughout the whole cycle. It is worth pointing out that the accuracy of a turbulent burning-speed model depends to a great extent on the unburned-gas turbulence evolution during combustion, which is normally thought to present a maximum near firing TDC (Baratta et al., 2008, Bozza et al, 2005, Heywood, 1988; Poulos & Heywood, 1983; Yoshiyama et al., 2001). Although the magnitude of such a maximum depends, among other factors, on the configuration of the combustion chamber and, to a lesser extent, on the intake system, its presence is generally modelled by means of heuristic approaches that add turbulence during compression and combustion. As an example, it can be seen that, in (Baratta et al., 2006, 2008; Poulos & Heywood, 1983), as soon as the combustion starts, Eq. (17) is no longer integrated, and the evolution of the turbulence intensity and of the spatial macro scale is calculated by assuming conservation of angular momentum for large scale eddies, which gives rise to an increase in the unburned-gas turbulence. Eqs. (16) and (17) have recently been modified by adding ad-hoc compression-related terms that account for the density variation in the cylinder during the compression and expansion strokes. If such terms are introduced into the K-k model equations, these also have to be integrated through the

The main advantage of the K-k model is that of its simplicity. However, since the energy cascade is governed by the dissipation constant c, such a constant has to be set by the user, taking the combustion chamber geometry into account. In (Baratta et al., 2006), a sensitivity

a. the 'prompt burning' approach: the mixture that is entrained through the flame front burns promptly. Thus, it follows that the increment in xb for each computational step of the combustion phase is:

$$d\mathbf{x}\_b = \frac{1}{m\_{cyl}} \rho\_u A\_{b\not\!\!/} S\_b \frac{d\theta}{\alpha} \tag{12}$$

Examples of models based on such an approach can be found in the literature (Baratta et al., 2008; Bozza et al., 2005; D'Errico et al., 2002; Wu et al., 1993; Yoshiyama et al., 2001).

b. the 'entrainment + burning' approach: the entrainment process is followed by a laminar burnout process in the region behind the flame front. The process is described by the following equations:

$$d\mathbf{x}\_e = \frac{1}{m\_{\rm cyl}} \rho\_u A\_{bf} S\_b \frac{d\theta}{\alpha} \tag{13}$$

$$d\mathbf{x}\_b = \frac{\mathbf{x}\_e - \mathbf{x}\_b}{\tau} \tag{14}$$

xe being the entrained mass fraction (Brown et al., 1996; Grill et al., 2006; Hattrel et al., 2006; Poulos & Heywood, 1983; Wahiduzzaman et al., 1993). The rate of laminar burnout in eq. (14) is assumed to be proportional to the unburned mass behind the flame front. Since the burnout is postulated to take place at the laminar flame speed and over a typical length scale of the turbulence microscale, the time constant in eq. (14) is given by the ratio of the Taylor turbulence microscale, , and of the laminar flame speed SL (Wahiduzzaman et al., 1993):

$$
\pi = \frac{\lambda}{S\_L} \tag{15}
$$

The calculation of the turbulent flame speed Sb, to be used in Eqs. (12) or (13), involves the evaluation of the laminar burning speed, SL, of the in-cylinder turbulence generation, as well as of the turbulence-flame interaction (Baratta et al., 2006). As will be shown later on, each of these aspects can exert a remarkable influence on the indicated cycle prediction. For this reason, the assessment of the accuracy of each of the adopted sub-models is often very difficult. A comparison between simulated and experimental or diagnosed combustionrelated quantities generally allows one to assess the overall simulation-model accuracy, which depends on how the different sub-model results are combined. A critical discussion of the main sub-models which are required is provided hereafter.

#### **3.2 In-cylinder turbulence sub-model**

In-cylinder turbulence modelling is required in a predictive combustion model, in order to quantify the increase in burning velocity, due to the turbulence, with respect to the case of laminar combustion. In the framework of thermo-dynamic models, a zero-dimensional model is usually applied, which provides a uniform value of the mean-flow velocity and of the root mean square (rms) turbulent fluctuation. The turbulence model results can also be

a. the 'prompt burning' approach: the mixture that is entrained through the flame front burns promptly. Thus, it follows that the increment in xb for each computational step of

b. the 'entrainment + burning' approach: the entrainment process is followed by a laminar burnout process in the region behind the flame front. The process is described by the

*e b*

xe being the entrained mass fraction (Brown et al., 1996; Grill et al., 2006; Hattrel et al., 2006; Poulos & Heywood, 1983; Wahiduzzaman et al., 1993). The rate of laminar burnout in eq. (14) is assumed to be proportional to the unburned mass behind the flame front. Since the burnout is postulated to take place at the laminar flame speed and over a typical length scale of the turbulence microscale, the time constant in eq. (14) is given by the ratio of the Taylor turbulence microscale, , and of the laminar flame speed

> *SL*

The calculation of the turbulent flame speed Sb, to be used in Eqs. (12) or (13), involves the evaluation of the laminar burning speed, SL, of the in-cylinder turbulence generation, as well as of the turbulence-flame interaction (Baratta et al., 2006). As will be shown later on, each of these aspects can exert a remarkable influence on the indicated cycle prediction. For this reason, the assessment of the accuracy of each of the adopted sub-models is often very difficult. A comparison between simulated and experimental or diagnosed combustionrelated quantities generally allows one to assess the overall simulation-model accuracy, which depends on how the different sub-model results are combined. A critical discussion

In-cylinder turbulence modelling is required in a predictive combustion model, in order to quantify the increase in burning velocity, due to the turbulence, with respect to the case of laminar combustion. In the framework of thermo-dynamic models, a zero-dimensional model is usually applied, which provides a uniform value of the mean-flow velocity and of the root mean square (rms) turbulent fluctuation. The turbulence model results can also be

of the main sub-models which are required is provided hereafter.

1 *e ubf b cyl <sup>d</sup> dx A S m*

> *b x x dx*

Examples of models based on such an approach can be found in the literature (Baratta et al., 2008; Bozza et al., 2005; D'Errico et al., 2002; Wu et al., 1993; Yoshiyama et al.,

(12)

(13)

(14)

(15)

1 *b ubf b cyl <sup>d</sup> dx A S m*

the combustion phase is:

following equations:

SL (Wahiduzzaman et al., 1993):

**3.2 In-cylinder turbulence sub-model** 

2001).

used by the in-cylinder heat-transfer sub-model, in order to derive a value for the convective heat-transfer coefficient.

A rather simple, but quite widely applied model is the K-k one (Poulos & Heywood, 1983), which is based on a zero-dimensional energy cascade from the mean flow to the viscous eddy dissipation. According to such a model, the rates of change in the mean flow kinetic energy ( <sup>2</sup> 1 2 *K mU cyl* ) and in the turbulent kinetic energy ( <sup>2</sup> 32 ' *cyl k mu* ) are, respectively:

$$\frac{dK}{dt} = \frac{1}{2}\dot{m}\_i v\_i^2 - P - K\frac{\dot{m}\_o}{m\_{cyl}}\tag{16}$$

$$\frac{dk}{dt} = P - m\_{cyl}\varepsilon - k\frac{\dot{m}\_o}{m\_{cyl}}\tag{17}$$

where: U is the average mean-flow velocity, u' is the rms turbulent velocity, '3 *i u L* is the

rate of turbulent kinetic energy dissipation per unit mass, Li is the characteristic size of the large-scale eddies, m is the mass in the cylinder, *mi* is the mass flow rate into the cylinder, *mo* is the mass flow rate out of the cylinder, *<sup>i</sup> v* is the jet velocity entering the chamber and P is the rate of turbulent kinetic energy production. Since the turbulence model does not spatially resolve the flow parameters, the production term P is empirically estimated from mean flow quantities:

$$P = 0.3307c\_{\beta} \frac{K}{L\_i} \left(\frac{k}{m}\right)^{\frac{1}{2}} \tag{18}$$

where the turbulent dissipation constant c is adjusted to give the expected profiles of u' and U throughout the whole cycle. It is worth pointing out that the accuracy of a turbulent burning-speed model depends to a great extent on the unburned-gas turbulence evolution during combustion, which is normally thought to present a maximum near firing TDC (Baratta et al., 2008, Bozza et al, 2005, Heywood, 1988; Poulos & Heywood, 1983; Yoshiyama et al., 2001). Although the magnitude of such a maximum depends, among other factors, on the configuration of the combustion chamber and, to a lesser extent, on the intake system, its presence is generally modelled by means of heuristic approaches that add turbulence during compression and combustion. As an example, it can be seen that, in (Baratta et al., 2006, 2008; Poulos & Heywood, 1983), as soon as the combustion starts, Eq. (17) is no longer integrated, and the evolution of the turbulence intensity and of the spatial macro scale is calculated by assuming conservation of angular momentum for large scale eddies, which gives rise to an increase in the unburned-gas turbulence. Eqs. (16) and (17) have recently been modified by adding ad-hoc compression-related terms that account for the density variation in the cylinder during the compression and expansion strokes. If such terms are introduced into the K-k model equations, these also have to be integrated through the combustion phase (Bozza et al., 2005).

The main advantage of the K-k model is that of its simplicity. However, since the energy cascade is governed by the dissipation constant c, such a constant has to be set by the user, taking the combustion chamber geometry into account. In (Baratta et al., 2006), a sensitivity

Numerical Simulation Techniques for the

affect the combustion velocity.

(Wahiduzzaman et al., 1993):

**3.3 Turbulent flame-speed model** 

Prediction of Fluid-Dynamics, Combustion and Performance in IC Engines Fuelled by CNG 279

few regions, rather than into a multitude of computational nodes. Thus, the turbulent kinetic energy actually includes a wide range of scales, which are partially superimposed onto the mean-flow spectrum. A turbulence production term, due to the unburned-gas compression, is also included (Morel et al., 1988). As far as the prediction of the convective heat-transfer coefficient is concerned, Morel & Keribar's model is expected to require fewer adjustments when it is applied to a new engine, due to its expanded physical basis compared to the previous models. However, when the accurate prediction of in-cylinder turbulence is required, as is the case of a predictive combustion model implementation, the proposed discretization of the k- model equations does not seem to be much more accurate than the very simple production term formulation given by Eq. (18). In particular, since a portion of the kinetic energy associated with the mean-flow scales is included in the turbulence terms, the simulated u' velocity may be different from the turbulent velocity which is presumed to

Turbulent burning-velocity models start from laminar burning speed (SL) data for the burning mixture, and are based on the evaluation of both the turbulence level and the turbulence-flame interaction. The latter can be modelled through either a non-fractal or a fractal approach, taking the turbulence level into account. In the first approach, the entrainment velocity Sb is the sum of a convective component and a diffusive one

<sup>1</sup> ' 1

rb being the burned-gas radius, Li the turbulence macroscale and Cs and Ck two fitting parameters. The term in brackets modulates the contribution of u', in order to avoid an overprediction of the flame speed at the initial stage of combustion, when the size of the

The second approach is based on the fractal and laminar flamelet concepts (Baratta et al., 2006; Gouldin & Miles, 1995), in order to represent the turbulent premixed flames. The laminar flamelet concept assumes that the combustion within a turbulent flame is confined to asymptotically thin moving laminar flamelets that are embedded in the turbulent flow. Since these thin layers behave like laminar flames, the turbulent burning velocity can be evaluated as the product of the surface area of the flamelets and the laminar burning speed, which is corrected to take stretch and flame curvature effects into account. The characterization and evaluation of the flamelet surface area, in relation to the turbulence properties, have been the subjects of many studies (Gülder and Smallwood, 1995; Gülder et al., 2000; North & Santavicca, 1990). According to the fractal theory, the general expression of the correlation that links the turbulent flame surface area and speed to the laminar ones is:

> *T b o LL i*

where AT and AL are the wrinkled and mean flame front surface areas, Sb and SL represent the turbulent burning speed and the laminar flame velocity, o and i are the so-called outer and inner turbulence cutoff length scales, and D is the fractal dimension of the flame front

*A S A S*

1 *bLs*

*S S Cu*

flame kernel is comparable with that of the turbulence eddies.

2 2

(19)

*kb i*

*D* 2

(20)

*Cr L* 

analysis of the computed mean velocity and turbulence intensity was carried out in a range of c values of between 0.75 and 3.0. The experimental turbulence data close to TDC, taken in a similar motored engine to that under study, were in good agreement with the calculated levels when the value c=3.0 was employed.

Fig. 8. Distributions of diagnosed and computed Sb/SL, versus (column on the left) and rb (column on the right), for different engine speeds – 2000cc 16v NA engine, CNG fuelling, bmep = 440 kPa, RAFR = 1.0, MBT timing (Baratta et al., 2006).

(Morel & Keribar, 1985) tried to develop a rather innovative flow model, with an expanded physical basis, in order to reduce the necessity of adjustments to account for engine-toengine differences. The model was initially developed with reference to convective heattransfer estimation in Diesel engines, and was then extended to SI engines and to the prediction of the turbulence level, for Sb prediction purposes (Morel et al., 1988). The cylinder volume is divided into multiple regions: the central core region, the squish region, the head recess region, and the piston cup region. Some axial and radial velocities are calculated, at each time step in each region, from mass conservation and piston kinematics, whereas the other radial and axial velocities, which are essentially three-dimensional, are included in the turbulence. A swirl equation, based on the conservation of angular momentum under the hypothesis of solid-body rotation, is solved for each region. The turbulent kinetic energy and the dissipation rate evolution are predicted by means of the k model equations (Wilcox, 1994), which were rewritten to account for a discretization into a few regions, rather than into a multitude of computational nodes. Thus, the turbulent kinetic energy actually includes a wide range of scales, which are partially superimposed onto the mean-flow spectrum. A turbulence production term, due to the unburned-gas compression, is also included (Morel et al., 1988). As far as the prediction of the convective heat-transfer coefficient is concerned, Morel & Keribar's model is expected to require fewer adjustments when it is applied to a new engine, due to its expanded physical basis compared to the previous models. However, when the accurate prediction of in-cylinder turbulence is required, as is the case of a predictive combustion model implementation, the proposed discretization of the k- model equations does not seem to be much more accurate than the very simple production term formulation given by Eq. (18). In particular, since a portion of the kinetic energy associated with the mean-flow scales is included in the turbulence terms, the simulated u' velocity may be different from the turbulent velocity which is presumed to affect the combustion velocity.

### **3.3 Turbulent flame-speed model**

278 Computational Simulations and Applications

analysis of the computed mean velocity and turbulence intensity was carried out in a range of c values of between 0.75 and 3.0. The experimental turbulence data close to TDC, taken in a similar motored engine to that under study, were in good agreement with the calculated

**<sup>0</sup>**

**<sup>0</sup>**

Fig. 8. Distributions of diagnosed and computed Sb/SL, versus (column on the left) and rb (column on the right), for different engine speeds – 2000cc 16v NA engine, CNG fuelling,

(Morel & Keribar, 1985) tried to develop a rather innovative flow model, with an expanded physical basis, in order to reduce the necessity of adjustments to account for engine-toengine differences. The model was initially developed with reference to convective heattransfer estimation in Diesel engines, and was then extended to SI engines and to the prediction of the turbulence level, for Sb prediction purposes (Morel et al., 1988). The cylinder volume is divided into multiple regions: the central core region, the squish region, the head recess region, and the piston cup region. Some axial and radial velocities are calculated, at each time step in each region, from mass conservation and piston kinematics, whereas the other radial and axial velocities, which are essentially three-dimensional, are included in the turbulence. A swirl equation, based on the conservation of angular momentum under the hypothesis of solid-body rotation, is solved for each region. The turbulent kinetic energy and the dissipation rate evolution are predicted by means of the k model equations (Wilcox, 1994), which were rewritten to account for a discretization into a

**Sb / SL [-]**

**Sb / SL [-]**

**3**

**3**

**0**

**3**

**Sb / SL [-]**

**6**

**9**

**12**

**6**

**9**

**12**

**6**

**9**

**CNG - N = 2000rpm exp. eq. (20)**

**CNG - N = 3300rpm exp. eq. (20)**

**0 10 20 30 40 50 rb [mm]**

**CNG - N = 4600rpm exp. eq. (20)**

**eq. (22)**

**eq. (22)**

**eq. (22)**

**12**

**CNG - N = 2000rpm exp. eq. (20)**

**CNG - N = 3300rpm exp. eq. (20)**

**CNG - N = 4600rpm exp. eq. (20)**

**330 345 360 375 390 405 420 [deg]**

bmep = 440 kPa, RAFR = 1.0, MBT timing (Baratta et al., 2006).

**eq. (22)**

**eq. (22)**

**eq. (22)**

levels when the value c=3.0 was employed.

**0**

**12**

**0**

**12**

**0**

**3**

**Sb / SL [-]**

**6**

**9**

**3**

**Sb / SL [-]**

**6**

**9**

**3**

**Sb / SL [-]**

**6**

**9**

**12**

Turbulent burning-velocity models start from laminar burning speed (SL) data for the burning mixture, and are based on the evaluation of both the turbulence level and the turbulence-flame interaction. The latter can be modelled through either a non-fractal or a fractal approach, taking the turbulence level into account. In the first approach, the entrainment velocity Sb is the sum of a convective component and a diffusive one (Wahiduzzaman et al., 1993):

$$S\_b = S\_L + C\_s \mu' \left(1 - \frac{1}{1 + C\_k \, r\_b^2 \sqrt{L\_i^2}}\right) \tag{19}$$

rb being the burned-gas radius, Li the turbulence macroscale and Cs and Ck two fitting parameters. The term in brackets modulates the contribution of u', in order to avoid an overprediction of the flame speed at the initial stage of combustion, when the size of the flame kernel is comparable with that of the turbulence eddies.

The second approach is based on the fractal and laminar flamelet concepts (Baratta et al., 2006; Gouldin & Miles, 1995), in order to represent the turbulent premixed flames. The laminar flamelet concept assumes that the combustion within a turbulent flame is confined to asymptotically thin moving laminar flamelets that are embedded in the turbulent flow. Since these thin layers behave like laminar flames, the turbulent burning velocity can be evaluated as the product of the surface area of the flamelets and the laminar burning speed, which is corrected to take stretch and flame curvature effects into account. The characterization and evaluation of the flamelet surface area, in relation to the turbulence properties, have been the subjects of many studies (Gülder and Smallwood, 1995; Gülder et al., 2000; North & Santavicca, 1990). According to the fractal theory, the general expression of the correlation that links the turbulent flame surface area and speed to the laminar ones is:

$$\begin{aligned} \text{A}\_{\text{T}} \bigvee\_{\text{A}\_{\text{L}}} &= \bigvee\_{\text{S}\_{\text{L}}} \bigvee\_{\text{S}\_{\text{L}}} = \left( \bigvee\_{\text{E}\_{\text{i}}} \bigvee\_{\text{E}\_{\text{i}}} \right)^{\text{D}-2} \end{aligned} \tag{20}$$

where AT and AL are the wrinkled and mean flame front surface areas, Sb and SL represent the turbulent burning speed and the laminar flame velocity, o and i are the so-called outer and inner turbulence cutoff length scales, and D is the fractal dimension of the flame front

Numerical Simulation Techniques for the

2005; D'Errico et al., 2002; Wu et al., 1993).

**0**

RAFR = 1.0, MBT timing.

**0.02**

**H**

**R**

**Rmax [deg-1]**

**0.04**

**0.06**

**exp. reference**

**pmax [bar]**

Prediction of Fluid-Dynamics, Combustion and Performance in IC Engines Fuelled by CNG 281

The results of the application of the Baratta et al. model are reported in Fig. 8, in terms of predicted Sb/SL profiles versus crank angle and versus burned-gas radius. In the figure, 'exp.' denotes the results of the multizone diagnostic model application (Catania et al., 2004). Two calculation results are plotted: one is obtained through Eq. (22) and the other with the original Eq. (20) in which Li is used as the outer cutoff scale. The general good agreement between the predicted and experimental results in Fig. 8 clearly shows the capability of the model to describe the flame-turbulence interaction in the overall flame propagation interval, from the development to the extinction of the flame. Furthermore, as can be inferred from this figure, Eq. (22) can extrapolate a value of Sb/SL that is almost equal to 1 as rb approaches zero, in agreement with the flame propagation theory (Baratta et al., 2006; Velherst & Sheppard, 2009). This supports the correctness of the introduction of the

reduced species and heat transfer across the flame front at low in-cylinder densities. If Eq. (20) is applied with o=Li, the modulation of Sb/SL is insufficient, and the model generally needs special measures for a correct end-of-combustion simulation (Bozza et al.,

The model developed in (Baratta et al., 2008) resulted to be in quite good agreement with the experimental or diagnosed combustion-related parameters. The main attractiveness of this model is given by its 'intrinsic' capability to reproduce, with good accuracy, the Sb/SL modulation through combustion, as well as its straightforward calibration. In fact, the crank angle at which the 0 term is evaluated in Eq. (22) is the only parameter that should be set for a specific engine, and it is virtually independent of fuelling and operating conditions. An example of model validation is given in the column on the left of Fig. 9: with the exception of a couple of points, the agreement of the maximum in-cylinder pressure pmax and the maximum heat-release rate HRRmax is within a few percent. Obviously, the model accuracy could be improved if an ad-hoc model tuning is carried out for each operating point, but even with a fixed calibration set for a given engine and fuelling, the degree of accuracy is

**<sup>0</sup>**

**pmax [bar]**

**0**

**0.02**

**H**

Fig. 9. Comparison between experimental/diagnosed and simulated quantities (left) and model sensitivity analysis (right) – 1200cc 8v TC engine, Methane fuelling, N = 3000 rpm,

**R**

**Rmax [deg-1]**

**0.04**

**0.06**

**reference turbulence +40% turbulence -40% convective coeff. +30% convective coeff. -30% D +0.01 D -0.01**

**0 400 800 1200 1600 bmep [kPa]**

*n*

to account for the

 

term *Abf* into the definition of o, as well as of the term <sup>0</sup>

acceptable, and the model can be used for predictive purposes.

**0 400 800 1200 1600 bmep [kPa]**

surface. The success and the applicability of these models depend to a great extent on the availability of reliable values for the inner and outer cutoff length scales and for the fractal dimension. Correlations to evaluate the cutoff length scales have been proposed, namely in (Gülder and Smallwood, 1995).

The inner cutoff length scale i is taken equal to the turbulence micro scale :

$$
\omega\_i = \eta = \mathrm{L}\_i \left( \mu' \bigvee\_{\nu} \bigvee\_{\nu} \right)^{-\sum\_{\substack{\ell}}} \tag{21}
$$

The outer cutoff length scale has been taken equal to the integral turbulence scale, Li, in many studies (D'Errico et al., 2002; Wu et al., 1993). However, although this choice could be acceptable for stationary flames issuing from burners (Gülder et al., 2000), in SI engine combustion, in which the mean flame front is assumed to be spherical, the wrinkling effect of turbulence on the flame front should be a function of the ratio between the characteristic flame-front and eddy dimensions, because the initially regular flame front surface progressively becomes more and more wrinkled as its dimension increases with respect to turbulent eddies. Therefore, (Baratta et al., 2006, 2008) replaced the integral length scale in Eq. (20) with a characteristic linear dimension of the flame front, i.e., the square root of its surface area. Furthermore, (Gülder et al., 2000) pointed out that, if the fractal geometry approach yields a true measure of the wrinkled surface area of the flame front, then Eq. (20) may not be a reasonable assumption for the turbulent premixed flames in the flamelet regime. In fact, small-scale turbulence can enhance the transfer of heat and species across the flame front. For this reason, (Baratta et al., 2006, 2008) introduced a pre-multiplying factor in Eq. (20), and argued that such a factor should be proportional to the average charge density, as an increased density means higher concentrations of reactive species near the flame front. Thus, the proposed formula for the turbulent burning speed was:

$$\frac{\mathcal{S}\_b}{\mathcal{S}\_L} = \left(\frac{\rho}{\rho\_0}\right)^n \left(\frac{\mathbb{C}\_L \sqrt{A\_{bf}}}{\eta}\right)^{D-2} = \left(\frac{\rho}{\rho\_0}\right)^n \left(\frac{\mathbb{C}\_L \sqrt{A\_{bf}}}{\mathbb{C}\_L \left(\mathbb{I}\_{\min} + \mathbb{S}\_p\right) \left(\frac{u' \mathbb{C}\_L \left(\mathbb{I}\_{\min} + \mathbb{S}\_p\right)}{\nu}\right)^{-\frac{\mathcal{N}}{\mathcal{N}}}}\right)^{D-2} \tag{22}$$

where hmin is the clearance height, Sp the instantaneous piston displacement from TDC, Abf the unwrinkled flame-front area, and CL a closure coefficient that depends on the engine fuelling and speed. n was found to be 1.25 for two different engines and 0 is the average density evaluated at a reference crank angle, which depends on the specific engine (Baratta et al., 2008).

Although predictive combustion models have been developed as phenomenological models, they can simply be considered as mathematical descriptions that agree with the S-shaped mass-fraction burned observations. It is reasonable to use either Eq. (19) or the fractal approach (Eqs. (20) or (22)), coupled to the 'prompt burning' or the 'entrainment + burning' approaches, provided appropriate choices are made for Abf, and the specific model's constants. However, though it is not a general rule, virtually all the fractal models used in the literature are coupled to the 'prompt burning' approach, whereas the non fractal ones are implemented with the 'entrainment + burning' formulation.

surface. The success and the applicability of these models depend to a great extent on the availability of reliable values for the inner and outer cutoff length scales and for the fractal dimension. Correlations to evaluate the cutoff length scales have been proposed, namely in

*u L*

The outer cutoff length scale has been taken equal to the integral turbulence scale, Li, in many studies (D'Errico et al., 2002; Wu et al., 1993). However, although this choice could be acceptable for stationary flames issuing from burners (Gülder et al., 2000), in SI engine combustion, in which the mean flame front is assumed to be spherical, the wrinkling effect of turbulence on the flame front should be a function of the ratio between the characteristic flame-front and eddy dimensions, because the initially regular flame front surface progressively becomes more and more wrinkled as its dimension increases with respect to turbulent eddies. Therefore, (Baratta et al., 2006, 2008) replaced the integral length scale in Eq. (20) with a characteristic linear dimension of the flame front, i.e., the square root of its surface area. Furthermore, (Gülder et al., 2000) pointed out that, if the fractal geometry approach yields a true measure of the wrinkled surface area of the flame front, then Eq. (20) may not be a reasonable assumption for the turbulent premixed flames in the flamelet regime. In fact, small-scale turbulence can enhance the transfer of heat and species across the flame front. For this reason, (Baratta et al., 2006, 2008) introduced a pre-multiplying factor in Eq. (20), and argued that such a factor should be proportional to the average charge density, as an increased density means higher concentrations of reactive species near the flame front. Thus,

<sup>3</sup> 0 0 <sup>4</sup>

where hmin is the clearance height, Sp the instantaneous piston displacement from TDC, Abf the unwrinkled flame-front area, and CL a closure coefficient that depends on the engine fuelling and speed. n was found to be 1.25 for two different engines and 0 is the average density evaluated at a reference crank angle, which depends on the specific engine (Baratta

Although predictive combustion models have been developed as phenomenological models, they can simply be considered as mathematical descriptions that agree with the S-shaped mass-fraction burned observations. It is reasonable to use either Eq. (19) or the fractal approach (Eqs. (20) or (22)), coupled to the 'prompt burning' or the 'entrainment + burning' approaches, provided appropriate choices are made for Abf, and the specific model's constants. However, though it is not a general rule, virtually all the fractal models used in the literature are coupled to the 'prompt burning' approach, whereas the non fractal ones

min

*Ch S*

*L p*

3 4 '*i*

'

 

min

  2

(22)

*D*

(21)

The inner cutoff length scale i is taken equal to the turbulence micro scale :

 *L*

the proposed formula for the turbulent burning speed was:

2

are implemented with the 'entrainment + burning' formulation.

*L bf L bf b*

 

 

*S C A C A*

*<sup>L</sup> L p*

*<sup>S</sup> uC h S*

*<sup>D</sup> n n*

et al., 2008).

*i i*

(Gülder and Smallwood, 1995).

The results of the application of the Baratta et al. model are reported in Fig. 8, in terms of predicted Sb/SL profiles versus crank angle and versus burned-gas radius. In the figure, 'exp.' denotes the results of the multizone diagnostic model application (Catania et al., 2004). Two calculation results are plotted: one is obtained through Eq. (22) and the other with the original Eq. (20) in which Li is used as the outer cutoff scale. The general good agreement between the predicted and experimental results in Fig. 8 clearly shows the capability of the model to describe the flame-turbulence interaction in the overall flame propagation interval, from the development to the extinction of the flame. Furthermore, as can be inferred from this figure, Eq. (22) can extrapolate a value of Sb/SL that is almost equal to 1 as rb approaches zero, in agreement with the flame propagation theory (Baratta et al., 2006; Velherst & Sheppard, 2009). This supports the correctness of the introduction of the

term *Abf* into the definition of o, as well as of the term <sup>0</sup> *n* to account for the reduced species and heat transfer across the flame front at low in-cylinder densities. If Eq. (20) is applied with o=Li, the modulation of Sb/SL is insufficient, and the model generally needs special measures for a correct end-of-combustion simulation (Bozza et al., 2005; D'Errico et al., 2002; Wu et al., 1993).

The model developed in (Baratta et al., 2008) resulted to be in quite good agreement with the experimental or diagnosed combustion-related parameters. The main attractiveness of this model is given by its 'intrinsic' capability to reproduce, with good accuracy, the Sb/SL modulation through combustion, as well as its straightforward calibration. In fact, the crank angle at which the 0 term is evaluated in Eq. (22) is the only parameter that should be set for a specific engine, and it is virtually independent of fuelling and operating conditions. An example of model validation is given in the column on the left of Fig. 9: with the exception of a couple of points, the agreement of the maximum in-cylinder pressure pmax and the maximum heat-release rate HRRmax is within a few percent. Obviously, the model accuracy could be improved if an ad-hoc model tuning is carried out for each operating point, but even with a fixed calibration set for a given engine and fuelling, the degree of accuracy is acceptable, and the model can be used for predictive purposes.

Fig. 9. Comparison between experimental/diagnosed and simulated quantities (left) and model sensitivity analysis (right) – 1200cc 8v TC engine, Methane fuelling, N = 3000 rpm, RAFR = 1.0, MBT timing.

Numerical Simulation Techniques for the

Vermont, USA, 2005.

April 3-6, 2006.

062805-1/11, 2008.

(3), pp. 289−306, 2010.

Engines, Vol.109, pp.1610-1620, 2001.

Series B, Vol. 46, No.1, pp. 75-85.

1999.

1974.

1996.

ON, Canada, September 11-14, 2005.

Prediction of Fluid-Dynamics, Combustion and Performance in IC Engines Fuelled by CNG 283

Badami, M., Millo, F., and Giaffreda, G., "Experimental and Computational Analysis of a

Baines, N.C., "Fundamentals of Turbocharging", Concepts NREC, White River Junction,

Baratta, M., d'Ambrosio, S., Spessa, E., and Vassallo, A., "Analysis of Cyclic Variability in a

Baratta, M., Catania, A.E., Spessa, E., and Vassallo, A., "Development and Assessment of a

Baratta, M., Catania, A.E., d'Ambrosio, S., and Spessa, E., "Prediction of Combustion

Baratta, M., and Spessa, E., "Turbocharged CNG Engines for Urban Transportation:

Baratta, M., Spessa, E., and Mairone, P., "Numerical Investigation of Turbolag Reduction In

Blair, G.P., "Design and Simulation of Four Stroke Engines", SAE R-186 ISBN 0-7680-0440-3,

Blizzard, N.C., and Keck, J.C., "Experimental and Theoretical Investigation of Turbulent

Bozza, F., Gimelli, A., Merola, S. S., and Vaglieco, B. M., "Validation of a Fractal Combustion Model Through Flame Imaging," SAE Paper No. 2005- 01-1120, 2005. Brown, A.G., Stone, C.R., and Beckwith, P. "Cycle-by-Cycle Variations in Spark Ignition

Catania, A.E., Misul, D., Mittica, A., and Spessa, E., "Unsteady Convection Model for Heat

Catania, A.E., Misul, D., Mittica, A., and Spessa, E., "A Refined Two-Zone Heat Release

Catania, A. E., Misul, D., Spessa, E., and Vassallo, A., 2004, "A New Quasi-Dimensional

the ASME ICED, Milwaukee, Wisconsin, USA, May 3-6, 2009.

Geometry Exhaust System", SAE Paper 2002-01-001, 2002.

High Performance Four-Stroke Motorcycle Engine Equipped with a Variable

Bi-Fuel Engine By Means of a 'Cycle-Resolved' Diagnostic Technique", ASME Paper ICEF2005-1214, 2005 Fall Technical Conference of the ASME ICED, Ottawa,

Multizone Combustion Simulation Code for SI Engines Based on a Novel Fractal Model", SAE Paper No. 2006-01-0048, SAE 2006 World Congress, Detroit, MI, USA,

Parameters, Performance and Emissions in CNG and Gasoline SI Engines", ASME Transactions, Journal of Engineering for Gas Turbines and Power, Vol. 130, pp.

Evaluation of Turbolag Reduction Strategies by means of Computational Analyses", ASME Paper No. ICES2009-76067, 2009 Spring Technical Conference of

HD CNG Engines By Means Of Exhaust Valve Variable Actuation and Spark Timing Control", KSAE International Journal of Automotive Technology, Vol. 11

Burning Model for Internal Combustion Engines," SAE Paper No. 740191,

Engine Combustion – Part I: Flame Speed and Combustion Measurement and a Simplified Turbulent Combustion Model", SAE Paper No. 960612,

Release Analysis of IC Engine Pressure Data," SAE Transactions, Journal of

Model for Combustion Analysis in SI Engines," JSME International Journal, 2003,

Multizone Combustion Diagnostic Model for the Analysis of Heat Release, Flame

The column on the right in Fig. 9 provides an example of the influence of the main submodels on the results of the overall model. In-cylinder turbulence and heat-transfer submodels are considered, and in addition, the influence of a slight change in the fractal dimension, D, is assessed. The 'reference' series is the same as in the first column, and represents the calibration in (Baratta et al., 2008). The deviation from the 'reference' calibration for each sub-model has been set on the basis of the uncertainty that can be expected in the adopted modelling framework. In particular, a deviation of 40% was considered for the turbulence level at the spark discharge, due to the significant approximations in the K-k model. An error of 30% is reasonable for the heat-transfer, especially if a diagnostic tool is not available for its calibration. Finally, the uncertainty on the fractal dimension D can be even higher than 0.01, since at present there is no agreement on its value (Baratta et al., 2006). As can be seen, for the considered deviation values, an increase in the heat-transfer coefficient has almost the same effect as a decrease in the fractal dimension. Both parameters can influence the model performance to a certain extent. As can be expected, the turbulence level exerts a remarkable influence on the overall model output. In particular, the bell-like shape of the u' profile versus crank angle, although obtained through empirical formulas, is very important to obtain an acceptable Sb/SL profile.

The above discussion confirms that the overall model accuracy depends on each specific sub-model formulation, as well as on the related calibration. A precise model prediction can be obtained by adopting very accurate sub-models, but also when the sub-models error cancel each other. A good predictive combustion model should be formulated and calibrated so as to be able to reproduce the engine indicated cycle with a reasonable accuracy over a wide range of operating conditions, and to capture the engine performance trends when a design or operation variable is modified.
