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

1. Introduction

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 compression-expansion 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 pre-mixed phase and the diffusive phase.

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

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:

$$m_{vc}={\displaystyle \sum m_e}-{\displaystyle \sum m_s}$$E1

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{d{m}_{vc}}{dt}=\frac{d{m}_e}{dt}-\frac{d{m}_s}{dt}$$E2

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

$f=\frac{m_f}{m}$E3

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

$$\stackrel{•}{f}=\frac{df}{dt}=\left(\frac{{\stackrel{•}{m}}_e-{\stackrel{•}{m}}_s}{m}\right)\left[\left(f_e-f_s\right)\right]$$E4

From equivalence ratio (mixture richness) definition:

$\varphi =\frac{\frac{m{}_f}{m_a}}{{\left(\frac{m{}_f}{m_a}\right)}_{sto}}=\frac{\frac{\stackrel{•}{m}{}_f}{{\stackrel{•}{m}}_a}}{{\left(\frac{\stackrel{•}{m}{}_f}{{\stackrel{•}{m}}_a}\right)}_{sto}}$E5

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

$\stackrel{•}{\varphi }=\frac{d\varphi }{dt}=\frac{1}{{\left(\frac{m_f}{m_a}\right)}_{sto}}\frac{\stackrel{•}{f}}{{\left(1-f\right)}^2}$E6

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{dQ}{dt}-\frac{dW}{dt}+\stackrel{•}{m_e}{h}_e-\stackrel{•}{m_s}{h}_s$E7

Since the work due to a volume change is:

$\frac{dW}{dt}=\stackrel{•}{W}=P\frac{dV}{dt}$E8

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

$\frac{dE}{dt}=\frac{d}{dt}\left(mu\right)={\left(m\frac{du}{dt}\right)}_{vc}+{\left(u\frac{dm}{dt}\right)}_{vc}$E9

or in terms of enthalpy:

$\frac{dE}{dt}=\frac{d}{dt}\left(mh\right)-\frac{d}{dt}(pV)$E10

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

${\left(m\frac{du}{dt}\right)}_{vc}+{\left(u\frac{dm}{dt}\right)}_{vc}=\stackrel{•}{Q}-\stackrel{•}{W}+{\stackrel{•}{m}}_{e}h{}_e-{\stackrel{•}{m}}_{s}h{}_s$E11

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 \varphi }\right)\frac{d\varphi }{dt}$E12

$\frac{dh}{dt}=\left(\frac{\partial h}{\partial T}\right)\frac{dT}{dt}+\left(\frac{\partial h}{\partial p}\right)\frac{dp}{dt}+\left(\frac{\partial h}{\partial \varphi }\right)\frac{d\varphi }{dt}$E13

$\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 \varphi }\right)\frac{d\varphi }{dt}$E14

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}$E15

$\frac{d\rho }{dt}=\frac{\frac{dp}{dt}-\rho R\frac{dT}{dt}-\rho T\frac{dR}{dt}}{RT}$E16

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 \varphi }\right)\frac{d\varphi }{dt}}{\left(\frac{\partial \rho }{\partial p}\right)}$$E17

substituting Eq. 17 in Eq. 16, rearranging and solving for $\frac{dp}{dT}$:

$\frac{dp}{dt}=\frac{-\frac{\rho }{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 \varphi }\right)\frac{d\varphi }{dt}}{\left(\frac{\partial \rho }{\partial p}\right)-\frac{1}{RT}}$E18

Solving Ec. 15 for $\frac{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 \varphi }\right)\frac{d\varphi }{dt}$$E19

Replacing ideal gas equation in Ec. 19 and solving for $\frac{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 \varphi }\right)\frac{d\varphi }{dt}\right]$$E20

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

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

$\left(\frac{\partial \rho }{\partial T}\right)=-\frac{p}{R{T}^2}$E22

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

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

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 $\frac{du}{dt}$:

$\frac{du}{dt}=\frac{\stackrel{•}{Q}}{m_{vc}}-\frac{p}{m_{vc}}\frac{dV}{dt}+\frac{1}{m_{vc}}\left(\stackrel{•}{m_e}{h}_e-\stackrel{•}{m_s}{h}_s-{\left(u\frac{dm}{dt}\right)}_{vc}\right)$E24

and defining:

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

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

$\frac{du}{dt}=B$E26

Substituting Ec. 26 in Ec. 12 and solving for $\frac{dp}{dt}$:

$\frac{dp}{dt}=\frac{B-\left(\frac{\partial u}{dt}\right)\frac{dT}{dt}-\left(\frac{\partial u}{d\varphi }\right)\frac{d\varphi }{dt}}{\left(\frac{\partial u}{\partial p}\right)}$E27

Replacing Ec. 27 in Ec. 20 and solving for $\frac{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 \varphi }\right)\frac{d\varphi }{dt}-B+\left(\frac{\partial u}{\partial \varphi }\right)\frac{d\varphi }{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)}$$E28

Now, considering:

$R=R(T,p,\varphi ,t)$E29

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 \varphi }\right)\frac{d\varphi }{dt}$E30

Differentiating ideal gas equation and solving for $\frac{dR}{dt}$

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

Replacing Ec. 30 in Eq. 31 and solving for$\frac{dp}{dt}$:

$\frac{dp}{dt}=\frac{\left(\frac{\partial R}{\partial T}\right)\frac{dT}{dt}+\left(\frac{\partial R}{\partial \varphi }\right)\frac{d\varphi }{dt}+\frac{dT}{dt}\frac{R}{T}+\frac{d\rho }{dt}\frac{R}{\rho }}{\frac{R}{p}-\left(\frac{\partial R}{\partial p}\right)}$E32

Replacing Ec. 32 in Ec. 27:

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

Defining:

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

Collecting terms containing $\frac{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 \varphi }\right)\frac{d\varphi }{dt}+\frac{d\varphi }{dt}\frac{R}{\rho }\right]+\left(\frac{\partial u}{\partial \varphi }\right)\frac{d\varphi }{dt}=B$E35

Defining:

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

Replacing Ec. 36 in Eq. 35 and solving for$\frac{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 \varphi }\right)\frac{d\varphi }{dt}+\frac{d\rho }{dt}\frac{1}{\rho }\right]-\left(\frac{\partial u}{\partial \varphi }\right)\frac{d\varphi }{dt}}{\left(\frac{\partial u}{\partial T}\right)+\frac{C}{D}\frac{p}{T}\left(\frac{\partial u}{\partial p}\right)}$E37

Since:

$m=\rho V$E38

Time differentiating and solving for $\frac{d\rho }{dt}$:

$\frac{d\rho }{dt}=\frac{\frac{dm}{dt}-\rho \frac{dV}{dt}}{V}$E39

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 \varphi }\right)\frac{d\varphi }{dt}\frac{1}{R}+\frac{dm}{dt}\frac{1}{m}-\frac{dV}{dt}\frac{1}{V}\right]-\left(\frac{\partial u}{\partial \varphi }\right)\frac{d\varphi }{dt}}{\left(\frac{\partial u}{\partial T}\right)+\frac{C}{D}\frac{p}{T}\left(\frac{\partial u}{\partial p}\right)}$E40

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{•}{Q}-{\stackrel{•}{m}}_{{}_{bb}}{h}_{bb}-\frac{dV}{dt}\left[\frac{mC}{V}+p\right]-\frac{dm}{dt}\left[\left(D\left(\frac{\partial u}{\partial \varphi }\right)-h_{cil}+u\right)-C\left(\frac{D}{R}\left(\frac{\partial R}{\partial \varphi }+1\right)\right)\right]}{m\left[C\left(\frac{1}{p}-\frac{1}{R}\left(\frac{\partial R}{\partial \varphi }\right)\right)+\left(\frac{\partial u}{\partial p}\right)\right]}$$E41

and:

$$\frac{dm}{dt}=\frac{\stackrel{•}{Q}-{\stackrel{•}{m}}_{{}_{bb}}{h}_{bb}-\frac{dV}{dt}\left[\frac{mC}{V}+p\right]-m\left[C\left(\frac{1}{p}-\frac{1}{R}\left(\frac{\partial R}{\partial \varphi }\right)\right)+\left(\frac{\partial u}{\partial p}\right)\right]\frac{dp}{dt}}{\left(D\left(\frac{\partial u}{\partial \varphi }\right)-h_{cil}+u\right)-C\left(\frac{D}{R}\left(\frac{\partial R}{\partial \varphi }+1\right)\right)}$$E42

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(\phi )=V_d\left[\frac{1}{r_c-1}+\frac{1}{2}\left[R_{LA}+1-cos\phi -{(R_{LA}^2-se{n}^2\phi )}^{\frac{1}{2}}\right]\right]$$E43

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

