Influence of Mesh and Combustion Parameters on a Spark Ignition Engine CFD Simulation

In the current paper, a CFD simulation study on a pent-roof, spark ignition engine was performed. The scope was to analyze the basic parameters that influence the performances of the simulation: mesh configuration, mesh refinement, turbulence model, and flame wrinkling model parameters. It was used as the OpenFOAM solver for spark ignition engines. The solver was modified by adding the layer addition/removal mechanism for the mesh. Further, a Zeldovich extended NOx model was implemented. Four different mesh configurations were tested with different refinement levels. A set of different turbulence models and initial parameters and different flame wrinkling model parameters were tested. It was observed how all these parameters influence the rate of heat release, cylinder pressure, temperature, the level of turbulent kinetic energy, the flame wrinkling factor, and finally the resulting concentration of nitride monoxide.


INTRODUCTION
The internal combustion engine has been the primary power plant that provides power to our transport means for more than a century. In the last decades, environmental problems, partially related to internal combustion engines are gaining attention. In the first place, this is the carbon dioxide emission that is responsible for causing the greenhouse effect and the related global warming and climate change issues. Further, there are other emissions, like nitrogen oxides, unburned hydrocarbons, carbon monoxide, and particulate matter. Hence, legislation pushes ever more stringent emission limits for internal combustion engines. Some governments announced the complete replacement of internal combustion engine-powered cars with other solutions, like electric or hydrogen-powered cars. However, there are still many limitations of these technologies and it seems that the internal combustion engine is going to be developed further parallel to these solutions.
Some of the combinations that are planned are hybrid systems where internal combustion engines are paired with electric motors. Other arrangements involve the use of alternative fuels which allow a lower carbon footprint, such as biofuels (bio-methanol, bio-ethanol, and biodiesel). Ammonia combustion is also in the focus as a possible, carbon-free fuel. Hydrogen combustion is also being tested in internal combustion engines, pure or blended with more conventional fuels. Another solution is the artificial synthesis of hydrocarbon fuels from captured carbon dioxide and green hydrogen known as e-fuels. In the development process, various blends among alternative and fossil fuels are being tested. Also, alternative combustion modes are being developed such as dual fuel engines, HCCI combustion, and similar. Gas after treatment is another field that is being developed in order to reduce combustion-related problems. However, it is hard to imagine a future without the internal combustion engine, especially in the marine and commercial environment. Hence, the development of IC engines remains an important topic of scientific investigations.
Engine development is performed via experimental studies where the engine is set on a test bench with the crankshaft connected to a brake that simulates the working torque. The engine is connected to the fuel, cooling, and exhaust gas systems which ensure its operation. Many parameters are measured, such as power, torque, fuel consumption, exhaust gas temperature, and composition. In-cylinder pressure is another important parameter together with components' mechanical and thermal stresses. Such laboratory tests are very complex and costly, but important.
A cheaper and faster alternative is the computer simulation which also gives additional insight into the IC engine's internal processes. There are various kinds of simulations, but the computational fluid dynamics (CFD) simulation can give insight into the engine cylinder processes, responsible for chemical reactions and pollutant formation. There are different commercial software that allows CFD analysis of engine processes. Some open-source projects offer CFD analysis.
OpenFOAM is an open-source, free-of-charge, CFD collection of solvers and libraries written in C++ [1][2][3]. Of its availability, flexibility, and collaboration between users, it is a good platform to use to approach uncommon CFD problems. The core of the spark ignition engine combustion calculation is the flame speed. Amirante [4] provides simple and workable expressions, suitable for spark-ignition engine simulations, that allow laminar flame speed calculations of some practical fuels. Flame propagation depends on the wrinkling process. Chatakonda et al. [5] proposed a new flame wrinkling model to apply to premixed flames where the turbulence scales can penetrate the flame. The model utilizes concepts of fractal geometry and Damköhler's limiting flame-speed scalings. Anetor et al. [6] performed a CFD investigation of a modern, high-performance 4-valve engine with a quasi-symmetric pent-roof combustion chamber using OpenFOAM. Sendyka and Noga [7] performed a combined theoretical, experimental, and numerical analysis of a dual injection, pentroof, spark ignition engine. They focused on the influence of the combination of port injection and direct injection on engine performance and emissions. Gundmalm [8] in his thesis performed a detailed CFD simulation using OpenFOAM, including the intake, compression, and combustion phases. The focus was on the investigation of the presence of an intake air tumble motion and later on the tendency to knocking phenomenon. Battistoni et al [9] simulated pure gasoline and pure ethanol combustion in a 4-valve singlecylinder SI engine. They assumed premixed combustion conditions and used detailed chemistry with the wellmixed combustion model approach. Simulation results are compared with data from optical engine experiments. Mashkour [10] and al. developed a 3D CFD simulation of a SI, pent-roof engine. They analyzed the influence of the equivalence ratio on cylinder pressure, temperature, heat transfer, and NOx emissions. Yin et al. [11] analyzed the combined effects of piston bowl geometry and a charge motion control valve (CMCV) on tumble flow and combustion features in GDI engines, using the commercial CFD software FIRE. In their paper, Mauro et al. [12] used a single-zone model and a CFD model. They were implemented for the simulation of the closed part of the cycle of a CFR engine in order to obtain an accurate heat release evaluation. The local turbulent flame development is modeled using the G-equation, which provides a strict correlation to the laminar flame speed which, in turn, is a chemical property of the gas mixture. Kutkan et al. [13] used the Algebraic Flame Surface Wrinkling model to simulate a high-pressure, turbulent premixed Bunsen flame experiment. They compared the results obtained with OpenFOAM and Ansys Fluent to the experiment. Yoshy et al. [14] used OpenFOAM to analyze the effects of the spark plug position on engine performance.
Here, an attempt to get a systematic insight into the influence of the mesh, turbulence, and combustion model parameters on calculated performance is done. The focus is on the heat release rate, cylinder pressure, and NO x emissions. The analysis can be used as guidance for further spark ignition engine simulations.
The paper is organized as follows. In section 2, some parts of the mathematical model behind the OpenFOAM software for spark ignition engine simulation are presented: the transport equations, combustion modeling, turbulence modeling, and NOx calculation. In section 3, some details about the simulated engine and the calculation boundary conditions are reported. Then the results with comments follow the selection of the mesh configuration, the analysis of the number of cells, the selection of the turbulence model and its parameters, and the influence of combustion model parameters. Finally, some conclusions are derived in the final section.

