Abstract - A code for computational simulation of internal combustion engines is presented. One- dimensional gas dynamics equations are used for model the flow through pipes and manifolds, and the remaining components in the engine (cylinders, valves, etc.) are modeled by using thermodynamic or 0D models. The numerical code developed is able to simulate sparkignition and compression-ignition, two-stroke and four- stroke, multi-cylinder and multi-valve engines, naturally aspirated or turbo-charged, and different geometries of the combustion chamber. The code was implemented in the scripting language Python, which is a dynamic object-oriented programming language that offers strong support for integration with other languages and tools. The numerical methods used in the discretization of the equations and implementation details are presented. Several test cases are included in order to show the performance of the code.

The modeling of reciprocating and rotary internal combustion (IC) engines is a multidisciplinary subject that involves thermodynamics, fluid mechanics, turbulence, heat transfer, combustion, chemical reactions, mathematical analysis, and numerical methods. Historically, different levels of approximation have been used to predict the performance of IC engines, from simple air standard cycles to complex 3D models including turbulence, chemical reactions, spray dynamics, etc. IC engine simulation can be classified into four categories, namely zero-dimensional single zone, 0D/1D single zone models, quasi-dimensional multizone models and multidimensional models. In 0D/1D models, the engine is represented as a network of pipes (intake and exhaust manifolds) interconnected among them with "devices" that simulate different parts of the machine (valves, cylinders, pipe junctions, etc.). One-dimensional CFD models are used for pipes and thermodynamic (or zero-dimensional) models for the above mentioned "devices".

For 0D models, most properties are averaged over the total volume and no spatial information is available. A survey of thermodynamic models for cylinders are presented by Blumberg et al. (1979), Mattavi et al. (1980), Heywood (1980), among others. These models rely on some understanding of the physics involved and try to capture the main features of the processes. By including the description of the most important aspects, the models have performed surprisingly well and are ideally suited for parametric studies. A zero-dimensional single-zone model is capable of predicting engine performance and fuel economy accurately with a high computational efficiency (Krieger and Borman, 1966; Foster, 1985; Assanis and Heywood, 1986). The major drawback of single-zone models is their inability to simulate the wave propagation into pipes and manifolds that strongly influence on volumetric efficiency. Also, these models are unable to account for fuel spray evolution and the spatial variation in mixture composition and temperature, both of which are essential in predicting harmful species formed during the combustion process.

On the other hand, multi-dimensional models, such as KIVA (Oran and Boris, 1981; Bracco, 1985; Amsden et al., 1985; Amsden et al., 1987; Varnavas and Assanis, 1996), resolve the cylinder space into fine grids, thus providing a considerable amount of spatial information. However, multi-dimensional models still employ phenomenological submodels describing fuel spray processes, and their simulation results may vary with assumed initial or boundary conditions. Consequently, the accuracy of the results cannot always be guaranteed. Furthermore, computational time and storage constraints still preclude these models from routine use for design purposes. Currently an intermediate step between zero-dimensional and multi-dimensional models has arisen, called quasi-dimensional. Multi-zone models (Hiroyasu and Kadota, 1976; Hiroyasu et al., 1983a,b; Kyriakides et al., 1986; Yoshizaki et al., 1993; Kouremenos et al., 1997; Rakopoulos and Hountalas, 1998; Jung and Assanis, 2005) can be effectively used to simulate new technology engine combustion systems, by combining the advantages of zero-dimensional models and multi-dimensional models. These models are able to provide the spatial information required to predict emission products with significantly less requirement on computing resources than for multidimensional models.

In this work we start to use 0D/1D models in order to obtain a computational tool that can predict with sufficient precision the performance of an IC engine at a relatively low computational cost, with the target of arriving to more sophisticated quasi-dimensional models in the future. Usually, the 0D/1D engine simulators employ explicit schemes for integration in time. Hence, the 0D models are generally formulated according to an explicit scheme. Because most of these 0D models are non-linear and also the need to reduce the time step of the simulation for stability reasons when explicit schemes are applied, choosing an implicit integration can provide more robustness to 0D/1D codes.

