Transported PDF Modeling of Ethanol Spray in Hot-Diluted Coflow Flame

This paper presents a numerical modeling study of one ethanol spray flame from the Delft Spray-in-Hot-Coflow (DSHC) database, which has been used to study Moderate or Intense Low-oxygen Dilution (MILD) combustion of liquid fuels (Correia Rodrigues et al. Combust. Flame 162(3), 759–773, 2015). A “Lagrangian-Lagrangian” approach is adopted where both the joint velocity-scalar Probability Density Function (PDF) for the continuous phase and the joint PDF of droplet properties are modeled and solved. The evolution of the gas phase composition is described by a Flamelet Generated Manifold (FGM) and the interaction by exchange with the mean (IEM) micro-mixing model. Effects of finite conductivity on droplet heating and evaporation are accounted for. The inlet boundary conditions starting in the dilute spray region are obtained from the available experimental data together with the results of a calculation of the spray including the dense region using ANSYS Fluent 15. A method is developed to determine a good estimation for the initial droplet temperature. The inclusion of the “1/3” rule for droplet evaporation and dispersion models is shown to be very important. The current modeling approach is capable of accurately predicting main properties, including mean velocity, droplet mean diameter and number density. The gas temperature is under-predicted in the region where the enthalpy loss due to droplet evaporation is important. The flame structure analysis reveals the existence of two heat release regions, respectively having the characteristics of a premixed and a diffusion flame. The experimental and modeled temperature PDFs are compared, highlighting the capabilities and limitations of the proposed model.


Introduction
Spray combustion is widely utilized in various engineering applications, such as industrial furnaces and propulsion systems. To achieve higher efficiency while minimizing the pollutant emission, novel combustion technologies are demanded. Among others, the MILD (Moderate or Intense Low-oxygen Dilution) combustion is demonstrated to be a promising technology [6,58]. By dilution of the reactants with the recirculated reaction products, the flame peak temperature is substantially reduced, resulting in a low production of NOx. Delft Spray-in-Hot-Coflow (DSHC) burner has been used to study fundamental aspects of MILD oxidation of bio-derived liquid fuels [10]. A first numerical study of an ethanol flame from DSHC database with the transported PDF method is reported in this paper.
Modeling of turbulent spray combustion is particularly challenging, because many physical and chemical processes including turbulence, atomization, evaporation, combustion and radiative heat transfer are involved and interact with each other [28]. These phenomena and processes have to be modeled in a proper way in the sense that the main physical characteristics have to be accounted for, but with a reasonable computational cost. For simplicity, many spray combustion studies have been carried out in the regime of dilute spray [9,24], and this approach is also deployed in this study.
The transported probability density function (PDF) method [46] has proven to be a powerful closure method for modeling turbulent reactive flow [23,47]. PDF method has been applied to spray combustion since the 1990s, and is still an active research area. Naud et al. [36][37][38][39] developed a hybrid finite-volume/transported PDF method, and systematically studied the modeling issues in the context of Lagrangian-Lagrangian approach. Beishuizen [4] studied the particle-turbulence interaction of turbulent spray flames. Ge and Gutheil [20] proposed a joint mixture fraction-enthalpy PDF method for modeling turbulent spray combustion, and recently developed a joint velocity-mixture fraction PDF model [21]. Bhattacharjee and Haworth [5] compared well-stirred reactor (WSR) model with PDF method for n-heptane and n-dodecane spray flames under engine conditions, concluding that the PDF method performs better due to the fact that the turbulent fluctuations have been taken into account. Pei et al. [44,45] simulated diesel engine combustion using composition PDF coupled with Reynolds-averaged k-ε model. Recently, attention has been paid to the FDF (filtered density function) method in conjunction with Large Eddy Simulation. Heye et al. [24] modeled the Sydney ethanol spray flame with a LES/FDF approach. Jones et al. [29] modeled a gas turbine combustor using LES with sub-grid probability density function to account for the sub-grid turbulence-chemistry interaction. Despite the contribution of these works, the micro-mixing model, small scale droplet models, as well as the combustion models for the application of PDF methods to spray combustion are still open issues.
The most outstanding advantage of transported PDF methods is that the mean reaction source term appears as a closed term. However, the direct use of detailed chemistry is computationally very expensive. Proper chemistry reduction is required for affordable yet accurate models. This is normally accomplished by either using reduced chemical mechanisms, or employing tabulated chemistry methods. In this study the Flamelet Generated Manifold model [42] is used, which falls in the second group of chemistry reduction methods. In the FGM model, the scalars, such as temperature, species mass fractions, density or progress variable source term, are stored in a lookup table as a function of a few independent variables-usually the mixture fraction and a progress variable. The scalars are then retrieved from the pre-built lookup table during turbulent combustion simulation according to the value of the modeled independent variables. The influence of turbulence fluctuations on the mean properties is accounted for through the joint PDF of the independent variables. For many applications, the shape of the PDFs of independent variables are simply assumed before simulation. For instance, mixture fraction is often presumed as a β-function with shape parameters determined by its mean and variance values. However, many studies [28,55] already pointed out that due to the presence of droplet evaporation, the β shape PDF is no longer valid for mixture fraction in spray combustion. For the PDF of progress variable, even more ambiguities exist, both β-function and δ-function have been reported in the literatures [9,14], and further studies are required. Alternatively, in this study, the transport equation of the joint PDF of gas phase properties is directly modeled and solved, such that the turbulence-chemistry interaction is considered in a more precise manner.
The application of the classical flamelet model to spray combustion was first made by Hollmann and Gutheil [25,26] to simulate a methanol/air diffusion flame and extended by them to a formulation using spray flamelets [20,22]. However, the high dimensionality of the spray flamelets makes them difficult to tabulate and use. A novel two dimensional spray flamelet, using mixture fraction and droplet evaporation rate as independent variables, was recently proposed by Olguin and Gutheil [43]. However, the shape of the PDF for droplet evaporation rate is still an open issue for the application of this model with presumed PDF methods. Chrigui et al. [8,9] applied a FGM model to Large Eddy Simulation of ethanol spray combustion, using presumed β-PDF for mixture fraction and δ-function for progress variable. To the authors' knowledge, application of FGM model for spray combustion in the context of transported PDF has not yet been reported.
The purpose of this paper is twofold: on the one hand, to validate the transported PDF modeling approach with FGM for spray combustion; and on the other, to increase the understanding of MILD spray combustion. This paper is structured as follows: some background information for the current study is firstly given in this section, followed by mathematical modeling approaches for continuous and dispersed phase respectively in Sections 2 and 3. The experimental configuration and numerical setup of the target flame are presented in Section 4. The results are discussed in Section 5, focussing on: the role of the "1/3" rule; the influence of the droplet initial temperature and of the evaporation model; the comparison with experimental data; the flame structure; and the temperature PDF. Major conclusions and future study are then emphasized in Section 6.

