Large Eddy Simulation of CO2 diluted oxy-fuel spray flames

We report results of a computational study of oxy-fuel spray jet ﬂames. An experimental database on ﬂames of ethanol burning in a coﬂow of a O 2 –CO 2 mixture, created at CORIA (Rouen, France), is used for model validation (Cléon et al., 2015). Depending on the coﬂow composition and velocity the ﬂames in these experiments start at nozzle (type A), just above the tip of the liquid sheet (type B) or are lifted (type C) and the challenge is to predict their structure and the transitions between them. The two-phase ﬂow ﬁeld is solved with an Eulerian–Lagrangian approach, with gas phase turbulence solved by Large Eddy Simulation (LES). The turbulence-chemistry interaction is accounted for using the Flamelet Generated Manifolds (FGM) method. The primary breakup process of the liquid fuel is neglected in the current study; instead droplets are directly injected at the location of the atomizer exit at the boundary of the simulation domain. It is found that for the type C ﬂame, which is stabilized far downstream the dense region, some major features are successfully captured, e.g. the gas phase velocity ﬁeld and ﬂame structure. The ﬂame lift-off height of type B ﬂame is over-predicted. The type A ﬂame, where the ﬂame stabilizes inside the liquid sheet, cannot be described well by the current simulation model. A detailed analysis of the droplet properties along Lagrangian tracks has been carried out in order to explain the predicted ﬂame structure and discuss the agreement with experiment. This analysis shows that differences in predicted ﬂame structure are well-explained by the combined effects of droplet heating, dispersion and evaporation as function of droplet size. It is concluded that a possible reason for the difﬁculty to predict the type A and B ﬂames is that strong atomization-combustion interaction exists in these ﬂames, modifying the droplet formation process. This suggests that atomization-combustion interaction should be taken into account in future study of these ﬂame types.


Introduction
Since in many combustion processes the main source for NOx formation is the oxidation at high temperature of the N 2 contained in air, a natural suggestion to reduce or eliminate the NOx emission, has been to separate N 2 and O 2 and use enriched air or pure O 2 as oxidiser. This is the concept of oxy-fuel combustion. This combustion technology has many advantages. In case of 100% pure oxygen and in absence of fuel bound nitrogen, NOx emission is no longer an issue. Second, the flue gas of this combustion process is predominantly CO 2 and H 2 O, by separating water vapor through cooling or compression, a CO 2 stream for carbon capture and sequestration (CCS) is available. Such a zero emission combustion system, is particularly appealing.
However, oxy-fuel combustion also faces several challenges. First of all, N 2 separation from air with current technologies is energy consuming and expensive. Second, switching to oxy-fuel combustion drastically changes the process conditions. Adiabatic flame temperature of combustion with O 2 is high and the resulting high local heat flux implies a heavy thermal load to the burner. Third, due to the high temperature, a small amount of N 2 remaining after incomplete separation has a large chance to be converted to NOx. Less severe conditions with moderate heat flux and lower emissions can be created by dilution of the O 2 with part of the produced CO 2 . The level of dilution appears as a process variable and research is still needed to find optimal the oxy-fuel combustion technology to be used in practical systems.
Local structure in spray flames can have a variety of types depending on the relative time scales of the process involved. This has been systematically reviewed recently by Sanchez et al. [2]. Detailed numerical simulations reveil the mechanisms leading to the different structures. Reveillon and Vervisch [3] did pioneering work by reveal the dilute spray flame structure using 2D DNS. They reviewed earlier spray flame regime diagrams and presented a new classification based on three dimensionless quantities: the fuel/air equivalence ratio within the core of the spray jet, the mean interdroplet distance to flame thickness ratio, and the evaporation time to flame time ratio. In jet-in-coflow flames these parameters can be influenced by changing fuel injection and coflow conditions. The influence of varying oxygen concentration in the coflow has been the subject of a limited number of studies in the literature. A number of references have addressed the range of oxygen concentrations lower than air. Reddy et al. [4] studied the variation in flame structure experimentally using kerosine as fuel and comparing flame structure as function of fuel injection pressure and coflow composition. Their study includes cases with oxygen percentage in the coflow varying from 21% down to 17%. The database of the Delft spray in coflow flames [5] covers cases with air as coflow and cases with hot diluted coflow with oxygen percentage around 10%. An extensive study on ethanol spray combustion in a coflow consisting of only O 2 and CO 2 , and covering a very wide range of oxygen concentration from 25% to 80% was done at CORIA (CNRS, University of Rouen and INSA of Rouen) and reported by Cléon et al. [1]. The goal of the present work is to report results of a computational study of the CORIA experiments. In the next sections we respectively describe the experimental setup and the simulation method Section 2, analysis of the results Section 3 and conclusions Section 4.

