A Phenomenological Model for Prediction Auto-ignition and Soot Formation of Turbulent Diffusion Combustion in a High Pressure Common Rail Diesel Engine

A new phenomenological model, the TP (Temperature Phase) model, is presented to carry out optimization calculations for turbulent diffusion combustion in a high-pressure common rail diesel engine. Temperature is the most important parameter in the TP model, which includes two parts: an auto-ignition and a soot model. In the auto-ignition phase, different reaction mechanisms are built for different zones. For the soot model, different methods are used for different temperatures. The TP model is then implemented in KIVA code instead of original model to carry out optimization. The results of cylinder pressures, the corresponding heat release rates, and soot with variation of injection time, variation of rail pressure and variation of speed among TP model, KIVA standard model and experimental data are analyzed. The results indicate that the TP model can carry out optimization and CFD (computational fluid dynamics) and can be a useful tool to study turbulent diffusion combustion.


Introduction
Today's strict emission level legislation limiting combustion-generated pollutants sets high demands on the development of cleaner engines.Since modern engines have already reached a high level of optimization through engineering measures a detailed understanding of the different processes leading to the formation of pollutants is required [1].Soot, which is a product of incomplete combustion reducing the combustion efficiency, is identified in flames and fires as yellow luminescence and is mostly an unwanted process in today's practical combustion devices.Soot deposits in engines or gas turbines are problematic and cause reduced efficiency, increasing the need for maintenance.Another severe problem is the effect of soot particles and their precursors on human health [2].
Due to the complexity of soot formation and the general lack of knowledge, a number of soot models have been proposed.The inception of soot particles is still not well-understood, but it is generally accepted that polycyclic aromatic hydrocarbons (PAHs) are the key soot-forming species.In the model it is assumed that the rate limiting step is the collision of two PAHs leading to the first three-dimensional structures held together by van der Waals forces.Mitchell and Frenklach [3] investigated the dynamics of agglomeration with simultaneous surface growth using a Monte Carlo method.It was found that the shape of the aggregates depends markedly on the surface growth, e.g., when surface growth is very fast the aggregates have spherical structures.To reduce the numerical demands of the method a reduction of the PAH polymerization model was proposed by Mauss [4].Assuming that the polymerization reactions are fast or in other words, assuming that the transition from PAHs to soot is a fast process, all PAHs can be assumed to be in a steady-state.The reduction leads to an infinite system of algebraic equations describing the local PAH concentrations from which the moments of the size distribution can be calculated.
Calculations of soot formation and oxidation in diesel engines can mainly be classified into different model approaches.Balthasar et al. [5] focused mainly on the fluid dynamics of the problem using simple models for the chemical processes.Soot formation is modeled with the help of empirical or semi-empirical models.The most recent approach is based on the laminar flamelet concept, solving the so-called Representative Interactive Flamelets (RIF) on line with the CFD code.A full description of the chemical processes as well as soot and NO x formation can be obtained without having to simplify the fluid dynamics.The latter approach is the most detailed, but time consuming, and might be more difficult to incorporate into existing CFD code.
Thus, there is a need to develop numerical calculations based on accurate modeling and analysis based on the latest research results.The research presented herein focuses on obtaining an understanding of the ignition phase and soot model in a real high pressure common rail diesel engine through numerical simulations.First, a new model called Temperature Phase (TP) model is proposed, the auto-ignition phenomenon is described in some detail and implementation processes are given.Second, the experimental set-up and boundary conditions required for the calculations, as well as the operating points investigated will be discussed.Next, the simulation results and the measured data are compared.In order to assess the performance of the matching results in terms of predicting cylinder pressure and heat release over a wide range of part-load conditions, an extensive parameter study varying injection timing and rail pressure has been conducted.