Transported PDF hybrid finite volume/particle method
The in-house RANS/transported PDF code "PDFD", which has already been successfully used in gaseous flames, evaporating sprays and coal combustion [4,14,36,52], is used in this study (where RANS refers here and in the following to averaged Navier-Stokes equations in the sense of Favre averages). The continuous phase is described by the joint velocity-scalar PDF, where the scalars considered are the FGM independent variables: mixture fraction Z and progress variable Y c . As detailed in the next section, the dispersed phase is described by the joint PDF of droplet position, velocity, temperature, diameter, and the gaseous properties "seen" by the droplets. To cope with the high-dimensionality, the joint PDFs are solved by Monte Carlo particle methods. In contrast with the Eulerian-Lagrangian approach, both the gas phase and the dispersed phase evolution are defined by Lagrangian equations, and we are therefore considering a "Lagrangian-Lagrangian" approach.
To overcome the bias error due to the limited number of computational particles in the Monte Carlo (MC) method, the continuous phase mean velocities and Reynolds stresses are calculated using a Finite-Volume (FV) method, in which the Favre-averaged Navier Stokes equations are solved [38]. Similar approaches were also used by Ge et al. [21], Bhattacharjee et al. [5] and Anand et al. [2]. Note that in our case, special attention is paid to the consistency between the particle velocity evolution and the Reynolds-stress and scalarflux second moment closures used in the RANS model [37,39]. Figure 1 illustrates the computational algorithm. The FV submodel provides the mean velocity and its gradient, mean pressure gradient, Reynolds stresses and mean turbulent dissipation rate to the MC part.
The fluctuating velocity increment of the gas phase particles is determined by the generalized Langevin model (GLM)-more specifically, by the variable C 0 formulation of the GLM presented in [39]-in correspondence with the LRR-IPM Reynolds-stress model used for the modeling of the pressure strain correlation in the FV part. The evolution of the gas phase composition is described by the FGM and the interaction by exchange with the mean (IEM) micro-mixing model: where the mixture fraction Z is not a conserved scalar in spray combustion and its source term S Z corresponds to the mass coming from droplet evaporation. On the other hand, by definition of the progress variable Y c , as we will see in Eq. 3, its source term does not include effects of the evaporated fuel andω FGM Y c corresponds to the chemical reaction source term, shown in Fig. 3, retrieved from the FGM lookup table as a function of the independent variables. The IEM mixing model reads: where ω φ = C φ /k is the modeled scalar variance decay frequency, with C φ the mixing model constant, set to C φ = 2. The mean turbulent kinetic energy (k) and turbulent dissipation rate (ε) are provided by the FV part. The two-way coupling source terms due to the drag force appearing in the momentum and Reynolds stress equations [4] are not included here since these effects proved to be small for the dilute spray considered. However, we do include the mean mass source terms in the mean continuity and momentum equations.

Combustion model
In flamelet-based models, the multi-dimensional turbulent flame is considered as a set of 1D flamelets. The 1D flamelets are characterized by different controlling parameters to describe the local variations of the real flame. For the FGM model, the controlling parameters are mixture fraction Z and a progress variable Y c . Different methods exist for the construction of the 2D FGM lookup table [41]. A commonly used one is to first calculate different steady flamelet equations with scalar dissipation rate increasing from a very small value to the extinguished value. These steady flamelets are then mapped in (Z, Y c )-space together with the unsteady extinguishing flamelet solution [8]. Another approach is to solve the unsteady process of a 1D diffusion flame in physical space from pure mixing until the steady flame is established. The flamelet solution at different time is then transformed into (Z, Y c )-space. Compared to the "extinguishing" FGM generated by the first method, the second method generates an "auto-igniting" FGM table, which is therefore more suitable to describe the auto-ignition process of the DSHC flame. Note that so-called unsteady flamelet/progress variable approaches have also been proposed where both ideas are combined. In that case, igniting and extinguishing flamelets are resolved for different scalar dissipation rates, as for instance presented recently in [27,40]. However, such approaches require one additional control parameter, and a 3D lookup table needs to be considered in (Z, χ , Y c )-space, with χ the scalar dissipation rate. Auto-ignition lookup tables can also be constructed by solving Perfectly Stirred Reactors (PSR) [17] or combining with premixed flamelets [16].
The 2D auto-igniting FGM table used in this study is generated with the CHEM1D code developed at the Eindhoven University of Technology [7]. The conterflow diffusion flame is solved in physical space at unity Lewis number where the boundary conditions are specified such that the specified strain rate is 100 s −1 and the fuel corresponds to pure C 2 H 5 OH vapor at its boiling temperature T boil = 351 K. Compared to PSR, this configuration takes into account the diffusion during ignition, therefore it is physically more representative of the reality. The chosen strain rate of 100 s −1 (corresponding to a rather low stoichiometric scalar dissipation rate of 1.85 s −1 ) is consistent with the observation made in [27], for a lifted methane/air flame with similar coflow conditions as the current case, that ignition happens in relatively low scalar dissipation rate regions. The detailed ethanol high temperature oxidation mechanism containing 57 species and 383 reactions by Marinov [33] is employed. The ignition process is illustrated by the temperature profiles in mixture fraction space with increasing time, as shown in Fig. 2. The progress variable in this study is defined as a weighted sum of species mass fractions as follows: where W k refers to the molecular weight of species k. The progress variable source terṁ ω FGM c (Z, Y c ) is shown in Fig. 3 in mixture fraction and normalized progress variable space. In principle, droplet evaporation influences the gas phase flamelet structure by consuming energy and adding fuel vapor. As explained in Section 2.1, the vaporized fuel is accounted for by the source term for mixture fraction. Many different approaches have been proposed to take into account the enthalpy loss effect of droplet evaporation. For example, spray flamelets [26], using total enthalpy as progress variable [3], partially premixed flamelet method [19] or generating FGM table by solving new spray flamelets equations as derived in [31] and [43]. In this study, as a first step of the model validation, an adiabatic gaseous FGM table generated by using pure fuel vapor as fuel stream is employed, and enthalpy loss is not included in this 2D FGM table. The influence of this scheme on the results will be further discussed in Section 5.

