Investigation of pollutants formation in a diesel engine using numerical simulation

The current study aims at simulating the in-cylinder combustion process in a diesel engine and investigating the engine performance and pollutant formation. The combustion simulation is performed on a 3D sector employing appropriate models for various physical and chemical processes contributing in the combustion phenomenon. The overall model includes Transition SST turbulence model, eddy dissipation model for turbulence chemistry interaction, Moss–Brookes model for soot calculation and Zeldovich mechanism for NO production other than the usual transport equations. The numerical solutions are based on the finite volume discretization of the governing partial differential equations. Engine performance has been studied in terms of pressure, temperature and heat release rate while the pollutants formation has been investigated in terms of soot and thermal NO production. The results show that the ignition delay is quite short and that the injection timing may be successfully employed to control the combustion behavior. The simulation results are quite consistent with the expected behavior of the target variables indicating that the CFD analysis can be successfully employed in the diesel engine design. The results validation may be acknowledged in view of the mesh independence test, literature comparison and justification of the models.


Introduction
The combustion in diesel engines is governed by complex processes, involving various inter-connected parameters influencing them such as the turbulence, swirl, fuel injection time, injection duration, spray orientation, size of the injector hole, etc. Computational fluid dynamics (CFD) has emerged as a powerful tool to resolve the inter-connectivity of these parameters and to have an understanding of the impact of each parameter separately on the combustion. The impact of these parameters on the combustion is quantified by observing their influence on the in-cylinder pressure, temperature, ignition delay, heat release rate, pollutant formation, etc.
In recent years, the demand of diesel engines is increasing due to its recognized thermal efficiency. Inspite of their higher thermal efficiency, the diesel engines are responsible for emitting more pollutants, e.g., soot and nitric oxide (NO) because of lack of the ample time required for mixing air and fuel uniformly before the ignition. These pollutants are extremely harmful to the lives and the environment and it poses a very serious problem. Improving engine design can reduce pollutants, but it is a challenging task. It requires the utmost attention and strenuous efforts of researchers and engineers to manufacture efficient and cheaper eco-friendly engines. The task of reducing pollution may be accomplished by improving the engine design to achieve complete and perfect combustion which is highly dependent on the uniformity of the air-fuel mixture. Improving engine design involves investigating the design of air inlet and outlet passages, injector nozzle, combustion chamber and fuel blends [1][2][3][4][5][6][7][8]. CFD techniques may be used to understand the in-cylinder combustion characteristics which lead to improving engine design to reduce the pollutants. Jafarmadar et al. [9] investigated the pollutants formation for various combustion chamber designs using CFD techniques. Sakata et al. [10] improved the engine design by introducing controlled fuel injection. Konno [11] noticed a significant reduction in pollutants formation by generating delayed turbulence. The late turbulence increased the air interaction with unburned fuel resulting in increased engine efficiency. Lin et al. [12] performed numerical simulations for various combustion chamber designs and investigated one which produces less amount of pollutants. Kim et al. [2] studied the impact of fuel spray direction on the combustion and emission characteristics. Wei et al. [13] examined the impact of nozzle angle on the uniformity of air-fuel mixture and its combustion by numerical simulations. Caroline et al. [14] and Shi and Reitz [15] investigated optimal piston bowl geometries by using a genetic algorithm. They used Bezier curves to design the bowl. Ikegami [16] pointed out that the spray and gas motions should be best adapted to ensure better air utilization. The turbulent mixing process is one of the essential factors that govern the progress of combustion and also pollutant formation. Enhanced turbulence is likely to promote mixing between fuel-rich and fuel-lean fluid particles thereby accelerating combustion.
In the current study, we perform combustion simulation to investigate and discuss the in-cylinder combustion characteristics, the pressure, temperature, heat release rate together with the soot and NO formation.

Problem statement
In this study, the engine under consideration is a 12 cylinder V-type engine. The bore and stroke of the engine are respectively 150 mm and 180 mm whose compression ratio is 13.1 and displacement volume of each cylinder is 3.18986 L. The fuel considered for the simulation is n-decane C 10 H 22 which appears as a major component of the diesel fuel and may be regarded as the representative of the diesel fuel. The ignition temperature of n-decane is 483.1 K. The injector of the engine has 6 holes distributed uniformly on the sac. It is assumed that the plumes emanating these holes are identical. Thus, we have a freedom of choice of one-sixth of the piston-cylinder assembly which leads to a reduction in the computational time and resources. The amount of fuel to be injected through each hole is 29.58 mg/cycle with injection pressure being 120 MPa. The orifice diameter of each hole is 0.254 mm. The spray emanates from the point with coordinates x = 0, y = 0, z = 2e −5 m and is directed in the chamber at an angle of 70 • with the cylinder axis. The injection starts at 691 • crank angle (CA) and ends up after 20 • CA. The combustion simulations on a sector expands from the closure of the intake valves to the opening of the exhaust valves. Thus this is a period during which all the valves are closed and so the intake and exhaust ports and the corresponding valves are not included in the computational domain. In the current study, this period expands from 570 • CA to 833 • CA. The in-cylinder flow is initialized by inducing a swirl of intensity 1.3 at the engine speed being 2000 rpm. Complete engine specifications are given in Table 1.