$R_{LA}=\frac{l}{a}$E44

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

$\frac{dV}{d\phi }=\frac{V_d}{2}\left[sen\phi +\frac{sen\phi cos\phi }{{\left(R_{LA}^2-se{n}^2\phi \right)}^{\frac{1}{2}}}\right]$E45

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

$t=\frac{\phi }{6rpm}$E46

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[sen\frac{\pi rpm}{30}t+\frac{sen\frac{\pi rpm}{30}tcos\frac{\pi rpm}{30}t}{{\left(R_{LA}^2-se{n}^2\frac{\pi rpm}{30}t\right)}^{\frac{1}{2}}}\right]$$E47

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{•}{m}}_{ar}}{{\stackrel{•}{m}}_{at}}=\frac{{\stackrel{•}{m}}_{ar}}{{\rho }_{0}i{V}_d\frac{rpm}{30j}}$$E48

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

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

$$\frac{p_{down}}{p_{up}}<1\stackrel{•}{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)$$E49

$$\frac{p_{down}}{p_{up}}\ge 1\stackrel{•}{m}=C_d{A}_{ref}{p}_{up}\sqrt{\left(\frac{\gamma }{R{T}_0}\right){\left(\frac{2}{\gamma -1}\right)}^{\frac{\gamma +1}{2(\gamma +1)}}}$$E50

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

$$A_{ref}=\pi d_v{L}_v$$E51

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

$$L_v(\phi )=L_{v\max }+C_2{\phi }^2+C_p{\phi }^p+C_q{\phi }^q+C_r{\phi }^r+C_s{\phi }^s$$E52

Coefficients C₂, C_p, C_q, C_r y C_s are determined with the following equations:

$$C_2=\frac{-pqrsh}{\left[\left(p-2\right)\left(q-2\right)\left(r-2\right)\left(s-2\right)c_{med}^2\right]}$$E53

$$C_p=\frac{2qrsh}{\left[\left(p-2\right)\left(q-p\right)\left(r-p\right)\left(s-p\right)c_{med}^p\right]}$$E54

$$C_q=\frac{-2prsh}{\left[\left(q-2\right)\left(q-p\right)\left(r-q\right)\left(s-q\right)c_{med}^q\right]}$$E55

$$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]}$$E56

$$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]}$$E57

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{\frac{dV}{dt}}{V}+\frac{\frac{dm}{dt}}{m}+\left(\frac{dT}{dt}\right)\frac{1}{T}\right]$E58

$\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)}$E59

To solve these equations heat release, heat transfer, blow by, ignition delay and chemical species formation models are required. In addition the terms$$\frac{\partial u}{\partial T}$$,$$\frac{\partial u}{\partial p}$$, $$\frac{\partial R}{\partial T}$$ y $$\frac{\partial R}{\partial 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

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:

$$I{D}_{[\mathrm{deg}]}=A{p}_c^n(EA)$$E60

$$I{D}_{\left[ms\right]}=\frac{I{D}_{[\mathrm{deg}]}}{0.006(rpm)}$$E61

where: $$\text{A=0}\text{.36+0}{\text{.22V}}_{\text{mp}},$$$${\text{V}}_{\text{mp}}\text{=}\frac{\text{c(rpm)}}{\text{30}},$$$$\text{n=0},$$ $$\text{EA=exp}\left[{\text{R}}_{\text{u}}\left(\frac{\text{1}}{{\text{RT}}_{\text{c}}}\text{-}\frac{\text{1}}{\text{17190}}\right)\left(\frac{\text{21}\text{.2}}{{\text{p}}_{\text{c}}\text{-12}\text{.4}}\right)\right],$$$$\text{E=}\frac{\text{618840}}{\text{NC+25}},$$$${\text{p}}_{\text{c}}{\text{=p}}_{\text{amb}}{\text{r}}_{\text{c}}{}^{{\text{n}}_{\text{c}}},$$$${\text{T}}_{\text{c}}{\text{=T}}_{\text{amb}}{\text{r}}_{\text{c}}{}^{{\text{n}}_{\text{c}}\text{-1}}\text{,}$$

$${\text{n}}_{\text{c}}\text{=1}\text{.30}\text{a}\text{1}\text{.37},{\text{R}}_{\text{u}}\text{=8}\text{.3143}\left[\text{J/mol}\text{K}\right].$$E62

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 follows: $$\text{A=3}\text{.45},$$$$\text{n=1}\text{.02},$$$$\text{EA=exp}\left[\frac{{\text{E}}_{\text{a}}}{{\text{R}}_{\text{u}}{\text{T}}_{\text{c}}}\right],$$$$\frac{\text{E}}{{\text{R}}_{\text{u}}}\text{=2100}.$$

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]: $$\text{A=2}\text{.4}{\varphi }^{\text{-0}\text{.2}}.$$

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}d{X}_b=\frac{dm}{dt}\left[\left(D\left(\frac{\partial u}{\partial \varphi }\right)-h_{cil}+u\right)-C\left(\frac{D}{R}\left(\frac{\partial R}{\partial \varphi }+1\right)\right)\right]$$E63

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{d{m}_c}{d\phi }\right)}_{Tot}={\left(\frac{d{m}_c}{d\phi }\right)}_{pre}+{\left(\frac{d{m}_c}{d\phi }\right)}_{dif}$$E64

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}}$$E65

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

