Phenomenological Modeling of Combustion Process in Diesel Engines Based on Stochastic Method

In order to satisfy the growing demand for the reduction of fuel consumption and pollutant emissions, various technologies have been employed in diesel engines. Consequently, to determine the optimal combustion control strategy, many parameters such as injection pressure, nozzle diameter, injection timing, injection quantity, and exhaust gas recirculation (EGR) rate should be selected properly corresponding to the engine operating conditions. It is difficult to obtain the appropriate strategies without understanding the change in combustion process when varying these parameters. To realize parametric studies on combustion control strategy of modern diesel engines, a phenomenological combustion model based on stochastic method was developed. In this model, the modeling of the spray tip and tail penetration after the end of injection, and interaction between the sprays of sequent injection stages were focused on to modify the stochastic combustion model for combustion simulation with multiple injection. The effects of swirl, wall impingement, and adjacent spray interaction are formulated simply to make the combustion model more accurate and computationally efficient. The simulation results were compared with experimental data from a singlecylinder test engine for pilot/main two-stage injection. The results reveal that the model has capability to accurately predict the combustion characteristics and emissions of diesel engine with pilot/main two-stage injection.


Introduction
Recently, continuous reduction of the oxides of nitrogen (NOx) and PM emissions for confronting increasingly stringent emission standards becomes a significant task for combus-tion improvement of diesel engines; especially the reduction of NOx emissions is quite difficult for aftertreatment for low-load operating conditions due to low exhaust gas temperature.Hence, in-cylinder combustion improvement is required for low-load operation, subsequently the advanced combustion concepts have been proposed for dealing with this situation, such as premixed charge compression ignition (PCCI) and low-temperature combustion (LTC).No matter whether it is PCCI or LTC, the main aim is to obtain more homogeneous mixture and suppress locally high combustion temperature.This aim is typically achieved with injection timing adjustment with high exhaust gas recirculation (EGR) rate.However, accompanying with reduction of NOx and PM, the CO and HC emissions are observed to increase and even exceed the emission standards [1][2][3][4].Meanwhile, long ignition delay promotes the fuel-air mixing before ignition; it leads to the strong dependence of ignition timing on chemical reaction, which makes the ignition timing control difficult.And fast combustion of the mixture premixed before ignition increases the combustion noise associated with rapid rise of incylinder pressure.Because of these barriers, the advanced combustion concepts are limited in low-load operation.
To extend the advanced combustion concepts to high-load operation and obtain the best reduction of emissions, many technologies, including high-injection pressure, small nozzle diameter, multistage injection, EGR, and high-intake boost, have to be integrated to realize incylinder combustion improvement.Consequently, to determine the optimal combustion control strategy, many parameters such as injection pressure, nozzle diameter, injection timing, injection quantity, and EGR rate should be selected properly corresponding to the engine operating conditions.It brings a heavy work and high cost on experiments, moreover it is difficult to obtain the appropriate strategies without understanding the change in combustion process when varying these parameters.Thus, an accurate and computationally efficient combustion model is needed for these parameters study.Considering the trade-off between accuracy and computational efficiency, the phenomenological combustion models are often adopted for the parametric study [5].
In phenomenological combustion models, the combustion process is mimicked mathematically based on the results of observed phenomena with focus on predominant and ratedetermining processes in combustion.To describe the heterogeneous distribution of fuel concentration and temperature and low-temperature combustion phenomena, some detailed phenomenological diesel combustion models are developed based on the spray behavior, which incorporate the calculation of chemical reaction in the spray region.According to the method for description of fuel concentration distribution, these models can be divided into two types.One type is to divide the spray into different zones based on the spray structure [6][7][8].In the other type, fuel concentration distribution is represented by probability density function (PDF) of fuel mass fraction [9,10].PDF method is dominant to couple the turbulent mixing and chemical kinetics to calculate the turbulent combustion process in diesel engines, and it is easy to blend with the emission models, such as NOx and soot models.Whence the pollutants production process can be predicted and analyzed with varying relevant parameters.
In modern diesel engines, multiple-injection has played a significant role in the latest engines to reduce NOx, smoke emissions, and combustion noise simultaneously [11][12][13].Proper selection of injection strategy from numerous injection parameters satisfying different operating conditions is a challenging work.Thus, to estimate the modern diesel engines combustion, the combustion process with multiple injections needs to be included in the phenomenological combustion model.Recently, some models had been proposed based on PDF method to calculate the combustion with multiple-injection [14][15][16].In these models, each injection is customarily treated to form a spray region, and relevant calculations including spray propagation, fuel evaporation, turbulent mixing, ignition, and chemical reaction are performed in each spray.According to the sprays propagation phenomena under multipleinjection situation, the combustion of the latter spray is significantly dependent on the thermodynamic states of the former spray.Therefore, the interaction between former spray and latter spray has to be considered in the model.Normally, the latter injection starts after the end of the former injection, so that the evolution of the spray after the end of injection (EOI) has to be imitated.In fact, the spray behaviors after EOI are indicated quite different from that during injection by the fundamental studies on turbulent gas and diesel jets [17][18][19][20][21].However, there is no simple model available in phenomenological models to describe the spray development and air entrainment after EOI.Currently, frequently used models of spray propagation are Wakuri's model, Siebers's model, and Hiroyasu's model that only satisfy the spray propagation during injection with steady state and do not involve any information for that after EOI.
According to the modeling concept from the developed models [14][15][16], the latter spray is assumed to be injected into the former spray region, so that the entrained gas of the latter spray only comes from the former spray zone.It should be noted that all of these interaction treatments assume that later spray entrains the gas from the former spray zone immediately when the later injection starts, whereas some fundamental studies discover that the mixture near the injector tip is diluted very fast after EOI [20,22].This implies that the gas near the injector tip should be very lean when the later injection starts, particularly for the early pilot injection case.Thereby, it is necessary to consider this phenomenon in appropriate form to imitate the interaction between sprays from the sequent injections.Meanwhile, the swirl flow deflects the spray from the initial injection direction [23].In multiple-injection case, the magnitudes of deflection are different for former and latter sprays, which results in that the swirl effects has to be involved for spray-to spray interaction modeling.
In addition, the phenomenological combustion models usually do not consider the spatial information.However, the early injection is a typical strategy for advanced combustion concepts, especially quite early pilot injection for multiple-injection strategy.The early injection causes the wall impingement, consequently the effects of different impinging positions and the spray flames propagation along the cylinder wall on fuel-air mixing need to be involved in the model.Thus, in this study, a phenomenological combustion model was developed based on a stochastic mixing model.The heterogeneity of fuel concentration and temperature is described using the PDF.Evolution of the mixing and combustion is represented by the temporal change of PDF due to the mixing, fuel injection, fuel evaporation, air entrainment into a spray, heat loss to the combustion chamber walls, and chemical reactions.To focus on the simulation of new combustion concepts, such as LTC and PCCI, the effects of EGR and multiple-injection was involved in the model.The model employs a quasi-global chemical kinetics model combined with stochastic mixing model to describe the low-combustion process caused by EGR effect.The key points of the study concentrated on are as follows: (1) Spray model developing to describe the spray evolution after end of injection (EOI).(2) Defining and calculating the presumed spray tail to consider the quite lean mixture near the nozzle exit after EOI.(3) Modeling spray-to-spray interaction between sequent injections and the swirl flow effects on it.(4) Introducing the effects of spray-wall interaction on fuel-air mixing.Finally, the simulations are performed in pilot/main two-stage injection cases, and are validated by experimental data from a single-cylinder diesel engine, including in-cylinder pressures, heat release rates, and NOx and PM emissions.

