Modelling of soot formation in laminar diffusion flames using a comprehensive CFD-PBE model with detailed gas-phase chemistry

A discretised population balance equation (PBE) is coupled with an in-house computational fluid dynamics (CFD) code in order to model soot formation in laminar diffusion flames. The unsteady Navier–Stokes, species and enthalpy transport equations and the spatially-distributed discretised PBE for the soot particles are solved in a coupled manner, together with comprehensive gas-phase chemistry and an optically thin radiation model, thus yielding the complete particle size distribution of the soot particles. Nucleation, surface growth and oxidation are incorporated into the PBE using an acetylene-based soot model. The potential of the proposed methodology is investigated by comparing with experimental results from the Santoro jet burner [Santoro, Semerjian and Dobbins, Soot particle measurements in diffusion flames, Combustion and Flame, Vol. 51 (1983), pp. 203–218; Santoro, Yeh, Horvath and Semerjian, The transport and growth of soot particles in laminar diffusion flames, Combustion Science and Technology, Vol. 53 (1987), pp. 89–115] for three laminar axisymmetric non-premixed ethylene flames: a non-smoking, an incipient smoking and a smoking flame. Overall, good agreement is observed between the numerical and the experimental results.


Introduction
Soot is particulate matter that is formed due to the incomplete combustion of hydrocarbon fuels. It is responsible for affecting the thermal efficiency of combustion devices, as well as Earth's climate and atmosphere due to its radiative properties; in addition, it has harmful effects on human health and its presence in emissions must be avoided [1]. Recently, it has been found that the contribution of soot to the greenhouse effect is far greater than previously thought (the second biggest contributor to the greenhouse effect behind CO 2 emissions) [2]. Furthermore, soot has a negative impact on human physiology as it is made up of polycyclic aromatic hydrocarbons (PAHs) that have been found to be carcinogenic and mutagenic [3]. At present soot formation is not well understood, largely due to a number of experimental limitations (such as measuring the size of incipient particles in the early stages of their formation). Computational simulations of soot formation are of vital importance for the design and optimisation of combustion equipment. However, modelling soot formation poses great challenges; apart from the limits in understanding and availability of experimental data, the computational problem of soot prediction is formidable, especially in the case of turbulent combustion where complex interactions between mixing, chemical reactions and soot formation/destruction processes require additional modelling.
It is particularly important, therefore, that modelling approaches are thoroughly evaluated in laminar flames before being ported to turbulent combustion, in order to gauge the effect of the various elements of the model and their importance before introducing further modelling errors pertaining to mixing-chemistry-particle formation interactions.
A considerable amount of research has been reported on soot formation in laminar diffusion flames. To begin with, an experimental and numerical investigation of a laminar co-flow methane-air diffusion flame was reported by Smooke et al. [4]. In that work, a quantitative description of the structure of the flame was obtained, without any soot computations being performed. Further experimental and numerical investigations of an ethylene flame, including soot prediction, were accomplished in [5] by employing a sectional soot model. Soot volume fraction profiles were obtained by varying the jet fuel dilution ratio.
In studies by Santoro et al. [6,7], several experiments were performed for the same configuration under different co-flow velocities. The following three sooting flame types appear by varying the inlet flow rates: a non-smoking, an incipient smoking and a smoking flame. In the non-smoking flame, the oxidation process completely destroys the soot particles within the flame area and no soot is emitted, whereas in the other two flame types soot escapes from the flame region. More studies on Santoro's flames can be found in [8][9][10], which use a simplified soot model involving two transport equations. These simplified soot models could predict only the mean properties of the particle size distribution (PSD), ignoring the polydispersity of soot particles. Moreover, these papers concentrated on the non-smoking and smoking flame experiments. Kennedy and Yam [9] simulated these two flame conditions and obtained relatively good agreement with the experimental data for the non-smoking flame, but not so good for the smoking one. They employed an OH oxidation collision efficiency and stated that a more accurate O 2 oxidation model is needed. The study of Liu et al. [8] obtained better results partly due to applying more sophisticated correction factors to the oxidation rates of both O 2 and OH. A later study of both flames can be found in [10], where the soot kinetics and correction factors employed are the same as in [8], but several radiation models are applied. In [11] all three flame conditions were studied with a polydisperse soot model, and the disagreement in the literature regarding the soot kinetics of O 2 and OH was investigated. Finally, we should mention the numerical investigation of [12] that focused on investigating the effect of gravity on the flame structure; the study is on a different flame and no comparison with experiments was attempted.
Although intensive research has been carried out over the last three decades, our knowledge of soot formation is still far from complete. No soot model exists that could be applicable to a wide range of fuel compositions and flame conditions [13]. Most of the experimental and numerical efforts focus on the average properties of the PSD, such as soot number density and soot volume fraction. Only recently has interest started to shift towards the nanoscale morphology and size distribution of particles [14]. To improve our insight into soot phenomena, more elaborate predictive soot models are needed to estimate the full PSD by solving a detailed population balance equation (PBE). The aim of this study is to apply a finite volume-discretised PBE model coupled with a detailed in-house computational fluid dynamics (CFD) code (BOFFIN) and a comprehensive gas-phase chemical mechanism, and to test the methodology by simulating all three Santoro flames (non-smoking, incipient smoking and smoking). The use of a discretised PBE allows for prediction of the complete PSD without any prior assumption of the shape distribution. Several discretisation methods have been proposed so far for solving the PBE (e.g. [15][16][17]); in this paper the PBE is solved with a finite volume technique, which is an extension of the original collocation finite element method to be found in Rigopoulos and Jones [16]. The paper thus explores the importance of combining comprehensive fluid dynamics, reactions and PBE modelling for the modelling of soot formation in laminar flames, a significant step before moving on to the modelling of turbulent flames where the additional uncertainties of modelling turbulence-chemistry particle formation interactions are introduced.