Governing equations
In the paper, the CFD code OpenFOAM is used. It uses the finite volume method in order to calculate the variation of the main variables in the fluid continuum. The main governing equations are described in the following part. The global continuity equation yields where ρ is the total mass density of the fluid, u is the gas velocity vector and t is the time.
The averaged momentum equation for the gas reads where p is the gas pressure, σ is the viscous stress tensor, g is the specific body force, which is assumed to be constant. The energy equation yields where e is the specific internal energy, excluding the chemical energy, and J is the heat flux vector which is the sum of heat conduction and enthalpy diffusion. Q c is the source term due to chemical heat release.

Spark ignition engine combustion modeling
There are several strategies to cope with the turbulent premixed combustion that takes place in spark ignition engines. Some of them are the Eddie Break Up [15], Thin Wrinkled Flame [16], and Flame Area Evolution [17]. The Weller flame wrinkling model implemented in OpenFOAM is based on the solution of two equations [18]. One is the transport equation for the spatial and temporal distribution of the progress variable b and the other solves the wrinkling factor Θ. The flame wrinkling is calculated by simple algebraic equations based on flame speed correlation: where  is the ratio of the 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, u' is the subgrid turbulence intensity, R η is the Kolmogorov Reynolds number. The flame front propagation is calculated by the following transport equation for the density-weighted mean reaction regress variable b: For a more detailed description of the combustion modeling [6] can be consulted.