Experimental setup & Simulation detail
In this study we simulate jet-in-coflow flames from the CORIA oxy-fuel spray combustion database [1]. Fig. 1 shows the dimension of the furnace in which the experiment was carried out and also shows the cross section of the computational domain, discussed below. The database concerns a series of flames with different combinations of coflow velocity and CO 2 dilution level of the oxidiser. A parameter a is used to characterize the degree of dilution of O 2 by CO 2 , and is defined as follows: where X denotes the mole fraction. In the experiment, the coflow velocity was changed by varying the coflow exit area with different insert units. In this way the coflow mass flow rate could be kept constant while varying the velocity [1]. Here we consider cases with two different coflow inserts, namely ''insert 95" and ''insert 200", respectively having coflow annulus outer diameter 95 mm and 200 mm and corresponding coflow mean velocity 0.51 m/s amd 0.11 m/s, respectively. For each insert we consider a case with a ¼ 40 and a case with a ¼ 60. An overview of the characteristics of the four case is given in Table 1.
In the experiments three types of flame structure have been observed, differing in the relative distance of the flame base to the atomization region [1]. The ''type A" and ''type B" flames are observed in cases with relatively small a (e.g. 40). The ''type A" flame is anchored at the nozzle by a small conical central flame, while the main flame stabilizes at the tip of the liquid sheet. The type B flame, found at higher coflow velocity, consists only of the main flame and anchors at the tip of the liquid sheet. Finally for larger a, e.g. a ¼ 60, also ''Type C" flame is observed, which stabilizes at far downstream of the dense region.
One of the flames in the database (case a60 À I95) has been simulated by Enjalbert [6] using massively parallel computing emplying the YALES2 solver and using LES with tabulated chemistry. In that simulation the computational domain covered the entire furnace interior, and the computation was done using a mesh with 27 M cells on 1024 processors and using a finer mesh with 215 M cells on 8192 processors. That study reached qualitative agreement of flame structure, but it also made clear that the modeling of the spray inlet conditions for this experiment is an important issue.
The simulation in this study is carried out using the open source CFD package -OpenFOAM [7]. New libraries have been created for the FGM storage and retrieval algorithms and are dynamically linked to a customized solver for spray combustion. The new solver is referred to as ''sprayFGMFoam". This new solver has been successfully applied earlier in the modeling of MILD spray flames from the DSHC dataset, created at Delft University of Technology [8,9]. We use LES with tabulated chemistry (FGM) and the simulations have been performed on 100 processors of Cartesius, the Dutch supercomputer. As a first step study, in this paper we are only interested in the near field structure of the spray flames, therefore a smaller computational domain is adopted, illustrated in Fig. 1. In order to study the influence of the computational domain and mesh resolution, three different meshes have been adopted; details are listed in Table 2. The computational domain in all cases is a 3D cylinder and a hexahedral structured mesh is used, see Fig. 2. The computational domain in the axial direction extends from Z ¼ À20 mm to Z ¼ 250 mm, where Z ¼ 0 mm is the location of atomizer exit. The reason for using this length is that it was observed that the axial gradients of all properties are already very small at Z ¼ 250 mm. Two different diameters of the computational domain have been considered, respectivelly called ''large" and ''small". The diameter of the ''large" domain equals 400 mm, the diameter of the furnace. In the case of the small domain it is 200 mm. The smallest cell size in all three meshes is 0.3 mm, appearing at the injector exit (first cell layer above the inlet). The difference between the ''small" and ''small-fine" meshes is mainly the growth ratio of cell size at downstream. In the ''small-fine" mesh, the cell grows slower streamwisely than that in the ''small" mesh, therefore the former one has a finer resolved region at the reaction zone compared to the latter one. A refinement in the entire domain may be more convincing, however, due to the limitation on the computation resources, this has not been done in the current study. The results of simulation respectively using these three cases will be discussed in Section 3.1.
The transport equations are spatially discretized with a Finite Volume Method (FVM). The convection and Laplacian terms are discretized respectively by second-order accuracy total variation diminishing (TVD) schemes Gauss vanLeer and Gauss vanLeer corrected. Implicit second-order method CrankNicholson is used for the temporal integration. It should be noted that these schemes are highly dissipative, and may leads to under-estimation of turbulent kinetic energy. However, similar schemes have been used in the simulation of Delft Spray-in-Hot-Coflow (DSHC) flames [9], one of which has similar flame structure with those in the current study. Comparison with the DSHC experiment shows that the current numerical approach is able to correctly predict main parameters and the flame structure, which are the focus of the current study. Therefore these highly dissipative numerical schemes are still considered acceptable in the current study. A fixed CFL number 0.5 was used during the simulation.