Geometric model
The three-dimensional geometric model of the engine is one-sixth of the engine cylinder which is a three dimensional sector as shown in Fig. 1 and may be defined as follows.
1. The origin is at the top of the cylinder centerline which is parallel to the y-axis so that the cylinder head lies in the zx-plane. This plane may also be regarded as the rθ-plane. Positive y-direction points upward being outward normal to the cylinder head and the negative y-direction points into the cylinder. The lengths y and r are measured in millimeters and the angle θ is measured in radians. 2. The three dimensional geometric model of the engine is a 3D sector of the pistoncylinder assembly with sector angle 60 • . This assembly comprises 60 • sector of the flat circular cylinder head, 60 • sector of the cylinder body, and 60 • sector of the piston crown lying at the bottom of the cylinder and the crevice region lying between the cylinder wall and the piston as shown in Fig. 1. The boundaries of this assembly may be expressed as follows.
a. Cylinder head:

Mathematical model
Mathematical model contains the following sub models.

Flow model
Mathematical model, that governs the flow, comprise In Eq. (2), it is assumed that S i = 0. i.e., there is no extra source contributing to govern the flow. In Eq. (3), k 1 stands for thermal conductivity. These equations are respectively known as continuity, momentum or momentum transport and energy equations. The general transport equation for any specie φ is given by

Turbulence model
It has been established that there is no single turbulence model suitable for every type of turbulent flow. The k − turbulence model works well in free-stream regions of flow which fails near the wall boundaries while k − ω model behaves quite well in rather opposite conditions like in the near-wall regions and is not suitable for free-stream flows. This situation leads to the design of a hybrid model, like the SST k − ω model which is switched to the k − model in the free stream regions and is switched to the k − ω model near the wall regions to capture the wall effects precisely [17]. In our study, firstly we tested the realizable k − model which failed. We did not test the k − ω model in view of the above arguments. Secondly we tested the SST k − ω model which solves the following two transport equations, Eqs. (5) and (6), ∂ρ ω ∂t This model worked well. Thirdly we tested the Transition SST model which comprises the above SST k − ω model and two additional equations one for the intermittency γ given by Eq. (7) and the other for the transition momentum thickness Reynolds number Rẽ θ t given by Eq. (8). This model was developed by Menter et al. [18]. ∂ρ This model did not give any significant improvement over the SST k − ω model and the results obtained by the SST k − ω and Transition SST models were quite comparable. Therefore, either of these two models could be used. However, we decided to use the 4-equation model namely the Transition SST model for more reliability, robustness and accuracy.

Chemical species transport
The transport equation for species produced due to chemical reaction is ∂ ∂t

Spray breakup model
For spray breakup, we use Kelvin-Helmholtz and Rayleigh-Taylor (KHRT) model designed by Beale and Reitz [19]. This model is a modified version of the KH and RT models and consists of two steps, the primary breakup and the secondary breakup. The KH model is used to predict the primary breakup of the intact liquid core of a diesel jet while to predict the secondary breakup of the individual drops, the KH model is used in conjunction with the RT model. The model equations are .

Turbulence-chemistry interaction
In non-premixed flames, turbulence is mainly responsible for the air-fuel mixture preparation. The rate at which this mixture is prepared dictates the rate at which combustion will take place. This rate of mixing of the reactants is determined by the characteristics of turbulence and is much less than the rate at which reactions take place. Therefore, it is quite reasonable to assume that the combustion is limited by the rate at which the air-fuel mixture is prepared and not by the rate of reactions determined kinetically. Similarly, in premixed flames, the rate of mixing of the cold reactants and hot products is very much less than the rate of reactions. These are the cases when the combustion is mixing-limited and so the chemical kinetic rates can be neglected. Eddy dissipation model, designed by Magnussen and Hjertager [20], is based on the eddy breakup model of Spalding [21]. The eddy dissipation model uses the rate of dissipation of eddies containing reactants and products to determine the reaction rate as the rate of eddy dissipation determines the rate of mixing. In this way, the model works under the assumption that the reactions are mixing controlled as only mixed fuel and oxidizer are burnt. Therefore, for the problem under consideration, this model is more suitable for the turbulence-chemistry interaction. In this model, the source term R i,r for i-th specie from the reaction r is where and In the above equations, N denotes the number of species, Y R and Y P denote respectively the mass fractions of reactant R and product P, A is a constant whose value is equal to 4 and B is a constant whose value is equal to 5.

