Comparing in Cylinder Pressure Modelling of a DI Diesel Engine Fuelled on Alternative Fuel Using Two Tabulated Chemistry Approaches.

The present work presents the comparative simulation of a diesel engine fuelled on diesel fuel and biodiesel fuel. Two models, based on tabulated chemistry, were implemented for the simulation purpose and results were compared with experimental data obtained from a single cylinder diesel engine. The first model is a single zone model based on the Krieger and Bormann combustion model while the second model is a two-zone model based on Olikara and Bormann combustion model. It was shown that both models can predict well the engine's in-cylinder pressure as well as its overall performances. The second model showed a better accuracy than the first, while the first model was easier to implement and faster to compute. It was found that the first method was better suited for real time engine control and monitoring while the second one was better suited for engine design and emission prediction.

1. Introduction

The modeling of internal combustion engines has been largely developed during past years. For this purpose, a multitude of industrial codes dedicated to the simulation of the engines is available on market (GT-Power, Diesel-RK, Ricardo-Wave, Fluent, etc.). These codes, despite the fact that they can be useful for predicting engine performances, are expensive for third world's university laboratory and their source codes are practically impossible to modify to implement new model or functions.

When modeling compression engine different approaches can be used with different level of complexity, such as thermodynamic 0D models, quasi-dimensional multizone models, and computational fluid dynamics (CFD) models [1, 2].

For this study a quasi-dimensional approach was selected; these models allows to compute efficient, economic, and fast calculations of engine performances as a function of different engine parameters. Using these types of models we can compute the different stages of a diesel engine cycle as in compression, injection, ignition delay, and combustion and exhaust stage.

The evaluated models in this study will be evaluated in terms of accuracy and speed of calculation.

2. Governing Equations

2.1. Fuel Injection and Vaporization Model

For both models the fuel spray characteristic is modeled using the phenomenological model of Razleytsev [3] and Lyshevsky [4]. That model has been used and implemented in the so called RK-model [5] and its results. The main aim of this part of the model is to determine the finesse of the atomized fuel from the injector nozzle. The finesse of atomization is characterized by the calculated Sauter mean diameter of fuel droplet. The model is a simple sequence of calculated parameters and is structured as follows.

First the average outflow velocity of fuel from the injector is calculated in m/s by where q is the cyclic fuel supply in Kg/cycle; RPM is the crankshaft rotational speed in rpm; ρ is the density of the fuel in kg/m3; d is the diameter of the injector hole in mm; φ inj is the duration of injection in degree of rotation of the crankshaft.

The criterion M characterizing the relationship between the surface tension force, inertia, and viscosity is calculated as follows: where μ is the coefficient of dynamic viscosity of the fuel at a temperature of 323 K in Pa·s; σ is the coefficient of surface tension of the fuel at a temperature of 323 K in N/m.

Weber number characterizing the relationship between surface tension force and inertia is determined as

The density of the charge at the end of the compression before the TDC is calculated as follows: where μ air = 28.9 Kg/Kmole is the molecular mass of air; ρ air is the density of air.

We then determine the fuel/air density ratio by

Finally the Sauter mean diameter of the atomized fuel is calculated in micron by where E 32 is an empirical factor depending on designs of the injector, whose recommended value is 1.7.

At this point the finesse of the pulverized fuel is determined and we can see it mainly depends on physical properties of the fuel such as viscosity and density. The next step of the model is to determine the kinetic of combustion of the atomized fuel, whose kinetic will be dependent of the earlier found Sauter mean diameter.

2.2. Model of the Kinetic of Combustion of the Atomized Fuel

This part of the elaborated model will permit us to determine parameters of the process such as the evaporation rate of the atomized fuel, its ignition delay, and the duration of combustion.

Pressure in the cylinder at the end of the compression before the TDC is where P 0 is the reference pressure which in our case corresponds to the inlet manifold pressure.

We then determine the theoretical constant of evaporation of the fuel by where K = 1/P · 106 is the evaporation constant related to cylinder pressure in m2/s where A is a constant characterizing the duration of evaporation of large drops in diesel engine in 1/s and its value is given as  2.4 in [3] but Kuleshov [5] used it as weight coefficient that can be varied in order to match experimental data; λ is the coefficient of excess air.

The full duration of combustion is then calculated as with ID being the ignition delay computed using the formula of Hardenberg and Hase [6] in both models