Stochastic combustion model
The stochastic combustion model in this study was first proposed in [10] for conventional combustion with single-stage injection in DI diesel engines.And then this model was modified in [24][25][26] to analyze the PCCI combustion process.Up to Ref. [26], this model has employed a stochastic turbulent mixing model coupled with a quasi-global chemical kinetic model to describe the heterogeneous combustion process in diesel sprays, and the specific submodels includes air entrainment, turbulent mixing, fuel vaporization, heat loss, and chemical kinetic model.

Model description
The essential of the stochastic combustion model focus on the turbulence controlled fuel-air mixing process.In this model, the heterogeneity in fuel concentration and temperature is described using probability density functions f(x, t), where x = (h, y 1 , y 2 , … , y n ), h denotes the specific enthalpy and y 1 , y 2 , … , y n the mass fraction of chemical species.Based on considering the in-cylinder pressure as uniform distribution over the entire cylinder, the time evolution of f(x, t) is expressed as Eq. ( 1) where (•) expresses the time derivative, p the in-cylinder pressure, v the specific volume, m(x, t) represents the mixing caused by turbulence, and ω the mixing rate.To calculate the f(x, t), the pressure Eq. ( 2) was derived in [5] to simultaneously solve with Eq. ( 1) where . is the average specific volume variation rate corresponding to engine speed.
However, these two differential equations are too complicated to obtain the analytical solution for f(x, t).Therefore, the Monte Carlo method is adopted to solve the equations.This is realized by treating numerous Monte Carlo elements as the fluid particles.The injected fuel and charged air in cylinder are assumed to be sufficiently small fluid particles compared to the microscale of turbulence, and the fluid particles are considered to have identical mass.In this process, the change of fluid particle's thermodynamic states due to turbulent mixing is described by the random collision and dispersion process occurred between fluid particles.It is calculated by modified Curl's model [27] in which a randomly selected pair of fluid particles in the spray collides and disperses with exchange of their thermodynamic states as expressed in Eq. ( 3) The x 1 and x 2 represent the thermodynamic states of two fluid particles, and X is a random number uniformly distributed in (0, 1).The mixing rate is represented by collision frequency that is estimated by turbulence kinetic energy in the spray.Based on this concept, the combustion process of diesel spray is described as shown in Figure 1.In this process, the combustion chamber is divided into spray zone and surrounding air zone after the injection start.The spray zone is formed by injected fuel particles and entrained air particles whose amount is determined by momentum theory [28].In the spray zone, the fuel-air mixing is taken into account with modified Curl's collision and dispersion model [27].The mixing rate is controlled by stochastic collision frequency that is dependent on turbulence kinetic energy in spray.Meanwhile, liquid fuel contained in each fluid particle is vaporized according to an average evaporation rate, sequentially chemical reaction calculation using Schreiber's model [29] is performed in each fluid particles for combustion calculation.In addition, the heat loss from wall is also considered using Woschni's equation [30].Thus, in the stochastic combustion model, fuel-air mixing and chemical reactions are simultaneously calculated.To simulate the combustion with multiple-injection, the situations shown in Figure 2 can be abstracted as that shown in Figure 3.The volume in the cylinder is divided into three zones to represent the sprays and ambient gas for two-stage injection situation, including the first spray zone, the second spray zone, and the ambient air zone.Before injection, the air and EGR gas in the cylinder is treated as ambient air zone that is presumed to be composed of a great number of fluid elements.Then after the start of the first injection, the fuel is injected into the cylinder as fluid elements, and entrained fluid elements from ambient air zone to form the first spray zone.After the end of the injection, the first spray tail is assumed to depart from nozzle exit to the downstream.When the second injection starts, the first and second spray zones exist in the cylinder and the fluid elements of ambient air zone are entrained into these two spray zones, respectively.Once the second spray tip overtakes the first spray tail, both of the fluid elements from the ambient air zone and the first spray zone are entrained into the second spray zone.Finally, these two spray zones are combined into one spray zone after all of the fluid elements of the first spray zone are entrained into the second spray zone.Obviously, the spray tip and tail penetration and air entrainment rate are primary parts to carry out the concept, and the spray tip penetration should include the spray evolution after EOI as the discussion in Section 1. Thereby, a new zero-dimensional spray penetration calculation was developed, and the air entrainment was calculated based on the developed model.

