**Thermodynamic Study of the Working Cycle of a Direct Injection Compression Ignition Engine**

Simón Fygueroa, Carlos Villamar and Olga Fygueroa

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/50028

## **1. Introduction**

74 Internal Combustion Engines

[10] C. Strozzi, J. Sotton, A. Mura, M. Bellenoue, Experimental and numerical study of the influence of temperature heterogeneities on self-ignition process of methane-air mixtures in a rapid compression machine, Combustion Science and Technology 180 (2008) 1829-1857. [11] J.F. Griffiths, J.P. MacNamara, C.G.W. Sheppard, D.A.Turton, B.J. Whitaker, The relationship of knock during controlled autoignition to temperature in homogeneities

[12] S.M. Walton, X. He, B.T. Zigler, M.S. Wooldridge, A. Atreya, An experimental investigation of iso-octane ignition phenomena, Combust. Flame 150 (3) (2007) 246-262. [13] G. Cho, G. Moon, D. Jeong, C. Bae, Effects of internal exhaust gas recirculation on controlled auto-ignition in a methane engine combustion, Fuel 88 (6) (2009) 1042-1048. [14] G. Mittall, C.J. Sung, A rapid compression machine for chemical kinetics studies at elevated pressures and temperatures, Combust. Sci. Technol. 179 (2007) 497-530. [15] C. Strozzi, Étude expérimentale de l'auto-inflammation de mélanges gazeux en milieux confines et sa modélisation avec une description cinétique chimique détaille, PhD

[16] Eliseu Monteiro, Combustion study of mixtures resulting from a gasification process of

[17] Lewis B., von Elbe G., Combustion, Flames and Explosions of Gases, 3rd Edition,

[18] Metghalchi, M., and Keck, J. C., Laminar Burning velocity of propane-air mixtures at

[20] Griffiths J.F, Q. Liao, A. Schreiber, J. Meyer, K.F Knoche, W. Kardylewski. Experimental and numerical studies of ditertiary butyl peroxide combustion at high pressures in a

[21] Clarkson J., Griffiths J.F., Macnamara J.P., Whitaker B.J. Temperature fields during the development of combustion in a rapid compression machine. Combustion and Flame

[22] Minetti R., Carlier M., Ribaucour M., E. Therssen, L.R. Sochet. Comparison of oxidation and autoignition of the two primary reference fuels by rapid compression. Proceedings

[23] Liou, T.M., Santavicca, D.A. Cycle resolved LDV measurements in a motored IC engine.

[24] Eliseu Monteiro, M. Bellenoue, J. Sotton, N.A. Moreira, S. Malheiro, Laminar burning velocities and Markstein numbers of syngas-air mixtures, Fuel 89 (2010) 1985-1991. [25] Eliseu Monteiro, J. Sotton, M. Bellenoue, N. A. Moreira, S. Malheiro. Experimental study of syngas combustion at engine-like conditions in a rapid compression machine.

[26] S. Verhelst, R. Sierens, A quasi-dimensional model for power cycle of a hydrogen-

[27] Alla A.G.H., Computer simulation of a four stroke spark ignition engine, Energy

[28] J. Hartman, How to Tune and Modify Engine Management Systems, Motorbooks, 2004.

Transitions of the ASME. Journal of Fluids Engineering 107 (1985) 232-240.

fuelled ICE, International Journal of Hydrogen Energy 32 (2007) 3545-3554.

Experimental Thermal and Fluid Sciences 35 (2011) 1473-1479.

Conversion and Management 43 (2002) 1043-1061.

high temperatures and pressure. Combustion and Flame 38 (1980) 143-154. [19] Ganesan V. Internal Combustion Engines. McGraw-Hill companies, 1995.

rapid compression machine. Combustion and Flame 93 (1993).303-315.

and fuel reactivity, Fuel 81 (7) (2002) 2219-2225.

Thesis, University of Poitiers, France, 2008.

Academic Press, 1987.

125 (2001) 1162-1175.

forest biomass. PhD Thesis, ENSMA, France, 2011.

of the Combustion Institute 26 (1996) 747-753.

Currently, one of the worldwide most used *energy sources* are fuels derived from oil, such as hydrocarbons, that burn with oxygen releasing large amount of thermal energy. This energy can be transformed into mechanical work by mean of internal combustion engines [1]. Internal combustion engine is a device that allows obtaining mechanical energy from the thermal energy stored in a fluid due to a combustion process [2].

It should be noted that in reciprocating internal combustion engines (RICE) the combustion products constitute the working fluid; this simplifies their design and produces high thermal efficiency. For this reason these engines are one of the lighter weight generating units known and thus actually are the most commonly used transport engines [3].

The RICE operate following a *mechanical cycle* consisting of two main parts: the first one is the closed cycle, where the compression, combustion and expansion processes are carried out, and the second one is the open cycle where the working fluid is renewed, known as the gas exchange process and constituted by the intake and exhaust processes [4]. When RICE is study, it is mandatory to determine the working fluid thermodynamic properties, as well as the mixture amount that enters and leaves the cylinder [5].

The *flow characteristics* in spark-ignition engines (SIE) or compression ignition engines (CIE) can be summarized according to [6] as follows: transient as a result of piston movement, fully turbulent for all cylinders due to engine velocities and admission duct dimensions, and three-dimensional due to the engine geometry that also varies during the cycle (contours varying with time) producing different local velocity fields.

During the *gas exchange* ondulatory and inertial phenomena processes occur, as well as instability in the processes that occur within the engine. The variation of in-cylinder pressure during the intake and the exhaust has a complex pattern, for this reason the

analytical calculation of gas exchange considering the above-mentioned phenomena is quite complicated and requires the use of specialized computer programs that use coefficients obtained experimentally [1].

The basis for calculation of non-stationary non-isentropic flow characteristic of RICE inlet and outlet ducts, and NO emissions was established in [7]. Various empirical correlations to take into account heat transfer during gas exchange process and to adjust the Reynolds number exponential factor in such a way that reduce to only one the adjustment coefficients were considered in[1,2,4 and 8].

A procedure widely used in both experimental and theoretical study of flow in engines, is to analyze the engine cycle in absence of combustion, simulating the process of compressionexpansion and making measurements to the engine operating at this condition [5].

In all RICE working cycle processes, there is heat transfer to the cylinder walls, which occurs with greater intensity during the combustion and expansion due to the high temperature gradients reached. Woschni [9] proposed equations to determine turbulent convective heat transfer considering average speed of in-cylinder gases and Annand [10] find correlations to calculate instant average coefficients for turbulent convection heat transfer using gas average temperature and proposed correlations to evaluate flame radiation emitted during combustion. Correlations for convection heat transfer taking into consideration surface change and cylinder enclosed volume as piston moves was established in [11]. A computer program to calculate the heat transfer in a RICE combustion chamber using models to consider turbulence was presented in [12]. A universal correlation for mixture flow in the admission and exhaust process, correcting the coefficients of Nusselt, Reynolds and Prandtl numbers was proposed in [13].

In present study the *compression* process is considered adiabatic and reversible, but in real engines there is heat transfer between working fluid, valves and cylinder walls. At the beginning of compression the fluid temperature is lower than the temperature of the surfaces that surround the cylinder volume, causing an increase in the fluid temperature, some instant later temperatures get equal and latter on, heat is transferred from the working fluid to the walls, therefore the politropic coefficient varies during the process [1].

The complexity of the combustion process in RICE because of untimely and incomplete combustion, dissociation and heat transfer, has encouraged the development of special techniques for carrying out studies. Adequate realization of this process is decisive in terms of the produced power and its efficiency having great influence on the engine life and reliability [14].

Various *models* have been proposed to study the combustion process such as the Wiebe burning law applicable to SIE and the Watson law, applicable to CIE [15]. These laws determine burned mass fraction and heat released depending on the crankshaft rotation angle. These models used physical constants obtained experimentally. The Rasselier and Withrow relationship, along with burning laws allow obtaining combustion pressure per crankshaft rotation degree. To quantify *ignition delay* there are many correlations as the one proposed by [16], [17], [18] or [19]. Models proposed by [15], [20] and [17] are used to calculate the burning factor, which is representative of the burned mass fraction in the premixed phase and the diffusive phase.

76 Internal Combustion Engines

obtained experimentally [1].

were considered in[1,2,4 and 8].

numbers was proposed in [13].

reliability [14].

analytical calculation of gas exchange considering the above-mentioned phenomena is quite complicated and requires the use of specialized computer programs that use coefficients

The basis for calculation of non-stationary non-isentropic flow characteristic of RICE inlet and outlet ducts, and NO emissions was established in [7]. Various empirical correlations to take into account heat transfer during gas exchange process and to adjust the Reynolds number exponential factor in such a way that reduce to only one the adjustment coefficients

A procedure widely used in both experimental and theoretical study of flow in engines, is to analyze the engine cycle in absence of combustion, simulating the process of compression-

In all RICE working cycle processes, there is heat transfer to the cylinder walls, which occurs with greater intensity during the combustion and expansion due to the high temperature gradients reached. Woschni [9] proposed equations to determine turbulent convective heat transfer considering average speed of in-cylinder gases and Annand [10] find correlations to calculate instant average coefficients for turbulent convection heat transfer using gas average temperature and proposed correlations to evaluate flame radiation emitted during combustion. Correlations for convection heat transfer taking into consideration surface change and cylinder enclosed volume as piston moves was established in [11]. A computer program to calculate the heat transfer in a RICE combustion chamber using models to consider turbulence was presented in [12]. A universal correlation for mixture flow in the admission and exhaust process, correcting the coefficients of Nusselt, Reynolds and Prandtl

In present study the *compression* process is considered adiabatic and reversible, but in real engines there is heat transfer between working fluid, valves and cylinder walls. At the beginning of compression the fluid temperature is lower than the temperature of the surfaces that surround the cylinder volume, causing an increase in the fluid temperature, some instant later temperatures get equal and latter on, heat is transferred from the working

The complexity of the combustion process in RICE because of untimely and incomplete combustion, dissociation and heat transfer, has encouraged the development of special techniques for carrying out studies. Adequate realization of this process is decisive in terms of the produced power and its efficiency having great influence on the engine life and

Various *models* have been proposed to study the combustion process such as the Wiebe burning law applicable to SIE and the Watson law, applicable to CIE [15]. These laws determine burned mass fraction and heat released depending on the crankshaft rotation angle. These models used physical constants obtained experimentally. The Rasselier and Withrow relationship, along with burning laws allow obtaining combustion pressure per crankshaft rotation degree. To quantify *ignition delay* there are many correlations as the one

fluid to the walls, therefore the politropic coefficient varies during the process [1].

expansion and making measurements to the engine operating at this condition [5].

Volume changes may be evaluated using an expression proposed by [15], which correlates the engine dimensions: compression ratio, displaced volume, combustion chamber volume, crank radius, connecting rod length, and crankshaft rotation angle. Average temperature during combustion process can be determined using in-cylinder pressure and ideal gas equation [4].

Calculation methods used to obtain equilibrium composition and final state of chemical species present in the combustion products of a fuel air mixture are well known and are referenced in the literature [21, 22, 23 and 24]. One of the most complete programs is perhaps the NASA-Lewis code CEC [25 and 26] which considers liquid and gaseous chemical species, is extremely versatile and can be used to calculate thermodynamic state, chemical equilibrium, rockets theoretical behavior and even Chapman-Jouguet detonation properties.

Computer programs for calculating CHO and CHON constant pressure combustion systems assuming that combustion products were composed of eight and ten chemical species were presented respectively in [27 and 28]. A code less general than the NASA code, limited to twelve chemical species CHON combustion systems specifically designed to be applied to the analysis of internal combustion engines processes was published in [29]. A program for calculating twelve species CHON constant volume combustion systems, applicable to temperatures up to 3400 K was presented in [30]. A program, valid for temperatures up to 6000 K, which can be calculate, both constant pressure and constant volume combustion for an eighteen chemical species CHON system is available in [31].

