Simulation Studies of Combustion in Spark Ignition Engine using Openfoam

The open source Field Operation and Manipulation (OpenFOAM) software was used to investigate the performance of a fully premixed, modern highperformance 4-valve, iso-octane, dual overhead cam (DOHC) engine with quasi-symmetric pent roof combustion chamber running at 1500 revolutions per minute. The peak pressure occurred at the TDC and had a value of about 30 bar. The results from this study show that the maximum combustion temperature occurred at approximately 95 degrees crank angle ATDC and has a volume averaged value of about 2700 K ° , whereas the actual computed peak temperature was found to be about 3000oK and it occurred at grid point 12630. The other temperatures which were found to be higher than the volume averaged temperature were found to be in the range 2968.81 2974.01 K to K ° ° and correspond to grid point positions 12630 to 12633.The flame-wrinkling factor, / t u S S = Ξ was found to be in the range 1.0 ≤ Ξ ≤ 3.8. The dynamics of the regress variable b was accurately predicted.


INTRODUCTION
There is presently increasing concern regarding environmental pollutions (greenhouse gas emissions, various oxides of nitrogen, particulate matters, etc.) originating from the operation of internal combustion engines. In view of this, stringent emission control laws are being promulgated by various Environmental Protections Agencies (EPAs). Consequently, there is an urgent need to control automotive contributions to urban air pollution and to achieve significant improvements in internal combustion engine (ICE) fuel economy. As a result of these stringent emission control laws and standards, research efforts have led to significant reductions in emissions from gasoline and diesel engines. In spite of the advent of new power sources such as hybrid electric vehicle (HEV), plug-in hybrid electric vehicle (PHEV), fuel cell, etc., for powering the automobile, it is worth mentioning that these prime movers still require augmentation with internal combustion engines with specific characteristics in order to improve their efficiencies.
Even though the electric vehicle is gradually making in-road as an alternative or possible replacement for the internal combustion engine, the complete realization of this prospect is presently not certain.
In view of the uncertainty surrounding having one hundred percent electric vehicle and the thought that after more than a century of development of the internal combustion engine a lot of prospect for further development still exists. It is worth emphasizing that, there are a lot of empirical evidences which show that as the gasoline and diesel engine technology evolves, they continue to show substantial improvements in efficiency, power-density, improved emission and operational capacity.
In order to comply with the aforementioned stringent regulatory controls, research efforts which are geared towards clean diesel and gasoline engines are fundamental in guiding the development of new generation internal combustion engines. Advances in cleaner and more efficient combustion are as vital as other parameters such as alternative fuels, improved catalyst technologies or electronically controlled fuel injectors. Research activities relating to clean combustion require a comprehensive understanding of fluid dynamics, reaction kinetics, heat transfer, turbulent combustion and their interactions within the combustion chamber of an internal combustion engine.
Recent developments in combustion modeling have greatly improved the quality of internal combustion engine simulations. Furthermore, the complex dynamics /physics of combustion kinetics, pollutant formation, interaction of turbulence phenomena with other physical processes can now be modeled and analyzed. Modeling of reaction kinetics in turbulent combustion has advanced remarkably with efficient techniques, such as, the intrinsic low dimensional manifold [1], the flame prolongation of low dimensional manifold [2], the flame generated manifold [3] or the in-situ adaptive tabulation [4]. With these techniques it is now possible to comprehensively describe reaction kinetics during the combustion process, where hundreds of species and thousands of reactions are the norm. In addition to the aforementioned techniques, sophisticated approaches 788 ▪ VOL. 48, No 4, 2020 FME Transactions such as the G-equations [5], flame surface density [6,7], artificial flame thickening [8,9,10] or conditional moment closure [11,12] are presently been employed to accurately predict flame propagation under the influence of turbulence. In spite of the attendant issues associated with in-cylinder flow and combustion, the Large Eddy Simulation (LES) of internal combustion engine is also being used to analyze a variety of complex issues like NOx formation [13], cycle-to-cycle variations [14,15] and flame-acoustics interactions [16,17]. As long as the demands for low emissions, higher efficiency and improved fuel economy of the internal combustion engine continues, research efforts into the mathematical and numerical techniques for modeling turbulent combustion and associated complex physical phenomena will continue to be of vital interest to engine manufacturers and the engine research community.
In general, two main approaches are usually used for conducting research aimed at exploring ways of optimizing engine design and operation. These methods are experimental techniques which entail building of a real prototype engine and testing it while the other involves numerical simulation of the various engine processes. Because neither of these approaches is able to fully capture all the relevant dynamics at play in an internal combustion engine, the symbiotic relationship existing between them is usually taken advantage of during engine optimization projects. The experimentalist is unaware of the micro details of the physical processes within the engine combustion chamber. The experimentalist only has a macro knowledge of important parameters like global averaged temperature, pressure, etc., whereas the modeler is uncertain about the accuracy of the predicted results. The dilemma that was just pointed out is usually resolved by using the experimental data to validate the modeling results. Once the modeling results have been successfully validated, they can then be used to explore conceptual engines and combustion modes prior to building any real prototypes. Moreover, validated numerical simulations can also be used for parametric studies with the concomitant benefit of savings in time and cost.
As a result of the foregoing, the open source Field Operation and Manipulation (OpenFOAM) was used to investigate the performance of a fully premixed modern, high-performance 4-valve, iso-octane, dual overhead cam (DOHC) engine with quasi-symmetric pent roof combustion chamber running at 1500 revolutions per minute.
The unique aspect of this study is the fact that, the OpenFOAM was used to exhaustively investigate the salient operating characteristics (pressure, temperature, flame-wrinkling factor, / t u S S = Ξ , turbulent kinetic energy, regress variable, etc.) of an internal combustion engine.
This study is structured as follows: section 2 provides a brief description of the Open source Field Operation And Manipulation (OpenFOAM), which is the numerical code that was used for this work. Section 3 presents a brief fundamental philosophy underpinning the design of the OpenFOAM. In section 4, the models characterizing the dynamics of turbulent and chemically reactive flows within the combustion chamber of internal combustion engine are presented. Section 5 describes the premixed turbulent combustion model that was employed in this study. The computational and numerical schemes where presented in sections 6 and 7 respectively. The computational domain, grid and boundary conditions were presented in section 8. The results and discussions of the in-cylinder flow parameters and the combustion process of an ICE were presented in section 9. The conclusion of the present study was given in section 10.