Experimental and simulation setup
The model presented here will be used to simulate three laminar axisymmetric diffusion flame configurations the experimental setup and geometry of which are described in [6,7]. As the flames are axisymmetric, the numerical domain can be represented in two dimensions. Figure 1 shows the domain and the boundary conditions. The inlet fuel is pure ethylene and the oxidizer is air. The central fuel tube has an internal diameter of 11.1 mm and the outer passage is 101.6 mm. Ambient temperature and atmospheric pressure conditions are applied on this system. The fuel and the oxidizer have uniform exit velocity profiles. In the non-smoking flame the fuel velocity is 3.98 cm/s and the oxidiser's velocity is 8.9 cm/s. In the incipient smoking flame the fuel velocity is increased to 4.75 cm/s, whereas the oxidiser velocity is kept the same as that of the non-smoking flame. In the smoking flame the fuel and oxidiser velocities are increased to 5.06 and 13.3 cm/s, respectively. Symmetry conditions are applied to the left and to the right planes, velocities and concentrations are specified at the bottom and zero gradient conditions are imposed at the top. The simulations are performed with a fine grid, consisting of 200 points in the axial and 100 points in the radial directions. A grid independence study was previously performed to ensure that this grid is sufficient. The grid is non-uniform and is finer towards the jet burner exit, as larger gradients exist in that area.

Gas-phase combustion
In this study we have sought to employ detailed gas-phase chemistry. The comprehensive and validated mechanism found in [18] is used for ethylene combustion, comprising 75 species and 529 reactions. NO x reactions are not included in this mechanism and N 2 is present as an inert species. While other authors have used gas-phase mechanisms that include N 2 chemistry, such as GRI3.0 [8], they removed the reactions and species related to NO x formation and focused on soot formation.