$$\beta =1-\frac{a_1{\varphi }_0^{b_1}}{I{D}_S^{c{c}_1}}$$E66

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

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

Table 1.

Empirical values for burning factor

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

$${\left(\frac{d{m}_c}{d\phi }\right)}_{Tot}=\beta {\left(\frac{d{m}_c}{d\phi }\right)}_{pre}+\left(1-\beta \right){\left(\frac{d{m}_c}{d\phi }\right)}_{dif}$$E67

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

$${\left(\frac{d{X}_b}{d\phi }\right)}_{pre}=C_1{C}_2{\left(\frac{\phi -{\phi }_0}{\Delta \phi }\right)}^{C_1-1}{\left(1-{\left(\frac{\phi -{\phi }_0}{\Delta \phi }\right)}^{C_1}\right)}^{C_2-1}$$E68

$${\left(\frac{d{X}_b}{d\phi }\right)}_{dif}=C_3{C}_4{\left(\frac{\phi -{\phi }_0}{\Delta \phi }\right)}^{C_4-1}\exp \left(C_3-{\left(\frac{\phi -{\phi }_0}{\Delta \phi }\right)}^{C_4}\right)$$E69

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

$$\frac{dQ}{d\phi }=a\left(\frac{Q_{pre}}{{\phi }_{pre}}\right)m_{pre}{\left(\frac{\phi }{{\phi }_{pre}}\right)}^{m_{pre}-1}\exp \left[-a{\left(\frac{\phi }{{\phi }_{pre}}\right)}^{m_{pre}}\right]+\left(\frac{Q_{dif}}{{\phi }_{dif}}\right)m_{dif}{\left(\frac{\phi }{{\phi }_{dif}}\right)}^{m_{dif}-1}\exp \left[-a{\left(\frac{\phi }{{\phi }_{dif}}\right)}^{m_{dif}}\right]$$E70

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

Constant	Heywood	Watson
Expresión	Valores sugeridos
C1	$$2+1.25E-8{\left(I{D}_{[ms]}rpm\right)}^{2.4}$$	$$2+0.002703*ID{a}^{2.4}$$	3
C2	5000	5000	5000
C3	$$14.2{\varphi }^{-0.644}$$	$$10{\varphi }^{1.505}$$	6.908
C4	$$0.79C_3^{0.25}$$	$$0.48C_3^{0.423}$$	1.4

Table 2.

Constants for Heywood and Watson expressions

Constant	Miyamoto
mpre	4
mdif	1.5 ID o 1.9 ID
a	0.69
φpre	+7^°

Table 3.

Constants for Miyamoto expression

3.7. Blow by model

Due to the in-cylinder pressure increase during compression and combustion processes, a part of the gasses (m_bb) 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\phi }+\frac{1}{V}\frac{dV}{d\phi }=\frac{1}{m}\frac{dm}{d\phi }+\frac{1}{T}\frac{dT}{d\phi }$E71

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

$m{C}_v\frac{dT}{d\phi }+C_{v}T\frac{dm}{d\phi }=\frac{dQ}{d\phi }-p\frac{dV}{d\phi }-\frac{{\stackrel{•}{m}}_{bb}{C}_{p}T}{\omega }$E72

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

$\frac{dp}{d\phi }=-\gamma \frac{p}{V}\frac{dV}{d\phi }+\frac{\gamma -1}{V}\frac{dQ}{d\phi }-\frac{\gamma {\stackrel{•}{m}}_{bb}}{\omega m}p$E73