Modified Model
Since auto-ignition is the first important event of diffusion combustion under diesel conditions, we began to construct a kinetic mechanism for n-heptane ignition based on the detailed model proposed by Barths et al. [6].However, Barths' mechanism (involving 324 chemical species and 1650 reactions) is too large to use in multidimensional turbulent diffusion computations.Therefore, only a few lumped or equivalent intermediate species (e.g., the first and second n-heptyl isomers were retained to represent all the n-heptyl isomers) and some relevant key reactions were selected on the basis of steady state assumptions.These reactions are crucial for reproducing the auto-ignition behavior of n-heptane/air mixtures from the low/intermediate regime to the high-temperature regime.The high temperature oxidation of n-heptane proceeds from the attack on the fuel by H, OH and HO 2 radicals to form n-heptyl radicals and through the break-up of these into C 2 H 4 , CH 3 and H radicals.These are oxidized by C 1 -C 2 -chemistry reactions.The low temperature chemistry of lower aliphatic hydrocarbons is characterized by the degenerate chain branching which may be illustrated by Table 1.This low temperature mechanism is no longer valid when the temperature increases beyond about 850 K.The competition of the reverse reaction of the first and second O 2 -addition with the subsequent internal H-abstraction reaction is the key to the understanding of the negative temperature dependence of the ignition delay.With increasing temperature, these reverse reactions become faster than their forward reactions, thereby stopping the reaction sequence.A transition to the high temperature mechanism must occur.The apparent negative temperature dependence of the intermediate branch shall be explained by discussing a simplified 2-step ad-hoc model with adjusted n-heptane ignition rate constants, which is shown in Table 2.A lot of parameters can be found in our previous paper [7][8][9][10].However, at that time we only deduced and discussed these equations in theory and did not use them in a real diesel engine.Now we have used them in three-dimension turbulent diffusion combustion with turbulence movement in a diesel engine.For high temperatures, we used the original KIVA-3V model [11], assuming a simplified single-step oxidation model that converts fuel and oxygen to carbon dioxide and water, because the temperature is high and reactions are very quick.We replaced the kinetic chemistry subroutine CHEM with our own specialized model.

Implementation in CFD Code
Two steps are needed to implement the TP model in CFD code (KIVA-3V).One is making the TP model subroutine and the other is coupling the TP subroutine to the main KIVA program.The following is the description of how to achieve the above two steps: (1) TP model subroutine The TP model subroutine structure is shown in Figure 1.First, we select different reaction mechanisms according to different grid temperatures.When the temperature T < 850 K, we the select 6-step n-heptane model from Table 1 to calculate the auto-ignition, which explains the diesel ignition delay "S" curve in the lower part.When the temperature is 850 K < T < 1100 K, we select the 2-step n-heptane model from Table 2 to calculate the "S" curve in the middle of some instability.When the temperature is 1100 K < T, we use a one-step n-heptane model to calculate the upper S curve in order to save computation time.Because each grid temperature is different, we select reaction mechanisms according to each grid temperature thus avoiding all grids having the same calculation.This calculation strategy not only satisfies the original design requirements, but also greatly improves the computational efficiency.At the same time, we add a real-time control program to monitor the simulation.After each elementary chemical reaction calculation, the corresponding pressure, heat release rate, and other related information are output.If the intermediate calculations exit, it shows some parameter calculation value is divergent.The output of the calculation information is very important and it gives the scientific basis for normal calculations.
(2) Coupling TP subroutine in KIVA main program To couple the TP model subroutine to the KIVA code, we need put it before the chem.fsubroutine.The relationship between the TP subroutine and the main program is shown in Figure 2. If we want to make the TP subroutine operate smoothly, we need to calculate KIVA's k-å turbulent movement model.Although the k-å model can provide unsatisfactory results, it can fit the calculation of the fuel injection engine.If the TP subroutine quits, we need to check the initial cylinder temperature, the cylinder pressure and other related initial data which are in the itape5 file.Because KIVA-3V can only recognize regular hexahedron shaped grids, we need to constantly adjust the establishment of x, y, z grid layers in three directions to ensure the normal calculation.

Soot Model
The modified soot model includes two parts: the PAH and the soot source term.The soot model used in the calculations relies on a detailed description of the chemical processes leading from the fuel molecules to complex soot particles based on the statistical momentum approach.It can be subdivided into the growth of polycyclic aromatic hydrocarbons (PAHs) in the gas phase reactions and the processes of particle inception, condensation, surface growth, and oxidation.