Spray penetration including the spray evolution after EOI
Indeed, because of the termination of injection, the spray is transformed from steady state into deceleration state that leads to the spray propagation different from that in the steady state.
Recently, entrainment wave theory was proposed by Musculus [18] based on analytical research on the turbulent jets behaviors during the deceleration state.This theory describes the propagation of the increased entrainment region from nozzle exit to jet downstream after EOI, which also presents the propagation of information of the fuel injection termination.Once the entrainment wave arrives at the tip of the jet, the jet penetration shifts from a relation of square root dependence on time to a relation of fourth root dependence on time, the former relation is widely proved in steady jet and the later relation is same as the result from the experimental research on unsteady turbulent jets [17].The entrainment wave was also observed in diesel jets, and a one-dimensional discrete model derived in [19] was successful to predict the entrainment wave in diesel jets and the penetration of diesel jets even after the entrainment wave reaches the jet tip.These indicate that the one-dimensional discrete spray model in [19] has potential to predict the diesel spray propagation after EOI.Thus, a zerodimensional model of diesel spray propagation was developed referring to this one-dimensional discrete spray model to involve the spray information of penetration and air entrainment rate after EOI in this study.To conduct this work, the modeling method of the one-dimensional discrete model is only used for the treatment of the spray tip part, and the same assumptions as in the one-dimensional discrete spray model were employed for the zero-dimensional spray model.To simplify the spray modeling, the spray is treated to be formed with a constant injection rate.
According to the entrainment wave phenomenon, the spray propagation process can be divided into two stages.As shown in Figure 4, the time when the entrainment wave front arrives at the spray tip is defined as transition timing, and the diesel spray tip penetration calculation is separated in two stages by the transition timing.During injection, the spray is in steady state; after the EOI, the entrainment wave propagates from nozzle exit to the spray tip, and the spray tip keeps in the steady state before the transition timing.After the transition timing, the total spray transforms into decelerating state.In this spray propagation process, the spray tip maintains in steady state before transition timing and then decelerates from the steady state after transition timing.Thus, the spray tip penetration is calculated in two stages, before and after transition timing.In the case of constant injection rate, the transition timing is twice of injection duration.The spray tip penetration (S tip ) can be expressed as where u tip is the spray tip velocity.Based on the momentum theory, the spray tip velocity is equal to the ratio of the momentum flux integrated over the tip cross-sectional area to the mixture mass flow rate integrated over the tip cross-sectional area as Eq. ( 19): is expressed as follows: Developments in Combustion Technology where ρ tip is the average density in the tip cross-section, A tip is the tip cross-sectional area, ū tip is the average velocity over the tip cross-section, and β is a factor that is introduced to consider the effects of axial velocity and fuel volume fraction profiles over the spray tip cross-sectional area on the momentum flux and fuel mass flow rate integrated over the spray tip cross-sectional area, the specific derivation of β can be found in the reference [19].The value of β is ranged from 1.0 for a uniform profile to 2.0 for a fully developed spray.
Meanwhile tip is obtained as Substituting Eqs. ( 19) and ( 20) into ( 21), the spray tip velocity is represented as follows: Before the transition timing, the spray tip is in steady state, so that the momentum flux and fuel mass flow rate integrated over the tip cross-sectional area are constant, and the ratio of them is equal to the fuel velocity at the nozzle exit (u f ) due to the constant injection rate.
where is equal to the injection rate of fuel mass .
Therefore, after substituting Eq. ( 6) into ( 9), the ū tip can be derived as where f, tip is the average fuel volume fraction over cross-section at the spray tip, as well as the average density in the tip cross-section (ρ tip ) can be calculated by the following equation according to [29], So far the spray tip penetration (S tip ) can be obtained by Eq. ( 12) cannot be solved as a continuous equation.Thus, the discrete method is used to calculate the spray tip penetration.The spray tip penetration at any time of t is obtained based on the spray tip penetration and the average velocity over the tip cross-section at the last time step (t - Δt) as Eq. ( 13) And the average velocity over the tip cross-section at time t is obtained as After the transition timing, the total momentum flux and fuel mass flow rate over the tip crosssection decelerate from the values of steady state because the information of fuel injection termination already arrived at the spray tip [18].However, Eq. ( 9) is valid also for the period after the transition timing, which had been demonstrated in [31].Thus, the same form as Eq. ( 10) are available for the average velocity over the tip cross-section, where f, tip and ρ tip are replaced by those after the transition timing, f, tip, atr and ρ tip, atr .
f,tip,atr tip,atr tip,atr , where ρ tip, atr is calculated by Eq. ( 11), in which f, tip is replaced by f, tip, atr .The spray tip penetration after transition timing can be calculated as Developments in Combustion Technology where t tr is the transition timing.
To obtain f, tip, atr , air volume and fuel volume in spray can be treated as separate regions as "air" region and "fuel" region shown in the figure, respectively.Before the transition timing, the fuel mass flow rate at the tip is constant, and the volumetric flow rate at the tip, is equal to that at the nozzle exit.However, after the transition timing, the fuel volumetric flow rate decreases from its initial value in the steady state because the information of fuel injection termination already arrived at the spray tip [19].To describe this process, the volumetric flow rate, was assumed to be shared by fuel and air for the period after the transition timing.Based on this assumption, part of air flows into "fuel" region from "air" region and mixes with fuel, so that the "fuel" region for the steady state is replaced by the "fuel + air" region after transition timing, and then the mixture of "fuel + air" region flows out through the tip with the volumetric flow rate of Considering the "fuel + air" region at time "t" as the control volume, the mixture flows out of the control volume from tip cross-section with the volumetric flow rate of Because of the fluid continuity, air will flow into the control volume with the same volumetric flow rate as compensation.Therefore, at any time "t" after transition timing, the total volume of air flow into "fuel + air" region is For simplicity, it is assumed that the air and fuel mix with each other immediately in the "fuel + air" region.Then, the fuel volume fraction in the fluid flowing out of the tip is calculated by the following equation where m f (t) is mass of fuel in the spray, t is the time from injection start, and C d is a model constant.Finally, X f,tip,atr is derived as (18) Whence the spray tip penetration after the transition timing is able to be obtained, and the method for solving Eq. ( 16) is same as that before transition timing.