Soot kinetics
Nucleation is a poorly understood mechanism in soot formation, as well as being the starting point of the whole process, and it generates the incipient and smallest particles. There are several proposals for describing the nucleation mechanism. They state that the precursors of soot formation could be polyacetylenes, ionic species or PAHs. In Leung et al. [19] the nucleation step is assumed to be a first-order reaction, dependent on the molar concentration of C 2 H 2 . In Lindstedt [20] this nucleation mechanism is modified and extended to include the concentration of C 6 H 6 . An even more detailed PAH nucleation step is used in Smooke et al. [21]. This inception model is employed as a function of several gas-phase species. The most widely approved PAH-based nucleation model involves the collision of two C 16 H 10 molecules according to gas kinetic theory [22]. In this study, the inception process is modelled as a first-order reaction based on C 2 H 2 concentrations whose reaction rate constant (k n ) can be found in Liu et al. [8]. According to the latter, the incipient soot particle is assumed to be approximately 2.4 nm in diameter with 700 carbon atoms (C min ). The nucleation rate (B 0 ) is shown in Equation (1), where N A is Avogadro's constant: Surface growth is a heterogeneous surface reaction and it takes place when C 2 H 2 molecules react with the soot particles, resulting in an increase of soot particle volume. To represent this process, some authors use simple first-order rate equations [19,20] and some others use the more detailed HACA mechanism [22,23]. Several growth models were tested in the context of our approach, but the best agreement was observed with a first-order reaction according to Harris et al. [24]; an expression due to [21] is used in this study. The growth rate expression is shown below: Oxidation is the process where gas-phase species react with the surface area of soot particles, resulting in a reduction of the soot particle size. The dominant oxidating species are O 2 and OH. The oxidation models that are used in this paper are the O 2 oxidation model found in [25] and the OH oxidation model described in [26]. Both oxidation models are also implemented in [27] and are shown below: In the above expressions, ρ s is the soot density. The latter is taken to be 1900 kg/m 3 in accordance with the nucleation step expression obtained from the Liu et al. study [8]. The values of the reaction rate constants k α , k B , k T , k z can be found in Nagle and Strickland-Constable [25]. The correction factors f O 2 and f OH are functions of temperature and are used on all three flame conditions as shown in Liu et al. [8].
Regarding coagulation, the experimental data of Megaridis and Dobbins [28] on the same flame show that the number density of soot primary particles remains almost constant throughout the growth region across the flame, thus indicating that neither coagulation nor aggregation mechanisms significantly affect their number. Therefore the coagulation rate is not sufficient to cause a reduction of the primary particle number density in the growth areas. Furthermore, the study of Liu et al. [8] indicates that coagulation might be significant only in the early stage of the inception process, but seems to be insignificant in the growth region. Based on these observations, no coagulation mechanism is employed in this study.

Coupled CFD-PBE model
In the in-house code BOFFIN, the transported quantities are cast into a general transport equation form and solved in cylindrical coordinates. An algorithm based on the pressure correction concept is used to handle the pressure-velocity coupling. The transport equations are discretised with a finite volume scheme, using an upwind formulation for the convection terms. The gravitational term is applied only to the axial momentum equation, in order to take into account buoyancy effects. To enhance the conservation properties of the finite volume technique, the total variation diminishing (TVD) scheme is applied using the van Leer method [29]. The algebraic equations resulting from the discretisation are solved using a conjugate gradient solver. An optically thin approximation radiation model is employed in the enthalpy equation to account for the radiation of the gas species (CO, CO 2 and H 2 O) and soot particles, as shown in [13]. The species' differential diffusion is included by using detailed thermal and transport coefficients of each species.
The PBE in the form employed here is shown below (Equation 7), and is a partial differential equation whose independent variables are space, time and an additional variable that is a measure of particle size, ξ -here this is chosen to be the particle diameter.
Here φ α are the species' concentrations, ρ is the mixture density and p is the particle diffusivity, set to a very small value as it has a negligible contribution. Thermophoresis is also applied in this study via the additional velocity u t i . The thermophoretic velocity drives the soot particles radially inward to the flame. It is independent of soot particle size and is only dependent on position and temperature gradients. Therefore it is expected to be stronger at the wings of the flame and act as diffusion on soot particles by lowering the maximum soot volume fraction. The thermorphoretic velocity is kept the same for all soot size classes within the same computational cell and is computed via Equation (8) as stated in [8,9]: The nucleation, growth and oxidation are collected at the right-hand of Equation (7), in anticipation of their treatment as source terms. The nucleation rate B 0 corresponds to the rate given by Equation (1), appropriately scaled in order to express the number of particles produced per unit volume of mixture and per unit diameter. The function G represents the combined effect of growth and oxidation, as both processes result in a continuous change of particle size. This equation is coupled with the species transport equations through the nucleation and growth terms, which are functions of the species' concentrations.
In this work, the PBE is discretised by a finite volume method and a TVD scheme. The method is an extension of the collocation finite element discretisation approach PBE that can be found in Rigopoulos and Jones [16]. The particle size coordinate is discretised into 400 elements using a uniform grid. The method results in eliminating the diameter from the list of independent coordinates, thus yielding a set of PDEs in space and time. These are brought into the form of the generic transport equation and coupled with the CFD code in a manner similar to the other scalars. The source terms of the gas species and population density are treated via a fractional step approach.