The working fluid properties function of its temperature, pressure and richness can be determined by applying thermodynamic basic equations for ideal gas mixtures, considering the mass fractions of each component in the mixture [32]. Also can be determined through routines such as FARG and ECP [33 and 34]. To complement the study of the combustion process, models to determine NO emissions as the extended Zeldovich mechanism, have been considered. The reason for the use of these models is because the specific constants of the reaction rate for NO are very small compared to the combustion rate, for this reason it is supposed that all the species present in the products with the exception of NO, are in chemical equilibrium.

The *expansion* process produces mechanical work from energy released during combustion and ends when exhaust valve opens. At this moment products are expelled from the cylinder initially at critical speed ranging between 500 and 700 m/s, and then are pushed by piston movement towards the upper dead center [4 and 15]. Towards the end of the exhaust during the valve overlap, part of the fresh mixture escapes contributing to the emission of unburned hydrocarbons and reducing the engine efficiency.

To investigate the gas exchange process, using gas dynamics to analyze the gas flow in transient processes with variable composition and variable specific heats, models such as [35] have been used.

To improve the gas exchange process we must advance inlet valve opening (AIVO) and delay exhaust valve closure (DEVC). Because of this, there is a period in which the two valves remain open simultaneously, this period is known as *valve overlap*, which helps to remove as much gas and admit as much air or fresh mixture. This is due to the depression originated in the inlet valve vicinity, due to the ejection effect produced by burned gas movement through the exhaust valve; this will contribute to increase efficiency and power produced by the RICE [1].

Two research methods are employed to study the working cycle of RICE. The first one is based on data acquisition from experimental tests and the second one is based only on mathematical simulation. The latter method is more versatile and reduces the required research empirical data depending on the employed calculating method and imposed simplifications. However, to validate mathematical simulation results, experimental parameters obtained in laboratory are required [5]. Use of numerical analysis methods has now greatly developed and increased due to increase of velocity and calculation capacity of modern computers. These methods provide faster performance, versatile and can handle more information than can be measured in an experimental test. However, results accuracy obtained by applying the models depends on made assumptions.

*Modeling* is a research technique employed in RICE, its use has grown in last two decades due to the cost decrease obtained by eliminating or reducing the laboratory tests, since these require a large amount of repetitive tests to obtain appropriate results, bringing time and money losses in preparing, calibrating, measuring, repairing and replacing testing engines. RICE designers must build more efficient engines due to higher fuel cost and new regulations on combustion emissions produced in the process that occurs inside engine. In order to optimize these designs, numerous trial and error tests are required. Implementing the tests implies expensive construction and testing of several prototypes. Modeling is a procedure that allows realizing numerous tests with relatively low cost.

To determine engine p vs. V diagram, working fluid is considered as an ideal gas, the mass entering the cylinder is calculated using a filling model that takes into account valve rise and discharge coefficient. Initial in-cylinder mass are residual gases, same quantity that was used as a reference value to control the expelled mass during exhaust. Instant volume was determined using an equation in terms of crankshaft rotation angle [15]. Final compression temperature was found from the first law of thermodynamic considering a uniform flow process and convection heat transfer. Power, mean indicated pressure and maximum pressure and temperature were calculated using the methods proposed in [1], [4] and [15]. Cyclic dispersion were studied using mean indicated pressure variation coefficient and pressure variation as a function of main combustion phase angle in the range of 10° before top dead center (TDC) and 10° after TDC, [1]. Calculations for exhaust process were similarly to those of the admissions process. A model to study the closed loop of a CIE limited pressure cycle, replacing the constant volume heat rejection process by an isentropic expansion process followed by a constant pressure heat rejection is proposed in [36].

There are commercial packages which constitute a very useful tool in the field of research and development of RICE being employed by different companies in the automotive sector. These include the ECARD (Engine Computer Aided Research & Development), developed by the IMST group, a global model that allows simulating engine operation throughout its working cycle, using similar complexity models for the different processes involved. The OpenWAM is a free, open source 1-dimensional gas-dynamics code produced by CMT group that can be used to predict the flow movement through the elements of an internal combustion engine. NEUROPART, uses neural networks to determine the product properties and composition influence on the exhaust emissions and particles formation. CHEMKIN, uses the chemical kinetic concepts to analyze fluids in gas phase through fluid dynamic simulation. EQUIL, calculates the composition at equilibrium of combustion products. PREMIX, calculates combustion speed for different fuels. SENKIN, allows to determine the time delay for different fuels and the combustion kinetic evolution depending on the species involved in a process.

## **2. Mathematical model**

78 Internal Combustion Engines

produced by the RICE [1].

To improve the gas exchange process we must advance inlet valve opening (AIVO) and delay exhaust valve closure (DEVC). Because of this, there is a period in which the two valves remain open simultaneously, this period is known as *valve overlap*, which helps to remove as much gas and admit as much air or fresh mixture. This is due to the depression originated in the inlet valve vicinity, due to the ejection effect produced by burned gas movement through the exhaust valve; this will contribute to increase efficiency and power

Two research methods are employed to study the working cycle of RICE. The first one is based on data acquisition from experimental tests and the second one is based only on mathematical simulation. The latter method is more versatile and reduces the required research empirical data depending on the employed calculating method and imposed simplifications. However, to validate mathematical simulation results, experimental parameters obtained in laboratory are required [5]. Use of numerical analysis methods has now greatly developed and increased due to increase of velocity and calculation capacity of modern computers. These methods provide faster performance, versatile and can handle more information than can be measured in an experimental test. However, results accuracy

*Modeling* is a research technique employed in RICE, its use has grown in last two decades due to the cost decrease obtained by eliminating or reducing the laboratory tests, since these require a large amount of repetitive tests to obtain appropriate results, bringing time and money losses in preparing, calibrating, measuring, repairing and replacing testing engines. RICE designers must build more efficient engines due to higher fuel cost and new regulations on combustion emissions produced in the process that occurs inside engine. In order to optimize these designs, numerous trial and error tests are required. Implementing the tests implies expensive construction and testing of several prototypes. Modeling is a

To determine engine p vs. V diagram, working fluid is considered as an ideal gas, the mass entering the cylinder is calculated using a filling model that takes into account valve rise and discharge coefficient. Initial in-cylinder mass are residual gases, same quantity that was used as a reference value to control the expelled mass during exhaust. Instant volume was determined using an equation in terms of crankshaft rotation angle [15]. Final compression temperature was found from the first law of thermodynamic considering a uniform flow process and convection heat transfer. Power, mean indicated pressure and maximum pressure and temperature were calculated using the methods proposed in [1], [4] and [15]. Cyclic dispersion were studied using mean indicated pressure variation coefficient and pressure variation as a function of main combustion phase angle in the range of 10° before top dead center (TDC) and 10° after TDC, [1]. Calculations for exhaust process were similarly to those of the admissions process. A model to study the closed loop of a CIE limited pressure cycle, replacing the constant volume heat rejection process by an isentropic

expansion process followed by a constant pressure heat rejection is proposed in [36].

obtained by applying the models depends on made assumptions.

procedure that allows realizing numerous tests with relatively low cost.

Present paragraph will develop fundamentals and mathematical equations that govern the phenomena occurring in CIE. For this purpose, volume control in Figure 1, which shows the mass and energy interactions with surroundings will be considered.

**Figure 1.** Engine control volume

It should be noted that the control volume during gas exchange processes works as an open system. During compression, combustion and expansion processes works as a closed system, that is why corrections should be made to take into account the exchanged mass due to leakage and the supply of fuel.

#### **2.1. Mass conservation**

The mass conservation principle establishes that the total mass change in a control volume is:

$$
\Delta m\_{\rm vec} = \sum m\_e - \sum m\_s \tag{1}
$$

The summation is used when there are several inputs and/or output flows. Expressing Ec. 1 in differential form and dividing by a time differential we obtain the mass time rate of change:

$$\frac{dm\_{\rm rec}}{dt} = \frac{dm\_e}{dt} - \frac{dm\_s}{dt} \tag{2}$$

To express last equation in terms of air and fuel mass entering the control volume, we define:

$$f = \frac{m\_f}{m} \tag{3}$$

Time differentiating and rearranging we obtain the fuel time rate of change:

$$\stackrel{\bullet}{f} = \frac{df}{dt} = \left(\frac{\stackrel{\bullet}{m} - m\_s}{m}\right) \left[ \left(f\_e - f\_s\right) \right] \tag{4}$$

From equivalence ratio (mixture richness) definition:

$$\phi = \frac{\frac{m\_f}{m\_a}}{\left(\frac{m\_f}{m\_a}\right)\_{\text{sto}}} = \frac{\frac{m\_f}{m\_a}}{\left(\frac{\bullet}{m\_f}\right)\_{\text{sto}}}\tag{5}$$

replacing Ec. 3 in Ec. 5 and time deriving:

$$\stackrel{\bullet}{\phi} = \frac{d\phi}{dt} = \frac{1}{\left(\frac{m\_f}{m\_a}\right)\_{sto}} \frac{\stackrel{\bullet}{f}}{\left(1 - f\right)^2} \tag{6}$$

#### **2.2. Energy conservation**

The first law of thermodynamics for an open system, disregarding changes in kinetic and potential energy can be written in differential form as:

$$\frac{dE}{dt} = \frac{d\mathbb{Q}}{dt} - \frac{d\mathbb{W}}{dt} + \stackrel{\bullet}{m}\_e h\_e - \stackrel{\bullet}{m}\_s h\_s \tag{7}$$

Since the work due to a volume change is:

$$\frac{dW}{dt} = \stackrel{\bullet}{W} = P \frac{dV}{dt} \tag{8}$$

and the first term on the left-hand side of Eq. 7 can be evaluated in terms of the internal energy:

Thermodynamic Study of the Working Cycle of a Direct Injection Compression Ignition Engine 81

$$\frac{d\mathbb{E}}{dt} = \frac{d}{dt}(mu) = \left(m\frac{du}{dt}\right)\_{\mathbb{U}^c} + \left(u\frac{dm}{dt}\right)\_{\mathbb{U}^c} \tag{9}$$

or in terms of enthalpy:

80 Internal Combustion Engines

define:

*vc e s dm dm dm*

To express last equation in terms of air and fuel mass entering the control volume, we

*mf <sup>f</sup>*

*df m m <sup>f</sup> f f dt m*

 

> *<sup>f</sup> <sup>f</sup> a a*

*<sup>m</sup> <sup>m</sup> m m*

 

1

The first law of thermodynamics for an open system, disregarding changes in kinetic and

*dE dQ dW mh mh*

*dW dV W P dt dt* 

and the first term on the left-hand side of Eq. 7 can be evaluated in terms of the internal energy:

*d f dt m f m*

 

*<sup>f</sup>* 1 *a sto*

*f*

*m*

*m*

*<sup>a</sup> sto*

<sup>2</sup>

*ee ss*

(7)

(8)

*f*

*m*

*m*

 

*dt dt dt*

*a sto*

 *e s e s*

Time differentiating and rearranging we obtain the fuel time rate of change:

From equivalence ratio (mixture richness) definition:

replacing Ec. 3 in Ec. 5 and time deriving:

potential energy can be written in differential form as:

Since the work due to a volume change is:

**2.2. Energy conservation** 

*dt dt dt* (2)

*<sup>m</sup>* (3)

(4)

(5)

(6)

$$\frac{dE}{dt} = \frac{d}{dt}(mh) - \frac{d}{dt}(pV) \tag{10}$$

Substituting Eqs. 8 and 10 in Eq. 7 we have:

$$
\left(m\frac{du}{dt}\right)\_{\rm vc} + \left(u\frac{dm}{dt}\right)\_{\rm vc} = \stackrel{\bullet}{Q} - \stackrel{\bullet}{W} + \stackrel{\bullet}{m}\_{e}h\_{e} - \stackrel{\bullet}{m}\_{s}h\_{s}\tag{11}
$$

Since internal energy, enthalpy and density are T, p and ϕ functions, their time rate of change are:

$$\frac{du}{dt} = \left(\frac{\partial u}{\partial T}\right)\frac{dT}{dt} + \left(\frac{\partial u}{\partial p}\right)\frac{dp}{dt} + \left(\frac{\partial u}{\partial \phi}\right)\frac{d\phi}{dt} \tag{12}$$