Model for the Dispersed Phase and Phase Interactions
Accurate prediction of the droplet dispersion and evaporation in turbulent flows are of crucial importance since they are considered rate limiting processes in modeling dilute spray combustion [28]. In the proposed Lagrangian modeling of the dispersed phase PDF, the evolution of the properties of the stochastic particles representing posible realizations of the turbulent spray is chosen to follow Lagrangian models for single droplets [36]. The dispersed phase stochastic particles will be denoted "parcels" in the following.

Droplet motion
For practical spray combustion, the droplet drag force and the gravitational force are dominant compared to other forces, for instance the buoyancy force and Basset force. Therefore the particle momentum equation, also known as BBO (Basset-Boussinesq-Ossen) [50], is greatly simplified: The droplet relaxation time τ p is determined by: where ρ p and ρ g respectively refer to the liquid droplet and gas-phase densities, and D p is the droplet diameter. The drag coefficient C D is given by the Schiller-Naumann semiempirical correlation: with the particle Reynolds number: In the above equations and in the following, the subscript "seen" refers to the properties seen by the droplets (the undisturbed fluid flow properties at the position of the droplet center, which modelling is detailed in Section 3.6). The subscripts "p" and "g" respectively refer to droplet and gas-phase properties.

Parabolic temperature profile
For droplet heating and evaporation processes, a variety of models with different levels of complexity exist [1,30,34,49,51]. Among them a widely used one is the "infinite conductivity model", in which the temperature distribution inside the droplet is assumed uniform. However the finite conductivity effects become important when the droplet heating process is fast as is the case in the hot-diluted coflow condition of this study. Fully resolving the heat conduction problem inside the droplets greatly increases the computational cost.
By checking droplet temperature distribution, T p (r, t), one observes that, except at the very beginning of heating, the shape of the curve T p (r) looks close to a parabola [51]. Hence, the finite rate heat conduction process is taken into account by assuming that the temperature profile between the droplet surface and its center is a parabola [15]: where T cntr is the temperature at the droplet center (at r = 0) and T surf is the droplet surface temperature (at r = R p , with R p the droplet radius). If we generalize the derivation of [15] for evaporating droplets as done in [13], we can consider the volume averaged droplet temperature T p , defined as: After considering the boundary condition at the droplet surface as in [15], its time evolution then reads: where Nu is the Nusselt number, given in Eq. 16, and where C p refers to the specific heat capacity and λ to the thermal conductivity. The subscript "liq" refers to liquid properties, while the subscript "m" refers to the gas-phase properties close to the droplet surface (considered in Section 3.5). The temperature T ∞ , far away from the droplet surface, is interpreted here as the seen temperature T seen . The surface temperature is obtained as: Note that the temperature at the center of the droplet can be obtained from Eq. 9 as: T cntr = 5T p − 3T surf /2.

Abramzon and Sirignano evaporation model
In the above equations, we introduced the notation for the temperature evolution due to evaporation such that: whereṁ p is the evaporation rate and L v (T surf ) is the latent heat of vaporization at droplet surface temperature T surf . Abramzon and Sirignano [1] proposed modified Nusselt and Sherwood numbers, Nu * , Sh * , deduced from "film theory", to account for the boundary layer thickening effect by the Stefan flow. They express the evaporation rate based on mass diffusion of vapor or based on heat transfer as: where ρ m and λ m are the average density and thermal conductivity of the film gas mixture (i.e. at temperature T m and composition Y m , as explained in Section 3.5), D vap is the binary diffusion coefficient of pure vapor in the film gas mixture and C p,vap is the specific heat of pure vapor at temperature T m . The modified Sherwood and Nusselt numbers read: with the correction factors, representing the relative change of the boundary layer thickness due to the Stefan flow, approximated as: B M and B T are the Spalding mass and heat transfer numbers. Sh and Nu are the Sherwood number and Nusselt number for flow around a solid sphere, and can be evaluated by the well known Ranz and Marshall correlation: where Sc m and Pr m are the Schmidt and Prandtl number respectively: In the Abramzon and Sirignano model, the evaporation rate is obtained using the second expression in Eq. 13, based on heat transfer. Using the equality with the first expression based on mass transfer, B T , Nu * and F T are obtained iteratively from: The Spalding mass transfer number B M is calculated from Eqs. 19 and 20: where Y vap refers to the mass fraction of fuel vapor, and W is the mean molecular weight of the seen mixture. The mole fraction of fuel vapor X surf vap at the droplet surface is calculated from the Clausius-Clapeyron equation (assuming that the equilibrium between gas and droplet surface has been reached at each time step):

Heat transfer
With the Abramzon and Sirignano evaporation model, the heat transfer to the droplet can be obtained from the definition of the Spalding heat transfer number B T as: where the first term on the right hand side represents the droplet temperature change due to the convective heat transfer, and the second term the droplet temperature decrease due to the evaporation. In the present model, we would rather keep the standard expression for convective heat transfer and express the heat transfer between the droplet and the gas-phase as:q as already used in the derivation of the parabolic temperature profile (implying the first term in the bracket in Eq. 10 and the last term in Eq. 11). However, as explained before, since the gas-phase modeling is based on an adiabatic igniting FGM, described by mixture fraction and progress variable only (no enthalpy heat loss included), the contribution ofq drop is not considered here.

Evaluation of film properties and influence of internal recirculation
Equations 10 to 20 completely describe the droplet heating and evaporation process under the assumptions of parabolic temperature profile, phase equilibrium and taking into account the effects of Stefan flow. The presence of evaporation creates large normal gradients of composition and temperature near the droplet surface. The gas-phase properties of the film mixture, denoted by subscript "m", are all evaluated at an intermediate temperature T m and composition Y m : Typically α = 1/3 is widely accepted for spray combustion simulation. It is the well known "1/3 rule" [18]. However, many other possibilities exist, for example α = 1 means to directly use the "seen" gas properties. The influence of whether employing the "1/3 rule" or not will be further analyzed in Section 5. The slip velocity between droplet and gas phase may induce internal circulation inside large droplets, which may enhance the droplet internal heat transfer. It is possible to take into account the inner recirculation by replacing the liquid thermal conductivity λ liq in Eq. 10 with the so-called "effective conductivity" [1]: where κ = 1.86 + 0.86 tanh 2.245 log 10 (Pe L /30) , with Pe L the Peclet number for droplet interior. However, the droplet size in this study is relative small, and the droplet internal recirculation is assumed negligible. Therefore, the physical thermal conductivity, λ liq , is used.

