The two-step model computes the net rate of soot generation, , in the same way as the one-step model, as a balance of soot formation and soot combustion:

(9–122)

In the two-step model, however, the rate of soot formation, , depends on the concentration of radical nuclei, :

(9–123)

The rate of soot combustion, , is computed in the same way as for the one-step model, using Equation 9–118 – Equation 9–120.

The default constants for the two-step model are applicable for the combustion of acetylene (C₂H₂). According to Ahmad et al. [10], these values should be modified for other fuels, as the sooting characteristics of acetylene are known to be different from those of saturated hydrocarbon fuels.

The net rate of nuclei generation in the two-step model is given by the balance of the nuclei formation rate and the nuclei combustion rate:

(9–124)

The rate of nuclei formation, , depends on a spontaneous formation and branching process, described by

(9–125)

(9–126)

Note that the branching term in Equation 9–125, , is included only when the kinetic rate, , is greater than the limiting formation rate ( number of particles/m³-s, by default).

The rate of nuclei combustion is assumed to be proportional to the rate of soot combustion:

(9–127)

where the soot combustion rate, , is given by Equation 9–118.

Since the Moss-Brookes model was mainly developed and validated for methane flames, a further extension for higher hydrocarbon fuels called the Moss-Brookes-Hall model is available in Ansys Fluent. Here, the extended version is a model reported by Wen et al. [701], based on model extensions proposed by Hall et al. [228] and an oxidation model proposed by Lee et al. [349]. The work of Hall [228] is based on a soot inception rate due to two-ringed and three-ringed aromatics, as opposed to the Moss-Brookes assumption of a soot inception due to acetylene or benzene (for higher hydrocarbons).

Hall et al. [228] proposed a soot inception rate based on the formation rates of two-ringed and three-ringed aromatics (C₁₀H₇ and C₁₄H₁₀), from acetylene (C₂H₂), benzene (C₆H₆), and the phenyl radical (C₆H₅) based on the following mechanisms:

(9–132)

Based on their laminar methane flame data, the inception rate of soot particles was given to be eight times the formation rate of species C₁₀H₇ and C₁₄H₁₀, as shown by

(9–133)

where = 127 x s^–1, = 178 x s^–1, = 4378 K, and = 6390 K as determined by Hall et al. [228]. In their model, the mass of an incipient soot particle was assumed to be 1200 kg/kmol (corresponding to 100 carbon atoms, as opposed to 12 carbon atoms used by Brookes and Moss [81]). The mass density of soot was assumed to be 2000 kg/m³, which is also slightly different from the value used by Brookes and Moss [81].

Both the coagulation term and the surface growth term were formulated similar to those used by Brookes and Moss [81], with a slight modification to the constant so that the value is 9000.6 kg m kmol^-1 s^-1 (based on the model developed by Lindstedt [379]).

For the soot oxidation term, oxidation due to O₂ (based on measurements and a model based on Lee et al. [349]) was added, in addition to the soot oxidation due to the hydroxyl radical. By assuming that the kinetics of surface reactions are the limiting mechanism and that the particles are small enough to neglect the diffusion effect on the soot oxidation, they derived the specific rate of soot oxidation by molecular oxygen. Therefore, the full soot oxidation term, including that due to hydroxyl radical, is of the form

(9–134)

Here, the collision efficiency is assumed to be 0.13 (compared to the value of 0.04 used by Brookes and Moss) and the oxidation rate scaling parameter is set to 0.015. The model constants used are as follows:

The kinetic mechanisms of soot formation and destruction for the Moss-Brookes model and the Hall extension are obtained from laboratory experiments, in a similar fashion to the NO_x model. In any practical combustion system, however, the flow is highly turbulent. The turbulent mixing process results in temporal fluctuations in temperature and species concentration that will influence the characteristics of the flame.

