3.3.2. Particle Transport Model

In this section, the equations describing the evolution of parcels of particles are presented. Models are presented for the particle motion and for heat and mass transfer (ie. evaporation and combustion).

3.3.2.1. Particle Acceleration and Trajectories

Particles with densities much greater than that of a the fluid phase (solid or liquid particles in gaseous flows) are primarily affected by drag forces and body forces. In this limit where ${\rho}_p \gg {\rho}_g$ , the particle acceleration is written [101]:

(3.622) $\dfrac {du_{p,i}}{dt} = \dfrac {3 {\rho}_g C_D {|{\bf u}_g - {\bf u}_p}|} {4 {\rho}_p d_p} (u_{g,i}-u_{p,i}) + \left(\dfrac {{\rho}_p - {\rho}_g} {{\rho}_p}\right) g_i$

where $u_{p,i}$ and $u_{g,i}$ are the $i^{th}$ component of the particle and gas velocities, respectively, $|{\bf u}_g - {\bf u}_p|$ is the vector magnitude of the velocity differences, ${\rho}_p$ and ${\rho}_g$ are the particle and gas densities, and $g_i$ is the $i^{th}$ component of the acceleration due to body forces. The particle diameter is $d_p$ and this should be considered to be the equivalent particle diameter corresponding to a spherical particle. The effects of non-sphericity on the acceleration can be accounted for through the drag coefficient, $C_D$ . The gas-velocity to be employed in (3.622) is taken from the Eulerian solution of the continuum field. For turbulent flows, the sum of the mean (resolved) velocity and a perturbation to that mean, accounting for the turbulent fluctuations, both contribute to the gas velocity. The effects of turbulent fluctuations are described in the following section, Particle Dispersion and Turbulence.

For a spherical particle, the drag coefficient is modeled using standard drag coefficient relations:

(3.623) $C_D = \begin{cases} 24(1 + Re_p) ^ {2/3}/Re_p \hfill \qquad for Re_p < 1000 \\ 0.424 \hfill \qquad for Re_p > 1000 \end{cases}$

The particle Reynolds number, $Re_p$ , is based on the slip velocity and the particle diameter

(3.624) $Re_p = \dfrac {{\rho}_g d_p {|{{\bf u}_g - {\bf u}_p}|}} {{\mu}_g}$

Yuen and Chen [116] recommend evaluating the viscosity in (3.624) based on the weighted average of the properties at the gas-phase side of the particle surface (weighted by two-thirds) and the gas-phase properties far from the particle (weighted by one-third), the so called ‘1/3 rule.’ So for properties at the surface and far field identified with a subscripted $f$ and $\infty$ , respectively, the viscosity would be

(3.625) ${\mu}_g = {\mu}_{\infty}/3 + 2{\mu}_f/3.$

A similar relationship is suggested for other transport coefficients (conductivity, diffusivity, etc.). Note that additional forces are relevant for particles with densities nearer to or less than the continuum phase (bubbles) and for particles with high Reynolds numbers. A comprehensive overview of the force on particles is available from Maxey and Riley [107]. Equation (3.622) can be linearized and written

(3.626) $\dfrac {du_{p,i}} {dt} = \dfrac {(u_{g,i}-u_{p,i})} {{\tau}_p} + \left(\dfrac {{\rho}_p - {\rho}_g} {{\rho}_p} \right) g_i$

(3.627) ${\tau}_p = \dfrac {4 {\rho}_p d_p} {3 {\rho}_g C_D {|{{\bf u}_g - {\bf u}_p}|}}$

where ${\tau}_p$ is the particle velocity response time. In the small Reynolds number limit where the drag coefficient is inversely proportional to the slip velocity, this linearization is particularly relevant. Given the particle acceleration from Eq. (3.622), the particle trajectories can be determined by integrating the simple ODE

(3.628) $\dfrac {dx_{p,i}} {dt} = u_{p,i}$

Since all the particles within a parcel are the same size, all parcel trajectories are determined by Eqn. (3.628), subject to turbulence effects described below in Particle Dispersion and Turbulence.

3.3.2.2. Particle Dispersion and Turbulence

In CFD modeling of turbulent flows, the full velocity spectrum is generally not resolved. Instead, certain velocity fluctuations are modeled, being represented through the turbulent kinetic energy, $k$ , which is half the sum of the squares of the velocity fluctuations. These velocity fluctuations tend to introduce random fluctuations in the particle velocities that result in real particles being dispersed relative to the mean continuum flow [117, 118]. For Lagrangian particle methods, this phenomenon is modeled in two ways: by perturbing the velocity of parcels of particles and by affecting the spatial extent of the parcel itself.

To account for the effects of the velocity fluctuations on the parcels or particles, the random walk model of Gosman and Ioannides [111], as modified by Shuen et al. [112], is employed. In this approach, the gas velocity employed in equations (3.622), (3.624), (3.626), and (3.627) is the sum of the mean gas velocity, $\langle u_{g,i} \rangle$ , and a fluctuating velocity that is sampled from a normal (Gaussian) velocity distribution with a standard deviation given by ${\sigma}_u = \sqrt{2k/3}$ . For LES, $k$ here would be replaced with the subgrid kinetic energy. Sampling from the inverted cumulative distribution function with a random number uniformly distributed between zero and unity, $RN$ , gives the appropriate fluctuating velocity. The total and fluctuating gas velocities are then

(3.629) $u_{g,i} = \langle u_{g,i} \rangle + u'_{g,i}$

(3.630) $u'_{g,i} = \sqrt{2} {\sigma}_u erf^{-1}(2 RN -1)$

where $erf^{-1}$ is the inverse error function. The time during which a given velocity fluctuation affects a given parcel is determined by the expected time that the particle takes to cross the eddy inducing the given velocity fluctuation. Small particles will tend to stay in an eddy for the duration of the eddy lifetime,