Turbulence-chemistry interaction
A tabulated chemistry method -Flamelet Generated Manifold (FGM) [10] -along with the Large Eddy Simulation (LES) technique have been employed for the Turbulence-Chemistry Interaction (TCI). The following equations have been solved: Table 1 Descriptions of experimental cases.

Case
Coflow Flame type where q is the density, u i the ith component of the velocity, Z the mixture fraction, l the dynamic viscosity and D the mass diffusivity. Subscript ''t" denotes the turbulent properties.
is the deviatoric part of the strain rate tensor S ij .
s ij is the sub-grid scale (SGS) stresses, and it is closed with the dynamic Smagorinsky model in the current study. S q ; S u i and S Z are respectively the source terms for continuity, momentum and mixture fraction due to existence of evaporating droplets. Their expressions are given in Table 3. In this table, V c is the volume of a computational cell, and N p is the number of droplets represented by a parcel. Y c is the progress variable, and is defined as follows in the present study: where W and Y are the molar mass and mass fraction, respectively. It is related to the scaled progress variable, C, by: where Y min c and Y max c are the minimum and maximum progress variable values, respectively.
The influence of turbulent fluctuations on the local flame structure is accounted for through the joint Probability Density Function (PDF) of the independent variables. In this study a presumed bfunction is used for the PDFs of both mixture fraction and progress variable. A transport equation (Eq. (9)) and an algebraic model (Eq. (10)) have been used for the SGS variances of mixture fraction and progress variable respectively, following the approach in [11].
The last term in Eq. (9) accounts for the creation of mixture fraction variance due to droplet evaporation as suggested by Pera et al. [12]. The model constant value a ¼ 0:5 is used in the current study, following the recommendation of Hollmann and Gutheil [13]. The model constant C v is set to 0:15 according to [14].
To build the FGM table, a laminar counterflow flame is first solved in physical space with the CHEM1D code [15], and then the results are mapped to the mixture fraction space, similar approach was also applied [16,17]. The fuel vapor at room temperature (300 K) is specified as fuel stream of this counterflow flame, and the corresponding coflow condition (details are given in Table 1) is specified as oxidizer stream. The chemical mechanism used for this calculation is the detailed ethanol oxidation mechanism developed by Marinov [18]. The steady state solution of this counterflow flame at a fixed strain rate a is considered as one steady flamelet. And the result of the unsteady evolution of this counterflow flame is considered as unsteady flamelet. The final FGM table used in the current study contains the steady flamelets at different strain rate (the red lines in Fig. 3), from very small to the extinguishing value, and the state of the unsteady flamelet at the extinguishing strain rate (the blue lines in Fig. 3). The flamelet data are tabulated as function of mixture fraction and progress variable. In Figs. 3c and d, the source term of progress variable, Y c , as a function of mixture fraction and scaled progress variable, C, for two a are given. These results clearly show that the dilution of the oxidiser by CO 2 makes the mixture less reactive and also reduces the flame peak temperature.

Dispersed phase modeling
The droplets are injected from the atomizer (Z ¼ 0 mm) using the Conditional Droplet Injection Model (CDIM) proposed by the authors [8]. The droplet initial size distribution is given asRosin-Rammler distribution, but the range of the injection angles and the velocity distribution within that range depends on droplet size. In [8] the model was developed for sprays generated using the Delavan SWB 0.75-30 hollow cone spray nozzle and here it is applied to a spray from the Delavan WDB 0.75-30 nozzle used in the CORIA experiments. No sub-grid dispersion model is used for droplets, due to the very fine grid resolution. Droplets are tracked in a Lagrangian manner, with governing equations for droplet evolution as follows: where U p;i ; T p and m p are the droplet velocity, temperature and mass, respectively. D vap denotes the mass diffusivity of fuel vapor, k the thermal diffusivity, C p;liq the heat capacity of the liquid, g i the gravitational force on ith direction, L v the latent heat for evaporation. The droplet relaxation time s p is determined by: where q p and q 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 semi-empirical correlation: The subscripts ''p" and ''g" respectively refer to droplet and gasphase properties. Subscript ''seen" denotes the gas phase properties ''seen" by the droplets. Subscript ''m" refers to the properties of the film gas mixture and is evaluated according to the ''1/3-rule".
B M is the Spalding mass transfer number and can be calculated as follows: Nusselt number Nu and Sherwood number Sh are used to consider the convective effect on heat and mass transfer, and are calculated according to the well known Ranz and Marshall correlation: where Sc m and Pr m are the Schmidt and Prandtl number respectively. Bird's correction [19] is applied for Nu to account for the reduction of heat transfer due to evaporation: 3. Results and discussion