Blow by mass can be found from mass conservation:

$\frac{dm}{d\phi }=-\frac{{\stackrel{•}{m}}_{bb}}{\omega }$E74

Defining:

$C=\frac{{\stackrel{•}{m}}_{bb}}{m}$E75

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

$dQ=Q_{e}dx-d{Q}_{lost}$E76

with:

$\frac{d{Q}_{lost}}{dt}=h_{c}A(T-T_{wall})$E77

and defining the following dimensionless variables:

$\tilde{p}=\frac{p}{p_1}$E78

$\tilde{V}=\frac{V}{V_1}$E79

$\tilde{T}=\frac{T}{T_1}$E80

$\tilde{Q}=\frac{Q_e}{P_1{V}_1}$E81

${\tilde{Q}}_{lost}=\frac{Q_{lost}}{P_1{V}_1}$E82

${\tilde{h}}_c=\frac{h{T}_1\left(A_{cc}-\frac{4V_{cc}}{D_p}\right)}{P_1{V}_1\omega }$E83

$\beta =\frac{4V_1}{D_p\left(A_{cc}-\frac{4V_{cc}}{D_p}\right)}$E84

dimensionless heat losses for unit crankshaft rotation angle are:

$\frac{d{\tilde{Q}}_{lost}}{d\phi }={\tilde{h}}_c\left(1+\beta \tilde{V}\right)\left(\tilde{p}\tilde{V}-{\tilde{T}}_{wall}\right)$E85

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\phi }=-\gamma \frac{\tilde{p}}{\tilde{V}}\frac{d\tilde{V}}{d\phi }+\frac{\gamma -1}{\tilde{V}}\left[\tilde{Q}\frac{dx}{d\phi }-{\tilde{h}}_c\left(1+\beta \tilde{V}\right)\left(\frac{\tilde{p}\tilde{V}}{\tilde{m}-{\tilde{T}}_{wall}}\right)\right]-\frac{\gamma C\tilde{p}}{\omega }$E86

$\frac{d{\tilde{Q}}_{lost}}{d\phi }={\tilde{h}}_c\left(1+\beta \tilde{V}\right)\left(\frac{\tilde{p}\tilde{V}}{\tilde{m}-{\tilde{T}}_{wall}}\right)$E87

$\frac{dm}{d\phi }=-\frac{Cm}{\omega }$E88

$\frac{d\tilde{W}}{d\phi }=\tilde{p}\frac{d\tilde{V}}{d\phi }$E89

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$$\overline{m}$$we obtain:

$$\overline{m}=\exp \left[\frac{-C_{bb}(\phi +\pi )}{\frac{\pi (rpm)}{30}}\right]$$E90

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:

$$Nu=a{\mathrm{Re}}^b$$E91

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

$$\frac{{\overline{h}}_{c}L}{k}=a{\left(\frac{\rho VL}{\mu }\right)}^b$$E92

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}{\mathrm{Re}}^b\Delta T+c\left(T_{gas}^4-T_{wall}^4\right)$$E93

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_{lost}={\overline{h}}_c*A*\left({\overline{T}}_{gas}-{\overline{T}}_{wall}\right)$$E94

where$${\overline{h}}_c$$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}_{gas}^{0.75-1.62m}{w}^m$$E95

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:

$$w=\left[C_{1w}{V}_{mp}+C_{2w}\frac{V_d{T}_{ref}}{p_{ref}{V}_{ref}}(p-p_{nocomb})\right]$$E96

Subscript ref is used for a reference state that may be the compression process beginning. The empirical constant values depend on the process are: C_1w = 6.18 during gas exchange process, C_1w = 2.28 during compression, combustion and expansion processes, C_2w = 0.0 during gas exchange and compression processes and C_2w = 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=C_{1H}{V}_d^{-0.06}{p}^{0.8}{T}^{-0.4}{(V_{mp}+C_{2H})}^{0.8}$$E97

where C_1H and C_2H are empirical constants which take into account local variations due to intake turbulence or combustion chamber geometry of the, their values are: C_1H = 130E-4 and C_2H = 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

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_{total}=m_{aire}+m_{gr}+m_f+m_{EGR}-m_{bb}$$E98

The theoretical incoming air mass is calculated by:

$$m_{aire}=\rho V_d$$E99

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

$$m_{aire}={{\displaystyle \sum \stackrel{•}{m}}}_{adm}t$$E100