Spray tail penetration
After EOI, the mixture near the nozzle becomes very lean due to the termination of the fuel supply.Therefore, it is assumed that the part of the spray near the nozzle acts as ambient air zone and a spray tail exists.Knowing the position of the spray tail, the start time of the interaction with the subsequent spray can be determined.
In this study, the spray tail position is determined as a cross-section of a spray where 10% of total fuel is contained up to the nozzle.It can be imagined that the tip of the spray formed by this 10% fuel equates the total spray tail.In other words, the tail penetration of the spray containing total injected fuel can be calculated by the tip penetration of the spray formed by an injection whose quantity is equal to 10% of total injection quantity.Figure 5 shows the examples of calculated tail penetrations compared to the tip penetrations for the injection quantity of 2.5 and 25 mg cases.

Spray-to-spray interaction modeling
Another primary respect of simulation for diesel combustion with multiple-injection is to describe the interaction between the sprays from sequent injections.Figure 6 shows the sprayto-spray interaction, the interaction occurs after the arrival of the second spray tip at the first spray tail, and the entrainment from the first spray to the second spray is involved to represent the interaction.
The entrainment behavior can be considered as that the spray entrains the ambient gas through the spray boundary that has a conical surface, therefore the entrainment rate can be represent as the product of air density, entrainment velocity over the spray boundary surface, and the area of spray boundary surface.After the second spray tip touches the first spray tail, the entrainment area of the second spray is divided into two parts by the first spray tail as shown in Figure 6.The ratio of the entrainment rate of these two parts (R e ) can be obtained as Eq. ( 19) with assumptions as the air and the first spray have the same density and the entrainment velocity uniformly distribute over the spray boundary surface in the ambient air zone and the first spray zone, respectively: where A air and A spray are the areas of the spray boundary surface in the ambient air zone and the first spray zone, respectively.C e is a coefficient given to describe the difference between the entrainment velocity over the spray boundary surface in the ambient air zone and the first spray zone.Based on Eq. ( 19) and the total entrainment rate of the second spray, it is able to obtain the amount of the first spray mixture entrained into the second spray.Meanwhile, the swirl effects on the spray-to-spray interaction cannot be ignored.The swirl flow does not only decrease the spray penetration [32] but also deviates the spray path [23].
As that shown in Figure 7, the swirl effects causes the second spray tip to overtake the first spray tail and tip earlier than the case without swirl flow effects, and reduce the overlap region between the first and second sprays.Thus, based on the assumption that the entrainment rate is proportional to the spray surface area, the ratio (R e ) of the entrainment rate from the ambient gas and that from the first spray have to be recalculated as follows: where A upt and A blt are the spray boundary surface areas up and below the first spray tail, respectively, and R SA is the ratio of surface area of the second spray in the first spray over A blt .
To calculate A upt and A blt , the reduction of spray penetration in injection direction by swirl flow are considered to determine the position at which the second spray tip touches the first spray tail.The ratio (R MJ ) between the momentum from injected fuel (M j ) and total momentum in the spray, which includes the momentums from the injected fuel (M j ) and the entrained gas (M s ), is introduced as a factor for the penetration.The specific expressions are as follows: where and u j are the mass flow rate and the velocity of fuel at the nozzle exit, is entrainment rate with swirl flow effect, which is calculated according to the equation proposed by Kau et al. [33], and ū s is the average swirl velocity over the total spray penetration.The spray penetration (S SE ) in the injection direction with swirl flow effect can be represented as R MJ S tip , where S tip is the penetration without swirl flow effect.When the S SE of the second spray tip is larger than that of the first spray tail, the second spray tip is treated to arrive at the first spray tail.
To calculate R SA , the spray deflection by swirl flow was considered in a simple manner as shown in Figure 8.The shadow part expresses the overlap between the two sprays.The deviation angles θ 1 and θ 2 , and the spray spreading angles α 1 and α 2 are assumed to have small values for this simplification.In this way, R SA is approximately proportional to the ratio of the angle between the first spray windward and the second spray leeward over the second spray spreading angle as in Eq. ( 24) Figure 8. Simplification of sprays interaction with swirl flow effect [34].
The ratio (R MS ) of the momentum from entrained gas over the total momentum in the spray, which is calculated by Eq. ( 25), is assumed to represent the degree of spray deflection instead of θ 1 and θ 2 in Eq. ( 24).