Turbulence modeling
The most common model used to simulate the effects of turbulence in the flows is the k-ε model which solves two transport equations, one equation to solve the turbulent kinetic energy k and the other to calculate its dissipation rate, ε. A more detailed description and the physical meaning of the various terms can be found in [6]. The RNG model was developed using Re-Normalization Group (RNG) methods to account for the effects of smaller scales of motion [19]. In the standard k- model, the eddy viscosity is determined from a single turbulence length scale, so the calculated turbulent diffusion is that which occurs only at the specified scale, whereas, in reality, all scales of motion will contribute to the turbulent diffusion. The RNG approach, which is a mathematical technique that can be used to derive a turbulence model similar to the k-, results in a modified form of the  equation which attempts to account for the different scales of motion through changes to the production term.
The k-ω is a similar two-equation model. Instead of  there is the ω equation, where ω is the specific rate of dissipation of the turbulent kinetic energy into internal thermal energy [20]. The SST k-ω model is another common turbulence model implemented in OpenFOAM that can be used for internal combustion engine simulations. This is a hybrid turbulence model that combines aspects of the k-ε and k-ω models and is designed to be more accurate and robust for a wider range of flow regimes. It is often used in complex geometries or high Reynolds number flows.

NO x model
Among the products of the combustion in spark ignition engines, Nitrogen Oxides are one of the most harmful, and the one that can be controlled to a certain extent by engine parameters and settings. It is indicated as NO x since nitrogen and oxygen combine in various atomic ratios, but all of them are more or less harmful to human health. If nitrogen is not present in the fuel, it comes to the reactions from the air. The quantification of the NO x is usually done using the "Zeldovich extended mechanism" proposed by Zeldovich in 1946 [21]. It is based on three chemical reactions: Initially, nitrogen monoxide NO is formed during combustion in high-temperature regions (starting at about 2000 K), mostly from the nitrogen contained in the air. Afterward, some nitrogen monoxide NO will convert to nitrogen dioxide (NO 2 ) and nitrous oxide (N 2 O). The mix of oxides of nitrogen is called NO x . The extended Zeldovich model was implemented in the code used for the paper. The lack of chemical species tracking of the simple combustion model used was solved by the assumption of a known concentration of relevant species in fresh air-fuel mixture and fully combusted products.
The transition from the concentrations in the fresh mixture and products was assumed as linearly proportional to the combustion progress variable b and local density .
The concentration of the remaining species needed for the extended Zeldovich model is calculated as follows [22]:

THE ENGINE
The engine used for the simulations in this article is an experimental single-cylinder, optical engine with a pentroof combustion chamber and four valves described in [9]. The basic data about the engine are reported in Table 1. It can feature both direct injection and port fuel injection and is equipped with a camless electrohydraulic valve lifting system that allows free control of valve timing. Optical access to the combustion chamber is possible through a 60 mm diameter quartz piston crown window. Details about the pent-roof, retrieved from images in [9], are reported in Figure 1., with dimensions in millimeters.

SIMULATIONS AND RESULTS
A series of simulations was performed to test the influence of the main calculation parameters like the discretization mesh type, number of cells, turbulence model and its initial values, and combustion model parameters. The validation was not carried out separately. All the simulation sets are compared to the experimental result. Specifically, the measured parameter was the cylinder pressure found in [9]. The main simu-lation initial and boundary conditions are presented in Table 2.

Mesh configuration
The mesh for internal combustion engines can be complex. It includes at least a moving piston and it can include a piston bowl, a pent-roof head, and valve recesses. In order to cope with the moving piston, the layer addition/removal strategy is implemented in the engine foam solver to reduce the number of cells during the calculation and to maintain the aspect factor of the cells [23]. As the piston travels toward the TDC (Top Dead Course), layers of cells are removed, and they are added as the piston travels to the BDC (Bottom Dead Course). For some simulations, the influence of the piston and head details can be neglected. Some papers used the geometry supplied with the OpenFOAM tutorial, kivaTest like [6] and [14]. In the first part of the investigation, 4 meshes for the same engine are built with various levels of refinement in order to check how the mesh influences the basic combustion performance. The meshes at the end of the compression are represented in Figure 2, with increasing refinement from a) to d). The first mesh is a 2D mesh that does not reproduce any of the geometric details seen in Figure 1.
In [24] it is shown that very simple 2D meshes can result in realistic results. The second mesh from the left is a 90° sector mesh, that represents 1/4 of the real domain, and also neglects many of the details. It is a 3D mesh with cyclic boundaries. The third mesh reproduces half of the domain with symmetry boundaries. Even if it is somewhat simplified compared to real geometry, it reproduces the basic shape of the pentroof. The last mesh is a full 3D mesh of the domain that reproduces all the minor details of the engine geometry, the pentroof, and the recession in the piston. The first three meshes were created with blockMesh, the OpenFoam basic meshes. The last one is created with Ansys ICEM CFD commercial software and then converted to OpenFOAM format. The meshes have a different number of cells.
The 2D mesh has 488 cells at the beginning of the compression and, due to the layer-removing mechanism, 228 at the end of the compression. The ¼ sector mesh has from 5856 cells at the beginning to 2736 cells at the end of the compression, the half mesh has 10560 to 4320 cells and the full mesh has 84325 to 16845 cells.