This project began with the development of a single cylinder four-stroke spark-ignition engine simulator. The mathematical model was based on a thermodynamic model for the cylinder and a one-dimensional gas dynamics description of the intake and exhaust systems (Nigro et al., 1999). Then, with the target focused on the prediction of more real situations, models for pipe junctions in multi-cylinder configurations were added. In order to check the suitability and reliability of this computational tool in industrial applications, this code was rewritten in the compiled language Fortran 90/95 and several test cases were solved validating the results with measurements (Alianak and Nigro, 2003). In the same reference cite an interesting optimization work was done, where the coupling between the 0D/1D code and other mechanical design tools for sizing camshafts are shown. This optimization task demanded an extensive usage of our IC engine simulator showing as a by-side product its inherent robustness. In addition of the overall performance prediction of IC engines, our goal is also to be able to make 3D simulations of in-cylinder flows. Due to the impossibility to simulate an IC engine completely with 3D models, it is common to use 0D/1D codes to simulate the whole engine except the component that will be studied in detail. Thus, in this application, the 0D/1D engine simulator could be used as a generator of dynamic boundary conditions. In this line of work, the Fortran 90/95 version of the code was coupled with the KIVA-3V (Amsden et al., 1993) code. The coupling was done applying a boundary condition for the 3D code based on the mass flow rate and the pressure obtained from the 0D/1D code, for which was adopted an absorbing boundary condition strategy to compute the flow state at the coupling interface (Bella et al., 2003). Generally, CFD-3D codes employ an implicit scheme of integration in time, therefore, the availability of a version of the 0D/1D simulator with implicit integration gives greater generality to the code as a generator of boundary conditions.

The new code was written using the language Python (van Rossum, 1990-2007) because it offers strong support for integration with other languages, such as C++, Fortran, etc. Also, Python offers the advantages of the object-oriented programming and has a clear and simple syntax. The numerical code developed is able to simulate spark-ignition and compression-ignition, two-stroke and four-stroke, multi-cylinder and multi-valve engines, naturally aspirated or turbo-charged. Also, different geometries of the combustion chamber are available.

As future work we will use this software tool as the simulation stage of an optimization package that allow to setup the engine for performance improvement. Moreover, our interest would be to use this kind of approximation in real time simulation software.

The paper is organized as follows. The next section presents the mathematical models used for each component. Next, the numerical methods applied and some details about the implementation are included. Finally, numerical results for different kind of engines are shown and some conclusions are presented.

II. MATHEMATICAL MODEL

A. Pipe model

One-dimensional unsteady flow equations are used for modeling pipes and manifolds. In order to include effects like variable cross-section, viscous friction, and wall heat transfer, some source terms are added to the inviscid gas dynamic model represented by the system of Euler equations. The resultant system of equations can be written as (Heywood, 1988)

(1)

where ? is the density; p is the pressure; u is the fluid velocity; F is the pipe cross-section area;

is the specific friction force, with the friction coefficient given by f = 2τw/?u2, τw being the viscous shear stress at the pipe wall and D the equivalent diameter of the pipe; E is the total specific energy of the fluid; and is the heat transfer per unit mass of fluid per unit time. The total specific energy is related to the internal energy per unit mass e and specific kinetic energy as

The equation of state used here corresponds to the ideal gas assumption with particular gas constant Rgas.

B. Cylinder model

A single-zone model is used to model the cylinder. In this model the charge is assumed to be a homogeneous mixture of ideal gases at all times. The equations of the model are the conservation of mass and the first law of thermodynamics

where m is the mass contained in the cylinder; is the instantaneous mass flow rate through the j-th inlet/outlet (for instance, the mass flow rate through intake and exhaust valves, fuel addition, leakages, etc.); e is the specific internal energy of the mixture; V is the cylinder volume; represents the heat release due to combustion; is the heat transfer rate; and represents the enthalpy fluxes through the j-th inlet/outlet.