Wall impingement effects on the turbulent mixing
Adjusting the injection time earlier than top dead center (TDC) is often used to realize the PCCI or LTC, especially for the multiple-injection case in which the pilot injection time is usually advanced to the middle even at the early stage of compression stroke.The early injection timing makes the spray flow into the squish region and impinges on the cylinder liner or piston top as shown in Figure 9. Due to the low temperature of the walls and/or the adherence of fuel on the piston top surface, oxidation reaction and mixing in the mixture are attenuated [35].To involve such effects in the stochastic combustion model, the reduction of fuel-air turbulent mixing rate was considered according to the volume ratio of the spray flowed into the squish area and total spray.The temperature effect mentioned above was not considered; however, the overall oxidation reaction rate is lowered as a result of reduced mixing rate.As shown in Figure 9, if the spray tip cross-section impinges on the bowl lip edge, the spray can be divided into two parts, squish part and bowl part.Thus the volume ratio between squish part and incremental spray (R sq, inc ) can be calculated as follows: up sq,inc up bl ( ) , ( ) where A up and A bl are the areas of the cross-sectional area at impinging timing over and below the bowl lip edge, respectively, and the C(θ) is a function of the angle between the spray central line and cylinder head, and it is used to describe the ratio between spray spreading velocity in squish region and bowl region.This function is selected as cot θ in this study, because it represents the ratio of the horizontal and vertical components of average spray tip velocity when the piston top is treated as a horizontal area.Sequentially, the volume ratio of spray flowed into squish region and total spray (R sq ) can be obtained as follows where t im is the impinging timing and V spray is the total spray volume.And the stochastic collision frequency (ω) that represents the turbulent mixing rate is calculated as follows: where C m is a constant, ω 0 is the collision frequency of free spray, and G j is the total turbulence energy generated by injection.G j is used with minus power (x), because the larger turbulence energy generated by injection causes the stronger mixing (larger ω).C m and x were selected as 6.2 and -0.2 respectively, which are calibrated by experimental heat release rates under different injection time cases.

