Diffusion flames are modeled on a simple one-dimensional domain of fixed size $L$. Only diffusive transport is considered, with no advective transport. The system pressure is assumed spatially uniform and constant in time. By convention, the lower boundary is taken to be oxidizer and the upper boundary is taken to be fuel. This results in the mixture fraction increasing monotonically from zero to one from the lower to the upper boundary. The mixture fraction is the local mass fraction of gas material originating in the fuel stream.

As noted, there is no advective flame strain, but the fixed domain size imposes a diffusive strain, and at steady state there is a balance between reaction and diffusion. As the domain size is decreased, the diffusive mixing rate increases (through the larger imposed gradients), and the flame eventually extinguishes when the reaction rates cannot keep up with the mixing rates.

This flame configuration was advocated by Pierce for application as a subgrid model for turbulent combustion [3], as it is neutral regarding opposed or anti-opposed flow, in contrast to, e.g., laminar flamelets that assume opposed flow.

Species equations

The equations solved for gas species mass fractions $y_k$ are

$$\prtl{y_k}{t} = -\frac{1}{\rho}\prtl{j_k}{x} + \frac{\doot{m}_k^\tprime}{\rho},$$

Where $\rho$ is density, $j_k$ is the mass diffusion flux, and $\doot{m}_k^\tprime$ is the species reaction rate per unit volume.

The unsteady term on the left-hand side of the species equation is retained to allow unsteady evolution to a desired steady state, that is, for computational convenience. As written, the equation neglects the time dependency of density. At steady state, the obtained solution is correct for the given configuration, so the time-independent assumption on density only affects the intermediate unsteady profiles.

Dirichlet boundary conditions are specified for species, that is, the species mass fractions are specified.

Ideal gases are assumed, and the density is computed as $\rho=MP/RT$, where $M$ is the mean molecular weight of the gas, $T$ is temperature, $P$ is pressure, and $R$ is the gas constant. The diffusion flux $j_k$ is given by $$j_k = -\rho D_k\prtl{y_k}{x} - \left(\frac{\rho D_ky_k}{M}\prtl{M}{x}\right),$$ where $D_k$ is the effective species diffusivity. Mixture average transport properties are assumed. An option for unity Lewis numbers is available, in which case the species diffusivities are all equal to the thermal diffusivity, and the term in parentheses for $j_k$ is ignored.

Energy equation

Energy can be solved in several forms, usually in terms of either enthalpy or temperature. The equation is simpler when solving for enthalpy, but this requires regular inversion to compute the temperature for evaluation of themochemical properties. This inversion requires solution of a nonlinear equation, converged to some tolerance. This is trivial, but when combined with coupled spatial solvers on stiff systems, such as done here, using, e.g., CVODE, the system is less stable than solving temperature directly. Here, the temperature equation is solved, and is given by $$\prtl{T}{t} = -\frac{1}{\rho c_p}\prtl{q}{x} - \frac{1}{c_p}\sum_kh_k\prtl{y_k}{t} + \frac{Q_r}{\rho c_p},$$ where $q$ is the heat flux, $h_k$ is the species enthalpy, $Q_r$ is the radiative source term (if radiation is used), and $c_p$ is heat capacity.

This equation for temperature is equivalent to the enthalpy equation given by $$\prtl{h}{t} = -\frac{1}{\rho}\prtl{q}{x} + \frac{Q_r}{\rho}.$$ This equivalence was verified in the code implementation by direct comparison.

The heat flux $q$ is given by $$q = -k\prtl{T}{x} + \sum_kh_kj_k.$$

Cantera is used for all thermochemical and transport properties, including diffusivities, thermal conductivity, viscosity, heat capacity, density, and chemical reaction rates, using available detailed kinetic mechanisms. SI units are used for all quantities (except kmoles are used instead of moles).

Radiation