DESCRIPTION OF OPENFOAM
The OpenFOAM numerical simulation code is a versatile tool which is gaining popularity and it is extensively being used for research and process improvements in the scientific and industrial environments. OpenFOAM implements operator-based implicit and explicit second-and fourth-order Finite Volume (FV) discretization in three-dimension and on curved surfaces, a second order Finite Element Method (FEM) solver and a particle tracking model. Flexibility in mesh management is realized by supporting un-structured polyhedral meshes and topological variations. Efficiency of execution is accomplished by the use of preconditioned Conjugate Gradient and Algebraic Multigrid solvers and the use of massively parallel computers in domain decomposition mode. Automatic mesh motion solver, where point motion is defined by only prescribing the boundary motion facilitates the setup of deforming mesh simulations.
The consistent object-oriented architecture of the code and the extensive number of available physical and numerical models recommends it as an excellent tool suitable for the simulation of complex problems.

EQUATION REPRESENTATION
The fundamental philosophy of the OpenFOAM design is that the solver applications are written using the OpenFOAM classes and have a syntax that closely looks like the partial differential equations being solved.
For example, the material or Lagrangian derivative, describing the rate of change of an intensive physical property (t) in time can be written as: which would have the OpenFOAM syntax given in Equation (2)  For compactness these equations are written in vector notations. The unit vectors in the x,y and z directions are denoted by i , j and k respectively.
The position vector r is defined as: the vector operator ∇ is given by: and the fluid velocity vector u is given by: u u x y z t i v x y z t j w x y z t k In order to simulate the dynamics of the processes prevailing inside the combustion chamber of a spark ignition engine, the conservative/transport equations of continuity, momentum, energy and other physical models such as combustion, heat transfer, equation of state, turbulence, etc. must be solved. These equations and the relevant sub-models or constitutive equations for Newtonian fluids are used in the present study and they are given below:

Continuity Equation
The continuity equation for species m can be written as: where ρ m is the mass density of species m, ρ is the total mass density and u is the fluid velocity. D is the coefficient of diffusion and it is assumed to be governed by Fick's Law of diffusion. The constitutive relations for c m ρ the source term due to chemistry and D are given below. By summing Equation (6) over all species the total fluid density equation becomes: since mass is conserved in chemical reactions.

Momentum equation
The momentum equation for the fluid mixture is given by: where p is the fluid pressure.
The viscous stress tensor σ is Newtonian in form, therefore; The first and second coefficients of viscosity μ eff and λ are defined below. The superscript τ denotes the transpose operator and I is the unit dyadic. The specific body force g is assumed constant.

Turbulent Kinetic Energy and Dissipation Models
The kinematic eddy viscosity μ eff can be evaluated in various ways, ranging from algebraic relations and local equilibrium assumptions to the solution of transport equations. The commonest method is to express μ eff as a function of the turbulent kinetic energy k and its dissipation rate ε leading to a two-equation turbulence model. Accordingly, the transport equations for the turbulence kinetic energy, k and its dissipation rate, ε are: and ( ) : Equations (10) and (11) are the standard k -ε model with some additional terms. The source term 3 1 represents the self-diffusion of turbulence with a diffusivity of μ eff /ρ. The term -ρε represents the decay of turbulence into thermal energy. The physical meaning of the various terms in the rate of turbulence dissipation, that is Equation (11) are analogously defined.

Internal energy
The internal energy equation is given by: where I is the specific internal energy, without the chemical energy. The heat flux vector J is the sum of contributions due to heat conduction and enthalpy diffusion, and it is given by: where T is the fluid temperature and h m the specific enthalpy of species m. The source term due to chemical heat release c Q is defined below.

Equations of state
The state constitutive equations are assumed to be those of an ideal gas mixture. Therefore:  Tables [18].

Chemical reaction
The chemical reactions occurring in the combustion chamber are modeled as: Where X m represents one mole of species m and a mr and b mr are integral stoichiometric coefficients for reaction r.
The stichiometric reaction coefficients must satisfy: So that mass is conserved in the chemical reactions. The chemical reactions are divided into two classes, namely, those that proceed kinetically and those that are assumed to be in equilibrium.
Kinetic reaction r proceeds at a rate r ω which is The reaction orders are a' mr and b' mr and they are not necessarily equal to a mr and b mr . Therefore, empirical reaction orders can be used. The coefficients k fr and k br are assumed to be of a generalized Arrhenius form: and,

( )
where E fr and E br are activation temperatures.
The rates of equilibrium reactions are implicitly determined by the constraint conditions: where ( ) c r K T the concentration equilibrium constant is assumed have the form: where T A = T/1000 K. From the foregoing, the reaction rates r ω is determined by Equation (20) or (21), the chemical source term in the species continuity equation is given by: and the chemical heat release term in the energy equation is given by: where Q r is the negative of the heat of reaction at absolute temperature zero.
and ( ) is the heat of formation of species m at absolute temperature zero. The transport coefficients in the equations above are: and D Sc c μ is an empirical constant with a standard value of 0.09. A Sutherland formula is used for μ iar as follows: where A 1 and A 2 are constants. The constant A 3 is taken to be (-2/3) in the calculation of turbulent flow. Pr and Sc are the Prandtl and Schmidt numbers and they are input constants.