The relationships among soot formation rate, temperature, and species concentration are highly nonlinear. Hence, if time-averaged composition and temperature are employed in any model to predict the mean soot formation rate, significant errors will result. Temperature and composition fluctuations must be taken into account by considering the probability density functions that describe the time variation.

A description of the modeling of soot-radiation interaction is provided in The Effect of Soot on the Absorption Coefficient.

During soot formation, the soot particles are present in a wide range of diameters. The particle population balance method uses the particle-size distribution (PSD) for modeling of the soot particle surface area involved in the calculations of important sub steps of soot formation process like soot particles surface growth and oxidation. The various stages of soot formation have an impact on the soot PSD. The evolution of the soot particles’ population can be represented by the population balance equations, also known as the Smoluchowski master equations:

(9–135)

where is the particle number density function that denotes the number of the soot particles of the ^th size class per unit volume, is the collision coefficient between particles of size classes and .

For the first particle size class (the first equation), the right-hand side contains only one negative source term. This term reflects the reduction in the first size class population due to coagulation between the particles of the current size class and any other particles resulting in formation of particles of higher size classes. For all other higher particle size classes, the positive source term appears in the right-hand side of the second equation indicating the growth of the particle number density function as a result of said coagulation between particles of lower size classes.

The principal difficulty in obtaining a solution for Equation 9–135 involves solving an infinite number of particle size classes. The two approaches commonly used for these problems are the sectional method and the method of moments.

The sectional method is based on discretization of the entire range of particle size classes into a predefined finite number of intervals. Considering that the particle mass and population will be spread over all intervals and will be also affected by other processes (such as, surface reactions and oxidation), a sufficiently large number of sections (or finite intervals) is required for obtaining an accurate solution. For these reasons, applying the sectional method to Equation 9–135 is computationally demanding. An efficient alternative to the sectional method is the method of moments. This method considers moments of the soot PSD functions.

The discrete form of the ^th moment of particle-size distribution function for the particle of the size class is defined as:

(9–136)

From Equation 9–136, the total particle number of the particle population is equal to zeroth moment:

The mass of a particle is proportional to its class; therefore, the total mass of the particle population is calculated as:

where is a constant that represents the mass of a bulk species molecule comprising the particle core.

Similarly, the total particles volume can be obtained by:

where is the bulk density of particle core (soot density).

Multiplying Equation 9–135 by the soot particle size class and using the expression of Equation 9–136, we obtain:

(9–137)

where is the right-hand side of Equation 9–135, which is a coagulation source term.

In soot formation, in addition to the coagulation process, the nucleation and the surface reactions are the two other important contributing factors. In the presence of nucleation and surface reactions, Equation 9–137 is written as

(9–138)

where is the nucleation source term, and is the source term due to surface reactions including surface growth and oxidation.

Equation 9–138 describes the change in the moments due to various sub-processes at a single point. In a flow system, the moments will also be affected by convection and diffusion. For flow systems, the transport equation for the moments of soot concentration can be written as:

(9–139)

The source term for computation of is described in the sections that follow.

In Ansys Fluent, the nucleation process is modeled as coagulation between two soot precursor species. The soot precursor is a user-defined gas phase species. Typically, the soot precursors are poly-cyclic aromatic hydrocarbons (PAH), and the soot nuclei formation is modeled as coagulation of two PAH molecules. The mean diameter of the PAH is calculated from the number of carbon atoms in it and the soot density. The concentration of the PAH is obtained from the gas phase mechanism. In many practical cases, where the chemical mechanism used in simulations is small and does not include PAH species, the smaller species, such as C₂H₂, which is a building block of PAH molecules, can be used as a precursor species.

The nucleation source term for the first moment is calculated as:

(9–140)

where is the molar concentration of the precursor species, and is the constant calculated by:

(9–141)

The nucleation rates calculated based on the kinetic theory using Equation 9–140 are generally very large. To properly scale them, Equation 9–140 is adjusted using a sticking coefficient :

(9–142)