$$\frac{d\hbar}{dt} = \left(\frac{\partial\hbar}{\partial T}\right)\frac{dT}{dt} + \left(\frac{\partial\hbar}{\partial p}\right)\frac{dp}{dt} + \left(\frac{\partial\hbar}{\partial\phi}\right)\frac{d\phi}{dt} \tag{13}$$

$$\frac{d\,\rho}{dt} = \left(\frac{\partial\rho}{\partial T}\right)\frac{dT}{dt} + \left(\frac{\partial\rho}{\partial p}\right)\frac{dp}{dt} + \left(\frac{\partial\rho}{\partial\phi}\right)\frac{d\phi}{dt} \tag{14}$$

Assuming the working fluid is an ideal gas, differentiating ideal gas equation and rearranging we have:

$$p\frac{dV}{dT} + V\frac{dp}{dT} = mR\frac{dT}{dt} + mT\frac{dR}{dt} + RT\frac{dm}{dt} \tag{15}$$

$$\frac{d\rho}{dt} = \frac{\frac{dp}{dt} - \rho \mathcal{R}\frac{dT}{dt} - \rho \Gamma \frac{dR}{dt}}{RT} \tag{16}$$

From Eq. 14:

$$\frac{dp}{dt} = \frac{\frac{d\rho}{dt} - \left(\frac{\partial\rho}{\partial T}\right)\frac{dT}{dt} - \left(\frac{\partial\rho}{\partial\phi}\right)\frac{d\phi}{dt}}{\left(\frac{\partial\rho}{\partial p}\right)}\tag{17}$$

substituting Eq. 17 in Eq. 16, rearranging and solving for *dp dT* :

$$\frac{dp}{dt} = \frac{-\rho \frac{dT}{T} \frac{dT}{dt} - \frac{\rho}{R} \frac{dR}{dt} - \left(\frac{\partial \rho}{\partial T}\right) \frac{dT}{dt} - \left(\frac{\partial \rho}{\partial \phi}\right) \frac{d\phi}{dt}}{\left(\frac{\partial \rho}{\partial p}\right) - \frac{1}{RT}}\tag{18}$$

Solving Ec. 15 for *dR dt* , simplifying and substituting in Ec. 18:

$$\frac{dp}{dt}\left|\left(\frac{\partial\rho}{\partial p}\right) - \frac{1}{RT}\right| = -\frac{\rho}{T}\frac{dT}{dt} - \frac{\rho}{R}\left|\frac{p}{mT}\frac{dV}{dt} + \frac{V}{mT}\frac{dp}{dt} - \frac{R}{T}\frac{dT}{dt} - \frac{R}{m}\frac{dm}{dt}\right| - \left(\frac{\partial\rho}{\partial T}\right)\frac{dT}{dt} - \left(\frac{\partial\rho}{\partial\phi}\right)\frac{d\phi}{dt} \tag{19}$$

Replacing ideal gas equation in Ec. 19 and solving for *dp dT* :

$$\frac{dp}{dt} = \frac{\rho}{\left(\frac{\partial \rho}{\partial p}\right)} \left[ \frac{\frac{dV}{dt}}{V} + \frac{\frac{dm}{dt}}{m} - \frac{1}{\rho} \left( \frac{\partial \rho}{\partial T} \right) \frac{dT}{dt} - \frac{1}{\rho} \left( \frac{\partial \rho}{\partial \phi} \right) \frac{d\phi}{dt} \right] \tag{20}$$

Differentiating ideal gas equation with respect to p and T, we obtain:

$$
\left(\frac{\partial \rho}{\partial p}\right) = \frac{1}{RT} \tag{21}
$$

$$
\left(\frac{\partial \rho}{\partial T}\right) = -\frac{p}{RT^2} \tag{22}
$$

and substituting Ecs. 21 and 22 in Ec. 20:

$$\frac{dp}{dt} = p\left[ -\frac{\overline{d}V}{V} + \frac{\frac{dm}{dt}}{m} + \left(\frac{dT}{dt}\right)\frac{1}{T} - \frac{RT}{p}\left(\frac{\partial\rho}{\partial\phi}\right)\frac{d\phi}{dt} \right] \tag{23}$$

Later equation, function of density, volume, mass and mixture richness will be used to obtain the in-cylinder pressure when time varies (indicator diagram).

The procedure to obtain a similar expression for temperature variation over time will be illustrated below. Solving Eq. 11 for *du dt* :

$$\frac{du}{dt} = \frac{\stackrel{\bullet}{Q}}{m\_{\rm rc}} - \frac{p}{m\_{\rm rc}} \frac{dV}{dt} + \frac{1}{m\_{\rm rc}} \left( m\_{\rm c} h\_{\rm c} - m\_{\rm s} h\_{\rm s} - \left( u \frac{dm}{dt} \right)\_{\rm rc} \right) \tag{24}$$

and defining:

Thermodynamic Study of the Working Cycle of a Direct Injection Compression Ignition Engine 83

$$B = -RT\frac{dV}{V} + \frac{1}{m} \left(\stackrel{\bullet}{Q} + \stackrel{\bullet}{m}h\_e \stackrel{\bullet}{h}\_e - \stackrel{\bullet}{m}s\_s h\_s - \left(\stackrel{\bullet}{u}\frac{dm}{dt}\right)\_{vc}\right) \tag{25}$$

On the other hand, introducing Ecs. 8 in Ec. 11, get the following expression:

$$\frac{d\mu}{dt} = B\tag{26}$$

Substituting Ec. 26 in Ec. 12 and solving for *dp dt* :

$$\frac{dp}{dt} = \frac{B - \left(\frac{\partial u}{\partial t}\right)\frac{dT}{dt} - \left(\frac{\partial u}{d\phi}\right)\frac{d\phi}{dt}}{\left(\frac{\partial u}{\partial p}\right)}\tag{27}$$

Replacing Ec. 27 in Ec. 20 and solving for *dT dt* gives:

$$\frac{dT}{dt} = \frac{\left(\frac{\partial u}{\partial p}\right) \left[ -\rho \frac{\frac{dV}{dt}}{V} + \rho \frac{\frac{dm}{dt}}{m} - \left(\frac{\partial \rho}{\partial \phi}\right) \frac{d\phi}{dt} - B + \left(\frac{\partial u}{\partial \phi}\right) \frac{d\phi}{dt} \right]}{\left(\frac{\partial u}{\partial p}\right) \left(\frac{\partial \rho}{\partial T}\right) - \left(\frac{\partial \rho}{\partial p}\right) \left(\frac{\partial u}{\partial T}\right)}\tag{28}$$

Now, considering:

82 Internal Combustion Engines

Solving Ec. 15 for *dR*

 1

 

(18)

 

(20)

(23)

(24)

(21)

(22)

*dT dR dT d*

 

*p RT*

(19)

*dT* :

> 

 

1 1

1 *p RT*

> 2 *p*

> > 1

*ee ss*

*vc vc vc vc*

 

Later equation, function of density, volume, mass and mixture richness will be used to

The procedure to obtain a similar expression for temperature variation over time will be

1

*du Q dV <sup>p</sup> dm mh mh u dt m m dt m dt*

 

*T RT*

*dp T dt R dt T dt dt*

*dp* 1 *dT dV V R dT R dm dT d p dp dt p RT T dt R mT dt mT dt T dt m dt T dt dt*

> *dp dt dt dT d dt V m T dt dt*

> >

 

*dV dm*

obtain the in-cylinder pressure when time varies (indicator diagram).

*dt* :

*dp dt dt dT RT d <sup>p</sup> dt V m dt T p dt*

*dt* , simplifying and substituting in Ec. 18:

*dV dm*

*dt*

and substituting Ecs. 21 and 22 in Ec. 20:

illustrated below. Solving Eq. 11 for *du*

and defining:

Replacing ideal gas equation in Ec. 19 and solving for *dp*

*p*

Differentiating ideal gas equation with respect to p and T, we obtain:

  

$$R = R(T, p, \phi, t) \tag{29}$$

And differentiating:

$$\frac{dR}{dt} = \left(\frac{\partial R}{\partial T}\right)\frac{dT}{dt} + \left(\frac{\partial R}{\partial p}\right)\frac{dp}{dt} + \left(\frac{\partial R}{\partial \phi}\right)\frac{d\phi}{dt} \tag{30}$$

Differentiating ideal gas equation and solving for *dR dt*

$$\frac{dR}{dt} = \frac{\frac{dp}{dt}}{p}R - \frac{dT}{T}R - \frac{d\rho}{\rho}R\tag{31}$$

Replacing Ec. 30 in Eq. 31 and solving for *dp dt* :

$$\frac{dp}{dt} = \frac{\left(\frac{\partial R}{\partial T}\right) \frac{dT}{dt} + \left(\frac{\partial R}{\partial \phi}\right) \frac{d\phi}{dt} + \frac{dT}{dt} \frac{R}{T} + \frac{d\rho}{dt} \frac{R}{\rho}}{T} \tag{32}$$

$$\frac{R}{p} - \left(\frac{\partial R}{\partial p}\right)$$

Replacing Ec. 32 in Ec. 27:

$$
\left[ \left( \frac{\partial u}{\partial T} \right) \frac{dT}{dt} + \left( \frac{\partial u}{\partial p} \right) \right] \frac{\left( \frac{\partial R}{\partial T} \right) \frac{dT}{dt} + \left( \frac{\partial R}{\partial \phi} \right) \frac{d\phi}{dt} + \frac{dT}{dt} \frac{R}{T} + \frac{d\rho}{dt} \frac{R}{\rho} \right]}{\frac{R}{p} - \left( \frac{\partial R}{\partial p} \right)} + \left( \frac{\partial u}{\partial \phi} \right) \frac{d\phi}{dt} = B \tag{33}
$$

Defining:

$$D = 1 - \frac{p}{R} \left(\frac{\partial R}{\partial p}\right) \tag{34}$$

Collecting terms containing *dT dt* and substituting Eq. 34 in Eq. 33:

$$
\left(\frac{dT}{dt}\right)\left[\left(\frac{\partial u}{\partial T}\right) + \left(\frac{\partial u}{\partial p}\right)\frac{p}{DR}\left\{\left(\frac{\partial R}{\partial T} + \frac{R}{T}\right)\right\}\right] + \left(\frac{\partial u}{\partial p}\right)\frac{1}{D}\left[\left(\frac{\partial R}{\partial \phi}\right)\frac{d\phi}{dt} + \frac{d\phi}{dt}\frac{R}{\rho}\right] + \left(\frac{\partial u}{\partial \phi}\right)\frac{d\phi}{dt} = B \tag{35}
$$

Defining:

$$C = 1 + \frac{T}{R} \left(\frac{\partial R}{\partial T}\right) \tag{36}$$

Replacing Ec. 36 in Eq. 35 and solving for *dT dt* :

$$\frac{dT}{dt} = \frac{B - \left(\frac{\partial u}{\partial p}\right) \frac{p}{D} \left[\frac{1}{R} \left(\frac{\partial R}{\partial \phi}\right) \frac{d\phi}{dt} + \frac{d\rho}{dt} \frac{1}{\rho}\right] - \left(\frac{\partial u}{\partial \phi}\right) \frac{d\phi}{dt}}{\left(\frac{\partial u}{\partial T}\right) + \frac{C}{D} \frac{p}{T} \left(\frac{\partial u}{\partial p}\right)}\tag{37}$$

Since:

$$
\rho m = \rho V \tag{38}
$$

Time differentiating and solving for *d dt* : Thermodynamic Study of the Working Cycle of a Direct Injection Compression Ignition Engine 85

$$\frac{d\rho}{dt} = \frac{dm}{dt} - \rho \frac{dV}{dt} \tag{39}$$

Replacing Ec. 39 in Ec. 37:

$$\frac{dT}{dt} = \frac{B - \left(\frac{\partial u}{\partial p}\right) \frac{p}{D} \left[ \left(\frac{\partial R}{\partial \phi}\right) \frac{d\phi}{dt} \frac{1}{R} + \frac{dm}{dt} \frac{1}{m} - \frac{dV}{dt} \frac{1}{V} \right] - \left(\frac{\partial u}{\partial \phi}\right) \frac{d\phi}{dt}}{\left(\frac{\partial u}{\partial T}\right) + \frac{C}{D} \frac{p}{T} \left(\frac{\partial u}{\partial p}\right)}\tag{40}$$