Results and discussion
The simulation results for the non-smoking, incipient smoking and smoking flame are presented in this section. The amount of experimental data available for the non-smoking experiment is substantially more than that for the two smoking flames. As such, the non-smoking flame was investigated first in order to validate the PBE and soot model. Subsequently the same sub-models were used for the incipient smoking and smoking flames.
We first present results for the temperature, axial velocities, OH, C 2 H 2 and soot volume fraction radial profiles at different heights from the burner exit, for the case of the nonsmoking flame. Figure 2 shows the radial profiles of axial velocity and temperature against the experimental results. The axial velocity profiles are shown in Figure 2   70 mm above the burner. Excellent agreement is achieved between the computed profiles and experimental results. Furthermore, in Figure 2(b) the radial temperature profiles are shown at 20, 50 and 70 mm above the burner. The profiles at 20 and 50 mm also show excellent agreement, but at 70 mm they exhibit a discrepancy. While the centreline temperature is still well predicted there, it is followed by an overprediction and later drops to the same levels as the measurements. This problem has been identified by other authors as well [8,10]. The most likely cause of this overprediction is the high oxidation in that region that reduces the local soot volume fraction, therefore emitting less radiation.
In Figure 3, the radial mole fraction of two species' profiles are shown at 7, 20 and 70 mm above the burner. Figure 3(a) shows that the mole fraction of C 2 H 2 at 20 mm is in very good agreement with the measurements. At 7 mm, however, there is a discrepancy between the simulation and experimental results. This occurs also in the study of Kennedy and Yam [9], although no attempt was made to explain it. A likely cause of this overprediction very close to the burner is fuel preheating; this is a point that will be elaborated further later when discussing the results for the soot volume fraction. Furthermore, in Figure 3  OH measurements could be as high as 50% [30]. These two species are of key importance for soot prediction, hence obtaining good agreement for them is vital.
Radial soot volume fraction profiles at 15 and 50 mm above the burner are shown in Figure 4(a). The soot volume fraction profiles are in a good agreement with the experimental results at both distances, apart from an overprediction of the peak at 15 mm. Figure 4(b) shows the axial profile of the soot volume fraction integrated along the radius. Agreement is good apart from an overprediction at 0.02-0.04 m; furthermore the peak value is located around 0.035 m, whereas in the experiments it was found to be at 0.04 m. This overprediction is likely to be related to the issue of fuel preheating mentioned before and to the radiation modelling. Regarding the former, the preheating was not actually measured in the experiments, and there is much uncertainty regarding what temperature values one should assign to the fuel inlet, if it were to be simulated. Although some studies have attempted to include it [31], the values were chosen on the basis of the results of their computations and are unlikely to match a different model such as the one used in this study. Regarding radiation, the alternative would be to use a discrete ordinates method; such a model was used in [8] and the results were compared with an optically thin model. It was observed that for the non-smoking flame the results of the centreline temperature and integrated soot volume fraction have not changed significantly, although both results are slightly improved. An investigation of the impact of the radiation model is beyond the scope of this paper.
Another set of data provided in the experiments is the pathline of the maximum soot volume fraction; this is computed by finding, at each axial location, the radial position that contains the maximum soot volume fraction. Comparisons with the experimental measurements at those points are shown in Figure 5; again, the overall agreement is reasonable apart from an underprediction of the number density near the inlet. It is worth mentioning that the Liu et al. study [8] had similar underprediction near the inlet but better agreement in the latter two points of the number density of Figure 5(a). Liu et al. [8] did not include the destruction of the number density of soot particles due to the soot oxidation mechanism, whereas in our study this has been taken into account. Thus the underpredicted number density suggests that possibly an additional nucleation mechanism is needed.
Contour plots of the nucleation rates and surface growth rates are shown in Figure 6. The nucleation rate is dominant very close to the burner sides (early stages of soot formation), with significant contribution at the centreline region in the later stages. On the other hand,    surface growth is dominant in the wings of the flame. Therefore soot predictions at the centreline are expected to be particularly sensitive to the nucleation kinetics.
In Figure 7 the complete PSD is shown for three different points in the axial direction. The PSDs are normalised with their respective moment values at each computational cell. In Figure 7(a) we can see that the number density is moving towards larger particle sizes    as the distance from the burner height is increased up to a certain point (tip of the flame), while in Figure 7(b) the normalised volumetric PSD shows a unimodal distribution. The particles are increasing their diameter size up to a maximum of approximately 100 nm. Unfortunately, there are no experimental PSD results for comparison.
We now consider the incipient smoking flame, where the integrated soot volume fraction profile is shown (Figures 8(a) and 8(b)). To observe the influence of the correction factors suggested in [8], we have conducted a simulation without using them; the result is that no soot is emitted from the tip of the flame (Figure 8(a)). By applying these factors, however, the oxidation strength is decreased and a small amount of soot volume fraction is emitted, as shown in Figure 8(b).
Moving on to the smoking flame, results are shown for the radial profiles of soot volume fraction at 40 and 70 mm above the burner (Figure 9(a)), while results for the integrated soot volume fraction are shown in Figure 9(b). As discussed before, the centreline predictions are particularly sensitive to the nucleation kinetics and could thus be improved by an improved nucleation model, as well as by the inclusion of preheating and detailed radiation modelling. The results shown in Figure 9(a) are very close to those in the study of Kennedy and Yam [9] and exhibit the same discrepancies. Furthermore, the pathline of the maximum soot volume fraction and number density in the smoking flame are shown in Figure 10. Overall, in these    comparisons the soot volume fraction is in reasonable agreement with the experimental results, although the centreline and peak values exhibit a significant underprediction, and generally the level of agreement is slightly worse than for the non-smoking flame. Finally, the PSD results for the smoking flame are shown in Figure 11. Again, the results of the volumetric PSD indicate that the mass distribution of the primary particles is unimodal, having almost an order of magnitude particle diameter span in different axial positions for both the non-smoking and smoking flame.