NO model
Oxides of nitrogen are mainly formed in the high-temperature burned gases, i.e., products of combustion, through chemical reactions involving nitrogen and oxygen atoms and molecules, which do not reach chemical equilibrium [22]. The main constituent of nitric oxides (NOx) emissions from diesel engines is nitric oxide NO or thermal NO. This thermal mechanism requires a temperature above 1800 K for its activation [23]. The species transport equation for thermal NO is given by In the above equation, Y NO denotes the gas phase NO mass fraction , D denotes the effective diffusion coefficient. The source term of thermal NO, the S NO , is determined from the extended Zeldovich mechanism [24]. The principal reactions governing the formation of thermal NO from molecular nitrogen are as follows:

Soot model
In the current study, we use Moss-Brookes soot model to determine the soot production during the combustion process. This model was designed by Brookes and Moss [25]. It has less empiricism and theoretically, it provides better accuracy as compared to the existing one-step Khan and Greeves [26] and two-step Tesner [27] models. This model solves transport equations for the soot mass fraction and the normalized radical nuclei concentration. These equations are In the above equations Y soot represents the mass fraction of soot, σ soot denotes the turbulent Prandtl number for soot, M denotes the mass concentration of soot, σ nuc denotes the turbulent Prandtl number for nuclei transport, b * nuc denotes the normalized radical nuclei concentration.
In our study, we tested the one-step and two-step models, but these did not work while the Moss-Brookes model worked well.

Solution procedure
The underlying combustion simulation is performed using the IC Engine tool which is included in ANSYS Workbench 15.0. All the tasks of simulation are completed using this tool. The solution procedure contains the steps described in the following sections.

Discretization of governing equations
The finite volume method is used for numerically solving the governing equations. The model equations are discretized on each finite volume turning them into a system of algebraic equations, which are then solved by the multigrid method.

Meshing
As the piston curve is highly non-linear, so to accommodate the piston boundary for appropriate meshing, the computational domain is divided into very few zones. Mesh views in the plane passing through the middle of the sector, at three different positions of the moving piston, are given in Fig. 2. This figure explores that structured mesh is occupying most of the region of the computational domain. In the zone near the piston, there are some cells other than hexahedrons for better adoption of the piston crown.

Mesh independence study
In order to investigate the optimum mesh resolution, we considered four mesh densities as shown in Table 2 and employed improvement in average values of the soot mass fraction and NO mass fraction in the combustion chamber as mesh resolution criteria. We note that the soot and NO mass fractions show converging behavior with convergence from below at every crank angle as the mesh is refined successively from Mesh 1 to Mesh 2 and then to Mesh 3 as may be seen in Fig. 3a, b. As the mesh is further refined to Mesh 4, the values of soot and NO decrease reflecting an increase in the error due to possibly cancelation error caused by the over-refinement of the mesh in Mesh 4. Therefore, Mesh 3

Initial and boundary conditions
In the current study, dry air is considered as the oxidizer. The ideal gas law is used for calculating its density. All the cells of the domain are initialized uniformly taking swirl number as 1.3, gauge pressure as 0.321259 MPa, x, y and z velocities as 0 m/s, temperature as 404 K, turbulent kinetic energy as 1 m 2 /s 2 , turbulent dissipation rate as 1 m 2 /s 3 , ndecane C 10 H 22 as 0, oxygen O 2 as 0.232 of the air, carbon dioxide CO 2 as 0.00046 of the air and water vapors H 2 O as 5e−7 of the air. There are five boundaries of the computational domain. The cylinder top face is assigned the temperature of 602 K, the wall of the cylinder is assigned the temperature of 567 K, and the face of the piston is assigned the temperature of 645 K. The remaining two vertical faces are the periodic faces.