The model is closed specifying the geometry of the combustion chamber, the heat release rate, the heat transfer rate through the cylinder walls, and the mass flow rate of air and fuel. The sub-models used here are presented in the following sections.

Geometry of the combustion chamber

The geometric data necessary for the cylinder model are the total surface area of the cylinder walls (A), the volume of the cylinder and its time derivative. By the time, it is possible to simulate conventional reciprocating engines, opposed-piston engines, and the MRCVC (Motor Rotativo de Combustión a Volumen Constante; Toth, 2004; Toth et al., 2000). For all these cases, we implemented analytical formulae to compute the geometric variables of the combustion chamber.

Heat transfer model

The instantaneous heat transfer rate that appears in Eq. (2) is calculated applying Nusselt-Reynolds-Prandtl (Nu-Re-Pr) numbers correlations as, for example, the one developed by Woschni (1967) or by Annand (1963). All of them allows to compute a film transfer coefficient hc with expressions like the following

(3)

where L is a characteristic length; ? is the gas thermal conductivity; and C, a and ß are constants.

Then, the heat transfer rate to the walls is

(4)

where T is the temperature of gas into the cylinder, and Twall is the temperature of the cylinder wall.

Heat release model

In order to modeling combustion we use several approaches and mathematical models, which have the goal to describe the actual heat release via combustion as exactly as possible by means of the so-called substitute heat release rates.

For spark-ignition engines, the mass fraction of burnt gases (xb) is computed by using a Wiebe function (Heywood, 1988)

(5)

In Eq. (5) mb represents the mass of burnt gases; T is the crank shaft angle; ?T is the duration of combustion; Tig denotes the angle at which burning starts; c and s are parameters, where s is designated the shape parameter and c accounts for combustion efficiency. The heat release rate can be computed as

(6)

Hc being the calorific heat content of the fuel and m f the total mass of fuel trapped into the cylinder.

where τ = (T - Tig)/?T, and O is the angular speed of the shaft. The coefficients proposed in the original model are

, where N is the engine speed, in rpm; and t d is the ignition delay, in ms.

c2 = 5000.

being the equivalence ratio at the time of ignition.

The double Wiebe function is an extension of the model used for spark-ignition engines, in order to describe the premixed and diffusive combustion periods observed in diesel engines (Ramos, 1989). The mass fraction of burnt gases can be written as

(8)

where xp is the mass fraction of fuel burnt in the premixed combustion period, xdi is the mass fraction of fuel burnt in the diffusive combustion period, ?Tp and ?Tdi are, respectively, the duration of premixed and diffusive combustion.

In diesel engines, the ignition delay time could be calculated as the difference between the time at which combustion starts (tig) and the time at which injection starts (tinj).

The time tig can be obtained from the following expression (Assanis, 1985):

(9)

which accounts for the pressure and temperature variations resulting from compression. The ignition delay time (td) as a function of T and p was correlated for a variety of fuels with the expression

being C, n and Ta constants (Assanis, 1985). Also, the empirical formula developed by Hardenberg and Hase (1979) for predicting the ignition delay time in direct-injection diesel engines was implemented in the code.

Scavenge model

The scavenge process in two-stroke engines is modeled via a semi-empirical model proposed by Blair (1990). The mass of delivered air trapped into the cylinder at exhaust closure (mat) is computed as mat = ?Smt, where mt is the total mass retained after the exhaust port closure, and the scavenging efficiency (?S) can be computed as

?S = 1 - exp (bSv + d) (10)

b and d being constants experimentally determined, which depend on the type of scavenge. Sv is the scavenge ratio by volume, defined as the ratio between the volume of air supplied during the scavenge period and the cylinder volume.

C. Valve model