The value of the sticking coefficient varies with the size of the precursor species (molecular weight). Table 9.4: Sticking Coefficient for Different PAH Species lists the sticking coefficient values for different precursors as proposed by Blanquart and Pitsch [64].

Table 9.4: Sticking Coefficient for Different PAH Species

Species Name	Formula	Molecular Weight (gr/mol)
naphthalene	C10H8	128	0.0010
acenaphthylene	C12H8	152	0.0030
biphenyl	C12H10	154	0.0085
phenathrene	C14H10	178	0.0150
acephenanthrylene	C16H10	202	0.0250
pyrene	C16H10	202	0.0250
fluoranthene	C16H10	202	0.0250
cyclo[cd]pyrene	C18H10	226	0.0390

Other suggestions found in the literature propose using the sticking coefficient that is proportional to the fourth power of the precursor molecular weight:

(9–143)

where is the constant that can be either calculated using curve fitting from Table 9.4: Sticking Coefficient for Different PAH Species or obtained from experimental data. Equation 9–143 could be used to approximate the sticking coefficient values for precursors that are not listed in Table 9.4: Sticking Coefficient for Different PAH Species.

Alternatively, the nucleation can also be specified as an irreversible kinetic reaction between two precursors:

(9–144)

(9–145)

where and are the power exponent and stoichiometric coefficients of the precursor species.

Equation 9–145 is analogues to Equation 9–140 having the same PAH species and with the following values for the constants:

The nucleation source terms for higher moments are calculated from the lower moment source terms using the following expression:

(9–146)

Once formed, the soot particles collide with each other affecting the size distribution of the soot particles population. The process of coagulation assumes that the resulting particle remains a sphere with an increased diameter. The coagulation process changes the number density but not the total mass of the particle. The source terms in the moment transport equation due to coagulation are calculated as:

(9–147)

The coagulation term for the second size class is equal to 0 because the coagulation does not change the total mass of soot formed.

The coagulation terms for higher size classes are calculated by:

(9–148)

where is the collision efficiency, which is dependent on the coagulation regimes.

The coagulation process can take place in different regimes:

The regimes of coagulation depends on the Knudsen number

(9–149)

where the gas mixture mean free path, and is the soot particle diameter. The mean free path of soot particles are calculated as follows:

(9–150)

where is the Boltzmann constant, is the pressure, and is the diameter of the gas molecules calculated by

(9–151)

where is the molecular weight of the gas.

The mean diameter of the soot particles is calculated using the definition of the soot moments in Equation 9–136 in the following manner:

(9–152)

where is the diameter of smallest particle, which is a single carbon atom.

For different regimes, different treatments are applied when solving Equation 9–148 as further described.

Continuum Coagulation ( < 0.1)

The collision coefficient in the continuum region is expressed as:

(9–153)

where is the Cunningham slip correction factor equal to 1+1.257 Kn, is the continuum collision factor calculated as:

where is the molecular viscosity of the mixture.

Substituting the collision coefficient (defined in Equation 9–153) into Equation 9–148 gives the source terms due to coagulation described below.

For the first moment, the moment source term due to continuum coagulation is:

(9–154)

where .

By introducing the concept of reduced moment , Equation 9–154 can be rewritten as:

(9–155)

Similarly, for = 2, 3 …. the moment source terms due to continuum coagulation are expressed as:

(9–156)

The coagulation terms in Equation 9–155 and Equation 9–156 consist of the fractional order moments as well as the negative order moments that need to be calculated for closure. Here, the process of interpolative closure proposed by Frenklach ([192]) has been used. In the interpolative closure, the positive fractional order moments are obtained using a polynomial interpolation. Since the numerical values of the subsequent moments differ by an order of magnitude, the logarithmic interpolation is used to minimize interpolation errors. The fractional moments are obtained using the following expression:

(9–157)

where, is the Lagrange interpolating polynomial.