PAH Calculation
According to recent literature [12] a "skeleton" reaction mechanism plays a major role in the PAH soot precursor whereby the material (C 2 H 2 and C 4 H 5 ) is the basis for reactions to form the first benzene ring through the "add-dehydrogenation-acetylene" mechanism of the C 10 H 8 representation of PAH components.The C 10 H 8 oxidation product is not CO 2 , but CO and HC.The PAH model calculation is a modified Frenklach-Mauss fuel pyrolysis model [13] of the formation of the first benzene ring, cyclization and oxidation, whose units are seconds (s), centimeters (cm), Kelvin (K), and mole (mol).
(1) Fuel pyrolysis First of all pyrolysis of hydrocarbon fuels occurs at high temperature to produce a variety of groups of small molecules, such as CH 3 , C 2 H 2 , C 4 H 2 and C 4 H 5 , etc.In this paper, using n-heptane as the diesel calculation model, the pyrolysis equation can be written as follows: (2) The formation of the first benzene ring Based on the pyrolysis of fuel through acetylene (C 2 H 2 ) and C 4 H 5 the equation for the first aromatic (benzene) ring formation reaction is: (3) Cyclization The formation of a single benzene ring will be followed by the so-called HACA mechanism for the formation of aromatic molecules.HACA (H-Abstraction, C 2 H 2 -Addition) means that the hydrogen atom desorption from the addition of phenyl acetylene molecules, referred to as "plus-acetylenedehydrogenation" mechanism, but this calculation is different from the FM in the detailed model, using literature [14] in the test results.Using naphthalene (C 10 H 8 , abbreviated as A2, i.e., an aromatic hydrocarbon with two benzenes) as the PAH the representation of the nuclear components during the chemical reaction is as follows: (4) PAH growth, rapid polymerization and oxidation PAH growth in the early stages mainly produces planar structures.Due to the special stability of certain aromatic components, it is difficult for the reverse reactions to occur, so these rings involving accumulation of hydrocarbon molecules will grow rapidly.Of course, this growth will not be unlimited, because at the same time, an opposite process occurs, namely the oxidation of hydrocarbons, so the growth of PAHs is inhibited after growth exceeds more than four PAH rings.Growth continues by relying on the so-called "linear heap together" approach (linear lumping technique), and its principles are still based on HACA.However, due to the involvement of OH and O and PAH-generated reactive species in chemical reactions, and the presence of some reverse reactions, the formation and consumption of PAH between the two processes reaches a certain "balance".PAHs inhibit the growth of lean combustion processes and taking into account the oxygen-rich state in the linear stack PAH accumulation process, the introduction of OH and O 2 causes oxidation.As the OH-based oxygen activation energy activation energy is only 2/3 that of O 2 , OH is a more effective oxidant than oxygen.The PAH oxidation product is not CO 2 , but rather CO and HCs, and the corresponding chemical equation is: (5)