Wall impingement effects on the air entrainment
Figure 10 shows abstracted diagram of the situation of spray impinging on the wall in the cylinder.The upper spray represents the initial stage of wall impingement.In this stage, it is reasonable to consider that the air entrainment is enhanced because the surface area of spray is enlarged by wall impingement [36][37][38].Thus, to improve the stochastic combustion model, a constant (C ETRM ) is given to multiply the entrainment rate of free spray after the wall impinging.In addition, the interaction between adjacent sprays is an important factor on air entrainment decrease, the primary reason can be considered that the adjacent sprays overlap after wall impingement as shown in Figure 10 (below sprays) decreases the entrainment area of spray, thereby the entrainment rate is suppressed.In order to involve these effects in a simple way, the ratio of the total spray volume and the chamber volume was used to represent the intensity of the interaction between adjacent sprays, and the interaction effect on entrainment rate is introduced as follows: where is the air entrainment rate affected by spray-volume increase, is the original air entrainment rate, n hole is the number of nozzle holes, V cyl is the volume in cylinder, and C vr is a constant to fit the experiment data.
C ETRM and C vr are set as 1.5 and 0.8, respectively, based on the comparison between calculated and experimental heat release rates.

Emission models 6.1. NOx model
Based on the assumption that most of the NOx is NO, the production of NO during combustion process is computed by the extended Zeldovich mechanism [39].Meanwhile, the NO normally produces at a high rate in the mixture with equivalence ratio around 1.0 accompanying with high temperature.Thus, the NO concentration is estimated in the fluid elements whose temperature is over 1200 K with the equilibrium species including C, CO, CO 2 , O 2 , O, OH, H 2 , H, H 2 O, N 2 , and NO.