PREMIXED TURBULENT COMBUSTION MODEL
The reactive processes taking place inside the combustion chamber of a spark ignition engine are highly complex and turbulent. Active research efforts aimed at understanding these physical processes are presently ongoing and one vital aspect of it is the search for numerical models which adequately and satisfactorily describe the processes involved. It is worth mentioning that the models should provide a better understanding of the processes involved, and provide the ability to predict the behavior of the combustion processes prevalent in the combustion chamber of spark ignition engines.
There are three categories of models which are often used to characterize the complex and turbulent nature of premixed combustion in spark ignition engines. They are the Eddy Break Up (EBU) [19], Thin Wrinkled Flame (TWF) [20,21] and Flame Area Evolution (FAE) [22]

models
In the Eddy Break Up model class the local burning rate is assumed to be controlled by small scale turbulent mixing with a characteristic time scale that is proportional to the turbulence dissipation scale. This proposition is incorrect and thus lead to concave flame front. For accurate near wall behavior the characteristic burn time should not be equated to the eddy turn over time τ t =κ/ε but derived from spectral descriptions of turbulence which takes into account the change in the spectrum of scales during flame kernel formation and on approaching a wall.
The Thin Wrinkled Flame models [23,24] do account accurately for the normally thin (of the order of the Kolmogorov microscale) and wrinkled flame structure. The wrinkling being responsible for the increased burning rate over a smooth flame burning at laminar burning velocity. The laminar flamelet assumption has the advantage of making it possible to include detailed chemistry and flame straining effects in a more computationally efficient manner. However, these models assume that the characteristic scale of wrinkling is proportional to the characteristic turbulence scales thereby leading to the same issues affecting those of the Eddy Break Up model.
The more recently developed Flame Area Evolution models [25,26] share some of the characteristics of the Thin Wrinkled Flame models such as the thin flame assumption and the laminar flamelet burning. The local wrinkled flame area is however, modeled by an evolution equation which allow incorporation of turbulence effects as well as effects generated by the flame itself and non-local effects [27]. This type of modeling has been found to produce realistic results in modeling combustion in spark ignition engines [28] where simulation of combustion was based on the Weller flame wrinkling model [27,29]. It provides a transport equation for the spatial and temporal evolution of the wrinkling factor Ξ. The flame wrinkle density Ξ is the flame area per unit area resolved in the mean direction of propagation. The use of Ξ provides physical justification for modeling the velocity field generated by the combustion, rather than the traditional density weighted ensemble averaging and gradient transport modeling. Moreover, Ξ may be categorized as a spectral density function in length scale phase space, thereby enabling detailed description of the turbulent/flame interaction processes.
In the present work, the 2-equation Weller model for flame distribution, that is, the transport equation for the spatial and temporal evolution of progress variable b and governing equation for the wrinkling factor Ξ is used to model combustion.
In premixed combustion the flow field consists of regions of unburnt reactants and regions of burnt products. The extent of combustion of the air/fuel mixture can be described in a variety of ways in terms of a progress variable, c such that 0 ≤ c ≤ 1, with the extremal values indicating the presence of unburnt or fully burnt mixture phases, and the transition between these extremal values marking the flame front. This can be linked directly to physical properties of the air/fuel mixture, for instance, by utilizing normalized temperature, T or product mass fraction, Y, as follows: Irrespective of how the progress variable is defined, the exact relationship is not important; however, it can simply be viewed as an artifice for tracking the quantity of air/fuel mixture burned. In the present work, we will use a progress or regress variable, b = 1 -c; with 0 ≤ b ≤ 1, where 0 and 1 represents unburnt and fully burnt air/fuel mixture respectively. Suffices B and u represents burned and unburned fuel-air charge.
In view of the foregoing, the models for the wrinkle density Ξ transport equation and the corresponding equation for the regress variable b which may be considered a normalized fuel concentration from its origins as a flame area model to a closed pair of coupled transport equations and simple one equation model derived therefrom are given below: The simplified transport equation for Ξ is given by: Since the solution of the spectral evolution equations coupled with Equation (36) is prohibitively expensive, simple algebraic equations are considered more appropriate. The approach used in the present study is based on flame speed correlation, Equation (37) of Gülder [30], which has proven to be very good when compared to the full spectral solutions is used in this study, and it is given in Equations (37a) and (37b).
where Ξ is the ratio of turbulent flame velocity to laminar flame velocity, Ξ = S t /S ub ; S t is the turbulent flame speed, S ub is the laminar flame speed, τ η is the Kolmogorov time-scale, u' is the subgrid turbulence intensity, s u is the surface filtered velocity of the flame, η ℜ is the Kolmogorov Reynolds number, σ s and σ t and are the resolved strain rates.
The flame front propagation is modeled by solving a transport equation for the density-weighted mean reaction regress variable b as follows: For the interested reader, the derivation of these models can be found in [2,12,15].
The final requirement for closure is to provide a model for the laminar flame speed S ub . This will largely be dependent on the characteristics of the type of flame under consideration. Under certain circumstances it might be reasonable to assume that the laminar flame speed is unaffected by strain and curvature, therefore, we set S ub = 0 ub S , the unstrained flame speed.
In the present study the laminar flame speed S ub is used and it is computed using the Gülder's correlation [30] and it is given by: where P and T are the pressure and temperature respectively, Φ is the equivalent ratio.
and subscript ub represents unburned state.
It is worth mentioning however, that this approach has a drawback. In regions of high Damköhler, Da (highly reactive species) and Reynolds, R e numbers, the flame-sheet reaction zone is very thin and the transition of the progress variable is extremely steep. In view of this, the progress variable cannot be explicitly resolved on the large eddy simulation mesh. Sometimes, the Thickened Flame approach is employed to resolve this issue [24].