Soot Source Term
The formation of soot can be subdivided into the processes of particle inception, surface growth, oxidation and coagulation. (

1) Particle inception
In the present study, we have followed the work of Tao and Golovitchev's, and assume that soot particle nucleation takes place at sites where diacetylene (C 4 H 2 ) or naphthalene (A2) exist.The choice of these two species is based on the argument that soot particle nucleation is not merely linked to the pyrolysis of the hydrocarbon fuel, as assumed by Lindstedt [15], but it also has its own, unique macro-molecular precursor formation chemistry.Strictly speaking, however, these two reaction steps do not represent the soot particle nucleation process correctly, because soot does not only contain carbon atoms in its incipient particles, but also considerable amounts of hydrogen.In addition, there are uncertainties regarding the determination of the rate constants of these two reactions.As will be shown below, the values of the effective activation temperature for soot nucleation used in Lindstedt's model are too high for use in the present study.Unless a very high pre-exponential factor is employed, computations with Lindstedt's nucleation model predict a rather low soot mass concentration.We have found that appropriate values for the effective activation temperature are below 5000 K, thus we assign a value of 1500 K for this parameter in the present study because it provides reasonable predictions compared with the available experimental data, especially for the short time available for soot formation immediately after ignition.Initiation of the process of soot and hydrocarbon fuels is not only related to the thermal solution, but also with a specific precursor-related chemical processes.Results show [16,17], that flame and carbon-related intermediate species, mainly PAH and polyacetylene, play an equally important role in soot particle nucleation, which would follow Equations ( 6) and ( 7): Both the total package reaction to remove a large number of intermediate species and reaction details will reduce the computation time, but also provides a determined carbon grain initial surface area and a simple way to follow-up on the surface growth and particle number density of the initial calculations.Based on these two reactions, it is easy to write the soot nucleation rate formula: (2) Surface growth We have adopted the idea of Frenklach and Wang's active site model [18] for the present study.However, the model was modified by the inclusion of the irreversible reaction to account for the activation of surface sites by OH radicals and the deactivation of surface radical sites by H 2 O.This reaction is found to be critical in the soot formation in diesel spray combustion and the rate constants are fully consistent with those of the HACA reactions of PAHs.Only part of the rate constant for reaction is adjusted, to achieve agreement with experimental data.
(3) Surface oxidation Soot oxidation is a heterogeneous process that takes place on the surface of soot particles and depletes the carbon atoms accumulated in the particles.However, since the mechanism of this process is still poorly understood, it is commonly described using empirical formulations such as the Nagle and Strickland-Constable semi-empirical model to account for soot oxidation caused by O 2 attack.This model is adopted in the present study.The soot oxidation is also sensitive to other oxidants such as OH radicals.We include Neoh et al.'s OH oxidation model and the rate, correlated to the collision efficiency of OH radicals with the soot surface.
(4) Particle coagulation Particle coagulation causes the number of soot particles to decrease.We assume a monodisperse size distribution of the particles, and then describe the coagulation by the Smoluchowski equation. (

5) Soot Source Terms
The formation of soot can be subdivided into the processes of particle inception, surface growth and oxidation, particle coagulation.For the PAHs a statistical approach is used to describe the size distribution function of the soot particles.The moments of the size distribution of the size are defined as: where N i is the number density of particle i with a mass , with m1 being the mass of the smallest unit occurring in a soot particle.The moment M 0 is equal to the total particle number density: The moment M 1 can be related to the volume fraction, which defines the ratio of the volume occupied by soot particles to the gas volume: where ρ s = 1800 kg/m 3 is the density of soot.
Particle inception can be modeled as the coagulation of two PAH-molecules.Coagulation of particles of the same type can be described by Smoluchowski's equation: This equation gives the change of the particle number in size-class i as a function of time.The first terms in Equation ( 13) describes the formation of new particles from smaller sized particles and the second terms the consumption of particles in the i-th size-class by coagulation with particles of all size classes.The frequency factor in the free molecular regime is given by: where k B is the Boitzmann constant, μ i,j is the reduced mass, r i is the radius of particles of class i and ε i,j is the size dependent coagulation enhancement factor due to attractive or repulsive forces between the particles.The Smoluchowiski equation can be formulated for the particle inception omitting the second terms: Equation ( 15) can be described in terms of the moments of the PAH size distribution by multiplying with r i and a summation over all size classes (note the change of the upper summation limit with respect to j): One can rewrite (7) with the help of where V i and m i are the volume and mass of particle of size class i respectively, and ε i,j the coagulation enhancement due to inter-particle forces.The enhancement factor due to van der Waals interaction is set to a constant value of 2.2.The mean soot particle diameter is given by: Tao's experiment involved a conventional heavy-duty DI diesel engine based on injection conditions.Sandia National Laboratories reproduced the experiment in the optical engine using similar operating conditions.There are many advantages than DEC [19][20][21] conceptual model, particularly in the liquid spray penetration, fuel vapor permeability, flame structure and soot concentration on the spatial correlation.Tao et al. knew, when the larger soot particles are in fuel-rich vortex downstream of the head region, there is a thin layer of a large concentration of soot particles in the very low quality flame of the external surface area.They gave many valuable details in soot formation (for example, particle size, and density, form the material before, and acetylene concentration) and response speed (for example, the formation of soot particle coagulation, surface generation, and oxidation).

Numerical Implementation
The KIVA-3V code is used to solve the time-dependent conservation equations of total mass, momentum, energy, and species mass concentrations in the reacting gas-phase mixture.Spray dynamics are modeled using the discrete-particle technique, in which each computational particle represents a number of droplets of identical size, velocity, etc.In addition to particle-turbulence interaction, spray atomization, droplet breakup, evaporation, collision and coalescence are treated using spray sub-models for a single-component, vaporizing fuel sprays are provided in the original KIVA-3V code.Turbulence modeling is based on the compressible k-ε model (accounting for dilatation effects) modified to include the effects of fuel droplets.A characteristic micro-mixing time scale is determined from the turbulence modeling and then further used for the chemistry-turbulence interaction.All these sub-models implemented in the KIVA-3V code have been proven relevant to the accurate simulation of diesel spray combustion.