Influence of Computational domain and grid resolution
As mentioned in Section 2.1, in order to retain a reasonable computational cost, and also in line with the focus of this study -the near field structure of the CORIA flame -a small simulation domain has been adopted. The comparison of the gas phase mean velocity predicted by simulations using different numerical meshes is displayed in Fig. 4. As can be seen from Fig. 4, the pre-dicted gas mean velocity profiles almost overlap with each other, and agree reasonably well with the experimental data.

Double flame structure
Fig. 5 displays the OH, temperature, mixture fraction and O 2 fields on a cross-section of case ''a60 À I95". From these properties, the ''double flame" structure is clear. In the near axis region, a wide region with high values of OH concentration and temperature is present. This region is referred to as the inner flame region. Going outwards, there is another thin region with high values of OH concentration and temperature, and it is called the outer flame region. These findings are consistent with the experimental observations reported in [1]. The mixture fraction field shows that Z reaches its maximum between these two flame regions. The O 2 , on the other hand, has reached its minimum in the same region. In our previous study of the A II case of DSHC flames, which has a similar ''double flame" structure as this case, we have discussed the mechanism of the formation of this inner and outer structure from the point view of combustion, and found that they are created by different species, main fuel or intermediate species, and are of different type, premixed or non-premixed. For more details, readers are referred to [20].