COMPUTATIONAL METHODOLOGY
In the present study, the reactive Reynolds Averaged Navier-Stokes (RANS) model (Equations (1) to (33)) and the b -Ξ two-equation model (Equations (34) to (39)) were used to investigate the performance of a fully premixed modern, high-performance 4-valve, iso-octane, dual overhead cam (DOHC) engine with quasisymmetric pent roof combustion chamber, running at 1500 revolutions per minute. The intake valves are canted -20° while the exhaust valves are canted +22°. The evolution of turbulence inside the combustion chamber of the spark ignition engine was modeled using the standard k -ε equation augmented with the wall functions.
The flame-wrinkling model that was used to model combustion in the present study was obtained by replacing Equation (35) with the equilibrium expression Equation (37) and using Gülder's correlation, that is, Equation (39) for calculating the laminar flame speed S ub . Equations (35) and (37a and 37b) are usually referred to as the one-equation and 2-equation Weller models respectively. The thermodynamical properties were calculated using the JANAF Tables [18]. The mixture was ignited at 15° Before Top Dead Center (BTDC).
The performance parameters investigated were the pressure p, temperature T, turbulence intensity u', and the regress variable b. These parameters were investigated over the range of crank angle -180 ≤ °CA ≤ 120 and air/fuel equivalence ratio, Φ = 1.0. work needs to be carried out on this chiefly because of the enhanced removal rate caused by flame compression. Analysis using results of the spectral model will be required. c. For cases involving very lean propane flames or other flames which are very strain-sensitive, a transport equation for the laminar flame speed is required. This equation is derived using empirical arguments involving the strain time scale and the strain-rate at extinction. The transport velocity is the same as that for the Ξ equation [31].

NUMERICAL SCHEME
The governing Equations (1) to (39) were discretized both in space and time using the appropriate discretization utilities in OpenFOAM. Each of the term in these Equations is represented individually in Open-FOAM code using the classes of static functions finite Volume Method (fvm) and finite Volume Calculus (fvc) methods respectively. fvm and fvc contain static functions, representing differential operators, for example, ∇ 2 , ∇ and ∂/∂t, etc. that discretize geomet ricField<Type>s. The purpose of defining these functions within the two classes, namely, fvm and fvc, rather than one, is to distinguish: • Functions of fvm that calculate implicit derivatives of and return an fvMatrix<Type> and, • Some functions of fvc that calculate explicit derivatives and other explicit calculations, returning a geometricField<Type>.
Finite volume discretization of each term is formulated by first integrating the term over a cell volume V. Most spatial derivative terms are then converted to integrals over the cell surface S bounding the volume using Gauss's theorem: where S is the surface area vector, can represent any tensor field and the star notation * is used to represent any tensor product, such as, inner, outer and cross and the respective derivatives: divergence φ ∇ ⋅ , gradient φ ∇ and φ ∇ × . Volume and surface integrals are then linearized using appropriate schemes. The schemes which were used for each of the terms in Equations (1) to (39) are tabulated in Table 1. Comprehensive details concerning these schemes can be found in [32].