In order to calculate the negative order fractional moments, an extrapolation from integer order moments is used

(9–158)

where is the Lagrange extrapolation. It has been observed that a quadratic interpolation of Equation 9–157 and a linear extrapolation of equation Equation 9–158 work reasonably well.

Coagulation in the Free Molecular Regime ( > 10)

In the free molecular regime, the collision coefficient is calculated as:

(9–159)

where is the coefficient for the free molecular collision efficiency calculation in the following form:

(9–160)

Using Equation 9–159 for the collision coefficient, Equation 9–155 and Equation 9–156 give the moment source terms due to coagulation in the free molecular regime.

Note that the collision coefficient written in the form of Equation 9–159 is non-additive. Therefore, Equation 9–155 and Equation 9–156 are now expressed as:

(9–161)

where is the grid function defined as:

(9–162)

Since the evaluation of this expression is not straightforward for =1/2, the grid functions are evaluated for integer values of . Then the grid functions for the fractional values are obtained using Lagrangian interpolation.

For example, for cases where three moments are solved, the following grid functions need to be calculated:

These functions are obtained by first calculating , , and and then applying Lagrangian interpolation:

Note that the grid function is symmetrical, that is .

Coagulation in Transition Regime (0.1 < < 10)

When the mean free path is of the order of the soot particle diameter, the coagulation term is calculated using the harmonic average from free molecular and continuum regimes in the following manner:

(9–163)

Once the soot is nucleated, carbon from the gas phase is deposited on the surface of soot particles, leading to a growth of soot particles’ size as well as an increase in the soot mass. The soot formed due to nucleation is typically very small. Therefore, carbon deposition at the soot surface is the dominant mode of soot formation. In addition, the soot particles lose mass due to oxidation of the soot surface by O₂, OH, or any other oxidation species. The process of surface growth and oxidation are kinetically controlled and quite complex. Resolving the whole kinetics of these sub-processes involves large mechanisms and therefore computationally expensive. However, the process of soot surface growth and oxidation can be approximated by using a small mechanism (like a single step global mechanism) or reduced mechanisms.

One of the most widely used mechanism to model the surface growth and oxidation is the Hydrogen Abstraction C₂H₂ Addition (HACA) mechanism. The HACA mechanism represents the process of surface growth and oxidation by using reduced chemistry involving only few reaction steps.

There are a few variants of HACA mechanism. The HACA mechanism used in Ansys Fluent for modeling surface growth and oxidation are presented in Table 9.5: Arrhenius rate parameters for HACA mechanism.

In Table 9.5: Arrhenius rate parameters for HACA mechanism, the first six reactions are related to the surface growth and the last two reactions are related to oxidation due to O₂ and OH, respectively. The reaction rates of these reactions have been taken from Appel et al. [22].

Table 9.5: Arrhenius rate parameters for HACA mechanism

Reaction
	(cm3mol-1s-1)		(kcal/mol)
	4.2 x 1013	0	13.0
	3.9 x 1012	0	11.0
	1.0 x 1010	0.734	1.43
	3.68 x 108	1.139	17.1
	2.0 x 1013	0
	8.0 x 107	1.56	3.8
	2.2 x 1012	0	7.5
	Neoh at al. (1981) model, =0.13

Since the oxidation by OH reaction (reaction 6 in Table 9.5: Arrhenius rate parameters for HACA mechanism) includes the collision between one gas-phase species and a surface species, the reaction rate constant is specified in terms of a sticking coefficient (reaction probability) . The unitless “Arrhenius-like” form of the sticking coefficient is defined as a probability where the computed values greater than one are set to the unity:

(9–164)

The sticking coefficient is converted to the mass-action kinetics rate constants using the collision frequency of the gas species with the bulk species:

(9–165)

In Equation 9–165, the term accounts for the conversion of unitless sticking coefficient to appropriate units of the rate constant. The gas/surface collision frequency is considered by the term . If the occupancy is unity for all species, then the production term is equal to one.