The radiative source term $Q_r$ is computed using an optically-thin approximation. The gas is spectral and several models are avalable, as implemented in RadLib [5]. For a spectral formulation with $n_{gg}$ gray gases, $Q_r$ is given by $$Q_r = -4\sigma\sum_{j=1}^{n_{gg}} k_j(a_jT^4-a_{j,0}T_0^4),$$ where $\sigma$ is the Stefan-Boltzmann constant, $k_j$ is the absorption coefficient of the $j^{th}$ gray gas, $a_j$ is the weight factor of the $j^{th}$ gray gas, and $T_0$ is the temperature of the lower boundary, corresponding nominally to the ambient surrounding air. $k_j$ and $a_j$ are computed by RadLib, which can use the Plank-mean assumption with $n_{gg}=1$, the Weighted Sum of Gray Gases model (WSGG) with $n_{gg}=4$, or the Rank Correlated Spectral Line WSGG (RCSLW) model with arbitrary $n_{gg}$. By default, a Planck-mean assumption is used considering species H $_2$O, CO $_2$, CO, and CH $_4$. In that case, $Q_r$ is given by $$Q_r = -4\sigma k(T^4-T_0^4).$$

Soot

The equation for soot transport is give by

$$\prtl{M_k}{t} = -\prtl{j_{M,k}}{x} + {S_{M_k}},$$

Where $M_k$ is the $k^{th}$ soot mass-moment. Alternatively, it can be defined as the soot section $k$ of size $M_k$ and corresponding to the number of soot particles per $\text{m}^3$. ${S_{M_k}}$ is the soot source term and is taken directly from SootLib [6]. It includes nucleation, growth, oxidation, and coagulation. $M_k$ is defined as

$$ M_k = \int m^k n(m) dm $$

where $n(m)$ is the number of soot particles per $\text{m}^3$ per kg of soot. Additionally, $m$ is the mass in kg per soot particle. The value of $M_k$ when $ k = 0 $ is taken to be the number of soot particles per $\text{m}^3$. Likewise, when $ k = 1 $, $ M_k $ is representative of the kg of soot per $\text{m}^3$, which is equal to the denisty of the soot times the mass fraction of the soot. Further values of $M_k$ can be shown, but do not have as interesting of physical interpretations. The average size of a soot particle, $\langle m \rangle $ is defined as $\frac{M_1}{M_0}$, or the average mass in kg per soot particle. Each soot mass-moment is spaced from its predecessor by about 20 orders of magnitude. Thus, to help solve this equation, the soot moments are multiplied by a scaling factor so their values are all about the same order of magnitude. This helps the solver converge, and after the solve, the mass moments are divided by the same scaling factor to revert them back to their original order of magnitude.

The flux of each soot mass moment, $j_{M_{k}}$ is given by

$$ j_{M_{k}} = -0.556 \nu M_k \frac{\nabla T}{T}, $$

where $\nu$ is the kinematic viscosity of the gas. This makes the soot transport equation a hyperbolic PDE, and therefore we solve it by upwinding.

Numerical solution

The partial differential equations given above are solved using a finite volume (FV) formulation on a nonuniform spatial grid. The species equation at grid cell $i$ is given by $$\prtl{y_{k,i}}{t} = -\frac{j_{k,i}^e-j_{k,i}^w}{\rho_i\Delta x_i} + \frac{\doot{m}_{k,i}^\tprime}{\rho_i},$$ where the superscripts $e$ and $w$ denote east and west faces of the grid cell, respectively.

The FV energy equation is given by $$\prtl{T_i}{t} = -\frac{q_i^e-q_i^w}{\rho_i c_{p,i}\Delta x_i} - \frac{1}{c_{p,i}}\sum_kh_{k,i}\prtl{y_{k,i}}{t} + \frac{Q_{r,i}}{\rho_i c_{p,i}}.$$

Default configuration

The default configuration corresponds to the composition of the TNF piloted nonpremixed jet flames (flames C, D, E, and F). The case is set up to use the GRI3.0 gas mechanism, which consists of 53 gas species and 325 reactions. The domain size is 4 cm, but this is meant to be varied.