Soot model
The soot model refers to the Moss's soot model [40].Moss's soot model is a semi-empirical soot model derived based on laminar diffusion flame.The soot particles inception and coagulation are considered for the calculation of soot particles number density, and the soot particles surface growth and oxidation are calculated to obtain the soot volume fraction.The soot oxidation rates per unit area by O 2 and hydroxyl radical (OH) are introduced to calculate the soot oxidation rate in this model.The soot oxidation rate per unit area by O 2 (R SO ) is calculated using the Nagle and Strickland-Constable (NSC) model [41].The soot oxidation rate per unit area by OH (R SOH ) is considered by referring Neoh's equation [42].

Test engine
The test engine is a water-cooled single-cylinder four-stroke-cycle direct-injection diesel engine equipped common-rail injection system.The standard specifications are given in Table 1.All of the experiments were performed at thermally steady states of the engine at a fixed speed of 1500 rpm, an inlet coolant temperature of 80°C, and a lubricating oil temperature of 80°C.The intake pressure was 0.1 MPa and the intake temperature was 35°C.Exhaust back pressure valve was fully open.The fuel was JIS No. 2 diesel fuel (density at 15°C = 820.2kg/m 3 and cetane index = 54.7).The averaged in-cylinder pressure of 50 cycles was used to calculate the heat release rate, which was measured using a piezoelectric pressure transducer (Kistler 6052A).

Calculation conditions
With the aim for validating the simulation of combustion with multistage injection, at first, the pilot/main two-stage injection strategies were conducted in experiment and calculation, and the calculation ran from intake valve closure (IVC) of -145°ATDC to exhaust valve open (EVO) of 125°ATDC.Since the start of main injection is normally set near TDC in the pilot/main injection strategy, the main injection timing (θ main ) was fixed at 1°ATDC, and pilot injection timing and quantity were varied.And the engine operating condition was selected as high load, for which the indicated mean effective pressure (IMEP) was set to 1.01 MPa, in order to observe the variation of combustion phases including premixed combustion and mixingcontrolled combustion affected by pilot injection condition.The main experimental conditions are listed in Table 2.

Pressure and heat release rate
To valid the combustion model, the in-cylinder pressures and the heat release rates were calculated and compared to the experimental data at first.

Emissions calculation
Based on the good agreement in the pressures and the heat release rates, the NOx and soot emissions were calculated.Figure 14 shows the NOx (left) and soot (right) emissions against the pilot injection timing.The results reveal that the NOx calculation obtains the similar emission level and the variation caused by pilot injection conditions changing with the measured data.Regarding to the soot emission, the soot emissions level obtained by the model is comparable with the measured data, and the calculated soot emissions are able to reproduce the increase with the increase in pilot injection quantity at fixed pilot injection timings, which is observed in the experimental data.Although the tendency of soot emissions variation with pilot injection timing retarded does not completely coincide with that of every experiment with different pilot injection quantity, the general tendency can be captured correctly.