$${\stackrel{•}{m}}_{adm}$$ 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{array}{l}{m}_c{C}_m{H}_n{O}_l+\frac{n+\frac{m}{4}+\frac{l}{2}}{\varphi }\left(O_2+3.7N_2\right)\to \\ \to y_1{H}_{2}O+y_2{H}_2+y_{3}OH+y_{4}H+y_5{N}_2+y_{6}NO+y_{7}C{O}_2+y_{8}CO+y_9{O}_2+y_{10}O\end{array}$$E101

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:

$$C:y_1+y_8=n$$E102

$$H:2y_1+2y_2+y_3+y_4=m$$E103

$$O:y_1+y_3+y_6+2y_7+y_8+2y_9+y_{10}=l+2\frac{\left(n+\raisebox{1ex}{$m$}\!\left/ \!\raisebox{-1ex}{$4$}\right.-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)}{\varphi }$$E104

$$N_2:2y_5+y_6=7.428\frac{\left(n+\raisebox{1ex}{$m$}\!\left/ \!\raisebox{-1ex}{$4$}\right.-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)}{\varphi }$$E105

$${\displaystyle \sum _{i=1}^{10}{y}_i=1}$$E106

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

$$H_2\leftrightarrow 2H{K}_1=\frac{y_4^2}{y_2}\left(\frac{p}{p_0}\right)$$E107

$$O_2\leftrightarrow 2O{K}_2=\frac{y_4^2}{y_9}\left(\frac{p}{p_0}\right)$$E108

$$H_2+O_2\leftrightarrow 2OH{K}_3=\frac{y_3^2}{y_2{y}_9}$$E109

$$\frac{1}{2}{O}_2+\frac{1}{2}{N}_2\leftrightarrow 2OH{K}_4=\frac{y_6^2}{y_9^{\frac{1}{2}}{y}_5^{\frac{1}{2}}}$$E110

$$H_2+\frac{1}{2}{O}_2\leftrightarrow H_{2}O{K}_5=\frac{y_1}{y_2{y}_9^{\frac{1}{2}}}{\left(\frac{p}{p_0}\right)}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$$E111

$$CO+\frac{1}{2}{O}_2\leftrightarrow C{O}_2{K}_6=\frac{y_7}{y_8{y}_8^{\frac{1}{2}}}{\left(\frac{p}{p_0}\right)}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$$E112

In Eqs. 101 to 111 p is the products pressure, y₁ to y₁₀ 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:

$$U_R=U_P$$E113

$$H_R=H_P$$E114

3.11. NOx formation model

The nitrous oxides (NOx) produced by a RICE are NO and NO₂ 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:O+N_2\iff NO+N$$E115

$$R_2:N+O_2\iff NO+O$$E116

$$R_3:N+OH\iff NO+H$$E117

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]$$E118

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[N\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[NO\right]\left[H\right]$$E119

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

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

$$\frac{d\left[NO\right]}{dt}=\frac{2R_1\left\{1-{\left(\left[NO\right]/{\left[NO\right]}_e\right)}^2\right\}}{1+\frac{\left(\left[NO\right]/{\left[NO\right]}_e\right)R_1}{\left(R_2+R_3\right)}}$$E121

and Eqs. 114 to 116 are transformed into:

$$R_1:k_1^{+}{\left[O\right]}_e{\left[N_2\right]}_e$$E122

$$R_2:k_2^{-}{\left[NO\right]}_e{\left[O\right]}_e$$E123

$$R_3:k_3^{-}{\left[NO\right]}_e{\left[H\right]}_e$$E124

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)$$E125

$$K_2^{-}=1.5E09\exp \left(-19500/T\right)$$E126

$$K_3^{-}=2.0E14\exp \left(-23650/T\right)$$E127

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:

$$Fuel+Air+\%EGR\stackrel{}{\to }Products$$E128

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\mathrm{........}{p}_{n-1}\Delta V+p_n\Delta V\right\}$$E129

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

$$W_{pum}=W_{exh}-W_{adm}$$E130

$$W_i=W_{\exp }-W_{com}$$E131

$$W_{net}=W_i-W_{pum}$$E132

Net power is obtained through the following expression:

$${\stackrel{•}{W}}_{net}=W_{net}\frac{rpm}{30j}$$E133

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

$$mip=W_i/V_d$$E134

$$mi{p}_{net}=W_{net}/V_d$$E135

$$mi{p}_{pum}=W_{pum}/V_d$$E136