Droplet behavior: a Lagrangian point of view
Since all the fuel that burns in the reaction zones eventually comes from evaporation of droplets, the understanding of the behavior of droplets can be beneficial to further unravel the mechanism of this ''double-flame" structure and the influence of the coflow conditions on it. Therefore, in this section we attempt to analyze the droplet behavior from a Lagrangian point of view.
In the simulation using OpenFOAM, a unique original ID is specified to every injected ''parcel", and it is carried by this parcel throughout its lifetime and is saved at each output time step. Through this original ID, the history of each parcel can be easily traced. Together with the original ID, the saved information for a parcel include: the current position, the original injection position, the current diameter, the original diameter, the current temperature, etc. Since, the current location of a parcel is available, the droplet ''seen" gas phase properties can be obtained by interpolating the gas phase information, e.g. resolved velocity, temperature, mixture fraction, etc., at the droplet location. With this information, a full Lagrangian track of each parcel can be drawn. Note that in the simulation each parcel represents a number of droplets that have identical properties, e.g. location, diameter, velocity, etc. When displaying the Lagrangian tracks', analysis, we represent the properties of any of these droplets.
In Fig. 6a the trajectories of some randomly chosen droplets are shown. Also displayed in this figure is the gas phase mean velocity field, and the mean position of iso-surface Y OH ¼ 0:001, indicating the flame front. It is clearly shown in this figure that the coflow is entrained by the spray, and is accelerated in the central region. Trajectories of some droplets have been quickly changed by the entrained coflow, and go vertically upwards following the gas phase in the central region. But others keep moving along their initial injection direction and more or less remain ballistic motions. The group of droplets that have been blown to the center survive longer than those keeping their initial direction of motion. By comparing the droplet trajectories with the flame front, indicated by the OH iso-surface, three groups of droplets can be identified. The first group contains the droplets that are blown to the center, and enter the inner flame region at small radial distance. Droplets in the second group reach the flame base and vanish there. In the last group, droplets have nearly straight trajectories, and pass through the flame base, and are then completely vaporized before the outer flame region. In order to have more insight in the behavior of these three different droplet groups, and their contributions to the flame structure, we will pick one representative droplet from each group, and analyze them in greater detail in the following. The three selected droplets are labeled as ''P 1 ", ''P 2 " and ''P 3 ", respectively, and are shown in Fig. 6b.
In Fig. 7, information along the trajectories of droplets ''P 1 ", ''P 2 " and ''P 3 " is shown. The droplet was injected at time zero, and information was sampled every 1ms until the droplet is compeletely vaporized. The coordinates of the circular symbols indicate the actual droplet positions at the sample times, other information has been radially shifted on the figure in order to have a clear view. The arrows on the most-left side are the vectors of the instantaneous ''seen" gas velocity, obtained by interpolating the gas phase velocity at the droplet current location. The vectors in the middle indicate the instantaneous droplet velocity. Besides indicating the actual spatial locations of the droplets, the circular symbols also show the droplet diameter (enlarged), via the size of the symbols, and the droplet temperature, via the color of the symbols. On the most-right side, the rectangular symbols give the instantaneous ''seen" gas temperature, with the color of the symbols. The legends with scales of each property are given in Figs. 7 b and 7c. Furthermore, in Fig. 8 quantitative information on droplet and ''seen" gas properties as a function of droplet age has been given. The angle between droplet and ''seen" gas velocity vectors, h U , is defined as: whereŨ p andŨ g are the droplet and ''seen" gas velocity vectors, respectively. kŨk is the length of vectorŨ. First of all, from Fig. 7, we see a clear difference between the histories of the three selected droplets. Trajectories of ''P1" and ''P2" bend within a short distance from the injector, while ''P3" moves along a nearly straight line. From Fig. 8 we see that the major differences between these three droplets are their original size -27 lm, 41 lm and 63 lm, respectively. This can be easily understood by consulting Eqs. (14)- (16), which describe the evolution of droplet velocity. The droplet relaxation time quickly increases with its size. So for droplet ''P 3 ", it takes much longer time to adapt to local gas velocity compared to droplet ''P 1 ". And this is confirmed by the history of magnitude of slip velocity, and angle between the droplet and ''seen" gas velocity vectors, shown in Fig. 8. Both these two properties of ''P 1 " rapidly decay to zero within 5ms, while this time is almost doubled for ''P 3 ". The difference in the relaxation time eventually results in different trajectories.
Besides the momentum exchange between the two phases, droplets also undergo mass and energy exchange with ''seen" gas. And again, quite different histories in these aspects are observed for these three droplets. Figs. 7 and 8 show that the size and temperature for ''P 1 " have kept almost unchanged values for a long time and distance when traveling in the central region. At the last period of ''P 1 ", it experienced a fast increase of temperature and decrease of size. The reason for this is that it entered the inner flame region at the end, and this is evidenced by the rise of the ''seen" gas temperature at the end, shown in Figs. 7 and 8. From Fig. 8, we can also see that the ''seen" gas mixture fraction for ''P 1 " also quickly increased after it had entered the inner flame region. We can infer from these observations that the inner reaction region is at least partially fed by the evaporation of small droplets in the central region. And this is consistent with the findings in [20]: the inner edge of the inner flame region is actually a premixed flame burning mixture created by the evaporation of small droplets in the central region. The premixed-like inner flame zone was also claimed in the experimental study in [1] based on the OH-PLIF results.
Droplet ''P 2 " has a larger original size compared to ''P 1 ", and therefore has longer relaxation time. However, its history is somehow similar to that of ''P 1 " in the sense that it only has very strong energy and mass change with the surrounding gas at the end of its lifetime. The difference with ''P 1 " is that it vanished at the flame base. Droplet ''P 3 " also kept its original size and temperature before entering the flame front, but it survived for about 8 ms in the hot  region behind the first flame front it met. And in this hot region, it remained at the boiling temperature (351 K), but its diameter quickly decreased, indicating a fast evaporation. It is also interesting to see that after the first flame front, the ''seen" gas temperature has decreased, and the mixture fraction has reached a very high value. This means that the fast evaporation of these large droplets has created locally a hot and rich region. This local ''fuel pool" supplies the combustion in both the outer and inner flame front. Indeed, studies in [20] showed that the outer edge of the inner flame region and the outer flame region are created by non-premixed combustion between cracked fuel and O 2 . OH-PLIF results in the experiment also revealed a fine and symmetric OH zone at the outer flame location, denoting a non-premixed flame front.
From the above discussions, we found that the double flame structure is strongly related to the behavior of spray polydispersity. The small droplets are mostly convected to the central region by the entrained coflow, and provide fuel for the combustion at the inner flame front. The droplets of intermediate size directly reach the flame base and vanish there. Large droplets can pass through the first flame front, and generate a local fuel pool behind the first flame front. This fuel pool is then responsible for the outer flame front and possibly also for the inner one. Fig. 9 gives the predicted OH field on a vertical cross-section. This is used as an indication of the flame structure. For comparison, the averaged OH-PLIF from experiment for cases ''a60 À I95" and ''a40 À I95" are shown in Fig. 10. The ''double-flame" structure is observed in both the predicted and experimental OH fields. The results of cases ''a60 À I95", ''a60 À I200" satisfy the description of the type ''C" flame, in which the flame is lifted far downstream the injector exit. However, compared to the experimental OH fields, the lift-off heights in both ''a60 À I95" and ''a40 À I95" have been over-predicted. The reason for the over-prediction is not clear yet, a first guess is that the recirculated hot gas in the experimental furnace may help to stabilize the flame at a lower axial location, but other causes such as chemical model may apply. To correctly take into account the influence of the hot gas recirculation, a full furnace simulation, as done in [6], is required. Since this kind of simulation demands enormous amount of computational resources, it has bot been carried out in this first study.