This equation will be used to determine the in-cylinder temperature when time varies. If Eqs. 25, 34, 36 and 37 are replaced in Eq. 23 and collecting terms we obtained:

$$\frac{dp}{dt} = \frac{\stackrel{\bullet}{Q} - m\_{bb} \, h\_{bb} - \frac{dV}{dt} \left[ \frac{m \, \mathrm{C}}{V} + p \right] - \frac{dm}{dt} \left[ \left( D \left( \frac{\partial u}{\partial \phi} \right) - h\_{cl} + u \right) - \mathrm{C} \left( \frac{D}{R} \left( \frac{\partial R}{\partial \phi} + 1 \right) \right) \right]}{m \left[ \mathrm{C} \left( \frac{1}{p} - \frac{1}{R} \left( \frac{\partial R}{\partial \phi} \right) \right) + \left( \frac{\partial u}{\partial p} \right) \right]} \tag{41}$$

and:

84 Internal Combustion Engines

Replacing Ec. 32 in Ec. 27:

Collecting terms containing

Replacing Ec. 36 in Eq. 35 and solving for

Time differentiating and solving for

*dT*

Defining:

Defining:

Since:

*R dT R d dT R d R*

 

> 

(32)

 

(37)

(33)

(34)

 

(36)

 

(38)

*p p*

*R dT R d dT R d R u dT u T dt dt dt T dt u d <sup>B</sup> T dt p R R dt p p*

<sup>1</sup> *<sup>p</sup> <sup>R</sup> <sup>D</sup> R p*

*dt* and substituting Eq. 34 in Eq. 33:

*dT u u R R u R d d R u d <sup>p</sup>* <sup>1</sup> *<sup>B</sup> dt T p DR T T p D dt dt dt*

(35)

<sup>1</sup> *T R <sup>C</sup> R T* 

*u Rd d ud <sup>p</sup>* 1 1 *<sup>B</sup> dT p D R dt dt dt*

*m V* 

 

*T DT p*

*dT dt* :

*dt uC u p*

*d dt* 

:

*dp T dt dt dt T dt*

*dt R R*

$$\frac{dm}{dt} = \frac{\stackrel{\bullet}{Q} - m\_{bb} \, h\_{bb} - \frac{dV}{dt} \left[ \frac{m\stackrel{\bullet}{C}}{V} + p \right] - m \left[ \mathcal{C} \left( \frac{1}{p} - \frac{1}{R} \left( \frac{\partial R}{\partial \phi} \right) \right) + \left( \frac{\partial u}{\partial p} \right) \right] \frac{dp}{dt}}{\left( D \left( \frac{\partial u}{\partial \phi} \right) - h\_{cl} + u \right) - \mathcal{C} \left( \frac{D}{R} \left( \frac{\partial R}{\partial \phi} + 1 \right) \right)} \tag{42}$$

Eqs. 41 and 42 will be used to obtain the indicator diagram (p vs. V or p vs. φ diagram) and burned mass fraction diagram (m vs. t diagram) respectively.

#### **2.3. Instant in-cylinder volume**

The instant volume inside the control volume in terms of the displaced volume, compression ratio, connecting rod length to crank radius ratio and crankshaft rotation angle, can be obtained through the following expression [15]:

$$V(\rho) = V\_d \left[ \frac{1}{r\_c - 1} + \frac{1}{2} \left[ R\_{LA} + 1 - c \cos \rho - (R\_{LA}^2 - \text{sen}^2 \rho)^2 \right] \right] \tag{43}$$

In this expression RLA is the connecting rod length (l) to crank radius (a) ratio.

$$R\_{LA} = \frac{l}{a} \tag{44}$$

Deriving Ec. 43 with respect to crankshaft rotation angle, we obtain:

$$\frac{dV}{d\rho} = \frac{V\_d}{2} \left[ \begin{array}{c} \kappa n \rho \cos \phi \\\\ \left( R\_{LA}^2 - \operatorname{sen}^2 \phi \right)^{\frac{1}{2}} \end{array} \right] \tag{45}$$

The time in seconds it takes to describe some crankshaft rotation angle can be calculated with the following expression:

$$t = \frac{\wp}{6 \text{ rpm}}\tag{46}$$

Solving previous expression for φ and replacing in Ec. 45 to make the corresponding conversion from degrees to radians we obtained:

$$\frac{dV}{dt} = 3V\_d (rpm) \left[ \text{sen} \frac{\pi \, rpm}{30} t + \frac{\text{sen} \frac{\pi \, rpm}{30} t \cos \frac{\pi \, rpm}{30} t}{\left( R\_{LA}^2 - \text{sen}^2 \frac{\pi \, rpm}{30} t \right)^{\frac{1}{2}}} \right] \tag{47}$$

Previous expression allows determining the in-cylinder volume variation with respect to time, while Ec. 45 will be used to calculate the volume variation with respect to crankshaft rotation angle.

### **3. Equations, models and calculations**

Models and assumptions used to analyze each of the thermodynamic processes that are carried out in a CIE will be presented in this paragraph. Routines commonly used in RICE are employed to calculate the thermodynamic properties of the chemical species formed during combustion. FARG and ECP routines [34] are used to determine properties depending on gas temperature. PER routine [29] is used to obtain same properties depending on mixture richness. DVERK routine [37] found in the International Mathematics and Statistics Library is used for solving differential equations systems by the Runge - Kutta Verner fifth and sixth order method.

#### **3.1. Admission process**

The parameter characterizing the admissions process is the volumetric efficiency defined as:

$$\eta\_v = \frac{\stackrel{\bullet}{m}\_{ar}}{\stackrel{\bullet}{m}\_{at}} = \frac{\stackrel{\bullet}{m}\_{ar}}{\rho\_0 i V\_d \frac{rpm}{300}}\tag{48}$$

It takes into account the losses in the inlet valve and all the admission system if atmospheric density value is used for ρ0.

Real mass air flow entering the cylinder is determined by the following equations [15] function of the ratio pdown/pup:

86 Internal Combustion Engines

rotation angle.

with the following expression:

 

*R sen*

*LA*

The time in seconds it takes to describe some crankshaft rotation angle can be calculated

6 *t*

Solving previous expression for φ and replacing in Ec. 45 to make the corresponding

30 30 3( ) <sup>30</sup>

Previous expression allows determining the in-cylinder volume variation with respect to time, while Ec. 45 will be used to calculate the volume variation with respect to crankshaft

Models and assumptions used to analyze each of the thermodynamic processes that are carried out in a CIE will be presented in this paragraph. Routines commonly used in RICE are employed to calculate the thermodynamic properties of the chemical species formed during combustion. FARG and ECP routines [34] are used to determine properties depending on gas temperature. PER routine [29] is used to obtain same properties depending on mixture richness. DVERK routine [37] found in the International Mathematics and Statistics Library is used for solving differential equations systems by the Runge - Kutta

The parameter characterizing the admissions process is the volumetric efficiency defined as:

*m m*

*m iV*

*ar ar*

*at d*

It takes into account the losses in the inlet valve and all the admission system if atmospheric

*v*

<sup>0</sup> 30

*rpm*

*j*

(48)

*rpm rpm sen t cos t dV rpm V rpm sen t dt*

*LA*

*rpm* 

2 *d*

*d*

conversion from degrees to radians we obtained:

**3. Equations, models and calculations** 

Verner fifth and sixth order method.

**3.1. Admission process** 

density value is used for ρ0.

*d*

*dV V sen cos sen*

2 2 2

1

(46)

1

<sup>2</sup> 2 2

*rpm R sen t*

30

 

(45)

(47)

$$\frac{p\_{down}}{p\_{up}} < 1 \qquad \qquad \qquad m = \frac{C\_d A\_{ref} p\_0}{\sqrt{R\_a T\_0}} \left( \frac{p\_{down}}{p\_{up}} \right) \left( \frac{2\gamma}{\gamma - 1} \left| 1 - \left( \frac{p\_{down}}{p\_{up}} \right)^{\frac{\gamma + 1}{\gamma}} \right| \right) \tag{49}$$

$$\frac{p\_{down}}{p\_{up}} \ge 1 \qquad \qquad \overset{\bullet}{m} = \mathbb{C}\_d A\_{ref} p\_{up} \sqrt{\left(\frac{\gamma}{RT\_0}\right) \left(\frac{2}{\gamma - 1}\right)^{\frac{\gamma + 1}{2(\gamma + 1)}}} \tag{50}$$

where pdown is the downstream pressure and pup is the upstream pressure. Although the discharge coefficient Cd varies during the process, in present study, we assume that it is constant and equal to its average value. The reference area Aref, usually called curtain area, since it depends on the valve lifting Lv, is taken as:

$$A\_{ref} = \pi \, d\_v \, \mathcal{L}\_v \tag{51}$$

A model proposed in [38], was used to theoretically determine the lifting valve profile which is function of the maximum lifting, Lv max and the crankshaft rotation angle φ:

$$L\_v(\boldsymbol{\wp}) = L\_{v\max} + C\_2 \boldsymbol{\wp}^2 + C\_p \boldsymbol{\wp}^p + C\_q \boldsymbol{\wp}^q + C\_r \boldsymbol{\wp}^r + C\_s \boldsymbol{\wp}^s \tag{52}$$

Coefficients C2, Cp, Cq, Cr y Cs are determined with the following equations:

$$C\_2 = \frac{-pqrsh}{\left[ (p-2)(q-2)(r-2)(s-2)c\_{med}^2 \right]} \tag{53}$$

$$C\_p = \frac{2qrsh}{\left[ \left( p - 2 \right) \left( q - p \right) \left( r - p \right) \left( s - p \right) c\_{mcl}^p \right]} \tag{54}$$

$$C\_q = \frac{-2\,\text{prsh}}{\left[ \left( q - 2 \right) \left( q - p \right) \left( r - q \right) \left( s - q \right) c\_{mcd}^q \right]} \tag{55}$$

$$C\_r = \frac{2pqsh}{\left[\left(r-2\right)\left(r-p\right)\left(r-q\right)\left(s-r\right)c\_{med}^r\right]}\tag{56}$$

$$C\_s = \frac{-2pqrh}{\left[\left(s-2\right)\left(s-p\right)\left(s-q\right)\left(s-r\right)c\_{med}^s\right]}\tag{57}$$

Recommended values for p, q, r, and s are: p = 6; q = 8; r = 10; s = 12.

Gas pressure and temperature variation over time in this process is calculated from Eqs. 23 and 40. Since a CIE only compresses air, the term corresponding to mixture richness variation with time is zero. For this reason the above equations are:

$$\frac{dp}{dt} = p\left[-\frac{dV}{V} + \frac{dm}{m} + \left(\frac{dT}{dt}\right)\frac{1}{T}\right] \tag{58}$$

$$\left(\frac{1}{2} + \frac{1}{2}\right)\left(\frac{1}{2} + \frac{1}{2} + \dots + \frac{1}{2}\right)$$

$$\frac{dT}{dt} = \frac{B - \left(\frac{\partial u}{\partial p}\right) \frac{p}{D} \left| \frac{dm}{dt} \frac{1}{m} - \frac{dV}{dt} \frac{1}{V} \right|}{\left(\frac{\partial u}{\partial T}\right) + \frac{C}{D} \frac{p}{T} \left(\frac{\partial u}{\partial p}\right)}\tag{59}$$

To solve these equations heat release, heat transfer, blow by, ignition delay and chemical species formation models are required. In addition the terms *<sup>u</sup> T* , *<sup>u</sup> p* , *<sup>R</sup> T* <sup>y</sup> *R p* must be determined by using routines FARG and ECP.

With Eqs. 43 and 47 we calculate the volume and its time derivative, respectively, while with Eqs. 49 or 50 depending on the case determine the flow mass. The in-cylinder accumulated mass is obtained by summation of the mass flows times the values obtained by Eq. 46.

## **3.2. Closed loop cycle**

Closed loop cycle corresponds to compression, combustion and expansion processes. Compression starts when inlet valve closes. Variation in pressure and temperature with time during this process is determined taking into account that only air is compressed (Ecs 58 and 59) and there are mass losses due to blow by. When fuel injection begins the mixture composition varies; therefore the expressions used to determine temperature and pressure time variation during the closed loop cycle are Eq. 23 and 40. When exhaust valve opens, begins the exhaust process.