/4 sector mesh, c) half mesh, d) full mesh
In Figure 3. the influence of mesh configuration on the cylinder pressure is presented. It can be seen that all 4 tested configurations catch the general trend of the pressure. The pressure gradually builds up after a certain period after the ignition (at -17 ºCA). The maximum pressure obtain differs a lot, from 17 bar for the 2D mesh to 48 bar for the 1/4 sector mesh. The halfmesh and full-mesh configurations give more reasonable maximal pressures that are closer to the 33 bar obtained experimentally in [9]. The pressure peak of these cases comes a few degrees of crank angle too late compared to the experiment. It could be further tuned by modifying surface temperatures or by tuning the heat transfer model, but it was not the scope of the paper. In Figure 4. the heat release rate for the different mesh configurations is shown. The heat release is intense and comes early for the 1/4 sector mesh. It comes later for the half-mesh and full-mesh configurations and the heat release is more spread. For the 2D mesh, the peak comes early, but the heat release is much lower than in other cases.

Figure 4. Mesh configuration influence on the rate of heat release
In Figure 5. the temperature for the different mesh configurations is displayed. The temperature reaches about 2700 K for all the cases, except for the 2D mesh for which the temperature reached is unrealistically low. The temperature peak is reached slightly earlier for the 1/4 sector mesh than for the other configurations.  Figure 6, the level of turbulent kinetic energy is presented. It varies from 1.4 at the beginning to 0.4 at the end of the simulated period with a characteristic trend. The case 1/4 sector mesh presents a jump of turbulent kinetic energy that follows the fast combustion. The case 2D mesh is not presented since its turbulent kinetic energy is out of the range and jumps to 15 m 2 /s 2 from the beginning.

Mesh refinement
In the second set of simulations, the mesh refinement, or cell number and size was analyzed. It selected the halfmesh configuration with different cell numbers. Since the layer addition/removal function is implemented in the code, the initial number of cells is reduced as compression is going on. The number of the cells in the meshes was, for the first, coarse mesh, from 5120 at the beginning (BDC) to 1920 cells at TDC, for the second, intermediate mesh it was from 15060 to 5220 cells and for the third, finest mesh it was from 31110 to 11730 cells. In the following figures, the influence of the number of cells on the results is presented. The meshes are named on the initial number of cells.
In Figure 7. it can be seen that the maximal cylinder pressure is higher for the coarser mesh and the pressure rise starts earlier. The cylinder pressure reaches 35 bar, which is slightly more than in the experimental case. In the same case, the pressure peak is reached a few crank angles later than in the experimental case. The reason behind the pressure curves can be seen in the next figure.

Figure 7. The number of cells influences cylinder pressure
In Figure 8. the heat release over the crank angle is presented. It can be seen that the coarser mesh gives origin to earlier combustion with a faster reach of the peak in heat release. For the other, finer meshes, the heat release starts later and is more gradual. The heat release is reflected in the cylinder pressure and in the cylinder temperature.
As can be seen in Figure 9. the temperatures are slightly higher in the early phase for the coarser mesh, while the temperatures are lower in the later phase. In Figure 10., the influence of cell size on the turbulent kinetic energy is shown. It can be seen that the influence is very weak, but generally smaller cells result in smaller Turbulent Kinetic Energy. In Figure 11., the influence of mesh refinement on the flame wrinkling factor is presented. The coarser mesh gives origin to an earlier increase of the flame wrinkling factor, which then causes the turbulent flame speed to increase. This is the reason for the earlier build-up of heat release speed and the higher cylinder pressure. In Figure 12., the influence of the mesh on the NO x concentration is shown. It can be seen, the characteristic pattern of NO x formation during the most intensive combustion period and a small decrease after which the concentration remains constant. It is surprising that the finest mesh gives origin to the greatest, and much higher concentration of NO x . The obtained NO x concentration cannot be used for absolute formed NO x quantities evaluation, but the trends in various influential factors can be applied.