Basis for the Numerical Simulation
The computations started at IVC and ended at EVO.The injection time unit from us to CAD (crank angle degree) can use the following form for a 4-stroke engine: where φ is crank angle degree, n is speed (r/min), t is time (μs).The wall temperature (450 K) was set such that calculated and measured pressures matched during the compression phase before injection started.The wall temperature was held constant during the computations.The swirl was set to 0.5 times the number of revolution of the engine (r/min) which corresponds to the measured swirl.The injection nozzle was located on the axis of symmetry.Since the shape of the piston is not on the axis of symmetry, 360° was modeled which corresponds to x-offset 4 mm and y-offset 1 mm.This made the model 360-degrees, offering exact calculation of the full geometry.Figure 3 shows the computational grid used in the simulations containing 204,000 cells.

Experiment Set-up
In the process of engine bench tests, the cylinder pressure must be strictly monitored for the various tanks, therefore, the cylinder head design needs to be modified with the four warm-plugs and the cylinder pressure sensor is installed on the transition casing.Figure 4 is the cylinder pressure sensor design and structure.The in-cylinder pressure was measured using a flush mounted quartz sensor from Kistler (type 6061B).The pressure transducer is very accurate on a relative basis but does not directly yield absolute values.In this measurement, the reference was taken to be atmospheric pressure and it had to be corrected for the applied boost pressure.The peak motored cylinder pressure is extremely sensitive to small changes in the pressure at intake valve closing (IVC).
Figure 5 shows a schematic of the test bench configuration.The test engine is a four cylinder, 2.8 L 4JB1 DI engine.Bore is 93 mm, and stroke is 102 mm.Swirl ratio is 2.4.Compression ratio is 18.2.Connecting rod is 168 mm.Injector protrusion is 1 mm.The injection system is a third generation Bosch Common-Rail featuring a maximum injection pressure of 2500 bar.The injector is equipped with a six-hole nozzle, where each hole has a diameter of 0.124 mm.Soot was measured with the ELPI (Electrical Low Pressure Impaction).The ELPI is a real-time particle size spectrometer designed at the Tampere University of Technology for real-time monitoring of aerosol particle size distribution.The ELPI measures airborne particle size distribution in the size range 0.03−10 µm with 12 channels.With a filter stage the size range can be extended down to 7 nm.The operating principle is based on particle charging, inertial classification in a cascade impaction, and electrical detection of the aerosol particles.The main components of the instrument are a corona charger, low-pressure cascade impaction and multi-channel electrometer.Euro III and Euro IV but here we discuss only one work condition due to the size limitations of this paper.The investigated load conditions are 2200 r/min, 25% load, 50% load and 75% load, and 2700 r/min, 75% load are from 13 work conditions which include different speeds, different torques and different rail pressures.Previous investigations have shown that the load predictions part represents a particularly challenging modeling task [22].The injector is slightly offset (by 1.0 mm) to the center of the cylinder axis to allow a better cooling of the narrow bridge between the exhaust valves.Whereas the engine experiments were performed using standard diesel fuel, a model fuel comprising n-heptane was used to represent diesel in the simulation.Rutland et al. showed that the model fuel displays very similar performance and emission characteristics to diesel.
An important input for the spray calculation is the injection rate profile.In this study, the fourth generations of EFS 8420 France products was used for consistency of measurement, reproducibility, responsiveness, the ability to spray twice for the assessment of its verification aspects.In order to accurately measure the volume of a single injection, we need to set the nozzle back pressure in the external environment, usually to 20 bar, when the back pressure built up and stabilized at around 20 bar rear pressure can be measured.Four generations of products in terms of injection rate measurement principle, monitor friendly interfaces, operability, and so the two jet measurement capabilities are greater than the three generations of product improvements, the fuel injection quantity per injection can be a data output, injection rate curves can be displayed in real-time monitoring