Seen properties and distribution of vaporized fuel
The seen gas velocity is described by the modified Generalized Langevin Model (GLM) proposed and implemented in PDFD by Naud [37]. This model generalizes the model of Minier [35] for the seen velocity based on the Simplified Langevin Model, ensuring that the modeling is consistent with a given Reynolds stress and scalar flux second moment closure in the limit of tracer particles by generalizing the derivation of Naud [39]. In this case, the chosen seen velocity GLM is consistent with the LRR-IPM Reynolds stress model used in the RANS submodel. As indicated by Eq. 23, the droplet evaporation rate is very sensitive to the way the seen gas temperature T seen and composition Y seen are evaluated. For most applications, these seen properties are evaluated using the mean values in RANS simulations [5] or filtered quantities in LES [8,9,11], neglecting the influence of turbulent fluctuations. Transported PDF methods potentially allow for more refined modeling of the seen composition and temperature. Jones and co-workers [29] developed a stochastic Markov model to account for the effect of the SGS fluctuations of gas-phase reactive scalars on droplet dispersion and evaporation in LES. In this study, as proposed by Naud and De Meester [13,14], for all stochastic droplets the seen composition and temperature are obtained from the FGM at given values (Z , Y c ). Every given characteristic time T L,seen (corresponding to a seen Lagrangian scalar correlation time based on the seen velocity GLM), these values are sampled from a given gas-phase stochastic particle present in the same computational cell. The gas-phase particle is chosen to be the one with enthalpy h(Z , Y c ) the closest to the enthalpy of the saturated mixture, evaluated for every droplet at T surf and Y surf . During the characteristic time T L,seen , this seen composition (Z , Y c ) evolves according to the IEM mixing model with the scalar variance decay frequency ω φ = 1/T L,seen . The seen properties are then obtained from the FGM table as Another important issue is how to distribute the fuel vapor generated by evaporation. Although more advanced methods exist [36], for simplicity, the vaporized fuel vapor is evenly distributed to all the gas phase particles present in the computational cell, similar scheme has been used in [5,24]. This approach, essentially reduces the variance of the gaseous properties.

Experimental setup
The schematic of DSHC burner is shown in Fig. 4. The liquid fuel (ethanol) is injected into the hot-diluted coflow by a pressure swirl atomizer. The hot-diluted coflow is generated by a secondary burner matrix working on Dutch Natural Gas (DNG) and air, to mimic the diluted air by recirculated combustion products in a large scale MILD combustion furnace. The air-to-DNG ratio together with the effects of the two perforated plates and the pipe length dictate the coflow temperature, oxygen concentration and turbulence levels. Comprehensive laser diagnostic measurements, including Laser Doppler Anemometry (LDA), Phase Doppler Anemometry (PDA) and Coherent Anti-Stokes Raman Scattering (CARS) have been conducted. Gas phase velocity components, temperature and O 2 volume fraction have been measured along the radial direction at coflow exit (Z = 0 mm). Droplets properties (velocity, diameter, number concentration and mass flux) have been measured along the radial direction at different axial locations (Z = 8, 10, 12, 15, 20, 30, 35, 40, 45 mm). Gas phase temperature has been measured with CARS technique along the radial direction at different axial locations (Z = 15, 20, 30, 40, 50, 60 mm). For further details about the DSHC burner and the database, the readers are referred to [10]. The experimental data will be compared with simulation results for validation purpose. In this paper we simulate one of the ethanol spray in hot-diluted coflow cases, namely the case "H II " in [10]. The main parameters for this case are described in Table 1. Subscript "cf " refers to the property of coflow, and the last column is the mass flow rate of the liquid fuel at the injector. One may notice that the coflow temperature and O 2 mole fraction shown in Table 1 are different from those reported in [10]. This is because the whole profile including the boundary layer was considered for the averaging in [10], while a representative condition at the plateau of the coflow profile is used in this study.

Computational domain
As mentioned in Section 1, in this study we restrict ourselves to the modeling of dilute spray combustion, no attempt is made on the modeling of film breakup and droplet formation during the atomization process. The droplet collisions, coalescence and agglomeration are also ignored. The inlet boundary is chosen such that it is sufficiently far from the atomizer tip to avoid the dense spray region but below the region where the ignition starts. In this case, the axial location Z = 8 mm is chosen as the inlet boundary. This is also the first axial location where the dispersed phase properties were measured. As the flame is statistically axisymmetrical, a 2D axisymmetrical simulation is conducted. The computational domain is indicated by the yellow rectangle in Fig. 5.

Gas phase boundary conditions
Due to the presence of droplets, LDA measurements for the gas phase velocity were only conducted at the coflow exit (Z = 0 mm), which can not be directly used for this study.
Although the PDA results at Z = 8 mm (inlet boundary of this simulation) for the small droplets (D < 6 μm) can be used as gas velocity, they are only available at limited points due to the availability of the small droplets, see the symbols in Fig. 6. These limited data points do not provide enough information for the accurate assignment of the inlet boundary of the dilute spray combustion. Furthermore, because the spray is issued into a hot coflow, some liquid fuel has already vaporized before Z = 8 mm, and possibly some reaction has already started before the computational inlet boundary. Therefore, accurate mixture fraction and progress variable radial profiles have to be provided as boundary conditions in order to correctly predict the dilute spray combustion behavior downstream. However, the mixture fraction and progress variable are not directly available from experimental measurements. Nevertheless, the necessary properties are available at the coflow exit (Z = 0 mm). A simulation of the entire spray flame, starting at Z = 0 mm, including the spray atomization process was conducted with ANSYS Fluent 15.0 to derive reliable boundary conditions at Z = 8 mm for the present study, following the approach reported in [32]. In the Fluent simulation, the pressureswirl atomizer is modeled with Linearized Instability Sheet Atomization (LISA) model. Turbulence is modeled by Reynolds Stress Model (RSM). And the turbulence-chemistry interaction, to be consistent with the current study, is also modeled by FGM model but with presumed shape PDF method. β-PDF is used for both mixture fraction and progress variable. To examine the reliability of the boundary conditions provided by the Fluent simulation, not only the results at Z = 8 mm but also at other axial locations are compared with