### **3.3. Exhaust process**

Equations used in this process are the same used during the admissions process, but noting the working fluid is a mixture of burned gases and heat transfer is higher than during admission because of the high temperature present.

During valve overlap we want to extract as much as possible of burned gases and taking advantage of the dynamic effects, increase the amount of fresh charge entering the cylinder. Equations used during this process are the same used during the intake and exhaust, but considering that there is simultaneously fresh charge entry and burned gases exit.

#### **3.4. Ignition delay model**

88 Internal Combustion Engines

Eq. 46.

**3.2. Closed loop cycle** 

begins the exhaust process.

**3.3. Exhaust process** 

Gas pressure and temperature variation over time in this process is calculated from Eqs. 23 and 40. Since a CIE only compresses air, the term corresponding to mixture richness

> *dV dm dp dt dt dT <sup>p</sup> dt V m dt T*

*u dm dV <sup>p</sup>* 1 1 *<sup>B</sup> dT p D dt m dt V dt uC u p*

To solve these equations heat release, heat transfer, blow by, ignition delay and chemical

With Eqs. 43 and 47 we calculate the volume and its time derivative, respectively, while with Eqs. 49 or 50 depending on the case determine the flow mass. The in-cylinder accumulated mass is obtained by summation of the mass flows times the values obtained by

Closed loop cycle corresponds to compression, combustion and expansion processes. Compression starts when inlet valve closes. Variation in pressure and temperature with time during this process is determined taking into account that only air is compressed (Ecs 58 and 59) and there are mass losses due to blow by. When fuel injection begins the mixture composition varies; therefore the expressions used to determine temperature and pressure time variation during the closed loop cycle are Eq. 23 and 40. When exhaust valve opens,

Equations used in this process are the same used during the admissions process, but noting the working fluid is a mixture of burned gases and heat transfer is higher than during

During valve overlap we want to extract as much as possible of burned gases and taking advantage of the dynamic effects, increase the amount of fresh charge entering the cylinder. Equations used during this process are the same used during the intake and exhaust, but

considering that there is simultaneously fresh charge entry and burned gases exit.

*T DT p*

1

*T* , *<sup>u</sup> p* , *<sup>R</sup> T* <sup>y</sup>

(58)

(59)

*R p* 

must be

variation with time is zero. For this reason the above equations are:

species formation models are required. In addition the terms *<sup>u</sup>*

determined by using routines FARG and ECP.

admission because of the high temperature present.

Ignition delay in CIE characterizes the heat quantity will release immediately occur the fuel auto ignition and has a direct influence on engine rumble and pollutants formation. The model presented in [19] points out that ignition delay depends on in-cylinder temperature and pressure, engine speed and accumulated fuel amount and may be calculated in degrees and in milliseconds with the following expressions:

$$ID\_{\text{[deg]}} = Ap\_c^\prime(EA) \tag{60}$$

$$ID\_{\lceil \frac{ms}{\lceil ms \rceil} \rceil} = \frac{ID\_{\lceil \deg \rceil}}{0.006(rpm)} \tag{61}$$

where: A=0.36+0.22V , mp mp c(rpm) V= , <sup>30</sup> n=0, <sup>u</sup> c c 1 1 21.2 EA=exp R - , RT 17190 p -12.4 