Conclusion
In this study, a stochastic combustion model was introduced to develop a phenomenological combustion model for modern diesel engines.In order to be able to describe the combustion process of advanced combustion mode, the spray propagations after EOI, spray-to-spray interaction, swirl effects, and wall impingement effects were modeled based on their physical phenomena in appropriate ways for the stochastic combustion model.Then, the developed combustion model was validated based on the experimental data from a single-cylinder diesel engine with pilot/main two-stage injection.The results revealed that the model is able to accurately predict the combustion of the diesel engine with pilot/main two-stage injection, and reasonable prediction of NOx and soot emissions can be obtained by this model.Specific conclusions are as follows: (1) A zero-dimensional spray propagation model was derived.The spray model is able to predict the spray evolution including spray tip penetration and overall air entrainment after EOI.It is capable to catch the spray propagation tendency after EOI as that well recognized by the fundamental study.Thus, based on introducing this model, the accuracy of the combustion model can be improved especially for the short injection duration case in which the ignition is later than the EOI.As for the multiple-injection case, the thermodynamic states of the gas in the cylinder before the later injection start can be predict more realistically.
(2) Presumed spray tail was proposed for zero-dimensional spray model, which allows to take account of the fast dilution of the mixture near the injector tip after EOI, which resulted by the terminated fuel supply and increased air entrainment rate near the injector tip after EOI.And thanks to the presumed spray tail, the interaction between sprays from sequent injections was formulated by the rate of entrainment into the later spray from the former spray and/or the surrounding gas.Meanwhile, spray deflection by swirl flow cannot be neglected when imitating the interaction between sprays from sequent injection.And the interaction with swirl flow effect was described based on a simple geometrical consideration and formulated using the momentum from injected fuel and entrained gas.In this way, the stochastic combustion model can be used for simulation of combustion with multistage injection.
(3) To consider the wall impingement effect in early injection timing case, after wall impinging, the volume ratio of the spray flowed into the squish area, and total spray was introduced to reduce the fuel-air mixing rate that can involve the combustion chamber shape effects in the phenomenological combustion model.
(4) Air entrainment rate enhancement caused by wall impingement was considered by a factor with air entrainment rate of free spray during initial stage of wall impingement.And the reduction of main spray entrainment rate by interaction between adjacent sprays was formulated by the ratio between spray volume and chamber volume to balance the effect of enhanced air entrainment rate in the late combustion period.According to the results, effects of wall impingement and adjacent spray interactions on the entrainment rate are helpful to well predict the heat release rate in initial and late combustion periods, respectively.

Figure 2
Figure2shows the shadowgraph images of two-stage injection in a constant volume combustion chamber.In Figure2, the part "A" is the first spray flame region.After the second injection starts, the first spray flame region and the second spray region are separated at first as shown in the left graph.Once the second spray tip penetration is long enough to overtake the boundary of the first spray flame region, the second spray subsequently penetrates into the first spray flame region and causes the interaction between the first spray flame and the second spray as shown in the right graph.

Figure 2 .
Figure 2. Shadowgraph images of two-stage injection in constant volume combustion chamber.

Figure 5 .
Figure 5. Calculated spray tip and tail penetrations for the injection quantity of 2.5 and 25 mg cases.

Figure 10 .
Figure 10.Spray propagation along the wall.

Figures 11 -
13 show the in-cylinder pressures and heat release rates for different pilot injection timings and quantities.The simulation results obtain the similar levels of pressures and heat release rates to those of the experiments.It is also observed that the model is able to capture the tendencies of the pressure and heat release rate when varying the pilot injection timing and injection quantity.

Figure 11 .
Figure 11.Effects of pilot injection quantity on the in-cylinder pressure and the heat release rate (θ pilot = -9°ATDC).

Figure 12 .
Figure 12.Effects of pilot injection quantity on the in-cylinder pressure and the heat release rate (θ pilot = -19°ATDC).

Figure 13 .
Figure 13.Effects of pilot injection quantity on the in-cylinder pressure and the heat release rate (θ pilot = -24°ATDC).

Table 1 .
Standard specifications of test engine.