Dispersed phase boundary conditions
The dispersed phase boundary conditions are assigned based on the experimental data for each droplet size class. Available dispersed phase boundary conditions include the droplet velocity components and their variance, the dispersed phase mass flow rate and the fraction of mass flow rate for each droplet class. The uncertainty exists in the droplet temperature at Z = 8 mm. As described in Section 3.2, the modeled droplet temperature is determined by its own initial state as well as the experienced surrounding gas phase conditions. Therefore, it is difficult to accurately calculate the droplet temperature at Z = 8 mm. However, single droplet simulations showed that due to the presence of high temperature coflow, the droplet temperature rises rapidly after injection. The Fluent simulation, in which the finite conductivity model is used, also confirms that the small droplets temperature at Z = 8 mm are close to the boiling temperature. Here, two sets of droplet temperature boundary conditions will be tested to examine the sensitivity of the results on the droplet initial temperature. As shown in Table 2, droplet temperatures are assigned differently depending on their size to account for their different thermal innertia.

Number of particles
With the time averaging and particle split/merge algorithms as described in [36], the required number of computational particles per cell is dramatically reduced. Two cases with 20 and 50 gas phase computational particles per cell are tested, no clear difference on the results is observed. According to the experiment, droplets with a diameter larger than 70 μm are rarely detected. In the simulation, the droplets have been divided into 7 classes ranging from 0 to 70 μm. For each droplet class, 10 nominal computational parcels per cell are used, which means a total of 70 dispersed phase parcels per cell. Within each hybrid outer iteration, 500 finite volume iterations, 10 gas particle Monte Carlo iterations and 10 droplet parcel Monte Carlo iterations are conducted respectively. More than 1000 hybrid iterations have been carried out for each case in this study to reach converged results. This is similar to the coal combustion modeling presented in [52], and in the same way, the use of a local time stepping algorithm also helped to increase the convergence rate since larger particle time steps can be used in regions with small velocities.

Cases
In the subsequent sections, the uncertainty of the boundary conditions, namely the droplet initial temperature at Z = 8 mm, as well as the influence of different sub-models will be discussed. Four cases with different boundary conditions and sub-models will be analyzed, see Table 3. Cases "C" and "D" impose relatively low temperatures ("T p 1" in Table 2) as droplet boundary condition, while in the other two cases, "A" and "B", droplets are set to temperatures that are closer to the boiling temperature. In cases "A", "B" and "C", the parabolic temperature profile model is used in contrast with the infinite conductivity model used in case "D". The influence of the "1/3" rule is studied by setting α in Eq. 23 to 1 in case "A" and to 1/3 in the other cases.

Role of the "1/3" rule
In Section 3, we saw that the film properties, λ m , ρ m , and D vap are widely involved in the droplet sub-models. Due to the large normal gradients of composition and temperature near the droplet surface created by droplet evaporation, it is, theoretically, not straightforward to define a proper condition at which these properties should be evaluated. As discussed in Section 3.5, the empirical "1/3" rule is widely accepted for spray combustion. Figures 7 and 8 respectively show the predicted droplet mean axial and radial velocity components for the four cases considered in addition to experimental data. The results are  plotted in a matrix of subplots with each subplot representing a certain droplet size class at a certain axial location. The droplet size increases from left to right in the matrix and the axial location increases from bottom to top. The difference between Case "A" and the other cases is only related to the way the gas phase properties used in the droplet evaporation and dispersion models are evaluated. For cases "B", "C" and "D" where the "1/3 rule" is applied, the gas phase mixture properties are evaluated at state (T m , Y m ) obtained from Eq. 23 with weighting factor α = 1/3. In case "A", α = 1, the "seen" gas properties at (T seen , Y seen ) are directly used. These properties eventually affect the dispersed phase behavior via the droplet dispersion and evaporation models as described in Section 3. For the sake of clarity, in this Section we only compare results of cases "A" and "B". It can be observed that case "A" considerably under-predicts the droplet mean velocity while results of case "B" are in better agreement with experimental data. It is especially clear for large droplets at high axial locations (the up-right part of the subplots matrix). This means that the direct usage of "seen" gas properties leads to an over-prediction of the droplet velocity decay rate. Droplet velocity quickly reduces when traveling downstream. In many spray applications, part of the droplets evaporate in a low temperature environment, where the difference between T surf and T seen is relative small, the averaging of gas phase properties may not make a significant difference. However, for conditions like the droplet-flame interaction and spray in hot-diluted coflow flame, where the conditions between droplet surface and surrounding gas are considerably different, this becomes very important. For example, in the current study, the gas temperature on droplet surface, T surf , is approximately equal to the droplet boiling temperature, T boil = 351 K. However, the "seen" gas temperature, T seen , could vary in a wide range from fuel vapor temperature, ∼351 K, to flame temperature, above 2000 K. The gas viscosity evaluated at (T seen , Y seen ) is in general higher than when obtained with the "1/3" rule, according to Sutherland's law [53]. This in turn results in a shorter droplet relaxation time, see Eqs. 5 to 7. Droplets in this case tend to relax to the gas phase velocity more quickly, as demonstrated by the results of case "A" in Figs. 7 and 8.
This example also illustrates the importance of droplet "seen" property model, since the "1/3 rule" is averaging the gas properties between the droplet surface and the "seen" gas. If the "seen" gas properties are not properly sampled, the results could also be different. Hereafter, we only show the simulation results obtained with the "1/3 rule", namely cases "B", "C" and "D".

Influence of droplet temperature boundary condition
Since the droplet temperature boundary condition is the main uncertainty for the modeling of this flame, cases with two different sets of droplet temperature boundary conditions will be analyzed. In this section, we focus on the results predicted by cases "B" and "C".
As already shown in Figs. 7 and 8, the mean droplet velocity predicted by cases "B" and "C" do not exhibit considerable difference. With both sets of droplet initial temperature, the droplet Reynolds stresses u 2 p are over-predicted in the near axis region for all the droplet classes, as depicted in Fig. 9. The reason for the over-prediction will be discussed later. It is fair to say that the different droplet temperature boundary conditions do not lead to significant differences in the droplet mean velocities and their higher moments (Figs. 7 to 10).
However, as expected, the droplet Sauter Mean Diameter (SMD), shown in Fig. 11, do unveil the differences caused by droplet temperature boundary condition. In general, both case "B" and "C" predict correct trend and magnitude of SMD, indicating the good performance of the droplet evaporation and dispersion model. Nonetheless, at the spray outer edge, the SMD is over-predicted by the "low" initial droplet temperature in case "C" and under-predicted by the "high" initial droplet temperature in case "B". These results mean that a better initial temperature for large droplet should be in between the value in "T p 1" and "T p 2" in Table 2. The predicted results for small value of SMD are almost identical. This is because the temperature of these small droplets rise very quickly to the so-called "wet bulb temperature" after injection. The initial temperature of the small droplets therefore has smaller influence on the results. The same trend is observed for the droplet number density in Fig. 12. The results for small droplets predicted by these two cases are quite similar to each other, except that in case "B", the droplets have a wider radial distribution. The number density for large droplets is lower in case "B" than that in case "C", indicating that less large droplets survive due to faster evaporation.
The influence of droplet initial temperature on gas phase velocity and Reynolds stresses is almost negligible, as shown in Figs. 13 and 14. A slight difference is observed on the gas phase temperature in Fig. 13 between case "B" and "C", related to different mixture fraction source terms in both cases. However, this difference does not include the effects of heat loss due to droplet evaporation, which is not considered here.