Turbulence model
Here, a preliminary test of the described turbulence models is performed. The initial parameters were: k=3 m 2 /s 2 for all the cases and ε = 250 m 2 /s 3 for both the k-ε variants and ω=5 s -1 for the SST k-ω case.  Figure 13. shows that the k- models (standard and RNG) predict cylinder pressure better while the k- SST tends to over-predict it since pressure rise starts earlier than for the measured case. This is the result of an earlier start and more intensive heat release obtained with the k- model, Figure 14. This behavior is reflected also by cylinder temperature, Figure 15. The turbulence model influences directly the level of turbulent kinetic energy, Figure 16. It can be seen that the used k- SST dissipates less the turbulent kinetic energy in the compression phase which leads to higher levels of turbulence during combustion. The flame wrinkling factor is directly influenced by turbulence, Figure 17. While the k- SST model causes an earlier rise of , the highest levels are obtained with the standard k- model. Very interesting is the influence of the used turbulence model on NO x concentration, Figure 18. While the k- SST and RNG k- result in similar NO x concentrations, the standard k- results in a roughly 30% higher final value. Generally, engine CFD simulations are performed for various sizes and types of engines. Different calculation meshes are used and also different sub-models. The level of turbulent kinetic energy in the cylinder directly influences the flame speed propagation, which is responsible for pressure and temperature build-up and hence also for pollutant formation. The level of cylinder turbulent kinetic energy is not readily measured, hence it is not simple to choose the initial values for the turbulence model. In the OpenFOAM engineFoam solver tutorial named KivaTest, initial k = 4, and initial ε = 450 are used while in [6] initial k = 2.661 and initial ε = 25.763 are used. [8] used k = 4 and ε = 450 for the engine at the speed of 8000 rpm. Here, an attempt to analyze the influence of the turbulence model's initial values is performed. The level of turbulent kinetic energy, k, was varied from 2 to 5 m 2 /s 2 , and the amount of the turbulent kinetic energy dissipation rate ε was varied from 25 to 500 m 2 /s 3 .
In Figure 19. the cylinder's turbulent kinetic energy is shown. It can be seen that both k and ε have an important influence on turbulence levels. The highest levels of turbulence are found for the case with an average initial level (k = 3), but a low TKE dissipation rate of ε = 25, and for the high initial turbulence (k = 5) case. Lower than average turbulence levels during the process are found for low initial turbulence (k = 2) or high turbulent energy dissipation rate cases (ε = 500). This is clearly reflected in the following figure, Figure  20. which shows the flame wrinkling factor . The cases with higher turbulence levels are from Figure 19. result in higher wrinkling factors. In Figure  21. it can be seen that the higher wrinkling factor leads to faster heat release rates dQ/dφ.
In Figure 22. the cy-linder pressure is presented. The calculated results can be compared with the experimentally obtained ones. It can be concluded that the case for average initial values of k = 3 and ε = 250 result in cylinder pressure levels that are close to those of the experimental pressure. Besides pressure, the heat release rate is reflected in the cylinder's average temperature, Figure 23. The maximal temperature reached is very similar for all the cases. The cases with higher turbulence level results in an earlier temperature build-up, and slightly higher temperatures. The temperature is the main influential factor for the NO x concentration. In Figure  24. it can be seen that the higher temperatures of the cases with higher cylinder turbulence levels result in higher NO x concentration.