Code validation
The numerical schemes implemented in the OpenFOAM code have been extensively validated in the error estimation studies [33,34] and comparisons with analytical solutions [31].

Computational domain and mesh
The computational domain over which solutions of the differential equations described in Equations (1) to (39) were solved are those lying between the cylinder surfaces, cylinder head and piston face. In computational fluid dynamics parlance, these surfaces are referred to as walls. In internal combustion engines two of these walls move -the piston face moves continuously in time and the valves move intermittently. If a finite-difference computational mesh which is fixed in space is used, much of the grid will not be utilized during parts of the piston cycle. A mesh which does not translate with the piston motion is therefore highly wasteful of computer resources, that is, both computer processing unit's (cpu's) time and storage memory. This problem was resolved by defining a mesh which translates with the motion of the piston such that the whole mesh is always confined between the cylinder and piston head. Various ways of accomplishing this have been proposed, but the Arbitrary Lagrangian-Eulerian (ALE) technique [35] has been found to be the most suitable. This is the method that was used in the present work. The hexahedral computational mesh employed in this study was generated from the same data that was used for the KIVA-3V 4-valve model [36]. The initial computing mesh was created by the OpenFOAM's kivaToFoam grid conversion utility. The mesh at the bottom dead center (180°CA BTDC) position consists of a total number of 30742 nodes, 85742 faces, 79522 internal faces and 27544 hexahedral cells. The minimum and maximum face areas were 2.74725·10 -7 and 4.11662·10 -5 m 2 , while the minimum, maximum and total volumes are 3.07069·10 -10 , 1.36137·10 -7 and 1.31297·10 -4 m 2 respectively. The maximum skewness of the concerned face was 3.73901 which was adjudged satisfactory. It is worth mentioning that the OpenFOAM has a lot of utilities for managing the computational mesh, thereby affording the opportunity to building a more appropriate mesh and thus reducing and improving upon the discretization errors and the accuracy of the numerical results. Figures 1 and 2 show the computational mesh for the piston at the bottom dead center (BDC) and top dead center (TDC) positions respectively. The automatic mesh motion solver was used to compute the piston and valve motions. The characteristics of the computational mesh and intake swirl specification when valve geometry was not modeled but specified are shown in Tables 2  and 3.

Boundary conditions
In order to completely specify a numerical model, it is necessary to supply the conditions at the boundaries of the solution domain for the variables. These conditions are usually obtained either by specifying the value of the dependent variable at the boundary, or the value of the associated flux, or a relation between the two. The specifications of the relevant boundary conditions employed in this work also follow these basic principles. Accordingly, the boundary conditions are specified in Table 4. u Velocity field Uniform internal field = (0,0,0) Ξ Flame-wrinkling factor, S t /S u Uniform internal field = 0.0

Initial conditions
The simulations were initialized at 180-degree BTDC and the air/fuel mixture is assumed to be at rest within the entire computation domain and calculations were terminated at 120-degree ATDC. The initial inlet temperature and pressure were set to 600°K and 1.9 bar and to a fixed value of 450°K at the walls respectively. The initial velocity distribution was assumed to vary linearly between the piston and the cylinder head, being zero at the cylinder head and equal to the piston speed at the piston head. Furthermore, the no-slip boundary conditions were imposed on the walls. At the boundaries of the computational domain, heat and momentum losses were calculated by assuming that boundary layers were turbulent and matching to a lawof-the-wall velocity profile.
The initial turbulent kinetic energy was assumed to be uniform and equal to 10% of the square of the mean piston speed, while the initial turbulent length scale was assumed to be 40% of the distance from each of the cells closest to the solid boundary. At engine speed of 1500 revolutions per minute (4.1125 m/s), the initial specific turbulent kinetic energy, k = 2.661 m 2 /s 2 , while the initial rate of dissipation of turbulence kinetic energy ε =25.763 m 2 /s 2 , respectively.
Since there was no ignition initially, that is, at time, θ start = 180°CA BTDC, we set the regress variable, b = 1.0 and the flame wrinkling factor, Ξ = S t /S u = 1.0 in the whole computational domain.
The complete initial conditions governing and controlling other parameters are summarized in Table 5 below:

Solution procedure
The PIMPLE algorithm was used to solve discretized Equations (1) to (39). The PIMPLE numerical scheme is a combination of the pressure-implicit split-operator (PISO) and the semi-implicit method for pressure-linked equations (SIMPLE) algorithm. The PIMPLE algorithm is an iterative scheme which is normally used for solving the momentum and mass conservation equations. Details concerning the operation of PIMPLE can be found in [34]. However, it is worth mentioning that PIMPLE has two iterative mechanisms namely, Outer Correctors and Inner Correctors. These iterative utilities ensure stability especially when dealing with large time steps where the maximum Courant number may consistently be greater than 1 or when the nature of the solution is inherently unstable.
In the present study, the Outer and Inner Correctors were set to 2, while the Non-Orthogonal Corrector was set to 1. The stability of the algorithm was further improved by setting the relaxation factor to 0.7. The tolerance and relative tolerance values which were used to signify solution convergence are tabulated in Table 6. Relaxation factor φ 0.7 -

RESULTS AND DISCUSSIONS
The details of the engine specification, fuel and ignition characteristics that was used in the present study are provided in Tables 7 and 8 below. The finite-difference approximations to the governing Equations (1) to (39) were solved with the OpenFOAM XiEngineFoam solver for internal combustion engines.    4 show the residuals of the velocity component in the x,y and z directions, turbulent kinetic energy and the regress variable respectively. The total number of iterations at which convergence was achieved for all the variables range from 5500 to 6000. In view of Figures 3 and 4, it safe to assume that the numerical scheme was well-behaved during the computation of the properties being investigated, that is, the mass continuity equation, temperature, pressure, regression variable, etc. 796 ▪ VOL. 48, No 4, 2020 FME Transactions Figure 5 shows the variation of pressure with crank angle. As can be seen from Figure 5, the peak pressure occurred a few crank angle degree ATDC and has a value of about 30 bar. The graph of volume averaged temperature versus crank angle is shown in Figure 6. Figure 6 shows that the maximum volume averaged combustion temperature occurred at approximately 95 degree crank angle ATDC and has a volume averaged value of about 2700°K, whereas the actual peak temperature was found to be 3000°K and it occurred at grid point 12630. The other temperatures within the combustion chamber which were higher than the volume averaged peak temperatures were found to be in the range 2968.81°K to 2974.01°K and corresponds to grid points in the range 12630 to 12633. It is worth mentioning that, either advancing (that is, taking the spark initiation further from the BTDC) or retarding the start of ignition (that is, taking the spark initiation closer to TDC) both have some benefits and disadvantages [37]. Advancing initiation of ignition increases the peak pressure within the engine combustion chamber, whereas retarding the start of ignition reduces the maximum pressure. Furthermore, advancing the spark timing increases the maximum temperature, while retarding initiation of reduces the peak temperature. The aforementioned characteristics show how retarding the initiation of spark can have a beneficial effect on NOx production while simultaneously reducing the tendency for the engine to knock. However, the life span of the exhaust valve could be adversely affected because retarding the ignition timing does increase the heat release temperature, but it could be useful in helping to achieve early catalyst light-off. Most importantly, retarding the ignition timing does decrease the power output with the concomitant reduction of cycle efficiency. Figures 5 and 6 show that the locations of the peak pressure and temperature are not correlated, that is, the maximum pressure occurred at a few crank angle degree ATDC position whereas the peak temperature occurred at about 95° ATDC which corresponds to about 5 ms apart. The reason for this could be due to the fact that the computed values of temperature were volume averaged. Figure 7 shows the variation of volume averaged turbulent kinetic energy/turbulent intensity with crank angle. It is worth noting that the graph shows two turning points. The maximum turbulent kinetic energy occurred at approximately 25 degree before the top dead center position and decreased to just below 1.2 m/s at the top dead center position. This could be due to the fact that turbulence is modified by flames (recall that ignition was initiated at 15-degree crank angle BTDC) through flow acceleration and changed kinematic viscosity. It is worth mentioning, however, that modified turbulence structure does alter the flame shape by wrinkling it thereby resulting in large flame area. The wrinkled and increased flame area has the overall effect of enhancing mixing and chemical reactions through temperature fluctuations. Furthermore, it is worth noting that as the power stroke continues, the turbulent kinetic energy continues to decrease in magnitude.  Figure 8 shows the variation of the regress variable with crank angle. As expected, the regress variable maintained a constant value of 1.0 (unburnt air/fuel mixture) until about 10 degrees before the top dead center position. Thereafter, the air/fuel mixture was rapidly consumed. Figures 9 to 12 show the combustion chamber volume distributions of pressure, temperature, turbulent kinetic energy and flame wrinkling fields at various grid points and crank angle positions stated in each of the figures. Figure 9 shows the spatial distribution of pressure within the combustion chamber at the top dead center position. From the figure we can see that the pressure was almost uniform with a value of about 30 bar within the entire combustion chamber. The spatial variation of temperature versus crank angle at the TDC position is shown in Figure 10. We can see from Figure 10 that even though ignition was initiated at 15-degree BTDC, the flame kernel was not spreading as rapidly as expected. This could be attributed to the dissipation of turbulent intensity as depicted in Figure 7. Plotted in Figure 11 is the volume evolution of the turbulent kinetic energy at the TDC position. It can be observed that  This is obviously due to the fact that piston motion is creating a lot of viscous shear, that is, large eddies of the order of integral length scale, l close to the combustion chamber walls. These large eddies in turn breaks into medium and smallest-scale turbulent eddies on the order of the Kolmogorov scale which then dissipates the turbulence energy via molecular viscosity. Figure 11 shows that the values of the volume averaged turbulent kinetic energy falls in the range 5.5 · 10 -3 ≤ k ≤ 3.8 · 10 0 .