To calculate the flow rates through the intake and exhaust valves, we use an analogy with the steady flow through convergent nozzles proposed by Benson (1982). The model assumes the passage area through the valve as the nozzle throat (whose state is represented by the subscript T in the equations below), and the nozzle connecting the cylinder (subscript C in the equations) and the end of the pipe (subscript P in the equations). Depending on the direction of the flow velocity with respect to the pipe end, the problem could be an inlet (from cylinder to pipe) or an outlet (from pipe to cylinder). In addition, the flow at the throat could be sonic or subsonic. The equations of the model are presented below.

Subsonic inlet:

(11)

where is the specific heat ratio, a is the sound speed, and

(12)

In system (11), the first equation accounts for the compatibility along the incoming Mach line ?±; the second equation is the mass conservation between T and P; the third and fifth equations represent the energy conservation between C and P, and C and T, respectively; the fourth equation represents an isentropic evolution between the cylinder and the nozzle throat; and the last equation is the condition on the pressure at the nozzle exit.

Subsonic outlet:

(13)

with

(14)

and D/Dt denoting the material derivative. From first equation to the last one in the system (13), they represent, respectively, the compatibility along the incoming Mach line, the mass conservation between T and P, the compatibility along the incoming path line ?0, the isentropic evolution between P and T, the energy conservation between P and T, and the condition on the pressure at the nozzle exit.

Sonic inlet: in this case, the system of equations is the same as Eq. (11) with the last equation replaced by the condition aT = uT.

Sonic outlet: in this case, the system of equations is the same as Eq. (13) with the last equation replaced by aT = uT.

D. Pipe junction model

The pipe junction model applied was proposed by Corberan (1992). If the junction is composed by r incoming pipes and s outgoing pipes, the model is expressed as

Mass conservation

where N = r + s is the total number of pipes at the junction, Fj is the cross-section area of the j-th pipe and nj its exterior normal.

Energy conservation

Compatibility equation along incoming Mach lines

Compatibility equation along incoming path lines

Equality of pressure at all branches in the junction

Equality of enthalpy at all outgoing branches in the junction

III. NUMERICAL IMPLEMENTATION

The development of the code presented aims, in addition to its use as a tool for simulation of IC engines, to be a generator of boundary conditions for CFD codes. Due to the computational cost of the CFD simulations on engines, its resolution is usually done using an implicit scheme of integration in time. Therefore, it is necessary that the 0D/1D engine simulator could work with an implicit scheme avoiding a severe restriction on the time step.

The use of an implicit scheme for the discretization of equations that model the several devices in the code do not represent a problem for its practical implementation, with the exception of models of valve and pipe junction. As presented above, the systems of equations that model valves and pipe junctions change with the flow regime and the direction of its velocity.

In the valve model, to determine whether the problem is an inlet or an outlet, we compare the cylinder pressure with the stagnation pressure (p0P) at the pipe end, which is given by (White, 1983)

(15)

The valve is considered open when the passage area is larger than a prefixed tolerance Then, if pC > p0P , the flow is established from the valve to the pipe, otherwise it is considered an outgoing flow from the pipe. To take into account the transition between subsonic and sonic regime flow, we use the following convex combination

(16)

where Rsubsonic and Rsonic represent the systems of equations that model the subsonic and sonic cases, respectively; and

(17)

a being a coefficient that adjusts the transition of χ between 0 and 1. In valve model, the compatibility equations along the characteristic curves are used to complete the system of equations. These must be solved according to an explicit scheme, implying that the states at the pipe end, the valves, and the cylinder, do not depend on the state at interior points in the pipe at time t + ?t. Thus, the system of equations for valve+cylinder is decoupled from the remaining equations. However, the resolution is done in a coupled way with the goal of being able to implement other valve models of full implicit type in the code. When the valve is either opening or closing, it is important to determine the direction of the flow due to the change in the system of equations that models each situation, which directly influences on the global convergence. To predict the flow direction it is assumed an isentropic flow through a convergent nozzle between the cylinder and the corresponding end of the pipe. The flow is established between the states U0 and Ue, which are assumed constant and are identified with either the state in the cylinder (UC) or the state at the end of the pipe (UP), depending on the relationship between the pressure in the cylinder and the stagnation pressure p0P. If pC ≥ p0P, then it is adopted U0 = UC and Ue = UP , the direction of the flow being from the cylinder to the pipe. When pC < p0P , then U0 = UP and Ue = UC, yielding to a flow from the pipe to the cylinder. The critical pressure