(3.631) ${\tau}_e = \sqrt{3/2}C_{\mu}^{3/4} k/{\epsilon}$

where $\epsilon$ is the turbulent energy dissipation rate and $C_{\mu} = 0.09$ . Larger particles will have sufficient slip velocity to cross the eddy. The eddy-crossing time is estimated as

(3.632) ${\tau}_C = -{\tau}_p \: ln \left[1 - \dfrac {L_e} {{\tau}_p {|{{\bf u}_g - {\bf u}_p }|}} \right]$

where the eddy length scale is $L_e = C_{\mu}^{3/4} k^{3/2}/{\epsilon}$ . The eddy interaction time, that is the time over which a given velocity perturbation affect $u_{g,i}$ in Equations (3.629) and (3.630), is the minimum of ${\tau}_e$ and ${\tau}_C$ . The particles comprising a parcel are presumed to be distributed about the center of the parcel, tracked by Eq. (3.628), in a normal (Gaussian) manner with the standard deviation in each direction given by ${\sigma}_i$ , where $i$ is either $x$ , $y$ , or $z$ . This distribution is written

(3.633) $f_{\sigma}(\vec{x}; \vec{x_o}, t) = \dfrac {N_p} {(2{\pi})^{3/2} {\sigma}_x {\sigma}_y {\sigma}_z} exp \left({ -\left[ \dfrac {(x-x_o)^2} {2{{\sigma}_x}^2} + \dfrac {(y-y_o)^2} {2{{\sigma}_y}^2} + \dfrac {(z-z_o)^2} {2{{\sigma}_z}^2} \right]}\right)$

where $N_p$ is the total number of particles in the parcel and $\vec{x_o}$ is the center of the parcel. The spatial extent of the parcel is thus determined by ${\sigma}_i$ . This term is affected by unresolved turbulent fluctuations that act differently on the particles across the parcel. Following Zhou and Yao [113], the spatial extent of the parcel is the mean square displacement over time