Conclusions
The ultimate challenge of soot modelling is the prediction of soot formation in turbulent flames, as the majority of industrial processes that produce soot operate under turbulent flow conditions. However, the prediction of soot in turbulent flames is a formidable challenge, as in addition to the multitude of uncertainties regarding its kinetics one has to tackle the issue of chemistry-particle-formation-turbulence interaction. In sooting flames fluid dynamics, chemical kinetics and soot formation processes acting on particles of various sizes are coupled. In this study we presented a methodology for modelling soot formation that involves a detailed PBE, solved with a discretisation approach, coupled with a CFD code and a detailed mechanism comprising 75 species and 529 reactions. Our aim was to obtain a comprehensive framework that can be used as a basis for extending to turbulent flames, after being tested and validated with laminar flames. Particular emphasis was put on soot polydispersity, which was modelled here by a detailed PBE; the complete PSD can thus be obtained without any prior assumption on the shape of the distribution. The method was subsequently applied to simulate soot formation in three laminar diffusion flames -a non-smoking, an incipient smoking and a smoking flame. Overall, very good agreement was accomplished along all of the three flames, apart for some discrepancies that can be primarily ascribed to uncertainties regarding the possible presence of fuel preheating in the experiments or certain choices in soot kinetics and radiation models. While these uncertainties cannot be fully addressed, the synergy of comprehensive fluid dynamics, detailed kinetics and PBE provides a sound basis for future extension of the methodology to turbulent sooting flames.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
Petros Akridis' work was sponsored by an Engineering and Physical Sciences Research Council (EP-SRC) Doctoral Training Award (DTA), while Stelios Rigopoulos gratefully acknowledges financial support from the Royal Society.