The value of the apparent activation energy E in this correlation is given by E = 618840/(CN + 25); however, this value was computed for the particular experiment reported in [6] and needs to be corrected in order to be used for different experimental conditions. In order to match experimental ignition delay values, especially for biodiesel fuel, the method proposed by Aghav et al. [7], where the value 618840 can be varied from the initial 618840 to higher values, was used.

2.3. Heat Release and Heat Transfer Laws

For the two models, we used a double Wiebe function [8] to model the heat release rate and the Woschni [9] law for heat transfer simulation.

3. Combustion Models

3.1. Method 1: Combustion Model Based on Borman and Krieger [10]

This method is a single zone model based on the first law of thermodynamics, mass conservation, and ideal gas laws. Following that we can write the governing equation for calculating variation of in-cylinder pressure with respect to crank angle as [11]

The Krieger and Bormann algorithm is based on polynomial fitting constant to compute the adiabatic index (ratio of specific heat) as function of crank angle. The fitting coefficients can be found in the Appendix.

The internal energy of combustion products at equilibrium for a reaction between air and a hydrocarbon CnH2n is given as where u corr = c exp⁡(D + E + F ) is the correction factor for the internal energy accounting for dissociation

The gas constant is computed as with R corr = c exp⁡(r 0log⁡λ + ((r 1 + r 2/T) + r 3log⁡(f 6 P))/λ).

The adiabatic index (ratio of specific heat) is computed as with c = ∂u/∂T.

3.2. Method 2: Combustion Model Based on Olikara and Bormann Method [12, 13]

This model is a two-zone thermodynamic model where the thermodynamics properties of combustion products and fuel are determined using polynomial curve fitted data from the Chemkin [14] or JANAF [15] tables and are given by

The in-cylinder pressure and temperature of burned and unburned gases derivative can be calculated following Ferguson's algorithm [13, 16] by where

The equations of state of the mixture are given by with (∂h/∂T) = c .

The combustion products are evaluated assuming that they are in equilibrium at a given temperature and pressure. Using Olikara and Borman [12] method we find combustion products mole fractions and we can then find the thermodynamic property of the mixture as in entropy, specific volume, enthalpy, and internal energy.

Considering the combustion reaction of a hydrocarbon we then have y is the mole fraction of the ith combustion specie at equilibrium and T the total number of mole. Following that we can write ∑ 10 y − 1 = 0.

The system of equation comprises 10 unknowns and 4 equations; we then need 6 complementary equations in order to be able to solve the system. For that we introduce six gas phase equilibrium reactions constant with p representing the pressure at which the reaction occurs in atmosphere. The reactions constants are computed using the following expression:

with the constants A, B, C, D, and E given in Table 1.

Table 1

Olikara and Borman constants.

	A	B	C	D	E
K p,1	0.432168	−0.112464 × 105	0.267269 × 101	−0.745744 × 10−4	0.242484 × 10−8
K p,2	0.310805	−0.129540 × 105	0.321779 × 101	−0.738336 × 10−4	0.344645 × 10−8
K p,3	−0.141784	−0.213308 × 104	0.853461	0.355015 × 10−4	−0.310227 × 10−8
K p,4	0.150879 × 10−1	−0.470959 × 104	0.646096	0.272805 × 10−5	−0.154444 × 10−8
K p,5	−0.752364	0.124210 × 105	−0.260286 × 101	0.259556 × 10−3	−0.162687 × 10−7
K p,6	−0.415302 × 10−2	0.148627 × 105	−0.475746 × 101	0.124699 × 10−3	−0.900227 × 10−8

We find ourselves with a set of 11 nonlinear equations with 11 unknowns which are solved using Newton Raphson iteration. A detailed description of the resolution algorithm can be found in [13]. The method was implemented using a modified Matlab script provided by [16] to insert new subroutines of the fuel injection and pulverization parameters as well as the biodiesel constant and the double Wiebe function of heat release rate.

4. Results and Discussion

For validation purposes the two models were implemented and compared with experimental results from [17], the engine characteristics are given in Table 2 where fuel injection timing refers to the actual crankshaft angle at which fuel starts coming out of the nozzle. Two sets of calculation were performed, one for diesel fuel and the second for biodiesel fuel. The characteristics of diesel and biodiesel fuel were taken from [17, 18] for method 1. To simulate biodiesel combustion for method 2 we used methyl butanoate thermodynamic data as surrogate [14].