Influence of the evaporation model
In this section, the two cases, "C" and "D", which differ only in the droplet evaporation model will be analyzed. As explained in Section 3.2, the parabolic temperature profile model assumes that the temperature distribution inside the droplet is a parabola from surface to center, while the infinite conductivity model assumes isothermal conditions inside droplet. The parabolic temperature profile model can be categorized as "conduction limit" model, while the infinite conductivity model is also called "fast mixing" model. Figure 15  shows the average difference between droplet surface and center temperature predicted by the parabolic temperature profile model as a function of axial location. It illustrates that after injection, the difference between the temperature at the droplet surface and at the center varies differently for the different droplet sizes. For small droplets, the difference continuously decreases as the droplets travel downstream. However, for droplets larger than 20 μm, the droplet surface temperature first quickly increases and the surface-center temperature difference initially becomes larger. After some time, the heat conducts to the center and the difference decreases. The surface-center temperature difference in the infinite conductivity model case of course remains zero during the droplet lifetime. It was demonstrated by Dombrovsky and Sazhin [15] that the temporal evolution of the droplet temperature predicted by the parabolic temperature profile model closely resembles the one obtained by solving the heat conduction problem inside the droplets except in a very short time period at the beginning. In Fig. 11, it is observed that the droplet SMD predicted using the parabolic temperature profile matches well the experimental data, while it is obviously under-predicted at the near axis region when using the infinite conductivity model. In case "D", where the infinite conductivity model is employed, the droplet number density of small droplets, Fig. 12, is over-predicted in the near axis region, and under-predicted in the spray outer edge. This is because the isothermal assumption in the infinite conductivity model results in a relatively lower droplet surface temperature, which determines the evaporation process. The droplet evaporation is therefore under-predicted by the infinite conductivity model compared to the parabolic temperature profile model. The small droplets in the near axis region do not vanish as fast as in Case "C". Similarly, the relatively slow evaporation of large droplets at the spray outer region does not generate as many small droplets as in Case "C". That also explains the different behavior of SMD predicted by these two models as described above.
The different evaporation models also have a noticeable impact on the gas phase mean velocity via two-way coupling. Since the droplet velocity in case "D" has been underpredicted, the gas phase "feels" less acceleration from the evaporated mass. Therefore, the  gas phase velocity in case "D" is lower than that in the case "C" as well as the experimental data. This is especially clear at high axial locations, see Fig. 13. Because of the fast evaporation of small droplets in the near axis region, combustible mixture has been quickly formed in this region. The ignition and combustion in this case therefore occur at a smaller radial location than that in case "C", as depicted in Fig. 16. The temperature fluctuation in case "D" is also stronger in the near axis region.

Comparison with experimental data
From the previous analysis, it is already clear that the "1/3" rule should be applied in the simulation; the two different droplet temperature boundary conditions do not produce significant difference in the results, but a correct initial temperature for large droplets is important; and the parabolic temperature profile model outperforms the infinite conductivity model. We now proceed to an overall comparison of the predicted results with the experimental data. In this section, the results from case "C", which has been shown to be the best among others, will be examined over experimental data to show the achievements and incapabilities of the current modeling approach.  First of all, as can be seen from Figs. 7 and 8, the droplet mean velocity components for all the droplet classes have been accurately reproduced at all the axial locations that have been checked. Note that the experimental data is only available within a certain radial range, this is related to the cone-shape spray structure and the number of sample data required to have a converged statistics in experiment. Out of the main spray region, the number density of the droplet becomes very small. In order to have a converged statistics, a sufficiently long measuring time is demanded, which is not always convenient for a spray combustion experiment. Therefore, only data in the main spray region, which is also considered of most interest, are available [10]. The same argument also holds for the simulation. Even though the iteration averaging algrithm is employed, the number of numerical parcels for droplets is only large enough to have converged statistics in the main spray region. Out of this region, the results include large statistical errors. For the droplet Reynolds stresses, u 2 p and v 2 p , shown in Figs. 9 and 10, satisfying agreement with experimental data has been achieved in the main spray region. Significant over-prediction is observed in the near axis region and at large radius. These discrepancies of droplet Reynolds stresses are believed to be caused by the un-converged statistics due to the insufficient number of droplet parcels. The predicted v 2 p is in general in better agreement with experiment data than u 2 p , see Fig. 10. However, the discrepancies for large droplets are obvious, and this may be related to the fact that the number density of large droplet is one to two orders of magnitude lower than that of the small droplets, as can be seen in Fig. 12.
Droplet Sauter Mean Diameter (SMD) is fairly well predicted by case "C", as shown in Fig. 11. As discussed before, the over-prediction at the outer edge of spray is attributed to the low temperature boundary condition for large droplets. The trend of the droplet number density is also well captured, Fig. 12. The slight over-prediction of the droplet number density may be due to the fact that some sample data were rejected during experiment when more than two droplets are present in the probe volume [10]. The above discussion demonstrates that the droplet sub-models, including the droplet dispersion model, evaporation model and seen gas model, are capable to accurately reproduce the dispersed phase behavior of the DSHC flame.  In Fig. 13 the modeled mean axial and radial gas-phase velocity components are compared to the experimental data obtained by PDA from the smallest droplets as tracers. Good agreement is obtained for the mean axial velocity. The discrepancy between the modeled and measured mean radial velocity could be explained by the fact that using small droplets as tracers in the near field corresponds to a non-uniform seeding of the flow. The small droplets, mainly moving outwards, do not represent in an unbiased manner the complete gas phase. A better, unbiased, mean velocity measurement would include more samples of coflow gas without droplets, with predominantly inwards velocity. From axial location Z = 45 mm, not enough tracer droplets are available due to the evaporation, therefore no experimental data are available thereafter. The predicted gas phase Reynolds stress u 2 is in reasonable agreement with measured data. The v 2 is under-predicted at low axial location, which can also be explained by the bias in the gas-phase radial velocity measurements based on a non-uniform small droplet seeding.
Gas phase temperature in the spray region has been measured with CARS technique [10]. A reasonable agreement with experimental data is obtained in case "C". The flame D [10,20]um D [20,30]um D [30,40]um D [40,50]um D[50,60]um Fig. 15 Mean difference between droplet surface and center temperature predicted by the "parabolic profile model (case "C") at varying axial locations peak temperature as well as the flame width are correctly predicted. The radial position of the peak temperature is slightly shifted towards the center. This may mainly be caused by the mixture fraction profile specified at the inlet boundary. As explained in [32], the Fluent simulation used for providing inlet boundary information for this study predicts a smaller spray angle compared to the experiment. Consequently the distribution of the free vapor is also narrower. Close to the spray axis, an opposite temperature trend is predicted. The simulation shows a small temperature peak in the center, while the temperature progressively decreases towards the center in the experiment. This is because near the spray axis many small droplets exist, and considerable gas phase enthalpy loss happens in this region due to the fast evaporation of small droplets. The enthalpy loss, however, can not be accounted for by the 2D adiabatic FGM table used in the current study. As a consequence, the temperature has been over-predicted in this region. This problem can be solved by including enthalpy loss as another independent variable of the FGM table, namely using a non-adiabatic FGM table.
It is noticeable in Fig. 16, that the gas phase RMS temperature is significantly underpredicted. It is also noticed that the predicted RMS temperature is somehow systematically about 100 K lower than the experimental counterpart. This systematical difference in RMS temperature is believed to be mainly caused by the exclusion of temperature fluctuation at the inlet boundary. From the experimental data at radial position 20 mm to 50 mm in the first subplot in Fig. 16, we clearly see that the 100 K deviation is equivalent to the temperature fluctuation in the coflow. The influence of temperature fluctuation at the inlet boundary will be further discussed in Section 5.6.