$$\mathbf{E} = \frac{618840}{\text{NC} + 25}, \quad \mathbf{P\_c = \mathbf{P\_{amb}} \mathbf{r\_c}^{n\_c}, \quad \mathbf{T\_c = \mathbf{T\_{amb}} \mathbf{r\_c}^{n\_c - 1}, \quad n\_c = 1.30 \quad \text{a } 1.37, \mathbf{R\_u = 8.3143 \left[ \text{J/mol} \,\text{K} \right]} \text{J}$$

Other models based on experimental data suggest correlations that use an Arrhenius expression similar to that proposed in [15], in which the constants estimated by [39] are as

$$\text{follows: A=3.45, }\ n = 1.02, \quad \text{EA=exp}\left\lfloor \frac{\text{E}\_{\text{a}}}{\text{R}\_{\text{u}}\text{T}\_{\text{c}}} \right\rfloor \stackrel{\text{E}}{\text{R}\_{\text{u}}} = 2100.1$$

Another model whose constants are the same as in previous case use an A term function of richness as shown by the following expression [16]: -0.2 A=2.4 . 

#### **3.5. Heat release model**

Considering the fourth term numerator in Eq. 41, which represents the heat released during the combustion process and applying the Watson relationship we obtain the following equation:

$$m\_c H\_i dX\_b = \frac{dm}{dt} \left[ \left( D \left( \frac{\partial u}{\partial \phi} \right) - h\_{cil} + u \right) - \mathbf{C} \left( \frac{D}{R} \left( \frac{\partial R}{\partial \phi} + 1 \right) \right) \right] \tag{62}$$

Model of fuel apparent burning will be used to represent the combustion process. It uses two empirical equations, one for the pre-mixed combustion phase and another for the diffusive combustion phase. The instantaneous total amount of heat released by crankshaft rotation degree is given by the sum of the two components:

$$
\left(\frac{dm\_c}{d\rho}\right)\_{Tot} = \left(\frac{dm\_c}{d\rho}\right)\_{pre} + \left(\frac{dm\_c}{d\rho}\right)\_{dif} \tag{63}
$$

#### **3.6. Burning factor**

Heat release model requires defining, depending on the process physical condition, the initial fuel amount burned during the pre-mixed phase. An initial fuel burning factor is used for this purpose [15] [17]. This factor estimated, depending on initial richness and delay period, what part of the injected fuel is burned during the pre-mixed phase. The difference is burned during the diffusive phase. The burning factor is defined as:

$$
\beta = \frac{\left(m\_c\right)\_{pre}}{\left(m\_c\right)\_{Tot}} \tag{64}
$$

and may be calculated by the following expression [15]:

$$\beta = 1 - \frac{a\_1 \phi\_0^{b\_1}}{ID\_S^{cc\_1}} \tag{65}$$

The a1, b1 and cc1 values which are shown in Table 1 [39, 19 and 15] depend on the used model.


**Table 1.** Empirical values for burning factor

Taking into account the heat released during each phase Ec. 63 becomes

$$
\left(\frac{dm\_c}{d\rho}\right)\_{Tot} = \beta \left(\frac{dm\_c}{d\rho}\right)\_{pre} + (1-\beta) \left(\frac{dm\_c}{d\rho}\right)\_{dif} \tag{66}
$$

Heat released during each phase is evaluated by the empirical expressions proposed in [39] and [40]. Equations proposed in [39] are:

$$\left(\frac{d\mathbf{X}\_b}{d\rho}\right)\_{pre} = \mathbf{C}\_1 \mathbf{C}\_2 \left(\frac{\rho - \rho\_0}{\Delta \rho}\right)^{\mathbf{C}\_1 - 1} \left(1 - \left(\frac{\rho - \rho\_0}{\Delta \rho}\right)^{\mathbf{C}\_1}\right)^{\mathbf{C}\_2 - 1} \tag{67}$$

$$\left(\frac{d\mathbf{X}\_b}{d\rho}\right)\_{\mathrm{dif}} = \mathbf{C}\_3 \mathbf{C}\_4 \left(\frac{\rho - \rho\_0}{\Delta \rho}\right)^{\mathbf{C}\_4 - 1} \exp\left(\mathbf{C}\_3 - \left(\frac{\rho - \rho\_0}{\Delta \rho}\right)^{\mathbf{C}\_4}\right) \tag{68}$$

Equation proposed in [40] that use the duration and heat released percentage in each phase, unlike [39] equations, is:

$$\frac{dQ}{d\rho} = a \left(\frac{Q\_{\rm pre}}{\rho\_{\rm pre}}\right) m\_{\rm pre} \left(\frac{\rho}{\rho\_{\rm pre}}\right)^{m\_{\rm pre}-1} \exp\left[-a \left(\frac{\rho}{\rho\_{\rm pre}}\right)^{m\_{\rm pre}}\right] + \left(\frac{Q\_{\rm dif}}{\rho\_{\rm dif}}\right) m\_{\rm dif} \left(\frac{\rho}{\rho\_{\rm dif}}\right)^{m\_{\rm df}-1} \exp\left[-a \left(\frac{\rho}{\rho\_{\rm dif}}\right)^{m\_{\rm df}}\right] \tag{69}$$

Table 2 shows the constants for Heywood [15] and Watson [39] equations and Table 3 shows the constants for Miyamoto [40] equation.


**Table 2.** Constants for Heywood and Watson expressions


**Table 3.** Constants for Miyamoto expression

### **3.7. Blow by model**

90 Internal Combustion Engines

**3.6. Burning factor** 

model.

Heat release model requires defining, depending on the process physical condition, the initial fuel amount burned during the pre-mixed phase. An initial fuel burning factor is used for this purpose [15] [17]. This factor estimated, depending on initial richness and delay period, what part of the injected fuel is burned during the pre-mixed phase. The difference

> *c pre c Tot*

1 1 1 0 1 *b cc S*

*a ID* 

The a1, b1 and cc1 values which are shown in Table 1 [39, 19 and 15] depend on the used

Value Hardenberg model Watson model Heywood model

a1 0.746 0.926 0.80 - 0.95

b1 0.35 0.37 0.25 - 0.45

cc1 0.35 0.26 0.25 – 0.50

 1 *cc c Tot pre dif*

Heat released during each phase is evaluated by the empirical expressions proposed in [39]

1 2 1

3 4 <sup>3</sup> exp

*dX C C <sup>C</sup>*

 

 

> 

<sup>2</sup> 1 1 <sup>1</sup> <sup>1</sup>

4 4 1 0 0

 

> 

*C C*

*<sup>C</sup> C C*

 

0 0

(66)

(67)

(68)

*dm dm dm dd d*

 

Taking into account the heat released during each phase Ec. 63 becomes

(64)

(65)

*m m*

is burned during the diffusive phase. The burning factor is defined as:

and may be calculated by the following expression [15]:

**Table 1.** Empirical values for burning factor

and [40]. Equations proposed in [39] are:

*b pre dX C C*

*d*

*b dif*

*d*

Due to the in-cylinder pressure increase during compression and combustion processes, a part of the gasses (mbb) is lost through the rings resulting in a produced power reduction. The model that will be described below takes into account such losses.

From logarithmic derivative of the ideal gas equation:

$$\frac{1}{p}\frac{dp}{d\rho} + \frac{1}{V}\frac{dV}{d\rho} = \frac{1}{m}\frac{dm}{d\rho} + \frac{1}{T}\frac{dT}{d\rho} \tag{70}$$

Applying the First Law of Thermodynamics in differential form to an open system, we obtain:

$$m\mathbf{C}\_v \frac{dT}{d\rho} + \mathbf{C}\_v T \frac{dm}{d\rho} = \frac{d\mathbf{Q}}{d\rho} - p\frac{dV}{d\rho} - \frac{m\_{lb}\mathbf{C}\_p T}{\alpha} \tag{71}$$

Replacing *dT/dφ* from Ec. 70 and rearranging:

$$\frac{dp}{d\rho} = -\gamma \frac{p}{V} \frac{dV}{d\rho} + \frac{\gamma - 1}{V} \frac{dQ}{d\rho} - \frac{\gamma}{\rho m} \frac{m\_{bb}}{\rho m} p \tag{72}$$

Blow by mass can be found from mass conservation:

$$\frac{dm}{d\rho} = -\frac{\stackrel{\bullet}{m}}{\stackrel{\bullet}{\rho}}\tag{73}$$

Defining:

$$C = \frac{\stackrel{\bullet}{m}\_{bb}}{m} \tag{74}$$

Net heat entering the system, considering that a part is lost by blow by, is:

$$dQ = Q\_{\epsilon}d\chi - dQ\_{\text{los}} \tag{75}$$

with:

$$\frac{dQ\_{\text{lost}}}{dt} = h\_c A (T - T\_{\text{wall}}) \tag{76}$$

and defining the following dimensionless variables:

$$
\tilde{p} = \frac{p}{p\_1} \tag{77}
$$

$$
\tilde{V} = \frac{V}{V\_1} \tag{78}
$$

$$
\tilde{T} = \frac{T}{T\_1} \tag{79}
$$

$$
\tilde{Q} = \frac{Q\_e}{P\_1 V\_1} \tag{80}
$$

$$
\tilde{Q}\_{\text{lost}} = \frac{Q\_{\text{lost}}}{P\_1 V\_1} \tag{81}
$$

Thermodynamic Study of the Working Cycle of a Direct Injection Compression Ignition Engine 93

$$\tilde{h}\_c = \frac{hT\_1 \left(A\_{cc} - \frac{4V\_{cc}}{D\_p}\right)}{P\_1V\_1\alpha} \tag{82}$$

$$\beta = \frac{4V\_1}{D\_p \left(A\_{cc} - \frac{4V\_{cc}}{D\_p}\right)}\tag{83}$$

dimensionless heat losses for unit crankshaft rotation angle are:

92 Internal Combustion Engines

Defining:

with:

*v v*

Net heat entering the system, considering that a part is lost by blow by, is:

Blow by mass can be found from mass conservation:

and defining the following dimensionless variables:

Replacing *dT/dφ* from Ec. 70 and rearranging:

*mC C T <sup>p</sup> d dd d*

 

 

*d*

*dp p dV dQ m* 1 *bb <sup>p</sup> d Vd V d m* 

*dm mbb*

*mbb <sup>C</sup> m*

( ) *lost c wall*

> 1 *<sup>p</sup> <sup>p</sup>*

1 *<sup>V</sup> <sup>V</sup>*

1 *<sup>T</sup> <sup>T</sup>*

1 1

*lost*

1 1 *lost*

*dQ hAT T*

*bb p*

(71)

(72)

(73)

(74)

*e lost dQ Q dx dQ* (75)

*dt* (76)

*<sup>p</sup>* (77)

*<sup>V</sup>* (78)

*<sup>T</sup>* (79)

*Qe <sup>Q</sup> P V* (80)

*<sup>Q</sup> <sup>Q</sup> P V* (81)

 

*dT dm dQ dV m CT*

 

> 

 

$$\frac{d\bar{Q}\_{\rm lost}}{d\varphi} = \tilde{h}\_c \left( 1 + \beta \tilde{V} \right) \left( \tilde{p} \tilde{V} - \tilde{T}\_{\rm wall} \right) \tag{84}$$

Replacing Ecs. 77 to 83 in Ecs. 72, 73 and 84 and applying expansion work definition, we have the system of equations:

$$\frac{d\tilde{p}}{d\rho} = -\gamma \frac{\tilde{p}}{\tilde{V}} \frac{d\tilde{V}}{d\rho} + \frac{\gamma - 1}{\tilde{V}} \left[ \tilde{Q} \frac{d\mathbf{x}}{d\rho} - \tilde{h}\_c \left( 1 + \beta \tilde{V} \right) \left( \frac{\tilde{p}\tilde{V}}{\tilde{m} - \tilde{T}\_{\text{wall}}} \right) \right] - \frac{\gamma C \tilde{p}}{\alpha} \tag{85}$$

$$\frac{d\tilde{\mathbf{Q}}\_{\rm losst}}{d\rho} = \tilde{h}\_c \left( 1 + \beta \tilde{V} \right) \left( \frac{\tilde{p}\tilde{V}}{\tilde{m} - \tilde{T}\_{\rm wall}} \right) \tag{86}$$

$$\frac{dm}{d\rho} = -\frac{Cm}{\rho} \tag{87}$$

$$\frac{d\tilde{W}}{d\rho} = \tilde{p}\frac{d\tilde{V}}{d\rho} \tag{88}$$

Solving this ordinary differential equations system one can calculate the mass lost by leaks, lost heat and power, and pressure variation.

In-cylinder mass varies with time and can be calculated depending on the speed of the engine. Replacing Ec. 73 in Eq. 87 and solving for *m* we obtain:

$$\overline{m} = \exp\left[\frac{-\mathcal{C}\_{bb}(\wp + \pi)}{\frac{\pi \text{(rpm)}}{\Im 0}}\right] \tag{89}$$

#### **3.8. Heat transfer models**

Calculation of instant heat transfer in RICE is a complex problem. Expressions for calculating the total flow heat by combined conduction, convection and radiation requires to use empirical correlations [41]. One of the more used is based on the relationship between the Nusselt number and the Reynolds number for forced convection:

$$\mathrm{Nu} = a \,\mathrm{Re}^b\tag{90}$$

Replacing the Nusselt and Reynolds numbers in Eq. 112, we have:

$$\frac{\overline{h}\_c L}{k} = a \left( \frac{\rho V L}{\mu} \right)^b \tag{91}$$

In this expression L represents the characteristic length, which in RICE is the piston diameter, a magnitude that does not vary.

Considering also separately the radiation heat transfer, gets the expression to evaluate the total heat flow [41]:

$$\frac{q}{A} = a \frac{k}{D\_p} \text{Re}^b \Delta T + c \left( T\_{\text{gas}}^4 - T\_{\text{wall}}^4 \right) \tag{92}$$

In these expressions a, b and c are constants that usually take the following values: a = 0.35 – 0.80, b = 0.70 – 0.80 and c = 0.576 σ, where σ is the Stefan – Boltzmann constant.

An expression used to calculate convection heat transferred is based on the model proposed in [18]:

$$Q\_{\rm lost} = \overline{h}\_c \, \, ^\ast A \, \, \ast \left( \overline{T}\_{\rm gas} - \overline{T}\_{\rm wall} \right) \tag{93}$$

where *<sup>c</sup> h* represents the overall convection heat transfer coefficient that can be determined from the expression:

$$
\overline{H}\_c = a^\* \, D\_p^{m-1} p^m T\_{\text{gas}}^{0.75-1.62m} w^m \tag{94}
$$

The constant values in this equation are: a = 0.13 and m = 0.8. Pressure must be in bar and the temperature in kelvin.

w is the average gas speed calculated by the following expression:

$$
\Delta w = \left[ \mathbf{C}\_{1w} V\_{mp} + \mathbf{C}\_{2w} \frac{V\_d T\_{ref}}{p\_{ref} V\_{ref}} (p - p\_{no\,\,\rm{comb}}) \right] \tag{95}
$$

Subscript ref is used for a reference state that may be the compression process beginning. The empirical constant values depend on the process are: C1w = 6.18 during gas exchange process, C1w = 2.28 during compression, combustion and expansion processes, C2w = 0.0 during gas exchange and compression processes and C2w = 3.24E-3 during combustion and expansion processes.

Another expression used to consider the convection heat transfer process of RICE proposed in [19] is:

$$\overline{h}\_c = \mathbb{C}\_{1H} V\_d^{-0.06} p^{0.8} T^{-0.4} (V\_{mp} + \mathbb{C}\_{2H})^{0.8} \tag{96}$$

where C1H and C2H are empirical constants which take into account local variations due to intake turbulence or combustion chamber geometry of the, their values are: C1H = 130E-4 and C2H = 1.40. Eq. 96 considers the increase in gas speed due to engine velocity and uses as characteristic length the volume rather than the piston diameter, as it is the case when using expressions proposed in [18].

#### **3.9. In-cylinder mass**

94 Internal Combustion Engines

total heat flow [41]:

from the expression:

the temperature in kelvin.

expansion processes.

in [18]:

use empirical correlations [41]. One of the more used is based on the relationship between

Re*<sup>b</sup> Nu a* (90)

*<sup>c</sup> h L VL a*

In this expression L represents the characteristic length, which in RICE is the piston

Considering also separately the radiation heat transfer, gets the expression to evaluate the

*<sup>q</sup> <sup>k</sup> a T cT T*

In these expressions a, b and c are constants that usually take the following values: a = 0.35 –

An expression used to calculate convection heat transferred is based on the model proposed

where *<sup>c</sup> h* represents the overall convection heat transfer coefficient that can be determined

The constant values in this equation are: a = 0.13 and m = 0.8. Pressure must be in bar and

1 0.75 1.62 \* *m m mm*

1 2 ( ) *d ref w mp w nocomb ref ref*

*V T*

Subscript ref is used for a reference state that may be the compression process beginning. The empirical constant values depend on the process are: C1w = 6.18 during gas exchange process, C1w = 2.28 during compression, combustion and expansion processes, C2w = 0.0 during gas exchange and compression processes and C2w = 3.24E-3 during combustion and

 

*w CV C pp p V*

4 4 Re*<sup>b</sup>*

*gas wall*

(92)

*Q hAT T lost c gas wall* \* \* (93)

*c p gas h a D pT w* (94)

(95)

*k*

*p*

0.80, b = 0.70 – 0.80 and c = 0.576 σ, where σ is the Stefan – Boltzmann constant.

*A D*

w is the average gas speed calculated by the following expression:

*b*

(91)

the Nusselt number and the Reynolds number for forced convection:

Replacing the Nusselt and Reynolds numbers in Eq. 112, we have:

diameter, a magnitude that does not vary.

Total mass that fills the cylinder is composed of air, fuel and residual gases. Fuel enters the cylinder when compression process end and fuel injection begins. Residual gases have two components: gases remaining in the cylinder at the end of the exhaust process and recirculated gases entering the cylinder with admitted air as a pollution control measure. There is also a mass lost by blow by. Therefore we have:

$$m\_{\text{total}} = m\_{\text{air}} + m\_{\text{gr}} + m\_f + m\_{\text{EGR}} - m\_{bb} \tag{97}$$

The theoretical incoming air mass is calculated by:

$$m\_{aire} = \rho V\_d \tag{98}$$

While the real mass entering the cylinder is obtained by the expression:

$$m\_{aire} = \sum \stackrel{\bullet}{m}\_{adm} t \tag{99}$$

*madm* is calculated with Eqs. 49 or 50 and time t is obtained from Eq. 46.

#### **3.10. Combustion products model**

This model allows determining theoretically the composition at equilibrium and combustion product thermodynamic properties of a fuel-air mixture whereas reactive products consist of ten chemical species. The combustion global equation for a ten chemical species system is:

$$\begin{aligned} m\_c \mathcal{C}\_m H\_n O\_l + \frac{n + \frac{m}{4} + \frac{l}{2}}{\phi} \mathcal{O}\_2 + 3.7 N\_2 \Big) \to \\ \rightarrow y\_1 H\_2 O + y\_2 H\_2 + y\_3 O H + y\_4 H + y\_5 N\_2 + y\_6 NO + y\_7 CO\_2 + y\_8 CO + y\_9 O\_2 + y\_{10} O \end{aligned} \tag{100}$$

Applying conditions of mass conservation to elements C, H, O and N and reactants we obtain five equations with eleven unknowns (product molar fractions and fuel mass) as shown in following expressions:

$$\mathbb{C} \text{ : } \ y\_1 + y\_8 = n \tag{101}$$

$$H \colon \ 2y\_1 + 2y\_2 + y\_3 + y\_4 = m \tag{102}$$

$$1O: \ y\_1 + y\_3 + y\_6 + 2y\_7 + y\_8 + 2y\_9 + y\_{10} = l + 2\frac{\left(n + \frac{m\zeta}{4} - \frac{1}{2}\zeta\right)}{\phi} \tag{103}$$

$$N\_2 \colon \ 2y\_5 + y\_6 = 7.428 \frac{\left(n + \frac{m}{4}\sqrt{\frac{1}{4}} - \frac{1}{2}\right)}{\phi} \tag{104}$$

$$\sum\_{i=1}^{10} y\_i = 1\tag{105}$$

Applying the chemical equilibrium to combustion reaction yields six additional algebraic equations to close the system:

$$H\_2 \leftrightarrow 2H \qquad\qquad K\_1 = \frac{y\_4^2}{y\_2} \left(\frac{p}{p\_0}\right) \tag{106}$$

$$CO\_2 \leftrightarrow 2O \qquad \qquad K\_2 = \frac{y\_4^2}{y\_9} \left(\frac{p}{p\_0}\right) \tag{107}$$

$$\rm{H}\_{2} + \rm{O}\_{2} \Leftrightarrow 2\rm{OH} \tag{108}$$

$$\rm{K}\_{3} = \frac{y\_{3}^{2}}{y\_{2}y\_{9}} \tag{108}$$

$$\frac{1}{2}O\_2 + \frac{1}{2}N\_2 \leftrightarrow 2OH \qquad \qquad \qquad K\_4 = \frac{y\_6^2}{\frac{1}{2}\frac{1}{y\_5^2}}\tag{109}$$

$$H\_2 + \frac{1}{2}O\_2 \leftrightarrow H\_2O \qquad\qquad\qquad K\_5 = \frac{y\_1}{\frac{1}{y\_2}\left(\frac{p}{p\_0}\right)^{-\frac{1}{2}}} \tag{110}$$

$$\text{CO} + \frac{1}{2}\text{O}\_2 \leftrightarrow \text{CO}\_2 \qquad\qquad \qquad \text{K}\_6 = \frac{y\_7}{\frac{1}{y\_8 y\_8^2}} \left(\frac{p}{p\_0}\right)^{-\frac{1}{2}}\tag{111}$$

In Eqs. 101 to 111 p is the products pressure, y1 to y10 are the species molar fractions and K1 to K6 are the equilibrium constants which are function of temperature. The system is solved by calculating the partial derivatives with respect to temperature, pressure, and mixture richness through FARG, ECP y PER routines. Applying the first law of thermodynamics to reactants and products depending on the type of combustion process we can find their thermodynamic properties. For constant volume combustion and constant pressure combustion respectively:

$$
\mathcal{U}\_R = \mathcal{U}\_P \tag{112}
$$

$$H\_R = H\_P \tag{113}$$

#### **3.11. NOx formation model**

96 Internal Combustion Engines

equations to close the system:

1 8 *Cy y n* : (101)

(104)

(105)

(106)

(107)

(108)

(110)

(111)

1

1 234 *H y yyym* :2 2 (102)

*n m*

(103)

1

2 4

2 0

9 0

 

*y p* 

> 2 3

2 9

2 6

(109)

1 2

1 2

2 2 9 5

1

7

*y y*

*y y*

<sup>0</sup> <sup>2</sup> 2 9

<sup>0</sup> <sup>2</sup> 8 8

*p*

 

*p*

 

*y y*

*y y*

 

*y p* 

> 2 4

*n m*

1 3 6 7 8 9 10

*Oy y y y y y y l*

2 56

*N yy*

4 2 : 22 2

4 2 : 2 7.428

10

1 1 *<sup>i</sup>*

Applying the chemical equilibrium to combustion reaction yields six additional algebraic

<sup>2</sup> *<sup>y</sup> <sup>p</sup> HH K*

<sup>2</sup> *<sup>y</sup> <sup>p</sup> OO K*

*i y* 

2 1

2 2

2 2 4 1 1

2 22 5 1

*<sup>y</sup> <sup>p</sup> H O HO K*

*<sup>y</sup> <sup>p</sup> CO O CO K*

22 6 1

In Eqs. 101 to 111 p is the products pressure, y1 to y10 are the species molar fractions and K1 to K6 are the equilibrium constants which are function of temperature. The system is solved by calculating the partial derivatives with respect to temperature, pressure, and mixture richness through FARG, ECP y PER routines. Applying the first law of thermodynamics to

*<sup>y</sup> O N OH K*

2 2 3

1 1 <sup>2</sup> 2 2

> 1 2

> > 1 2

<sup>2</sup> *<sup>y</sup> H O OH <sup>K</sup>*

The nitrous oxides (NOx) produced by a RICE are NO and NO2 and can be calculated by Zeldovich mechanism. This theory postulates that the production of oxides during the combustion process can be explained through the following chemical reactions:

$$R\_1 \colon \ O + N\_2 \Leftrightarrow NO + N \tag{114}$$

$$R\_2 \colon N + O\_2 \Leftrightarrow NO + O \tag{115}$$

$$R\_3 \colon N + \mathcal{O}H \Leftrightarrow \mathcal{NO} + H \tag{116}$$

A differential equation which allows finding the NO concentration by unit time during combustion, using Zeldovich reactions and considering the possibility of occurrence in both directions can be obtained using the basic theory of chemical kinetics:

$$\frac{d\left[NO\right]}{dt} = k\_1^+ \left[O\right] \left[N\_2\right] + k\_2^+ \left[N\right] \left[O\_2\right] + k\_3^+ \left[N\right] \left[OH\right] - k\_1^- \left[NO\right] \left[N\right] - k\_2^- \left[NO\right] \left[O\right] - k\_3^- \left[N\right] \left[OH\right] \text{ (117)}$$

Considering the steady state condition for [N], which is equivalent to assume its change is small compared to other interesting species and proceeding in similar mode to the previous case, we obtained:

$$\frac{d\left[\left[N\right]{}\right]}{dt} = k\_1^+ \left[O\left[\left[N\_2\right]\right] - k\_2^+ \left[N\right]\left[O\_2\right] - k\_3^+ \left[N\right]\left[OH\right] - k\_1^- \left[NO\right]\left[N\right] + k\_2^- \left[NO\right]\left[O\right] + k\_3^- \left[NO\right]\left[H\right] \tag{118}$$

$$\frac{d\left[N\right]}{dt} = 0 \tag{119}$$

Substituting Ec. 118 in Ec.117 and using the equilibrium reactions (Ecs. 114 to 116) gets the following expression:

*dt*

$$\frac{d\left[NO\right]}{dt} = \frac{2R\_1\left\{1 - \left(\left[NO\right]\right\|\left[NO\right]\_\epsilon\right)^2\right\}}{1 + \frac{\left(\left[NO\right]\right\|\left[NO\right]\_\epsilon\right)R\_1}{\left(R\_2 + R\_3\right)}}\tag{120}$$

and Eqs. 114 to 116 are transformed into:

$$R\_1 \colon \begin{array}{c} k\_1^+ \end{array} \Big[ \bigcirc \big]\_\varepsilon \Big[ N\_2 \big]\_\varepsilon \Big] \tag{121}$$

$$R\_2: \ k\_2^- \left[ \text{NO} \right]\_e \left[ O \right]\_e \tag{122}$$

$$R\_3: \ k\_3^- \left[NO\right] \Big\|\ f\big\|\tag{123}$$

Terms enclosed in square brackets represent the chemical equilibrium concentration for corresponding species obtained using routines FARG and ECP. Constants for Zeldovich reactions, depending on temperature, are expressed in the following way.

$$K\_1^+ = 7.6E13 \exp\left(-38000/T\right) \tag{124}$$

$$K\_2^- = 1.5E09 \exp\left(-19500/T\right) \tag{125}$$

$$K\_3^- = 2.0E14 \exp\left(-23650/T\right) \tag{126}$$

Positive sign indicates reaction tendency to form products, while negative sign indicates tendency to form reactants.

#### **3.12. Exhaust gas recirculation (EGR) model**

Exhaust gas recirculation (EGR) is a very effective technique employed to reduce nitrogen oxide emissions (Lapuerta, 2000). The method involves replacing a part of the intake air with exhaust gases during the admissions process. In this way a decrease in the amount of air available for combustion reduces the final temperature of combustion, and therefore lowers the production of nitrous oxides emissions. This process can be expressed in the following way:

$$\text{Fuel} + \text{Air} + \text{ }\%\text{EGR} \longrightarrow \text{Products}$$

Recirculated gases will be taken at engine exit and will be introduced into the cylinder at a temperature greater than air admitted. Composition and properties calculation of involved species will be made with routines FARG and ECP.

#### **3.13. Work, mean effective pressure and efficiency**

To determine the work produced by closed loop cycle, we consider that it is represented by the loop enclosed area in a p-V diagram. The area was calculated with the trapezoidal rule [42]:

$$A = W = \left\{ p\_0 \Delta V + p\_1 \Delta V + p\_2 \Delta V + p\_3 \Delta V \dots \dots p\_{n-1} \Delta V + p\_n \Delta V \right\} \tag{127}$$

Gas exchange work, indicated work and net work are obtained with the following equations:

$$\mathcal{W}\_{pump} = \mathcal{W}\_{exh} - \mathcal{W}\_{adm} \tag{128}$$

$$\mathcal{W}\_i = \mathcal{W}\_{\text{exp}} - \mathcal{W}\_{\text{com}} \tag{129}$$

$$\mathcal{W}\_{net} = \mathcal{W}\_i - \mathcal{W}\_{pump} \tag{130}$$

Net power is obtained through the following expression:

$$\stackrel{\bullet}{W}\_{net} = W\_{net} \frac{rpm}{30j} \tag{131}$$

Mean indicated pressure, mean net indicated pressure and mean indicated gas exchange pressure are calculated with equations:

$$\text{mip} = \mathcal{W}\_i \not\!\!\!/ V\_d \tag{132}$$

$$\text{mip}\_{net} = \text{W}\_{net} \not\!\!/ \text{V}\_d \tag{133}$$

$$\text{mip}\_{pum} = \text{W}\_{pum} \not\!\!/ V\_d \tag{134}$$

Efficiencies are determined with:

98 Internal Combustion Engines

tendency to form reactants.

following way:

[42]:

**3.12. Exhaust gas recirculation (EGR) model** 

species will be made with routines FARG and ECP.

**3.13. Work, mean effective pressure and efficiency** 

and Eqs. 114 to 116 are transformed into:

11 2 : *e e R kO N* (121)

2 2 : *e e R k NO O* (122)

3 3 : *e e R k NO H* (123)

<sup>1</sup> *KE T* 7.6 13exp 38000 (124)

<sup>2</sup> *KE T* 1.5 09exp 19500 (125)

<sup>3</sup> *KE T* 2.0 14exp 23650 (126)

Terms enclosed in square brackets represent the chemical equilibrium concentration for corresponding species obtained using routines FARG and ECP. Constants for Zeldovich

Positive sign indicates reaction tendency to form products, while negative sign indicates

Exhaust gas recirculation (EGR) is a very effective technique employed to reduce nitrogen oxide emissions (Lapuerta, 2000). The method involves replacing a part of the intake air with exhaust gases during the admissions process. In this way a decrease in the amount of air available for combustion reduces the final temperature of combustion, and therefore lowers the production of nitrous oxides emissions. This process can be expressed in the

*Fuel Air EGR Products* % 

Recirculated gases will be taken at engine exit and will be introduced into the cylinder at a temperature greater than air admitted. Composition and properties calculation of involved

To determine the work produced by closed loop cycle, we consider that it is represented by the loop enclosed area in a p-V diagram. The area was calculated with the trapezoidal rule

*A W pV pV pV pV p V pV* 0123 1 ........ *n n* (127)

reactions, depending on temperature, are expressed in the following way.

$$\eta\_i = \frac{W\_i}{m\_c \ H\_i} \tag{135}$$

$$
\eta\_{int} = \frac{\mathcal{W}\_{net}}{m\_c \ H\_i} \tag{136}
$$

#### **3.14. MECID computer program**

MECID program consist of a routines group for determining: optimum lifting of intake and exhaust valves, gas leakage across piston rings (blow by), composition variation of chemical species that constitute the working fluid, effect of heat transfer, ignition delay, heat release, residual gases and recirculated gases, as well as nitrous oxides formation and pressure drop during intake and exhaust processes due to valve contraction. To calculate all parameters it uses the mathematical relationships and models proposed in preceding paragraphs. Figure 2 shows the MECID program flow diagram.

Using the ENGINE routine the program requested information on engine works site conditions (temperature and pressure) and on engine geometric characteristics (piston bore and stroke; valves diameter, maximum lift and average discharge coefficient, compression ratio, connecting rod length to crank radius ratio and engine speed) With VALVE routine requested information on calculation conditions (calculation angular interval, advance and delay of intake and exhaust valve closure, richness, residual gases ratio, recirculated gas percentage and cylinder wall average temperature). Selected heat transfer, ignition delay and heat release models which will be used in the application, as well as the angles when combustion begins and ends, the program allows to perform studies to determine what should be the working conditions to ensure maximum fuel energy utilization the influence of different variables on engine operation. Results of more important studies will be presented in next paragraph.

**Figure 2.** Block diagram and routines used by MECID computer program

## **4. Results and discussion**

## **4.1. Effect of burning process beginning on CIE operating parameters**

For combustion duration constant and varying combustion start we obtained the p and T vs. φ diagrams shown in Figure 3. We note that if the combustion process is advanced too much, from 360° to 310°, it reached a high maximum pressure while the piston is at some point in the compression stroke. This brings as a consequence, greater work consumption and reduced cycle efficiency. If combustion beginning delays too much, 360° to 410°, the maximum pressure is very low. If combustion beginning advance is such that the process occurs around TDC we obtain intermediate maximum pressures during expansion process and a higher maximum temperature, representing greater thermal energy available to do work. In view of the above, for subsequent studies we will take a 345° angle as combustion beginning angle. Maximum temperature is obtained for a greater angle than the maximum pressure one because the burning model uses constant combustion duration.

**Figure 3.** p vs φ and T vs φ diagrams for various combustion beginning angles

100 Internal Combustion Engines

connecting rod length to crank radius ratio and engine speed) With VALVE routine requested information on calculation conditions (calculation angular interval, advance and delay of intake and exhaust valve closure, richness, residual gases ratio, recirculated gas percentage and cylinder wall average temperature). Selected heat transfer, ignition delay and heat release models which will be used in the application, as well as the angles when combustion begins and ends, the program allows to perform studies to determine what should be the working conditions to ensure maximum fuel energy utilization the influence of different variables on engine operation. Results of more important studies will be presented in next paragraph.