Table 2

Engine specifications.

Number	Particular	Specifications
1	Make	Kirloskar oil engine
2	Model	DAF 8
3	Rated brake power (kW)	6
4	Rated speed (rpm)	1500
5	Number of cylinder	1
6	Bore × stroke (mm)	95 × 110
7	Compression ratio	17.5 : 1
8	Fuel injection timing	23°

The experimental imep (indicative mean effective pressure) was not given in the experimental work reported by [17], it was therefore evaluated using the experimental pressure trace by integrating the area under the P-V diagram using the trapezoidal rule, and the formula used was then with V being the displaced volume of the engine, in m3.

Simulations were computed at rated engine speed and the pressure plot was made and compared with the experimental results (Figures 1, 2, 3, and 4). Both models evaluated the combustion duration at 72 CAD and 87 CAD, respectively, for diesel and biodiesel fuel. This can be explained by the fact the viscosity and density of biodiesel are significantly higher than that of conventional diesel, which tends to generate pulverized droplets with higher mean diameter that tend to burn slower.

Figure 1

Comparison of computed and experimental pressure trace for model 1 for biodiesel fuel.

Figure 2

Comparison of computed and experimental pressure trace for model 2 for biodiesel fuel.

Figure 3

Comparison of computed and experimental pressure trace for model 1 for diesel fuel.

Figure 4

Comparison of computed and experimental pressure trace for model 2 for diesel fuel.

The pressure traces obtained from the two model show that they reproduce well the in cylinder pressure for both type of fuels in the given engine configuration. It can be seen that the maximum pressure decreases when biodiesel fuel is used for the two simulations, which is equally the trend observed experimentally and in other works such as [19-22]. This was expected since biodiesel fuel has higher oxygen content; the combustion coefficient tends to be higher in that case. Another trend observed in the pressure trace is that the cylinder pressure tends to be under predicted during the exhaust phase, especially for biodiesel fuel; this could be explained by the experimental conditions uncertainty.

The values of computed parameters are given in Table 3 in terms of maximum pressure, imep, maximum pressure occurrence, and ignition delay. In order to have an evaluation of the accuracy of both model, the mean relative error was calculated according to the four parameters given above, it was then observed that the first model gave a mean relative error of 18% while the second one had given a mean relative error of 11%.

Table 3

Comparison of simulated engine's combustion performances with experimental data.

	Diesel			Biodiesel
	Model 1	Model 2	Experiment	Model 1	Model 2	Experiment
Max pressure (MPa)	7.948	7.87	7.87	8.37	8.73	8.53
Imep (MPa)	0.52	0.47	0.6	0.522	0.6	0.7
Maximum pressure occurrence (CA after TDC)	4	3	5.8	3	3	5
Ignition delay (CAD)	8.69	8.69	8.5	4.37	4.37	4.3

5. Computation Speed

Another assessment made for the two models was the speed of computation; the first model was computed in 6 seconds, while the second took about 40 to 45 second to run completely on 2 Ghz dual core Pentium personal computer.

With this assessment, we can say that the first model could be suitable for real time engine control while the second could be used for engine's design and also emission evaluation since it computes in cylinder pressure by determining combustion species states at each time step. Another development that could be made for the second model can be the implementation of the Zeldovich mechanism [1, 11] to compute nitric oxide emissions.

6. Conclusion

The aim of the present study was to develop and implement two different models for diesel engine fuelled on biodiesel based on tabulated chemistry. The two models showed good predictability of engine performance, with the second approach giving a better accuracy. The first approach was faster to compute the cylinder pressure than the second and was found more suitable for engine monitoring and control while the second one was found better suited for engine design and emission prediction.

Table 4

Krieger and Borman coefficients.

		0	1	2	3	4	5	6
a			0.692	3.917e − 5	5.29e − 8	−2.29e − 11	2.7758e − 17
b		3049.33	−0.057	−9.5e − 5	2.153e − 8	−2e − 12
c u	2.32584
c r	0.004186
d		10.41066	7.85125		−3.71257
e		−15001	−15838		9613
f		−0.10329	−0.38656		0.154226	−14.763	118.27	14.503
r		−0.2977	11.98	−25442	−0.4354