The surface reactions are fast. The concentration of active radical sites are in steady state. The concentration of the active radical sites can be represented in terms of a fraction of the active sites:

(9–166)

The concentration of active sites on the soot surface is . Total active sites at soot surface are approximated as:

(9–167)

The active radical sites are calculated as:

(9–168)

where is the carbon atom diameter.

The moment source terms due to the surface growth reaction and oxidation are calculated as follows:

In Equation 9–169 through Equation 9–171, is the number of carbon atoms removed or added:

For the -th moment, the total source term due to surface growth and oxidation is

(9–172)

Note that with the HACA mechanism, acetylene is consumed during the surface growth of soot. Since Ansys Fluent generates the flamelets before the solution, the flamelets have no knowledge of the soot process, which can cause an over-prediction of the soot surface growth and eventually higher soot yields with non-premixed/partially premixed combustion models.

In Coagulation, soot particle coagulation assumes coalescent soot growth, whereby the larger particle resulting from collision of two smaller particles is always spherical. However, in many practical combustion devices, the growth of soot may not be limited to the coalescent regime. The colliding particles can stick to each other and form chain-like structures called aggregates. The aggregate has a larger surface area than a coalesced single sphere of the same mass. The formation of the aggregates, therefore, enhances the rate of coagulation process. Ignoring the aggregate formation could lead to errors in calculations of the coagulation and the soot surface growth rates.

In the regime of aggregate collision, the soot particle size distribution becomes bivariate. The number of primary particles within the aggregate represents the internal coordinate, while the aggregate mass and size represents the external coordinate. Resolving a two-dimensional PDF is not trivial and computationally expensive. Ansys Fluent has adopted the approach proposed by Kazakov and Frenklach [298] in which the method of moments for modeling the soot formation is extended to account for soot aggregation. The current formulation separates the two coagulation regimes, namely, the coalescence and aggregation, switching between the two based on a user-specified critical diameter. The coagulation is assumed to be coalescent if the average particle diameter is less than the critical diameter. When the soot particle diameter reaches the critical value, the coagulation is automatically switched to the aggregation regime.

In the presence of the aggregation process, Equation 9–136 through Equation 9–139 represent the size distribution and the transport of the soot aggregate mass (that is, the concentration moments). In addition to the aggregate moments, similar to Equation 9–136, the discrete form of concentration moment for the distribution of the primary particles in the aggregates is defined as:

(9–173)

where is the particle ^th moment, and is the number of the primary particles in ^th class. Similar to the aggregate mass (or global) moments, here also the concept of the reduced particle moment can be used:

(9–174)

where is the ^th reduced moment of the primary particles within the aggregate.

By definitions in Equation 9–136 and Equation 9–173, the particle zeroth moment and the global zeroth moment are identical, that is:

(9–175)

Also, the average number of the particles in an aggregate, , is given by the first two moments as follows:

(9–176)

Similar to the aggregate moments, particle moments will also evolve due to the various soot formation sub-steps, such as, nucleation and coagulation, as well as particle transport in the physical space. The governing equations for the change in the particle moments can be written as:

(9–177)

Because Equation 9–175 already yields the solution for = 0, Equation 9–177 ) is only solved for > 0.

The first particle moment gives the average number of the primary particles (Equation 9–176) from which the particle mass can be calculated. = 0 because the particles total mass does not change due to coagulation process.

For higher moments, in equation Equation 9–177 needs to be computed. Using Equation 9–148 along with the definition of particle moments in Equation 9–173, can be written as

(9–178)

where is the aggregate collision coefficient. Its value is dependent on the collision regime, which is calculated using the Knudsen number (Equation 9–149).

Aggregate Coagulation in Continuum Regime ( < 0.1)

The aggregate collision coefficient in the continuum regime is expressed as:

(9–179)