Flame structure
The fuel vapor released from the droplets makes the spray combustion show characteristics of both premixed and diffusion flames. The Flame index is a commonly used parameter to identify the premixed and non-premixed flame in the context of DNS or LES [3,48]. Although the applicability of the flame index in the RANS simulation is still questionable, the "averaged" flame index can still reveal some major characteristics of the flame structure of the DSHC flame. The flame index is normally defined as the product of the spatial gradients of fuel and oxidizer mass fraction: In most literature, only the fuel from the inlet or from the droplet vapor is used in the calculation of the flame index. This is acceptable in the cases where the global one step reaction is used or the chemical reaction is infinitely fast [55]. In cases where the intermediate reactions are important, the intermediate species such as H 2 and CO should also been considered as fuel. Therefore, in this study, we define the following three flame indices: Their contour plots are shown in Fig. 17. The contour plots of mean mass fractions of C 2 H 5 OH, O 2 , H 2 , CO and H 2 O are shown Fig. 18 to better illustrate the structure of this flame. In all these plots two iso-surfaces of gas temperature 1400 K and 1600 K are imposed to indicate the heat release region. According to the flame indices and the species mass fraction contour plots, the DSHC flame can be divided into the following five regions: the spray core region, the inner flame region, the center region, the outer flame region and the coflow region, they are indicated on the contour plot of flame indices Fig. 17. The spray core region is the region near the axis, where the small droplet accumulate and quickly evaporate. Next to it is the inner flame region, further outward is the center region, and the outer flame region. The coflow region is sitting most outside, and is not involved in any combustion.
In the spray core region, the O 2 entrained from the coflow before the lift-off of the flame coexists with the fuel vaporized from small droplets. The O 2 has a negative radial gradient while the fuel vapor has a positive gradient due to the cone shape spray generated by the pressure-swirl atomizer. Therefore in this region, the F I C 2 H 5 OH is negative, but it does not necessarily mean that a diffusion flame exists here, because the temperature in this region is relatively low. This negative F I C 2 H 5 OH is a sign of local production of fuel vapor, which is quite different from a gaseous flame case.
In the inner flame region, presence of a secondary flame front is revealed by a local peak in the RMS temperature profile in Fig. 16 (e.g. r = 8 mm at Z = 30 mm) and the positive F I C 2 H 5 OH in Fig. 17. This inner flame front is produced by the reaction of the premixed C 2 H 5 OH and O 2 coming from the spray core region. However, due to the dilution of the coflow, the O 2 in the spray core region is not enough to fully convert the C 2 H 5 OH to the final reaction products, CO 2 and H 2 O. The temperature in this diluted rich region is relatively low and this is particularly suitable for pyrolytic or reforming stages leading to H 2 and CO [6], see Fig. 18. Therefore, we see also a negative F I CO and F I H 2 in the inner flame region. Again, the negative F I CO and F I H 2 do not correspond to a diffusion flame here, but to a local generation of CO and H 2 . Supporting this observation, the high concentration of CO and H 2 have also been found in the furnace or well-stirred reactor (WSR) under  [56] and numerical [12] studies. In the DSHC experimental results, the inner flame is only detectable at axial locations from 30 mm onward, because the heat release from this inner flame front, has been partially compensated by the heat loss due to droplet evaporation. In the present simulations in which the heat loss is not considered, the inner flame front is clearly revealed by the shoulder in the mean temperature profile at the inner branch of the spray flame. The phenomenon that the inner flame front locates at more or less the same radial location with varying axial locations is also correctly captured by the simulation. The reason for this is that the inner flame front is mainly formed by the vaporized fuel vapor, and therefore behaves like a gaseous jet flame. The spreading rate is much lower than that of the outer flame region which is aligned with the spray outer edge.
Further outward, there is the center region. The peak of the temperature radial profile appears in this region, but it is mainly due to the heat released from the inner and outer flame regions rather than local combustion. Reaction barely happens in this region, because of the lack of O 2 . The intermediate species will continuously react with the O 2 coming from the coflow at the outer flame region, which is also the main heat release region. In the outer flame region, a diffusion-like flame is formed, see the negative F I CO and F I H 2 in this region in Fig. 17. The final combustion products are also formed in this region, see the contour plot of mean H 2 O mass fraction in Fig. 18. Finally, almost no droplets survive beyond the outer flame region, and the flow composition in this region remains the same as that of the coflow.
As discussed above, under the hot-diluted coflow conditions the incomplete oxidation of fuel and the further reaction of the intermediate species are spatially separated. This results in a more distributed heat release region and a lower peak temperature, which are precisely desired in the MILD combustion technology for the reduction of NOx formation [57,58]. Figure 19 compares the gas phase temperature PDF obtained from experiment and simulation respectively. The solid red line represents the predicted temperature PDF and the solid black line is the temperature PDF obtained from the CARS measurements. To guarantee the convergence of the PDF, at least 1000 CARS samples are considered for each point in the experiment. The Monte Carlo particles are sampled over at least 5000 Langrangian time steps. The PDF comparison are carried out for 9 different locations as illustrated by the black dots on the temperature contour plot in Fig. 18.