Influence of coflow conditions
Cases ''a60 À I95" and ''a60 À I200" show that with the same degree of dilution by CO 2 , the flame lift-off height decreases with the reduction of coflow velocity. This trend is in agreement with the experimental observation, but the reduction of the flame liftoff height is less significant in the simulation than in reality. This is probably due to the absence of the flame-atomization interaction in the simulation. In the experiment, heat is emitted from the hot flame zone to the liquid at the injector exit. When the flame gets closer to the injector, a larger amount of heat is received by the liquid. This may followed by a enhanced atomization if considerable temperature rise has been caused by the radiative heating, since the liquid surface tension decreases with rising temperature. As a consequence, smaller droplets are produced by the atomization process, and this in turn results in a even smaller lift-off height because small droplets decay quickly to gas phase velocity and evaporate much faster than large ones. Simulation results (not shown here) demonstrate that indeed the smaller the droplets injected at the atomizer location, the lower the flame lift-off height.
The influence of the flame-atomization interaction is even more significant in the cases ''a40 À I95" and ''a40 À I200", for which the Fig. 8. Droplet and ''seen" gas information as function of droplet age. Droplets were injected at Time ¼ 0ms. Top: droplet diameter (Dp) and temperature (Tp); middle: ''seen" gas mixture fraction Z and temperature Tg ; bottom: magnitude of slip velocity (U slip ), and angle between droplet and ''seen" gas velocity vectors (hU ).  flame type ''B" and ''A" were respectively observed in experiments, while flame type ''C" and ''B" are predicted by the simulation. As revealed by the experimental OH fields, the flames of these two cases anchor at the tip and inside the liquid sheet, respectively. A very strong flame-atomization interaction can therefore be expected. However, due to the absence of experimental information on the dense region, in the simulation the atomization process has not been directly modeled, rather individual droplets of the same initial distribution have been injected in all cases. Therefore, it is not surprising that the flame lift-off heights in cases ''a60 À I95" ''a40 À I95" and ''a40 À I200" have been overpredicted, and the flame type ''A" cannot be reproduced.

Conclusion
In the present study we investigated the oxy-fuel spray jet flames with a LES/FGM approach. The two-phase flow field was solved with an Eulerian-Lagrangian approach. Results on different computational domains and with different grid resolutions have been compared, and based on the results a proper numerical mesh was chosen. Comparison of the predicted OH mass fraction field with measured OH-PLIF showed that the simulation is able to capture the double flame structure observed in the experiment. A Lagrangian analysis of the droplet properties provided relevant information on the relation between dispersion, evaporation and combustion. Small droplets can very quickly adapt to local gas phase velocity, and are mostly convected to the central region by the entrained coflow. Before reaching the inner flame region, they only experience small changes in both the size and the temperature. And they quickly vanish after entering the flame front, and hence they continuously supply fuel to the inner flame region. Intermediate droplets follow a trajectory that ends at the flame base. Large droplets are able to pass through the first flame front they meet and survive even for a while in the hot postcombustion region. Once these large droplets have entered the hot region, they remain at boiling temperature, and experience a very fast evaporation, and this create a hot and rich fuel pool locally, which eventually supports the combustion at the inner and outer flame regions. For the type ''C" flame, which is stabilized far downstream the dense region, some major features were successfully captured, e.g. the gas phase velocity field and flame structure. The flame lift-off height of type ''B" flame was over-predicted. The type ''A" flame, where the flame stabilizes inside the liquid sheet, was not described well by the simulation, indicating the limitation of the current modeling approach. A possible reason is that strong atomization-combustion interaction exists in the ''A" and ''B" types of flame, modifying the droplet formation process. This suggests that atomization-combustion interaction should be taken into account in future study of these flame types.