CONCLUSION
The open source Field Operation and Manipulation (OpenFOAM) software was used to investigate the performance of a fully premixed modern, high-performance 4-valve, iso-octane dual overhead cam (DOHC) engine with quasi-symmetric pent roof combustion chamber running at 1500 revolutions per minute. The numerical results were found to qualitatively reproduce the basic features that are usually observed in internal combustion engine test bed experiments. In the XiEngineFoam solver, the flame front propagation was modeled by regress variable b, with the relation b = 1c, where c is the progress variable. The flame propagation was modeled by solving transport equation of the regression variable, Equations (35) to (39).
In specific terms, the salient conclusions from this study can be summarized as follows: 1. The results of the present work were found to qualitatively reproduce the general characteristics of the pressure distribution which are similar to those usually occurring within combustion chamber of internal combustion engines. The peak pressure occurred at a few crank angle degree ATDC and has a value of about 30 bar.
2. The present study shows that the peak volume averaged temperature of about 2700°K was within the general ballpark value of temperatures that are usually prevalent within the combustion chamber of a spark ignition engine. It is worth mentioning that the combustion duration was also well within typical industry standard, however, the location of the peak temperature lagged behind that of the maximum pressure. The simulation also shows that the actual maximum temperature within the combustion chamber was 3000°K and it is located within a very small region inside the combustion chamber. 3. The flame-wrinkling factor, Ξ = S t /S u was found to be in the range 1.0 ≤ Ξ ≤ 3.8 which translates to a magnification of nearly four times that of the quiescent value. 4. The dynamics of the regress variable b was accurately predicted as shown in Figure 8.