(18)

determines if the nozzle throat is choked or not. If pe > pcrit, the predicted state at valve (UT), i.e. the nozzle throat, is taken as

(19)

If pe = pcrit, then it is assumed

(20)

To solve the system of equations (1), the Finite Element Method stabilized with the Streamline Upwind/Petrov- Galerkin (SUPG) technique (Hughes and Mallet, 1986; Tezduyar and Hughes, 1983) was used. Time derivatives were discretized applying the trapezoidal finite difference scheme.

The code was implemented in the scripting language Python (van Rossum, 1990-2007). Python is a dynamic object-oriented programming language that can be used for many kinds of software development. It offers strong support for integration with other languages and tools, and comes with extensive standard libraries. Object-oriented programming allows to develop the code in an organized manner, and the possibility of integration with other languages make it suitable for solving the coupling between the 0D/1D simulator and CFD codes. Another feature that Python offers is the possibility of writing higher-level parts of large-scale scientific applications and driving simulations in parallel architectures like clusters of PCs. This can be done by using, for instance, the package mpi4py (Dalcín et al., 2005) which provides bindings of the Message Passing Interface (MPI) standard (MPI-Forum, 1994- 2008). This feature can be useful to include the engine simulator in a parametric optimization strategy or to couple with other codes that run in parallel.

Regarding the computational implementation, a class was defined for each device with methods that perform the computation of sub-models and the contributions to both global residue and global jacobian matrix. For instance, the class Cylinder contains the methods combustion, geometry, heat_transfer, etc. All of these classes derive from the classes that define an specific set of equations for each device. For example, the class Flow involve the equations of mass, energy and linear momentum conservation. This allows to add new features to the code as, for example, transport of chemical species with very few changes. Currently, only the class Flow is implemented in the code. The constructor of each class defines attributes for the object, some of them are required data and the remaining ones are optional (defined through default values). The nonlinear global system of equations is solved via the Newton-Raphson method. The linear system arising at each nonlinear iteration can be solved either using functions contained into the package NumPy (Oliphant, 2006), or with the package petsc4py (Dalcín, 2005-2007). petsc4py is a set of Python bindings for PETSc, the Portable, Extensible Toolkit for Scientific Computation (Balay et al., 1995-2007). Currently, we are working to improve the efficiency of the code using Cython (Behnel et al., 2007-2009).

The code also allows to define the parameters of calibration and the operational variables as functions of the engine speed, the cycles, and the time.

IV. RESULTS

In the next sections we will present the results obtained from the application of the code to the simulation of some IC engines. The cases were selected to show the spectrum of IC engines that can be simulated.

A. Four-stroke spark-ignition engine test

The first example of application is a 8.4 l. V10 four-stroke spark-ignition engine. Tables 1 and 2 contain the main data of the engine.

Table 1: Main cylinder data of the V10 engine.

Table 2: Valve data of the V10 engine.

The equivalence ratio, the coefficients of heat transfer through the cylinder walls, and the crank angle where the combustion starts, are functions of the engine speed. These functions were obtained from the experimental data available. The engine was tested at speeds ranging from 1600 rpm to 6000 rpm.

Figures 1 to 3 show the computed indicated power, the torque and the average mass flow rate of air versus the real curves obtained experimentally.

Figure 1: Indicated power as a function of rpm.