Results and discussion
As stated earlier, our objective is to study the mechanical and combustion performance of the engine under consideration. Mechanical performance is determined by the average in-cylinder pressure variations while the combustion performance may be determined by the in-cylinder average temperature variation, heat release rate, soot formation, NO production, etc. The simulation of in-cylinder combustion starts at 570 • CA which is 150 • CA before TDC and ends at 833 • CA which is 113 • CA after TDC. Fuel is injected at 691 • CA, which is 29 • CA before TDC, for a duration of 20 • CA. The injection profiles of the mass flow rate and the velocity of the fuel are shown in Fig. 4. Figure 5a shows volume averaged static pressure plotted against the crank angle. We note that the pressure starting with the initial value of 0.3 MPa at 570 • CA increases slowly with the compression till the crank angle 650 • is reached. This is followed by a sharp rise of the pressure making the curves highly steep at the crank angle 691 • at which the fuel is to be injected. At this crank angle, the pressure becomes around 5 MPa which is almost 17 times the initial pressure. After the fuel injection is started, the pressure curve becomes steeper and steeper and within a duration of 30 • CA, it becomes 16.8243 MPa at the TDC. As a general trend, the pressure curve has a single high peak at the TDC with an asymmetric steep variation on both sides of the peak indicating a high burning rate before TDC and high thrust after TDC.  Figure 5b shows the mass averaged temperature curve for the whole simulation duration while Fig. 5c shows the maximum static temperature. From Fig. 5b, we note that the temperature rises relatively rapidly with compression as compared to the behavior of the pressure at the start of the simulation which may be attributed to the high pressure-volume work done by the piston on the gas resulting into a rapid rise in the temperature in view of the ideal gas law. As the simulation proceeds, the temperature keeps on rising so that it becomes 850 K at the crank angle 691 • where the fuel is to be injected. This temperature is high enough to initiate spontaneous ignition as the fuel is injected and the air-fuel mixture is ready because the ignition temperature of the decane fuel has been reported to be 483.1 K. The highest average temperature is achieved at the TDC which is slightly above 1200 K. Although the maximum average temperature is around 1200 K, the maximum local temperature attains a value as high as 3004 K as shown in Fig. 5c which is an indicator of NO production locally in the cells of high temperature. Figure 5c also reflects that the maximum temperature does not show a significant rise until the ignition starts. The steep rise at the crank angle 694 • indicates the onset of ignition reflecting ignition delay of only 3 • CA. The peak of the local maximum temperature occurs at 11.45 • CA before TDC and has the value 3004.2161 K that may be the CA at which the combustion is taking place at its highest rate. Almost all of the heat is released within the duration of around 23 • CA starting with the ignition at 694 • CA till 717 • CA that is 3 • CA before TDC as depicted in Fig. 5d. The maximum heat release occurs at the same crank angle at which the peak of the local maximum temperature occurred. Figure 6a presents mass fraction of soot. We note from the figure that there is no significant soot formation in the 8 • CA duration after the start of ignition occurring at 694 • CA. However, at 702 • CA, soot formation starts and its rate of production increases quite rapidly at 706 • CA. This rate continues increasing till the point of inflection of the soot curve occurring at 714 • CA. After this point, the production rate starts decreasing quite insignificantly while it comes to zero quite rapidly at around 760 • CA. In this way, soot formation continues for quite a long duration after the heat release rate has become zero at 717 • CA indicating the completion of the combustion process. By comparing Figs. 5d and 6a we see that during the interval [717 • CA, 769 • CA] soot production occurs through endothermic reactions as the heat release rate is negative in this duration. This production could be stopped by the soot oxidation process whose study is not part of the present work.
In Fig. 6b, the NO mass fraction is plotted against the crank angle. A comparison of Figs. 5d and 6b shows that the major NO production starts when the heat release rate becomes significantly high and its production continues till the combustion stops when the heat release rate becomes zero, i.e., at 717 • CA. A comparison of Figs. 5c and 6b indicates that the minimum temperature at which NO production starts is above 1800 K which is achieved at 1.4 • CA after the start of ignition. The same temperature requirement for NO production has been reported in the literature by [23]. As the local high-temperature zone increases, the NO production also increases. The NO mass fraction at the end of the combustion is 0.000477528 and the same has been reported by [9] for one of his piston designs. The above findings also validate our simulation results. Now, we examine the local concentrations of soot and NO at various locations in the cylinder. In Figs. 7 and 8 , mass fractions of soot and NO are given in different horizontal planes in the cylinder at 833 • CA. The figures show that the concentrations of soot and NO are higher in the piston bowl. This is because, in view of the spray angle taken in the present design, the fuel is sprayed in the piston bowl and most of the air-fuel mixture lies in the piston bowl and travels with the piston. Moreover, post-combustion production of soots also takes place in the piston bowl.