Efficiencies are determined with:

$${\eta }_i=\frac{W_i}{m_c{H}_i}$$E137

$${\eta }_{inet}=\frac{W_{net}}{m_c{H}_i}$$E138

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

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.

Combustion 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
$${\stackrel{\text{•}}{\text{W}}}_{\text{net}}$$ [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

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

Engine speed	1500 rpm	2000 rpm	2500 rpm
ηv	0.92	0.69	0.55
β	50.57	45.33	40.89
Wi [kJ]	1.42	1.18	1.03
$${\stackrel{\text{•}}{\text{W}}}_{\text{net}}$$ [kW]	14.63	16.01	17.28
mip [kPa]	2314.47	1922.25	1678.68
ηi	0.71	0.74	0.77
pmax [kPa]	6703.23	6034.92	5506.50

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.

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
$${\stackrel{\text{•}}{\text{W}}}_{\text{net}}$$ [kW]	10.99	11.36	11.68
mip [kPa]	1200.22	1233.52	1258.36
ηi	0.45	0.46	0.47
Pmax [kPa]	7134.93	7791.11	8349.81
Tmax [K]	2400.71	2399.34	2392.72

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

Atmospheric pressure	100 kPa	90 kPa	80 kPa
ηv	0.72	0.72	0.72
β	45.25	45.93	46.59
Wi [kJ]	0.74	0.65	0.58
$${\stackrel{\text{•}}{\text{W}}}_{\text{net}}$$ [kW]	11.13	9.90	8.80
mip [kPa]	1217.03	1073.99	948.11
ηi	0.45	0.44	0.44
pmax [kPa]	7231.67	6406.81	5669.70
Tmax [K]	2401.91	2391.06	2380.09

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

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:

• As the engine speed increases the start of combustion should anticipate maintaining nearly constant maximum cycle pressure as well as the angle at which this occurs.

• At higher compression ratios, the greater the potency and mean effective pressure, this coincides with the theoretical internal combustion engines behavior.

• At higher atmospheric pressure the engine volumetric efficiency will increase which ensures a greater air intake and thus an increase in produced power.

Notation

A	area
c	piston stroke
Cd	discharge coefficient
Cp	constant pressure specific heat
Cv	constant volume specific heat
d( )	preceding term ordinary derivative
Dp	piston diameter
dv	valve diameter
E	energy
e	specific energy
f	fuel air ratio
h	specific enthalpy
H	enthalpy
hc	convection heat transfer coefficient
$h_c$	mean convection heat transfer coefficient
Hi	lower heating value
ID	ignition delay
j	engine strokes
k	conduction heat transfer coefficient
Lv	valve lift
m	mass
$\dot{m}$	mass flow
$\tilde{m}$	adimensional mass
mip	mean indicated pressure
Nu	Nusselt number
p	pressure
$\tilde{p}$	adimensional pressure
pdown	downstream pressure
pup	upstream pressure
Q	heat
$\tilde{Q}$	adimensional heat
$\dot{Q}$	heat flow
q	heat per unit mass
R	gas constant
rc	compression ratio
Re	Reynolds number
RLA	connecting rod length to crank radius ratio
rpm	engine speed rpm
T	temperature
$\tilde{T}$	adimensional temperature
t	time
u	specific internal energy
U	internal energy
V	volume
$\tilde{V}$	adimensional volume
v	specific volume
Vmp	piston mean velocity
W	work
$\dot{\text{W}}$	power
w	specific work
Xb	burned mass fraction
yi	i gas fraction

Greek. letters

β	premixed phase burning factor
γ	polytropic coefficient
∂( )	preceding term partial derivative
η	efficiency
μ	viscosity
π	3.1416
ρ	density
φ	crankshaft rotation angle
ω	angular velocity
ϕ	equivalence ratio (mixture richness)

Subindex

a	air
adm	admission
ar	real air
at	theoric air
bb	blow by
cc	combustion chamber
cil	cylinder
com	combustion
com	compression
d	displaced
dif	diffusive phase
e	input
EGR	exhaust gas recirculated
exh	exhaust
exp	expansion
f	fuel
gr	residual gasses
i	indicated
lost	lost
max	maximum
net	net
no comb	without combustion
P	products
pre	premixed phase
pum	pumping
R	reactants
ref	reference condition
s	output
sto	stoichiometric
Tot	total
v	volumetric
vc	control volume
wall	cylinder wall