where and are the collision diameters of the aggregates of different size classes and , respectively. The rest of notation is the same as in equation Equation 9–153. For coalescent collision, the collision diameter and the particle diameter are assumed to be the same. However, due to aggregation, the collision diameter can differ significantly from the particle diameter depending on the fractal structure of the aggregate. An aggregate is composed of a number of primary particles, which are assumed to be spherical and of equal size. The collision diameter of the aggregate is calculated from the following relation:

(9–180)

where is the diameter of the primary particle, and is the fractal dimension that describes the fractal structure of the aggregate. A value of between 1.7 and 2.0 is found to be reasonably good and suggested in the literature for aggregate coagulation. In the pure coalescent coagulation, = 3. The mass of the aggregate can be specified in the following manner:

(9–181)

From Equation 9–180 and Equation 9–181, the aggregate collision diameter can be written as:

(9–182)

Substituting the value of the collision diameter from Equation 9–182 into Equation 9–179 and replacing the value of the collision coefficient in Equation 9–178, we obtain the aggregate coagulation source term for the ^th moment:

(9–183)

where,

(9–184)

Equation 9–183 involves binary moments of the soot aggregate size and the number of primary particles. The calculations for binary moments are not trivial and require a multi-dimensional PDF of the particle size distribution. Kazakov and Frenklach [298] suggested approximating the binary moment using two one-dimensional moments as follows:

(9–185)

Thus, Equation 9–183 can be written as:

(9–186)

Equation 9–186 contains the fractional and intermediate primary particle moments, which are calculated using the same interpolation method employed to aggregate moments (see Equation 9–157 and Equation 9–158).

Aggregate Coagulation in Free Molecular Regime ( > 10)

The collision coefficient in the free molecular regime can be written as:

(9–187)

where and are the mass and collision diameter of the aggregate, respectively. The collision diameter is calculated from Equation 9–182. Similar to the coalescent coagulation in the free molecular regime calculations, difficulties arise when calculating the source due to the coagulation of the aggregates (Equation 9–178) because of the non-additive collision coefficient term from Equation 9–187. In a manner similar to Equation 9–161, can be expressed as:

(9–188)

where is the aggregate grid function. Equating the right-hand sides of Equation 9–178 and Equation 9–188 and substituting the collision coefficient from equation Equation 9–187, the aggregate grid function can be written as:

(9–189)

The collision diameter in Equation 9–189 is given by equation Equation 9–182. Therefore, Equation 9–189 can be rewritten as

(9–190)

Using Equation 9–136, Equation 9–173, and Equation 9–185, the aggregate grid function can be expressed as follows:

(9–191)

Aggregate Coagulation in Transition Regime (0.1 < < 10)

Similar to Equation 9–163, the coagulation term in the transition regime is calculated using the harmonic average from the free molecular and the continuum regimes.

The coagulation and surface growth rate source terms for global moments are calculated using Equation 9–155, Equation 9–156, Equation 9–161, and Equation 9–172. These source terms are dependent on the diameter and the surface area of the coagulated particle. When the aggregation model is enabled, the resulting diameter and surface area of the aggregate will be larger than those of the spherical particle. The coagulation and surface growth rate terms for aggregate moments are modified as further described.

Modified coagulation term

The modified continuum regime equation Equation 9–155 that includes the change in collision diameter, is follows:

(9–192)

Similarly, Equation 9–156, which is the coagulation source for ^th moment, is calculated as:

(9–193)

In the free molecular regime, the grid function in Equation 9–162 is modified to account for aggregation as follows:

(9–194)

Equation 9–194 can be written in the alternative form of reduced moments as follows:

(9–195)

Equation 9–192 through Equation 9–195 are obtained by substituting the modified value of the collision diameter from Equation 9–182 into Equation 9–155, Equation 9–156, and Equation 9–161.

Modified surface growth and the oxidation rate terms

Similarly, the surface growth and the oxidation rate for different moments are enhanced by a factor and calculated as follows:

(9–196)