(3.634) ${{\sigma}_x}^2 = \sum \left[{{u'_{p,i}}^2 ({\Delta}t_i)^2}\right]$

where $u'_{p,i}$ satisfies the ODE

(3.635) $\dfrac {d{u'}_{p,i}} {dt} = \dfrac {({u'}_{g,i}-{u'}_{p,i})} {{\tau}_p}$

and ${\Delta}t_i$ is the time over which ${{u'}_{p,i}}$ acts. Equation (3.635) is obtained by subtracting the instantaneous particle equations from the mean particle equations. Note that large particles with large ${\tau}_p$ are not dispersed appreciably.

3.3.2.3. Mass and Energy Exchange between Particles and the Gas Phase

Particles may exchange mass and energy with the gas phase according to conservation principles across the interface. The physics of mass and energy transfer are described in detail here because the anticipated applications include a wide range of physical phenomena that have not been described together. Included in the phenomena of interest are metal and hydrocarbon particle combustion as well as particle condensation and evaporation. Each of these is anticipated to be strongly energetic in the sense that the product of the evaporation rate with the enthalpy of evaporation and combustion is expected to be substantial. Further, for evaporation and condensation applications, it is expected that vapor pressures will range over a sufficiently wide range that the transition from evaporation to condensation should be correctly described as this often limits evaporation and condensation rates. For metal combustion, the associated temperatures are sufficiently high so that radiative heat transfer must be considered.

The theory for droplet vaporization and combustion has generally been developed based on a large number of simplifications [95, 96] including spherical-symmetry, unity Lewis numbers, droplet surfaces at the boiling temperature, and infinitely fast conduction through the droplet. Recent work provides guidance as to how these assumptions can be relaxed [119].

3.3.2.3.1. Theory for spherically symmetric flow

In general, the theory of heat and mass transfer is developed for spherically symmetric systems and then corrected to account for increased transfer associated with advection and asymmetry. In this section, relations are developed based on the assumption of spherical symmetry and corrections is provided in Extension to multiple oxidizers.

(3.636) $\dfrac {{\partial} \rho} {\partial t} + \dfrac {1} {r^2} \dfrac {{\partial}} {{\partial}r}(\rho r^2 \nu) = 0$

can be integrated to give

(3.637) $\dot{m_o} = 4 \pi \rho r^2 \nu$

where $\dot{m_o}$ is the rate of mass evaporation from the particle, $\nu$ is the Stefan velocity, directed normally away from the particle, and $\rho$ is the local vapor density considering both the particle vapor and the continuum gas concentrations. A coordinate transformation to the variable

(3.638) ${\xi}_T = \int_{r_o}^{\infty} {\dfrac{1} {4 \pi r^2 ({\lambda}/{c_p})} } dr$

or

(3.639) ${\xi}_F = \dot{m_o} \int_{r_o}^{\infty} {\dfrac {1} {4 \pi r^2 (\rho D_F)}} dr$

which represents the ratio of the Stefan velocity to the thermal diffusion velocity, greatly simplifying the species and energy conservation equations. In Equations (3.638) and (3.639), $\lambda$ is the thermal conductivity of the vapor, $c_p$ is its specific heat at constant pressure, and $D_F$ is the evaporating species diffusion coefficient. The subscript 0 on the evaporation rate indicates that this evaporation corresponds to that for the spherically-symmetric case; corrections relating the evaporation rate for the spherically-symmetric case to that with finite-slip velocities are provided in Extension to multiple oxidizers. In the spherically-symmetric case, conservation equations for conserved scalars, ${\beta}_k$ , can be written

(3.640) $\dfrac {\partial {\beta}_k} {\partial {\xi}_k} + \dfrac {{\partial}^2 {\beta}_k} {\partial {{\xi}_k}^2} = 0$

for which an analytic solution

(3.641) ${\beta}_k = C_1 + C_2 e ^ {-{\xi}_k}$

is readily obtained. In non-reacting, non-radiating flows, any mass fraction or temperature can be a conserved scalar. Other conserved scalars are provided below. The appropriate from of (3.638) and (3.639) to be used depends on the variable represented by (3.640). Since the diffusion coefficients appearing in (3.639) are most important near the particle surface, the diffusion coefficient to be employed is that most relevant at the surface. The choice will be clearly identified below. The application of the Dirichlet boundary conditions at both the surface and far from the droplet and the application of a Neumann boundary condition at the surface relate the boundary conditions for the conserved scalar to $\dot{m_o}$ through ${\xi}_k$

(3.642) ${\xi}_{k,f} = \dot{m_o} \int_{r_o}^{\infty} {\dfrac {1} {4 \pi r^2 (\rho D)}} dr = \ln \left(1 + \dfrac {{\beta}_{\infty} - {\beta}_f} {- {\dfrac {d \beta} {d{\xi}_k}}|_f }\right)$

where the subscript $f$ indicates quantities evaluated at the droplet surface, the so-called film state. For reacting flows we will consider two such conserved scalars in the present work. The general configuration considered is a fuel droplet reacting in an oxidizing atmosphere; in the event that, for example, an oxidizer droplet is reacting in a fuel atmosphere, then the “oxidizer” and “fuel” described would be switched. Also, if there is no reaction, the results generalize to droplet evaporation where the “fuel” is the evaporating component. Allowing variable properties, the temperature oxidizer coupling function is

(3.643) ${\beta}_{T-O} = \int_{T_o}^T {c_P} dT + \dfrac {Y_O W_F q_{comb}} {{\nu}_O W_O}$

where $W_i$ is the molecular weight and ${\nu}_i$ is the stoichiometric coefficient of species $i$ . The standard enthalpy of combustion for fuel and oxidizer, per unit mass of fuel evaluated at the film temperature, $T_f$ , is

(3.644) $q_{comb} = h_{F,f} + \dfrac {h_{O,f} W_O {\nu}_O} {W_F {\nu}_F} - \sum \left[ \dfrac {h_{p,f} W_p {\nu}_p} {W_F {\nu}_F} \right]$

where the last summation is taken over the produces of reaction, $p$ . Note that when there is no combustion, the second term in Eqn. (3.644) is ignored. For ${\beta}_{T-O}$ , the thermal diffusivity is the relevant diffusivity so that Eqn. (3.638) is employed in conjunction with ${\beta}_{T-O}$ .

The constants $C_1$ and $C_2$ in (3.641) are evaluated using the boundary conditions for the temperature and oxidizer at the droplet surface and in the far field. The boundary conditions employed for temperature are that the heat flux into the particle is balanced by the sum of the enthalpy of vaporization, the heating of the particle and any radiative losses

(3.645) $4 \pi r^2 \lambda \left(\left.\dfrac {dT} {dr} \right\vert_{f} \right) = -\dot{m_o} c_{p,f} \left(\left.\dfrac {dT} {{\xi}_T} \right\vert_{f} \right) = \dot{m} h_{vap} + Q_{rad} + m_p c_{v,p} \dfrac {dT_p} {dt}$

Here, the enthalpy of vaporization is $h_{vap}$ , the particle specific heat is $c_{v,p}$ , and the radiative heat loss over the droplet surface is

(3.646) $Q_{rad} = 4 \pi {r_p}^2 \alpha (\sigma {T_p}^4 - G_{in}/4)$

where $\alpha$ is the particle absorptivity, $\sigma$ is the Stefan-Boltzmann constant, and $G_{in}$ is the incident radiation. The incident radiation is the radiation intensity integrated over all directions, that is, the entire $4 \pi$ steradian solid angle. If the radiation intensity is $I$ , then $G_{in} = \int_{4 \pi} I d\Omega$ where $d\Omega$ is the differential solid angle. If the particle is not opaque, the particle absorptivity will be a function of the particle size [120]. Note that, in Eqn (3.645), the evaporation rate appears without the subscript 0, indicating that this is the evaporation rate corrected for finite-slip velocities as prescribed in (3.661) below. This is appropriate because the heat and mass flux to the surface are both increased by the relative droplet motion, while the other terms in (3.645) are not affected by that. The oxidizer is presumed to not be absorbed by the surface so that a no-flux boundary condition is employed

(3.647) $\dfrac {dY_O} {d{\xi}_T} = -Le_{O,f} Y_{O,f}.$

Because the thermal diffusivity is used in definition ${\xi}_T$ for ${\beta}_{T-O}$ , the ratio of the thermal to mass diffusivity in the form of the oxidizer Lewis number appears in (3.647). Taking the derivative of Eqn. (3.643) evaluated at the droplet surface and using (3.645) and (3.647), we obtain

(3.648) $\left(\left.\dfrac {d{\beta}_{T-O}} {d{\xi}_T} \right\vert_{f}\right) = - h_{vap} - \dfrac {Q_{rad}} {\dot{m}} - \dfrac {m_p c_{v,p}} {\dot{m}} \dfrac {dT_p} {dt} + \int_{T_o}^{T} {\dfrac {dc_{p,f}} {d{\xi}_T}} {dT} - \dfrac {Y_{O,f} W_F q_{comb}} {Le_{O,f} {\nu}_O W_O}$

that will appear in the denominator of Eqn. (3.642). Also the temperature and oxidizer mass fraction must approach their far field values at large radii. Additional assumptions are required to identify the temperature and the oxidizer mass fraction at the surface; these will be discussed later. Applying these boundary conditions to (3.641) for ${\beta}_{T-O}$ provides an expression for the rate of evaporation in terms of ${\xi}_T$ evaluated at the surface of the particle

(3.649) ${\xi}_{T,f} = \dot{m_O} \int_{r_o}^{\infty} {\dfrac {1} {4 \pi r^2 (\lambda/ c_p)}} {dr} = \ln{[1+B_{T-O}]}$

where the Spalding transfer number associated with ${\beta}_{T-O}$ is

(3.650) $B_{T-O} = \dfrac {\int_{T_O}^{T} {c_{p,\infty}} dT - \int_{T_O}^{T} {c_{p,f}} dT + q_{comb} \left[\dfrac {(Y_{O,\infty} - Y_{O,f} ) W_F} {{\nu}_O W_O} \right] } {h_{vap} + \dfrac {Q_{rad}} {\dot{m}} + \left(\dfrac {m_p c_{v,p}} {\dot{m}} \dfrac {dT_p} {dt} \right) - \int_{T_O}^{T} {\dfrac {dc_{p,f}} {d{\xi}_T} } dT + \dfrac {Y_{O,f} W_F q_{comb}} {Le_{O,f} {\nu}_O W_O} }$

Here the subscripts $f$ and $\infty$ indicate the states at the particle surface (on the gas-phase side of the interface) and the ambient far-field environment, respectively. The denominator in (3.650) represents a variety of potential sinks for the enthalpy at the droplet surface. These sinks include, in the order in which they are written, the enthalpy associated with vaporizing the particle, the radiative losses from the surface, the enthalpy conducted into the particle, the enthalpy flux to the gas-phase due to variable specific heats, and the enthalpy of combustion lost as a consequence of oxidizer leakage. Note that the radiative absorption and emission from the particle surface, $Q_{rad}$ , are included here in the surface boundary condition, but the radiative losses from a flame around the droplet must be accounted for by modifying $q_{comb}$ to provide an effective heat release, the heat release decremented by the near flame radiative losses. As far as the droplet is concerned, these radiative losses are far enough away to not affect the film state except to the extent that radiative flux to the surface is affected.

The above expressions comprise a relatively complete definition of the physics of particle evaporation, combustion, and interaction with a radiative field accounting for variable thermophysical properties. These expressions are simplified by making the assumption that no oxidizer penetrates a flame to reach the surface, $Y_{O,f} = 0$ , and by setting $T_f = T_O$ leading to

(3.651) $B_{T-O} = \dfrac {\int_{T_f}^{T_{\infty}} {c_{p,\infty}} dT + \dfrac {Y_{O,\infty} W_F q_{comb}} {{\nu}_O W_O} } {h_{vap} \dfrac {Q_{rad}} {\dot{m}} + \left( \dfrac {m_p c_{v,p}} {\dot{m}} \dfrac {dT_p} {dt} \right) }$

where the denominator is referred to as the effective enthalpy

(3.652) $h_{eff} = h_{vap} + \dfrac {Q_{rad}} {\dot{m}} + \left(\dfrac {m_p c_{v,p}} {\dot{m}} \dfrac {dT_p} {dt} \right).$

Equation (3.651) is employed in the numerical models. A coupling function for fuel and oxidizer is similar written

(3.653) ${\beta}_{T-O} = \dfrac {Y_F} {W_F} - \dfrac {Y_O} {{\nu}_O W_O}$

For a species that is evaporating, the flux boundary condition is

(3.654) $\left(\left.\dfrac {dY_F} {d{\xi}_F} \right\vert_{f} \right) = Y_{F,p} - Y_{F, f}$

and the fuel-oxidizer Spalding mass transfer number is

(3.655) $B_{F-O} = \dfrac {Y_{F,f} - Y_{F,\infty} + \dfrac {Y_{O,\infty} W_F} {{\nu}_O W_O} } {Y_{F,p} - Y_{F,f}}$

where $Y_{O,f} = 0$ has been assumed as in Eq. (3.651). Note that only one of $Y_{O,\infty}$ or $Y_{F,\infty}$ will be non-zero based on the current assumption of zero leakage through flames; if combustion is occurring then it will be $Y_{F,\infty}$ that is zero. The diffusion coefficient appropriate for the fuel-oxidizer coupling function is that for the fuel so that the second of Eqn (3.638) and (3.639) is used with the diffusion coefficient specifically that of the fuel, and the equivalent of Eqn. (3.649) for the fuel-oxidizer system

(3.656) ${\xi}_{F,f} = \dot{m_o} \int_{r_o}^{\infty} {\dfrac {1} {4 \pi r^2 (\rho D_F)}} = \ln [1+B_{F-O}]$

3.3.2.3.2. Extension to multiple oxidizers

Multiple oxidizers is discussed in detail later, in Multiple Oxidizers.

3.3.2.3.3. Correlations for finite slip velocities

The above relationships for the heat and mass transfer are derived for a spherically symmetric field around the droplet and are valid for droplets with zero slip velocity in the absence of buoyancy. Empirical correlations are available in terms of the Nusselt and Sherwood numbers parameterized by Reynolds, Schmidt, and Prandtl numbers to describe the effect of finite slip velocities in modifying the heat and mass transfer by reducing the boundary layer thickness. The Schmidt and Prandtl numbers are defined:

(3.657) $Pr = \dfrac {c_{p,g}{\mu}_g} {{\lambda}_g}$

(3.658) $Sc = \dfrac {{\mu}_g} {{\rho}_g D_g}$

These quantities with the subscript $g$ represent appropriate averages for transport properties in the gas phase boundary layer around the particle.These gas-phase quantities are evaluated using an appropriate averaging process that will generally be analogous to Eqn. (3.625). The alternative is to tabulate these quantities using Equations (3.665) and (3.666). The Nusselt number describes a dimensionless heat transfer rate to the droplet for a given difference between the ambient and surface temperature,

(3.659) $Nu_f = \left(\left. 2 r_p c_{p,g} \left.\dfrac {dT} {dr} \right\vert_{f} \right) \middle/ \left(\int_{T_f}^{T_{\infty}} c_{p,\infty} dT + \dfrac {Y_{O,\infty} W_F q_{comb}} {{\nu}_O W_O} \right.\right)$

The correction to the evaporation rate employed in the present work is based on measurements by Ranz and Marshall (1952)

(3.660) $\dfrac {Nu_f} {Nu_{f, Re=0}} = (1+0.3 Re ^ {1/2} Pr ^ {1/3})$

This can be introduced into (3.650), and the evaporation rate can be written

(3.661) $\dot{m} = \dfrac {Nu_f} {Nu_{f, Re = 0}} \dot{m_o}$

The Nusselt number for zero slip velocity, $Nu_{f, Re=0}$ , is included in the relations of the previous section, specifically in Eq. (3.650). A similar correlation,

(3.662) $\dfrac {Sh_f} {Sh_{f, Re=0}} = (1+0.3 Re ^ {1/2} Sc ^ {1/3})$

can be used for the Sherwood number,

(3.663) $Sh_f = \left(\left. -2 r_p \left.\dfrac {dY_F} {dr}\right\vert_{f} \right) \middle/ \left(Y_{F,f} - Y{F,\infty} + \dfrac {Y_{O,\infty} W_F} {{\nu}_O W _O} \right.\right)$

which is the dimensionless mass transfer coefficient, so that the evaporation rate for finite slip velocities can also be written

(3.664) $\dot{m} = \dfrac {Sh_f} {Sh_{f,Re=0}} \dot{m_o}$

where ${\dot{m_o}}$ should be taken from Eqn. (3.656).

3.3.2.3.4. Evaporation rates and effective diffusivities

The evaporation rate is seen in Eqn (3.649) and (3.656) to depend linearly on an area-weighted gas-phase diffusivity. Since accurate evaluation of this weighting is not feasible for the particle transport model, effect diffusivities,

(3.665) $\left(\dfrac {\lambda} {c_p} \right)_{eff} = \left[4 \pi r_p \int_{r_o}^{\infty} {\dfrac {1} {4 \pi r^2 (\lambda/c_p)}} {dr} \right]^{-1}$

(3.666) $\left(\rho D_F \right)_{eff} = \left[4 \pi r_p \int_{r_o}^{\infty} {\dfrac {1} {4 \pi r^2 (\rho D_F)}} {dr} \right]^{-1}$

are defined. The one-third rule defined in Eq. (3.625) provides a rough guideline for evaluating these diffusion coefficients in some cases. The effective thermal diffusion coefficient is combined with Equations (3.649) and (3.661) to give an evaporating rate based on thermal diffusion

(3.667) $\dot{m} = 4 \pi r_p \left(\rho D_F \right)_{eff} \dfrac {Sh_f} {Sh_{f,Re=0}} \ln \left[1+ B_{T-O} \right].$

The strong dependence of the evaporating rate on the diffusion coefficient, coupled with the fact that the diffusion coefficients depend strongly on variations in the compositions and temperature of the gases around the droplet, mean that the reasonable but judicious choice of diffusion coefficients can often match observed experimental measurements. Similarly, the mass transfer driven evaporation rate can be written with the use of Equations (3.656), (3.662), and (3.664)

(3.668) $\dot{m} = 4 \pi r_p \left(\rho D_F \right)_{eff} \dfrac {Sh_f} {Sh_{f,Re=0}} \ln \left[1+B_{T-O} \right]$

It is necessary that the evaporation rate predicted by Eqn. (3.668) be equal to that predicted by (3.667). These equations show that thermal and mass diffusion vary both through their diffusion coefficients and through difference in boundary layer thickness attributable to finite slip velocities. Both these effects are combined in an effective Lewis number, ${Le_{eff}}$ , to give

(3.669) $Le_{eff} = \dfrac {\int_{r_o}^{\infty} {\dfrac {dr} {4 \pi r^2 (\rho D_f )}} {\dfrac {Nu_f} {Nu_{f,Re=0}}}} {\int_{r_o}^{\infty} {\dfrac {dr} {4 \pi r^2 (\lambda/c_p)}} {\dfrac {Sh_f} {Sh_{f,Re=0}}}}$

Equations (3.667), (3.668), and (3.669) are used in the computational model for the evaporation rate. Note that energy and mass conservation must vie the same evaporation rate; this requirement determines the film state as described in the following section, Closure of film state with effective heat transfer coefficient.

3.3.2.3.5. Closure of film state with effective heat transfer coefficient

The system described by Equations (3.667) and (3.668) is closed with two additional assumptions. First, the film conditions, $T_f$ and $Y_f$ , are related through the Clausius-Clapeyron relationship

(3.670) $P_{F,f} = P_{ref} exp \left[ - \dfrac {h_{vap}} {R} \left( \dfrac {1} {T_f} - \dfrac {1} {T_{ref}} \right) \right]$

where the partial pressure gives the mole fraction through $X_F = P_{F,f}/P$ that can subsequently be converted to the mass fraction with the relationship $Y_F = X_F W_F / {\Sigma}_k X_k W_k$ . As provided by Lefebvre [121], such a relationship is

(3.671) $h_{vap} = h_{vap,ref} \left(\dfrac {T_{crit}-T_f} {T_{crit}-T_{ref}} \right)^{0.38}$

for $T_f < T_{crit}$ and zero otherwise. If the critical point temperature is not provided, the code sets $T_{crit}$ to a very large value, essentially making $h_{vap}$ independent of temperature. Second, the droplet heating is related to the difference between the film temperature and the droplet temperature by assuming an effective internal heat transfer coefficient in the form of an internal Nusselt number, $Nu_p$ , for the particle so that

(3.672) $m_p c_{v,p} \dfrac {dT_p} {dt} = 2 \pi r_p Nu_p {\lambda}_p (T_f - T_p).$

This internal Nusselt number, which is different from that for the external heat transfer indicated in Eqn. (3.661), can be estimated based on the results of numerical studies where the internal droplet was resolved [119, 122]. There the Nusselt numbers for no circulation and for rapid circulation were identified as 6.58 and 17.9, respectively, and a transition region was also identified based on the liquid Peclet number

(3.673) $Pe_p = \dfrac {2 {\rho}_p c_{v,p} U_{surface} r_p} {{\lambda}_p}$

where

(3.674) $U_{surface} = \dfrac {12.69 {|{ {\bf u}_p - {\bf u}_g }|} Re_p ^ {1/3}} {16} \left(\dfrac {{\mu}_g} {{\mu}_p} \right)$

which is based on the maximum surface velocity, $U_{surface}$ . This transition was empirically fitted [119] to

(3.675) $Nu_p = 6.58 \left[1.86 + 0.86 tanh \left[2.245 log_{10} \left(Pe_p/30 \right) \right] \right].$

In [119], both the evaporation rate and the surface temperature were reproduces using Eq. (3.675) with little error compared to simulations incorporating a detailed internal droplet convection model. Naturally, the liquid Peclet number should be zero if the particle is below the freezing temperature of the particle constituent.

The evaporation rates indicated in Equations (3.667) and (3.668) must be equal, subject to the constraints of Equations (3.670) and (3.672) based on the closure approximation employed in this section. Equating Equations (3.667) and (3.668) leads to a nonlinear equation for the surface temperature that is to be solved. This is readily solved using Newton’s method as described here. Newton’s method is an iterative root-finding method written in the form

(3.676) $T_f^{n+1} = T_f^n - g(T_f^n)/g'(T_f^n)$

where $g(T_f) = 0$ is the equation for which the root will be found and the superscript $n$ refers to the iteration number. Equating Equations (3.667) and (3.668) and using (3.669) gives

(3.677) $g(T_f) = (1+B_{T-O})^{Le_{eff}}-1-B_{F-O}$

Differentiation gives

(3.678) $g'(T_f) = Le_{eff}(1+B_{T-O})^{(Le_{eff}-1)} \dfrac {dB_{T-O}} {dT_f} - \dfrac {dB_{F-O}} {dT_f}$

(3.679) $\dfrac {dB_{T-O}} {dT_f} = \dfrac {- \left[c_{p,\infty}(T_{\infty}-T_f) + \dfrac {Y_{O,\infty} W_F q_{comb}} {{\nu}_O W_O} \right] \dfrac {dh_{eff}} {dT_f} - c_{p,\infty} h_{eff}} {h_{eff}}$

(3.680) $\dfrac {dB_{F-O}} {dT_f} = \dfrac {\left(Y_{F,f}-Y_{F,\infty}+\dfrac {Y_{O,\infty}W_F} {{\nu}_O W_O} \right)} {(Y_{F,p}-Y_{F,f})^2} \dfrac {dY_{F,f}} {dT_f}$

(3.681) $\dfrac {dY_{F,f}} {dT_f} = \dfrac {-W_F W_g \dfrac {P_{atm}} {P_{ref}} \dfrac {h_{vap}} {R} \left[ \dfrac {1} {{T_f}^2} + \dfrac {0.38} {(T_{crit}-T_f)} \left(\dfrac {1} {T_f} - \dfrac {1} {T_{ref}} \right) \right] exp \left[- \dfrac {h_{vap}} {R} \left(\dfrac {1} {T_f} - \dfrac {1} {T_{ref}} \right) \right] } {W_F - W_g + W_g \dfrac {P_{atm}}{P_{ref}} exp \left[- \dfrac {h_{vap}} {R} \left(\dfrac {1} {T_f} - \dfrac {1} {T_{ref}} \right) \right]}$

(3.682) $\dfrac {dh_{eff}} {dT_f} = \dfrac {2 \pi r_p Nu_p {\lambda}_p} {\dot{m}} - \dfrac {Le_f \left[Q_{rad} + 2 \pi r_p Nu_p {\lambda}_p (T_f - T_p) \right]} {\dot{m} (1+B_{F-O}) \left[\ln (1+B_{F-O}) \right]} \dfrac {dB_{F-O}} {dT_f}$

In writing these expressions for $g'$ , it is assumed that $Y_{O,f}=0$ and $T_O = T_f$ .

Because of the strong nonlinearities in $g(T_f)$ , care must be taken in providing initial conditions to solve these relations. This is conducted through a two-stage procedure. First, the boiling point temperature is identified, then a temperature just under the boiling temperature is used as an initial guess for the iterative solution that determines the film temperature. This procedure arises from a consideration of the shape of the $g(T_f)$ . For realistic temperatures, those for which $0 < Y_f < 1$ , both $g(T_f)$ and $g'(T_f)$ , are monotonically strongly increasing in magnitude. An initial guess with a temperature that is too low (less than $T_f$ ) results in a prediction, with the Newton’s method, of a very high temperature on the successive iteration due to the small magnitude of the derivative, $g'$ , for a small T. Typically, this second iteration will result in a temperature for which $Y_f > 1$ that leads to negative values of $B_m$ and the iteration diverges into non-physical regimes. However, an initial guess that is within the physically reasonable regime ( $0 < Y_f < 1$ ) and yet above the final $T_f$ will reliably converge. Therefore, the initial guess of $T = 0.99999 T_{boil}$ is used as an initial guess in determining the film temperature. This method has been tested for a wide range of conditions and appears robust except for those scenarios where the denominator of $B_T$ takes on negative values. (Since $T_f$ can be less than $T_p$ and this can result in the denominator of $B_T$ taking on negative values, in which case the iteration may fail.)

Because of the cost and potential stability issues in this method, the fast conduction limit described in Closure for surface state assuming fast conduction is used in the code instead.

3.3.2.3.6. Closure for surface state assuming fast conduction

The surface state described in the previous section is the most physically realistic state that can be obtained without solving a differential equation for the heat transfer with the droplet. However, determining this state involves the iterative solution of a system of nonlinear equations. In the previous section a robust method of solving these equations is provided, but a simpler approximation may provide suitable results under certain conditions. This simpler approach is to assume that the heat transfer within the droplet is fast relative to the heating of the droplet. This is equivalent to taking the zero Biot number limit, $Bi \approx Nu_f {\lambda}_f / Nu_p {\lambda}_p \to 0$ , in Eqn. (3.672), in which case $T_f = T_p$ . With the film state determined by the droplet temperature, the film mass fraction is directly obtained from Equations (3.670) and (3.671), and the mass transfer number is obtained from Eq. (3.655) with the evaporation rate following from Eqn. (3.668). The thermal transfer number is obtained from the mass transfer number by equating Equations (3.668) and (3.667), and the droplet heating is obtained y solving for the enthalpy change of the particle in Eq. (3.650) to obtain

(3.683) $\dfrac {m_p c_{v,p}} {\dot{m}} \dfrac {dT_p} {dt} = -h_{vap} - \dfrac {Q_{rad}} {\dot{m}} + \dfrac {\int_{T_f}^{T_{\infty}} c_{p,\infty} dT + \dfrac {Y_{O,\infty} W_F q_{comb}} {{\nu}_O W_O} } {B_{T-O}}$

Sirignano [123] and coworkers have demonstrated that for many conditions, this particular limit is a poor approximation for at least some of the particle lifetime. In the present models, this limit is employed in two situations: (a) if the particle temperature is within 1% of the wet bulb temperature, and (b) if the particle temperature exceeds the wet bulb temperature. As the droplet temperature approaches the wet bulb temperature, employing this limit is inconsequential. In the latter case, the rate of droplet cooling is likely to be over predicted, but this is a scenario for which the convergence of the film state otherwise is not guaranteed. In the interest of creating a more robust model, and because this particular situation is not anticipated to be predominant, we employ this simpler limit.

3.3.2.4. Conserved scalars and transfer numbers for various applications

In Correlations for finite slip velocities, expressions are provided for the particle evaporation rate, (3.668) and (3.667), using a model system of fuel evaporating from the particle and reactive with an oxidizer that diffuses from the ambient gas. Of significance in those expressions are the transfer numbers defined in Theory for spherically symmetric flow in Equations (3.651) and (3.650), and these in turn are based upon conserved scalars defined in Equations (3.643) and (3.653). In this appendix, expressions for alternate transfer numbers and conserved scalars are provided for two additional systems: a simpler system in which evaporating or condensation occurs without any reaction in the boundary layer (ie. water droplet evaporation or condensation) and a more complicated system in which multiple oxidizers are involved in the oxidation of the evaporated fuel (relevant to metal oxidation).

3.3.2.4.1. Simple evaporation and condensation

When species evaporate or condense but do not otherwise react in the boundary layer surrounding the particle, the species mass fraction of the evaporating/condensing species itself is a conserved scalar. In this case, the mass transfer number corresponding to Eq. (3.650) is simply

(3.684) $B_{F-O} = \dfrac {Y_{F,f}-Y_{F,\infty}} {Y_{F,p} - Y_{F,f}}$

and the heat transfer number corresponding to (3.651) is

(3.685) $B_{T-O} = \dfrac {\int_{T_f}^{T_{\infty}} c_{p,\infty} dT } {h_{vap} + \dfrac {Q_{rad}} {\dot{m}} + \left(\dfrac {m_p c_{v,p}} {\dot{m}} \dfrac {dT_p} {dt} \right)}$

3.3.2.4.2. Multiple Oxidizers

For metal oxidation, it is possible to have multiple oxidizers that simultaneously react with the evaporating metal. These can be expressed with a series of parallel single-oxidizer reactions of the form $F + {\nu}_{O_i}O_i \to \Sigma {\nu}_{P_{i,j}} P_{j,i}$ . Useful conserved scalars that can be formed with this set of reactions include

(3.686) ${\beta}_{F-O} = \dfrac {Y_F} {W_F} - {\sum}_i \dfrac {Y_{O_i}} {{\nu}_O W_{O_i}}$

(3.687) ${\beta}_{T-O} = \int_{T_f}^{T} c_p dT + {\sum}_i \dfrac {W_F q_i} {{\nu}_O W_{O_i}} Y_{O_i}$

where the enthalpy of reaction of the $i^{\rm{th}}$ oxidizer, $O_i$ , with $F$ is given by $q_i$ . These are analogous to the conserved scalars defined in Equations (3.643) and (3.653). Using these conserved scalars, a mass transfer number analogous to Eqn. (3.650) is found to be

(3.688) $B_{F-O} = \dfrac {Y_{F,p}-Y_{F,\infty} + {\sum}_i \dfrac {W_F q_i} {{\nu}_{O_i} W_{O_i}}} {Y_{F,p}-Y_{F,f}}$

and the heat transfer number is

(3.689) $B_{T-O} = \dfrac {\int_{T_f}^{T_{\infty}} c_{p,\infty} dT + {\sum}_i \dfrac {Y_{O_i, \infty} W_F} {{\nu}_{O_i} W_{O_i}} } {h_{vap} + \dfrac {Q_{rad}} {\dot{m}} + \left( \dfrac {m_p c_{v,p}} {\dot{m}} \dfrac {dT_p} {dt} \right)}$

These can be employed in Eqn. (3.667) and (3.668) to provide expressions for the particle mass burning rate as a function of the various oxidizer mass fractions far from the particle.

3.3.2.5. Energy Exchange between Particles and the Gas Phase Without Mass Transfer

In this section, a simpler scenario is considered where evaporation from and condensation onto particles is presumed to be negligible. For example, metal particles in dry air at ambient temperatures are unlikely to participate in evaporation or condensation. Models are presented in this section to treat these scenarios.

The models described in the previous sections are ill posed to solve these problems because the formulation is based on a balance between diffusion and the Stefan convective velocity, which is proportional to the negligible $\dot{m}$ . The models in this section are used for the heated particles model.

In the event that there is no mass transfer, the particle heating rate is determined based on the balance between conductive and radiative transfer. The conductive heat transfer can be expressed using an effective heat transfer coefficient, thereby taking advantage of available Nusselt number correlations indicated in Simple evaporation and condensation. Then the relationship for the droplet heating is

(3.690) $m_p c_{v,p} \dfrac {dT_p} {dt} = 2 \pi Nu_f r_p {\lambda}_f (T_g - T_f) + 4 \pi \alpha r_p^2 (G_{in}/4 - \sigma T_f^4).$

For the particle with no mass transfer, the closure of the film (or surface) temperature is obtained by equating the external heat flux, described by the right hand side of Eqn. (3.690), with the internal heat flux, as indicated on the right hand side of Eqn. (3.672) to obtain the following quartic constraint, which is solved for $T_f$

(3.691) $Nu_f {\lambda}_f (T_g - T_f) + 2 \alpha r_p (G_{in}/4 - \sigma T_f^4) = Nu_p {\lambda}_p (T_f - T_p)$

3.3.2.6. Further Notes on Radiative Heat Transfer

The expression employed for radiative droplet heating in Eqn. (3.646) is appropriate for relatively large and opaque particles, referred to as the geometric optics limit. For different particles, the expression for the absorptivity will change, but otherwise the expressions remain appropriate. To identify the appropriate absorptivity, it is useful to compare absorption coefficients for particle clouds. The absorption coefficient, $\kappa$ , in units of inverse length, is an effective cross-sectional area per volume. For large particles with a condense phase absorptivity give by $\alpha$ , the absorption coefficient is the summation of the cross-sectional areas time the absorptivity

(3.692) $\kappa = \dfrac {1} {V_c} {\sum}_p \left[\pi r_p^2 \alpha \int_{V_c} f_{\sigma}({\bf x}; {\bf x_o}, t) d{\bf x} \right].$

For comparison, in the small-particle or Rayleigh limit, the absorption coefficient is proportional to the volume fraction with the proportionality coming from the complex index of refraction. For intermediate particles where the particle optical depth is comparable to the particle radius, intermediate limits are appropriate, and the absorption can range from being proportional to the particle area to being proportional to the particle volume. The appropriate absorption coefficient, and from it the particle emissivity, cases can be determined as described in the available texts [120, 124].

3.3.2.7. Input Parameters for Particle Evolution

A large number of parameters are required to specify the evolution of the particles. The parameters provided in the input file are described in table Input parameters related to the particle evolution provided through the input file. In addition, table Variables passed from the gas-phase Eulerian solver required to evolve the particles specifies those variables that must be obtained from the gas-phase continuum flow. Because parcels of particles have finite extent and may span several control volumes, it is sometimes necessary to interpolate values of these gas-phase variables from several control volumes.

Table 3.12 Input parameters related to the particle evolution provided through the input file
Variable	Input name	Units	Description
${\rho}_p$	INJP_DENp	$kg/ m^3$	Particle density
$T_p$	INJP_Tp	$K$	Particle temperature
$W_F$	INJP_MWp	$kmol/kg$	Molecular weight of fuel or particle component
$c_{v,p}$	INJP_Clp	$J/kg/K$	Particle specific heat
$\alpha$	INJP_Emp	none	Particle absorptivity,emissivity
$Pr$	INJP_Prf	none	Film Prandtl number
$Sc$	INJP_Scf	none	Film Schmidt number
$P_{ref}$	INJP_Prefvapp	$Pa$	Reference pressure for vaporization for particle
$T_{ref}$	INJP_Trefvapp	$K$	Reference temperature for vaporization for particle
$h_{vap,ref}$	INJP_Hvaprefp	$J/kg$	Reference enthalpy of vaporization for particle
$T_{crit}$	INJP_T1cp	$K$	Critical temperature for particle
$T_{fr}$	INJP_T1fp	$K$	Freezing temperature for particle
${\sigma}_p$	INJP_STp	$N/m$	Particle surface tension
${\mu}_p$	INJP_mup	$kg/m/s$	Particle viscosity
$Pr_p$	INJP_Prl	none	Particle Prandtl number
$q_{comb}$		$J/kg$	Enthalpy of combustion for vapor species at $T_f$
${\nu}_i$		mol/ mol Fuel	Stoichiometric coefficients (relative to evaporating species)

Table 3.13 Variables passed from the gas-phase Eulerian solver required to evolve the particles
Variable	Units	Description
${\rho}_{\infty}$	$kg/m^3$	Gas-phase density
$P$	$Pa$	Pressure
$u_g, v_g, w_g$	$m/s$	Mean gas velocity
$k$	$m^2/s^2$	Turbulent kinetic energy
$\epsilon$	$m^2/s^3$	Turbulent kinetic energy dissipation rate
$T_{\infty}$	$K$	Gas temperature
$W$	$kg/mol$	Gas-phase mean molecular weight
$Y_{F,\infty}$	N/A	Mass fraction of fuel in gas phase
$Y_{O,\infty}$	N/A	Mass fraction of oxidizer in gas phase
$c_{p,\infty}$	$J/kg/K$	Gas-phase specific heat
${\mu}_g$	$kg/m/s$	Gas-phase viscosity
$G_{in}$	$W/m^2$	Incident radiation