**Figure 2.** Block diagram and routines used by MECID computer program

**4.1. Effect of burning process beginning on CIE operating parameters** 

For combustion duration constant and varying combustion start we obtained the p and T vs. φ diagrams shown in Figure 3. We note that if the combustion process is advanced too much, from 360° to 310°, it reached a high maximum pressure while the piston is at some

**4. Results and discussion** 

Figure 4 shows how varies the accumulated burned mass fraction and the instant burned mass fraction depending on combustion start angle. Burned mass is responsible for pressure and temperature rapid increase during combustion process. It is observed that peak values are the same for all three cases considered. This is because richness and volumetric efficiency remain unchanged.

**Figure 4.** Accumulated and instant burned mass fraction for various combustion beginning angles

#### **4.2. Effect of combustion process duration on CIE operating parameters**

In present case, the start of combustion process is kept constant. Figure 4 shows that at lower combustion duration the greater the maximum pressure reached since a great amount of heat is released in a short time period, as seen by the burned curves slope. In these curves can be noted that in considered cases, the maximum burned fractions do not change because the richness was kept constant and the obtained volumetric efficiency remained almost constant, therefore the total in-cylinder mixture mass is the same for all referred cases.

**Figure 5.** Pressure and accumulated burned fraction for different combustion duration

Table 4 presents the principal results obtained for operation parameters. From Table 4 it can be concluded that the more efficient cycle producing higher power, is one in which combustion takes 45° and therefore this duration will be used in future studies.