Figure 2: Torque as a function of rpm.

Figure 3: Average mass flow rate of air at intake port as a function of rpm.

B. Two-stroke spark-ignition engine test

This test case was taken from the literature (Blair, 1990) and consists in a two-stroke spark-ignition single-cylinder research engine denominated QUB 400. The engine speed is 3000 rpm at full throttle. To model the crankcase compression we use a cylinder without combustion, then the computational model is composed by two cylinders as can be observed in Fig. 4. The engine has six transfer ducts between the cylinder and the crankcase, which were modeled as a unique pipe with the same total cross-sectional area (see Fig. 4). Ports were assumed as poppet valves with

Figure 4: Computational model of the QUB 400 engine.

appropriate discharge coefficients, and where the valve lift was defined in such a way to make the passage area computed by the code the same as the instantaneous passage area of the port. The complete set of data can be found in the literature by Blair (1990). Figure 5 shows the pressure in the cylinder, the crankcase, the transfer duct, and the intake and exhaust ports as a function of the crank angle during a cycle. Mass flows through transfer and exhaust ports are shown in Fig. 6, where positives values represent incoming flow to cylinder. Table 3 shows some performance characteristics: power, indicated mean pressure (IMEP), scavenge efficiency (?S), peak cylinder pressure (Pmax) and the crank angle at which it occurs. In order to compare, results reported by Blair (1990) and experimental data are included. These results were in a very good agreement with those presented by Blair (1990) coming from his numerical simulations and also with his experiments.

Figure 5: Variation of pressure during a cycle.

Figure 6: Mass flow rates through transfer and exhaust ports.

Table 3: Performance characteristics of the QUB 400 en- gine.

C. Four-stroke diesel engine test

The test case selected to show the application of the code to a diesel engine was the truck diesel KamAZ-7405, taken from the web-page http://www.diesel-rk.bmstu.ru/. The main engine data are the following:

Table 4 contains the operational data of the three cases solved. In that table, pi and Ti are the intake manifold mean pressure and temperature, respectively; pe is the exhaust manifold mean pressure; and minj is the mass of fuel injected per cycle.

Table 4: Operational variables of the KamAZ-7405 diesel engine.

The combustion was modeled using two Wiebe functions, and the ignition delay time was calculated with the correlation proposed by Hardenberg and Hase.

The calibration of the code was done at 2200 rpm, defining the parameters for the two Wiebe functions, the coefficients of the heat transfer model, and the duration of the combustion. The results are presented in Tab. 5 together with the relevant experimental results, where SFC is the specific fuel consumption and is the mean mass flow rate of air through the intake system. There is in general a good agreement with experimental data, being the error within typical error margins for the type of code employed.

Table 5: Performance characteristics of the KamAZ-7405 diesel engine.

V. CONCLUSIONS

A code for computational simulation of internal combustion engines was presented. The mathematical models and the numerical methods used were described briefly. Some models and the solution strategies were reformulated to use an implicit scheme of integration in time. The code was written in the language Python in order to take advantage of the object-oriented programming and the possibility of integration with other languages. Some test cases were solved, including spark- and compression-ignition, two- and four-stroke IC engines. After code calibration, results were in good agreement with experimental data.

As future works we propose to use the code presented as a generator of boundary conditions for CFD-3D codes through an appropriate strategy of coupling. In addition, we intend to incorporate the 0D/1D simulator within a code of parametric optimization for IC engines.

ACKNOWLEDGMENTS This work has received financial support from Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina, grants PIP 5271/05), Universidad Nacional del Litoral (UNL, Argentina, grant CAI+D 2005-10-64) and Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT, Argentina, grants PICT 12-14573/2003, PME 209/2003). Authors made extensive use of freely distributed software as GNU/Linux OS, Python, G95 compiler, Octave, among many others. Special thanks are due to Lisandro Dalcín for his support with the Python language, and also to Mario Storti for his contribution into the discussion of the ideas.