Cylinder Pressure and Heat Release Rate
For each of the simulations cylinder pressure and heat release rate are compared with measured values.Figure 7 shows the cylinder pressure traces and corresponding heat release rates for the measured and simulated parameter study varying the work conditions.As can be seen in Figure 4, the overall agreement in both ignition delay and peak cylinder pressure are excellent.For 2200 r/min 64 Nm, the cylinder pressure is changed quickly when it is 20 °CA (Crank angle) BTDC (before top dead center).So the ignition delay time is 12°CA.The biggest experimental pressure is 5800 KPa, which is near to TP model (5820 KPa) and far from KIVA original model (6020 KPa).The heat release rates are similar between TP model (45 J/°CA) and experiment (42 J/°CA), but in the original KIVA model it is 55 J/°CA.This is due to fact the TP model is a multi-step chemistry model and KIVA is a one step model.But the difference between the original model and the improved model decreases with the fuel injection quantity, because as the fuel injection quantity increases, power increases, the cylinder temperature and pressure increase and the difference between the multi-step and one step models decreases.The cylinder pressures are well fitted between simulation and experiment for all other work conditions.The simulation heat release rates are higher than the experiments because of different calculation steps, but the errors are small and the TP model can be used for simulating a real diesel engine.fuel by the original model.The experimental value is 0.09 g/kg fuel.This shows the improved model is nearer to the experimental result than the original model.If we use the original model, the maximum soot value is 2.57 g/kg fuel at 3°CA TDC.There is 1°CA difference between the original model and the improved model, which is in line with numerous papers [23][24][25].When the engine is at 2200 r/min and 135 Nm, the maximum soot occurs after 7°CA TDC with 2.02 g/kg fuel by the improved model and the final soot is 0.0577 g/kg fuel.The peak soot is also at this point according to the original model, but the peak value is 2.42 g/kg fuel with final generation 0.0692 g/kg fuel.The experimental value is 0.041 g/kg fuel.The final soot by the improved model is lower than by the original model but closer to the experimental.The improved accuracy was about 25%.This behavior has been observed under other conditions.

Summary of Investigated Load Conditions
Figure 9 summarizes the different parameter variations performed.It is seen that a variation of rail pressure leads to soot variation.When the rail pressure is 980 bar, soot source is the least (about 0.022 g/kg fuel).When the rail pressure is 1190 bar, soot source is the biggest (about 0.141 g/kg fuel).It means soot is related to capacity.When the capacity is increased, soot is increased.But it is not a proportional relationship.This is due to different working conditions caused by different injection times.For the initial injection, there is less fuel, and the amount of soot produced is smaller.For the second injection, the jet occurs at a high cylinder temperature.The chemical reaction is very quick and the ignition delay is very short.The sprayed fuel is quickly ignited and the generated soot is oxidized continuously, so the final generated soot is relatively small.The original KIVA model and the improved TP model have the same tendency, but the improved model is closer to the experimental than the original model.In particular when the engine is at 2200 r/min (211 Nm) and 2700 r/min (195 Nm), the gaps are about 50%.This shows that the improved model is more realistic than the original model.

Conclusions
(1) Our improved combustion model can accurately predict cylinder pressure.Since the original KIVA-3V model reaction operates under the assumption that combustion is complete, it overestimates the heat release rate, in fact after all the fuel is oxidized, CO oxidation is still underway.Since the original KIVA-3V computational model of the activation energy has no real meaning, any calculation of cylinder pressure curves does not reflect the diesel ignition delay, while the use of this improved model can be used to calculate the ignition delay of diesel fuel.
(2) KIVA-3V calculates the heat release rate of the original model on the assumption that the ratio of combustion products CO and CO 2 remains the same, but this does reflect the actual situation.The instantaneous heat release rate by the original model difference from the test value is 30%.The instantaneous heat release rate by the improved model is in good agreement with the experiment.
(3) The original KIVA-3V model did not calculate soot precursor-PAH, and did not consider the PAH growth and the initial particle density, so the soot formation by the original KIVA-3V model cannot reflect the actual situation.However, the modified soot model considers these factors, and the soot determination accuracy increases about 50%.

Figure 2 .
Figure 2. Relationship between TP subroutine and KIVA main program.

Figure 3 .
Figure 3. Computational grid used in the simulations.

Figure 4 .
Figure 4. Cylinder pressure sensor structure and design.

Figure 5 .
Figure 5. Schematic of the experimental set-up.

Figure 9 .
Figure 9. Summary of parameters variation study.

Table 1 .
Low temperature mechanism

Table 2 .
Middle temperature mechanism