**Table 4.** Summary of main results obtained when varying combustion process duration.

### **4.3. Effect of angular velocity on CIE operating parameters**

Pressure variation depending on crankshaft rotation angle when engine speed change, is shown in Figure 6. It can be noted from this figure that higher pressures are reached when the engine operates at lower speeds. This is because while the engine works at lower speeds, greater will be its volumetric efficiency. If the richness is kept constant, the greater the mass which is burned during combustion process and thus more energy will be released.

**Figure 6.** In-cylinder pressure for various engine speeds

102 Internal Combustion Engines

of heat is released in a short time period, as seen by the burned curves slope. In these curves can be noted that in considered cases, the maximum burned fractions do not change because the richness was kept constant and the obtained volumetric efficiency remained almost constant, therefore the total in-cylinder mixture mass is the same for all referred cases.

**Figure 5.** Pressure and accumulated burned fraction for different combustion duration

combustion takes 45° and therefore this duration will be used in future studies.

Combustion

•

Table 4 presents the principal results obtained for operation parameters. From Table 4 it can be concluded that the more efficient cycle producing higher power, is one in which

process duration 30º 45º 60º

ηv 0.72 0.72 0.72

Wi [kJ] 0.60 0.73 0.75

Wnet [kJ] 0.57 0.65 0.64

Wnet [kW] 9.56 10.98 10.72

mip [kPa] 994.75 1200.21 1232.88

ηi net 0.35 0.40 0.39

Pmax [kPa] 7689.17 7134.93 6614.27

Tmax [K] 2235.64 2400.71 2496.92

Pressure variation depending on crankshaft rotation angle when engine speed change, is shown in Figure 6. It can be noted from this figure that higher pressures are reached when

**Table 4.** Summary of main results obtained when varying combustion process duration.

**4.3. Effect of angular velocity on CIE operating parameters** 


**Table 5.** Summary of main results obtained varying engine speed.

Table 5 presents the principal results obtained for operating parameters. From this table it can be concluded that at higher rpm volumetric efficiency is reduced due to increased air velocity in the intake system and increased pressure friction losses. Additionally, as there is less time to fill the cylinder, which results in lower air intake, the maximum pressure reached is reduced as well as the indicated work. Net power increases due to the influence of increasing rpm. Indicated efficiency tends to increase with increasing rpm, because burning the same fuel amount at higher rpm and net power is higher. Figure 6 shows that maximum pressure is obtained in all cases close to the TDC. One can also appreciate that the higher the rpm, the mass fraction burned during the premixed phase and the mean indicated pressure are reduced

## **4.4. Effect of compression ratio on CIE operating parameters**

Now we will consider the effect of varying the compression ratio on the main engine operating parameters assuming valve inlet pressure equal to atmospheric pressure. Figure 7 shows the in-cylinder pressure variation and the accumulated burned mass fraction for various compression ratios and in Table 6 we can see the summary of main results obtained for different compression ratios.

**Figure 7.** Pressure and accumulated burned fraction for different compression ratios

By analyzing Figure 7 and Table 6 we can conclude that: i) volumetric efficiency is unaffected by compression ratio since is not dependent on it, ii) at higher compression ratios there is a reduction of burned mass fraction in the premixed phase, iii) there is an indicated work increase because maximum pressure and temperature increase, resulting in increased efficiency and mean pressures, iv) pumping work increases when compression ratio increases as pressures are greatest during exhaust, v) cycle net work increases because net work and power increase at higher rates than pumping work.


**Table 6.** Summary of main results obtained varying compression ratio.

### **4.5. Effect of atmospheric pressure on CIE operating parameters**

It is well known that increasing the engine speed the start of combustion must anticipate in order that maximum pressures remain nearly constant and occur near TDC. For this reason in present study the combustion duration will remain constant and, its beginning will be advanced. Figure 8 shows the indicator diagram for various engine speeds. In Figure 8 which shows the indicator diagram for various engine speeds, it is observed that the maximum pressure in all studied cases occurs near TDC and its value is approximately constant. Since optimum injection angle for each speed is not known it is not easy to superimpose the curves which is what is desired in an engine in order to obtain maximum power and efficiency in a wide speeds range.

**Figure 8.** Pressure and accumulated burned fraction for various atmospheric pressures


**Table 7.** Summary of main results obtained varying engine speed.

The summary of results in Table 7 shows that increasing the engine speed i) volumetric efficiency is reduced which is the expected behavior for CIE, ii) mass fraction burned during the premixed phase is reduced and iii) mean pressure decreases possibly due to a combustion beginning advance not enough for the increase in speed considered.

## **5. Conclusions**

104 Internal Combustion Engines

•

for different compression ratios.

**4.4. Effect of compression ratio on CIE operating parameters** 

**Figure 7.** Pressure and accumulated burned fraction for different compression ratios

work and power increase at higher rates than pumping work.

**Table 6.** Summary of main results obtained varying compression ratio.

**4.5. Effect of atmospheric pressure on CIE operating parameters** 

By analyzing Figure 7 and Table 6 we can conclude that: i) volumetric efficiency is unaffected by compression ratio since is not dependent on it, ii) at higher compression ratios there is a reduction of burned mass fraction in the premixed phase, iii) there is an indicated work increase because maximum pressure and temperature increase, resulting in increased efficiency and mean pressures, iv) pumping work increases when compression ratio increases as pressures are greatest during exhaust, v) cycle net work increases because net

Compression ratio 18 20 22 ηv 0.72 0.72 0.72 β 45.33 42.68 40.24 Wi [kJ] 0.73 0.75 0.77

Wnet [kW] 10.99 11.36 11.68 mip [kPa] 1200.22 1233.52 1258.36 η<sup>i</sup> 0.45 0.46 0.47 Pmax [kPa] 7134.93 7791.11 8349.81 Tmax [K] 2400.71 2399.34 2392.72

It is well known that increasing the engine speed the start of combustion must anticipate in order that maximum pressures remain nearly constant and occur near TDC. For this reason

Now we will consider the effect of varying the compression ratio on the main engine operating parameters assuming valve inlet pressure equal to atmospheric pressure. Figure 7 shows the in-cylinder pressure variation and the accumulated burned mass fraction for various compression ratios and in Table 6 we can see the summary of main results obtained

> A methodology to theoretically determine the thermodynamic working cycle of direct injection compression ignition engines was presented.

A computer program based on the application of the above methodology was introduced.

Results of several studies performed with the program was presented and discussed. From some results we can conclude:


## **Notation**



### **Greek letters**

106 Internal Combustion Engines

**Notation** 

A area

c piston stroke

Dp piston diameter dv valve diameter

e specific energy f fuel air ratio h specific enthalpy

Hi lower heating value ID ignition delay j engine strokes

�� adimensional mass mip mean indicated pressure

�� adimensional pressure pdown downstream pressure pup upstream pressure

�� adimensional heat

q heat per unit mass

Nu Nusselt number

E energy

H enthalpy

Lv valve lift m mass �� mass flow

p pressure

Q heat

�� heat flow

Cd discharge coefficient

Cp constant pressure specific heat Cv constant volume specific heat d( ) preceding term ordinary derivative

hc convection heat transfer coefficient ℎ� mean convection heat transfer coefficient

k conduction heat transfer coefficient

some results we can conclude:

Results of several studies performed with the program was presented and discussed. From



coincides with the theoretical internal combustion engines behavior.

ensures a greater air intake and thus an increase in produced power.


## **Subindex**



## **Author details**

Simón Fygueroa

*National University of Colombia, Colombia, Turin Polytechnic Institute, Italy, Polytechnic University of Valencia, Spain, Mechanical Engineering Department, University of the Andes, Mérida, Venezuela, Pamplona University, Mechanical Engineering Department, Pamplona, Colombia* 

Carlos Villamar *University of the Andes, Mérida, Venezuela,* 

Olga Fygueroa *University of the Andes, Mérida, Venezuela, Polytechnic University of Valencia, Spain Altran-Nissan Technical Center, Barcelona, Spain* 

### **6. References**


[3] Rajput RK. Internal combustion engines. New Delhi: Laxmi Publications Ltd; 2007.

108 Internal Combustion Engines

e input

exh exhaust exp expansion f fuel

P products pre premixed phase pum pumping R reactants

s output sto stoichiometric

Tot total v volumetric vc control volume wall cylinder wall

**Author details** 

Simón Fygueroa

Carlos Villamar

Olga Fygueroa

**6. References** 

dif diffusive phase

gr residual gasses i indicated lost lost max maximum net net

no comb without combustion

ref reference condition

*National University of Colombia, Colombia,* 

*Polytechnic University of Valencia, Spain,* 

*University of the Andes, Mérida, Venezuela,* 

*University of the Andes, Mérida, Venezuela, Polytechnic University of Valencia, Spain* 

*Altran-Nissan Technical Center, Barcelona, Spain* 

*Mechanical Engineering Department, University of the Andes, Mérida, Venezuela, Pamplona University, Mechanical Engineering Department, Pamplona, Colombia* 

[1] Jovaj MS. Automotive vehicle engines. Moscu: Editorial MIR; 1982.

[2] Payri F. and Desantes JM. Motores de combustión interna alternativos. Valencia, (Spain): Servicio de Publicaciones Universidad Politécnica de Valencia; 2011.

*Turin Polytechnic Institute, Italy,* 

EGR exhaust gas recirculated