Combustion model parameters u' and 
Engine power, efficiency, and emissions performance are highly dependent on the combustion rate. The combustion rate depends on temperature, pressure, mixture composition, and other parameters. In the simulation, the combustion rate is modeled as described in the mathematical model. Among other parameters, u' and  coef are the most influential model parameters. Here, an attempt to gain an insight into the nature of their influence is performed. u' was varied from 0.5 to 2.4 and  coef from 0.1 to 1.  On the other hand, the  value shows a very low and late increase from the initial value of 1 for the combination of small values of u' and  coef . In Figure 26. the influence of the flame wrinkling factor  on the heat release rate is presented. It can be seen that the combustion speed is greatest for the case with u' = 2.4 and  coef = 1.6. The combination of u'=1 and  coef = 1 gives also a relatively fast combustion, while other combinations of parameters lead to a late and slow combustion result. The same trends are reflected also in Figure 27. which shows the cylinder pressure with the experimentally obtained result. It can be observed that the high value of u '= 2.4 leads to an overestimation of cylinder pressure. Very good results are obtained for the intermediate values of u' = 1.6 with  coef = 0.3 and u' = 1 with  coef =1. For lower values of the model coefficients, the cylinder pressure results are too low. In Figure 28., the influence of the combustion parameters on the turbulent kinetic energy is presented. It can be seen that it rises proportionally more for the combinations of parameters that lead to more intensive combustion. In Figure 29. the influence of the flame wrinkling parameters on the cylinder temperature are shown. It can be seen that the parameters that increase the flame wrinkling mostly lead to higher temperatures and earlier temperature peaks. The calculated concentration of NO x is reported in Figure 30. It can be seen that the highest concentration by far is obtained for the case with u' = 1 and  coef = 1. As expected, for the cases that resulted in low cylinder temperatures, NO x concentrations are the lowest. The final NO x concentration is lower for the highest temperature case (u' = 2.4,  coef = 0.3) than for the case that underestimated the cylinder pressure and temperature (u '= 1,  coef = 0.3).  Figure 31. a graphical presentation of the flame propagation in the cylinder with time is presented. Speci-fically, the regression variable b is presented for the moments of -10, 0, and 10 CA.

CONCLUSION
In the present paper, the basic parameters of a CFD simulation of a spark ignition, internal combustion engine are analyzed.
The simulation mesh that reproduces half of the domain resulted in the best compromise for accuracy, simplicity, and calculation time. The simpler meshes, a 2D and a 1/4 sector mesh showed excessively fast or slow combustion compared to the experimental results. The more complex, full 3D mesh did not behave better than the half mesh, it was more complex to build and lead to a slower calculation.
By increasing the mesh refinement (number of cells in the domain, smaller cells) with fixed all other parameters, the flame propagation speed is dampened. Hence, heat release rate and maximum pressure, and temperature are decreased. The influence on the turbulence kinetic energy and the flame wrinkling factor is smaller and not so direct. The implemented NOx model is sensitive to mesh refinement and resulted in a big increase in concentration for the finest mesh tested. The mesh that yields 15060 cells at the beginning and 5220 cells around TDC was selected as the best compromise for this size of the engine.
The RNG k-ε model showed to be the most appropriate. The tested k-ω over estimated cylinder peak pressure while the standard k-ε resulted in excessive NOx production. Clearly, it is possible that the results could be tuned by choosing different initial values for turbulence and turbulence dissipation, or different model parameters. The initial parameters influence was investigated for the RNG k-ε model. The turbulence level is recalculated after the internal velocity conditions, but initial values do influence all the other observed parameters. The turbulence kinetic energy can be reduced by choosing a low initial value or a high initial value of turbulent kinetic energy dissipation. The best results are obtained for the case with initial values of k =3 and ε=250. The influence of the turbulence model's initial parameters on NOx emissions is low. Finally, the influence of the flame wrinkling model parameters u' and Xi coef are tested. They can be used to tune the flame speed to a certain level. The heat release rate can be boosted by increasing any of the two coefficients. Best results are obtained for the combination of u' = 1.6 and Xi coef = 0.3. The case with 1 for both the parameters gave also very good results, but unexpectedly raised a lot the concentration of calculated NOx.
The simple extended Zeldovich-based NOx model results in concentrations of the same order of magnitude as in other similar investigations. It is known that it can be tuned with different model coefficients, but that is not the point. It can be used to observe the influence that some working parameters of the engine have on the final NOx emissions.