Temperature PDF
The under-prediction of the temperature variance discussed in Section 5.4 is clearly shown here by means of narrower temperature PDF distributions. The reason for this can be explained by the PDF at point "C", which is located in the coflow region and is only influenced by the coflow inlet boundary condition. The predicted temperature PDF of point "C" is a δ-function at the experimental mean temperature, whereas the experimental data show a Gaussian distribution within the range of 1000 K-1700 K. This means that the mean temperature of the coflow has been exactly represented by the FGM table, but no temperature fluctuation at the inlet is considered. The exclusion of the temperature fluctuation at the flame inlet boundary consequently reduces the temperature fluctuations in the whole simulation domain. The temperature PDF prediction could be improved if the temperature fluctuations at the inlet would be included. The implementation of the inlet boundary temperature fluctuations can be done provided that the non-adiabatic FGM table is applied. Temperature fluctuations can be imposed by supplying gas phase Monte Carlo particles with fluctuating enthaly loss/gain values. This will be done in a future study.
The adiabatic FGM table also has a significant influence on the temperature PDF. The gas phase energy loss due to droplet evaporation can not be properly considered by the current adiabatic FGM table. This is the direct reason for the absence of the lower temperature tail of the PDF for points A, E, F, H, I, J , and the shift of the PDF towards the high temperature at these points. For points B and F, however, the whole PDF have been shifted towards the low temperature. This may be explained by the differential diffusion of H 2 . In reality, the H 2 diffuses faster than other larger molecular species. In this case, the H 2 formed in the inner flame region diffuses faster outward than CO, therefore at the edge of the outer flame region, some pure H 2 combustion may occur. This can be proved by the very high temperature samples (>2200 K), which is higher than the adiabatic C 2 H 5 OH flame temperature under the diluted coflow condition, at points B and F. Correct account of differential diffusion of H 2 itself is still a big challenge for combustion modeling, and is out of the scope of current study. Readers interested in this topic are referred to [54]. Besides the two discrepancies mentioned above, comparison of the predicted temperature PDF with experimental data demonstrated the ability of transported PDF method.

Conclusion
In this paper, we reported a first numerical investigation of the Delft Spray-in-Hot-Coflow flame with transported PDF method and FGM model. The in-house hybrid finitevolume/particle transported PDF code "PDFD" is used for the simulation. A Lagrangian-Lagrangian approach is employed to describe the two-phase turbulent flow field in the dilute spray region. The continuous phase is described by a joint velocity-scalar PDF, and the dispersed phase is described by a joint PDF of droplet position, velocity, temperature, diameter, and gaseous properties "seen" by the droplet. The inlet boundary conditions for the dilute spray simulation was fulfilled by the results from a complete spray simulation with ANSYS Fluent 15.0 and the available experimental data. An uncertainty exists on the boundary condition for the droplet initial temperature. Two sets of different droplet temperature boundary conditions were tested. It was shown that the initial temperature of small droplets has a negligible effect on the results, while the initial temperature of large droplets considerably influences the droplet SMD and number density downstream. For future study, the temperature boundary condition is suggested to take values in between the two sets of temperature tested in this study. A parabolic temperature profile model was used to describe the droplet heating and evaporation process. Its performance was compared with the widely used infinite conductivity model. By being able to take into account the finite heat conduction process inside droplets, the former shows superior performance over the latter in terms of better agreement with experimental data on both gas and dispersed phase properties. The influence of the "1/3" averaging rule was examined in detail. The results showed that the "1/3" averaging rule has a considerable influence on the droplet dispersion behavior, especially when the difference between the temperature at the droplet surface and in the surrounding environment is large: for example for droplet-flame interaction and in the spray-in-hot-coflow situations. A direct use of the "seen" gas properties for the evolution of the gas phase properties in the evaporation and dispersion models under these conditions leads to a too fast decay of droplet velocity.
The current modeling approach was further validated by comparing the results predicted by case "C", where the parabolic temperature profile model and the "1/3" rule were applied, with available experimental data. Droplet velocity, Sauter Mean Diameter and number density are all in good agreement with measured data, showing that the spray sub-models including the evaporation and dispersion models used in this study are suitable for modeling the DSHC flame. Gas phase velocity also matches well the available experimental data. The mean temperature was predicted with reasonable agreement with experimental data. However it was over-predicted in the central part of the spray where intensive droplet evaporation happens. The reason for this is that the heat loss due to evaporation can not the properly considered by the adiabatic FGM table used in this study. The temperature variance has been significantly under-predicted due to the exclusion of the temperature fluctuations at the inlet boundary. The use of an adiabatic FGM table and leaving out the H 2 differential diffusion effects both contribute to the large discrepancy between the experimental and modeled temperature PDF.
The structure of the modeled flame was analyzed. It was found that two heat release regions exist. The inner flame region is formed by the reaction between the fuel vapor and the coflow oxidizer entrained from the flame base below lift-off height, and mainly behaves like a premixed flame. The outer flame region, which is the main heat release region, is mainly created by the further oxidation of the intermediate species formed at the inner flame region, and shows characteristics of diffusion flame. The incomplete oxidation of fuel and the further reaction of the intermediate species are spatially separated, resulting in a more distributed heat release region and a lower peak temperature, which are desired conditions in the MILD combustion technology for the reduction of NOx formation.
The transported PDF method together with the FGM model presented in this study show promising performance on modeling dilute spray combustion. The two phase flow field as well as the flame structure could be properly reproduced. However, further improvements are required in order to have a more precise prediction of the temperature PDF. These improvements include taking into account the enthalpy effects, considering differential diffusion effects, and including the temperature fluctuations at the inlet boundary.