3.3.3. Coupling the Lagrangian and Eulerian Fields

The previous section provides a description of the particle evolution given a gas-phase environment. In this section the means by which the gas-phase environment for a particle is determined from an Eulerian solution presumed to use a control-volume or similar approach. Following this, the effect of the particle field on the Eulerian field is described. Finally, to ensure that the coupled evolution proceeds in a physically realistic manner, it is necessary to identify limits on the time step size.

3.3.3.1. Gas-phase environment for parcels

When a parcel of particles spans more than a single control volume, it is appropriate to employ a weighted average of the gas properties in the control volumes spanned by the parcel. The average is weighted by the number of particles in a given control volume, which is obtained from $f_{\sigma}({\bf x} ; {\bf x_o},t )$ defined in Eqn. (3.633). Thus, the average value of a gas-phase variable, $\phi$ , for a parcel is

(3.693) $\bar{\phi} = \dfrac {\int_{V} {\phi ({\bf x}, t) f_{\sigma}({\bf x}; {\bf x_o}, t)} d{\bf x}} {\int_{V} {f_{\sigma}({\bf x}; {\bf x_o}, t)} d{\bf x}}$

where the integral volume may span more than one control volume. This procedure is employed for all gas phase variables that appear in the droplet evolution equations in Particle Transport Model. For transport properties (i.e. various diffusion coefficients), the gas phase properties employed are evaluated using the ‘1/3 rule’ as indicated for viscosity in Eqn. (3.625).

3.3.3.2. Effect of the particle phase on the gas phase

The source terms provided in the previous section show how the gas-phase properties affect the conservation of mass, momentum, and energy for the particles. From Newton’s third law, the action of the gas phase on the particles must be balanced by an action of the particles on the gas phase. This action is determined directly from the change in the state of the particle phase as described in the following.

The source of the mass for a given control volume, $V_c$ , located at $({\bf x}, t)$ is determined by summing the changes in the masses of all the particles in that control volume over the gas-phase time step $\delta t$

(3.694) $S_{mass}({\bf x},t) = - \dfrac {1} {V_c} {\sum}_P \left[\dfrac {m_p (t+ \delta t) - m_p(t)} {\delta t} \int_{V_c} {f_{\sigma}({\bf x}; {\bf x_o}, t)} d{\bf x} \right].$

Here, the summation is presumed to occur over the $P$ parcels that contribute to the control volume at $({\bf x}, t)$ , and the addend in the square brackets corresponds to the mass change and particle number distribution in each of the $P$ parcels. The similar equation for species concentration incorporates the conversion of the mass of the particle component to the mass of the gas-phase component. For pure evaporation, the conversion is trivial, but for combustion where the fuel evaporation corresponds to oxidizer consumption and the product formation, the expressions become complicated. The general expression for the species mass source term is

(3.695) $S_{mass,i}({\bf x}, t) = - \dfrac {1} {V_c} {\sum}_P \left[\dfrac {{\nu}_i W_i} {W_F} \dfrac {m_p (t+ \delta t) - m_p(t)} {\delta t} \int_{V_c} {f_{\sigma}({\bf x}; {\bf x_o}, t)} d{\bf x} \right].$

Here, ${\nu}_i$ is the number of moles of species $i$ produced when a mole of the particle component evaporates, and $W_i$ is the molecular weight of species $i$ . The particle component is presumed to be the fuel, denoted with the subscript $F$ , and for pure evaporation reduces to equation (3.694).

Similarly, the source term for the $j$ -momentum equations is determined by summing the changes in the particle momentum over all of the particles in a control volume.

(3.696) $S_{mom, j}({\bf x},t) = - \dfrac {1} {V_c} {\sum}_P \left[ \left(\dfrac {m_p(t+\delta t)u_{p,j}(t+\delta t) - m_p(t)u_{p,j}(t) } {\delta t} - m_p(t) g_i \right) \int_{V_c} {f_{\sigma}({\bf x}; {\bf x_o}, t)} d{\bf x} \right]$

Note that the last term in Eqn. (3.696) describes the effect of gravitational acceleration on the particle field.

The energy source term must account for the enthalpy of vaporization and heat of combustion associated with the particles as well as the particle heating. It must also account for the change in the enthalpy of the gas associated with any gases that evaporated or condensed. In addition, it is necessary to separate out the contributions associated with radiative transport since these do not necessarily affect the local control volume, but are transported over length scales determined by the radiative transport equation. The enthalpy source term is

(3.697) $S_{enthalpy}({\bf x},t) = -S_{rad} - \dfrac {1} {V_c} {\sum}_P \left[ \dfrac {m_p(t+\delta t) \left[h_p(t+\delta t) - h_p(t) \right]} {\delta t} \int_{V_c} {f_{\sigma}({\bf x}; {\bf x_o}, t)} d{\bf x} \right] - \dfrac {1} {V_c} {\sum}_P \left[ \dfrac {\left[m_p(t+\delta t) - m_p(t)\right] {\sum}_i {\nu}_i W_i h_i (T_g)/W_F} {\delta t} \int_{V_c} {f_{\sigma}({\bf x}; {\bf x_o}, t)} d{\bf x} \right] - \dfrac {1} {V_c} {\sum}_P \left[ \dfrac {\left[m_p(t+\delta t) - m_p(t)\right] (q_{comb}-h_{vap})} {\delta t} \int_{V_c} {f_{\sigma}({\bf x}; {\bf x_o}, t)} d{\bf x} \right]$

where the radiative transport source term is

(3.698) $S_{rad}({\bf x},t) = \dfrac {1} {V_c} {\sum}_P \left[Q_{rad} \int_{V_c} {f_{\sigma}({\bf x}; {\bf x_o}, t)} d{\bf x} \right],$

which provides the interface between the particle field and the radiation transport equation. In general, the solution of the radiant transport equation requires an absorption coefficient and a contribution to the emission. These must be defined in such a way that they agree with the radiation absorbed and emitted by the particle field. To do this, an energy balance for the particle field, but in an Eulerian frame, is employed. Within the energy conservation equation, the radiant source term appears as a divergence of the radiation heat flux $\nabla \cdot q_R$ . This divergence of the radiation heat flux can be related to the radiation intensity and the radiation emitted by [119]

(3.699) $\nabla \cdot q_R = \langle 4 {\kappa}_{net} \sigma T_{net}^4 \rangle - G_{in} \langle {\kappa}_{net} \rangle$

where the subscript $net$ indicates that contributions from the particles must be combined with those from gases, smoke, and anything else that participates with the radiative field. The angular brackets indicate that the quantities on the right hand side must be defined based on the appropriate summation or averaging process over all of these participating media, which may have varying temperatures. In order to balance the radiant energy flux in and out of the particle with their contribution to Eqn. (3.699), it is necessary to separate these particle contributions from the net radiant source term. To leading order, this can be done in an additive manner so that the first term in Eqn. (3.699) is split as $\langle 4{\kappa}_{net} \sigma T_{net}^4 \rangle = \langle 4 {\kappa}_p \sigma T_p^4 \rangle + \langle 4 {\kappa}_{others} \sigma T_{others}^4 \rangle$ and the second term is split with $\langle{\kappa}_{net}\rangle = \langle{\kappa}_p\rangle + \langle{\kappa}_{others}\rangle$ . Here the subscript $p$ indicates the contribution associated with the particle field and the subscript $others$ indicates all other contributions (gases, soot, etc.). Corrections to this leading order approximation are related to the optical thickness of the control volumes over which the radiation solve occurs. These corrections arise because a portion of the intensity is absorbed within the control volume. As long as the control volumes are all optically thin, then this correction is not important, but it must be accounted for in the radiation solve if that term is important.

Separating out only the contribution of the particle field to the radiation solve gives

(3.700) $\left(\nabla \cdot q_R \right)_p = \langle 4{\kappa}_p \sigma T_p^4 \rangle - G_{in} \langle{\kappa}_p \rangle$

where the angular brackets indicate that the quantities on the right hand side must be defined based on the appropriate summation over all of the particles in the control volume. In order to ensure energy conservation between the Lagrangian solution of the particle field and the eulerian solution of the radiation transfer and energy equations, these appropriate sums are obtained by integrating over the particle field in a given control volume as in Eqn. (3.698). Here, the emission and absorption contributions of $Q_{rad}$ are separated

(3.701) $S_{rad}({\bf x}, t) = \dfrac {1} {V_c} {\sum}_P \left[ 4 \pi r_p^2 \alpha \sigma T_p^4 \int_{V_c} {f_{\sigma}({\bf x}; {\bf x_o}, t)} d{\bf x} \right] - \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] G_{in}$

so that, equating the first and second terms on the right-hand side of Equations (3.700) and (3.701) gives

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

and

(3.703) $\langle {\kappa}_p \rangle = \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] G_{in}$

This information is crucial because it defines an energy conservation interface between the particle field and the radiation field.

The source terms defined in Equations (3.694) through Effect of the particle phase on the gas phase should be added directly to the gas-phase conservation equations in the manner appropriate for the chosen numerical method. In the above relations, it was presumed that the dimensions of the conservation equations solved are those of density, density time velocity, and density time enthalpy. The radiant flux source term in Eqn. (3.698) should similarly be added to the radiative transport equation, but in the consistent manner indicated in the above paragraph.

3.3.3.3. Time step control

The Lagrangian and Eulerian fields are advanced using an explicit operator splitting approach. The particles are presumed to be advanced using an ordinary differential equation solver capable of handling a stiff system. The system can be stiff because the magnitude of the forcing function, the right hand side, of the various particle evolution equations can vary over orders of magnitude. While the particle state is evolved, the gas-phase state is presumed to be constant, except as described in the following paragraph. The source terms indicated in the previous section are accumulated during the particle evolution, and then they are applied to the gas-phase conservation equations as it is advanced.

There are certain requirements that limit the particle evolution time step. Particle should not move more than the length of the control volume side without re-evaluating the gas-phase state using the methods described in Gas-phase environment for parcels. Particles should not evolve for longer than the eddy interaction time defined as the lesser of times defined in Eqn. (3.631) and (3.632) without determining a new value for $u_{g,i}$ in Eqn. (3.629) and (3.630). These requirements do not necessarily limit the gas-phase evolution time step.

There are also limitations to the gas-phase evolution time step imposed by the particle evolution. The source of this limitation is the requirement that the treatment of the gas-phase state as constant during the particle evolution not lead to nonphysical or inaccurate results. For example, if the particles transferred almost all of their momentum to the gas phase in a single time step (because they were small, for example), then a problem could arise. Specifically, if the subsequent momentum source term were large, then the gas could be accelerated to velocities that exceed the initial particle velocities. In the subsequent time step, the particles would accelerate and the solution procedure could thereby destabilize. There are two means of avoiding such problems. One approach is to employ an iterative implicit coupled advance of the Eulerian and Lagrangian state. Such an approach would depend on the specific algorithms employed for the Eulerian solver. The approach that will be discussed here is a limitation on the time step size for the coupled Eulerian-Lagrangian system. This limitation is based on the idea that the change in the Eulerian state relative to the Lagrangian state should not be too dramatic. This not only provides stability, but helps limit numerical errors.

For mass conservation, the limitation on the time step is based on the idea that the mass added to the cell should not dramatically affect the pressure (or whatever variable changes to allow additional mass to entry the control volume gas phase). With $S_{mass}$ in units of density per unit time, the appropriate limitation on the time step is a time step that leads to no more than a change of ${\delta}_{\rho} \rho$ in $\rho$ .

(3.704) $dt_{max} < {\delta}_{\rho} \rho / S_{mass}$

That is, the fractional change in the density is limited to ${\delta}_{\rho}$ . It is expected that ${\delta}_{\rho}$ is substantially smaller than 0.1 is appropriate, but this is subject to evaluation and will depend somewhat on the details of a given simulation. for the conservation of species, similar expressions apply, but with the second term that provides a sort of absolute tolerance in addition to the relative tolerance indicated above

(3.705) $dt_{max} < {\delta}_{Y_i} \rho Y_i / S_{mass,i} + {\eta}_{Y_i} \rho / S_{mass,i}.$

Intuition suggests that the limitation of ${\delta}_{Y_i}$ do not need to be as severe as those on ${\delta}_{\rho}$ , for stability and that ${\eta}_{Y_i}$ should be substantially smaller than typical magnitudes for $Y_i$ for accuracy. Again, the specific values will depend on how sensitive the system is on $Y_i$ . For momentum, a similar expression is provided with relative and absolute tolerances

(3.706) $dt_{max} < {\delta}_u \rho u_i / S_{mom, i} + {\eta}_u \rho / S_{mom, i}$

where ${\delta}_u u_i$ and ${\eta}_u$ are indicative of the acceptable uncertainties in the velocity field. Time-step control based on the enthalpy exchange is easiest to think of in terms of temperature

(3.707) $dt_{max} < {\delta}_h \rho c_{p,g} T_g / S_{enthalpy} + {\eta}_h \rho c_{p,g}/ S_{enthalpy}$

where ${\delta}_hT_g$ and ${\eta}_h$ are indicative of the acceptable uncertainties in the temperature field footnote{ Instead of using $S_{enthalpy}$ that includes both the chemical and sensible enthalpy changes, it is better to use a measure of the change in the sensible enthalpy source term (see Eqn Effect of the particle phase on the gas phase) where the next to the last set of brackets includes only the sensible contribution of the species heating and not the chemical enthalpy contribution.

3.3.3.4. Particle-surface interactions

Particles interacting with a surface can (1) bounce off of the surface and return to the flow with a different trajectory, (2) stick to the surface and remain as a deposit or (3) shatter so that smaller droplets are formed that leave the surface with various trajectories. The appropriate behavior depends primarily on the ratio of the droplet kinetic energy to the surface energy as determined by the Weber number

(3.708) $We = \dfrac {{\rho}_p d_p {|{\bf u}_{p,i}|^2}} {\sigma}$

where $\sigma$ is the surface tension of the condensed phase. A complete model description is available in [124] with criteria for droplet sticking and bouncing. That paper did not address droplet shattering in detail and the droplet shattering model is provided in the following paragraph.

Droplets will shatter if the criteria

(3.709) $We ^ {0.5} Re_p ^ {0.25} > K_{crit} = 57.7$

is satisfied, which occurs for relatively large droplets traveling at relatively high velocities. Satellite droplets are presumed to form with uniform sizes given by

(3.710) $d_{p,shatter} = max \left(\dfrac {7.9 \cdot 10^{10} \sigma We^{1.4}/Re^{2.8}} {{\rho}_p {|{\bf u}_{p,i}|^2}}, d_p \right)$

where the second term in the max function ensures that the empirical relationship given as the first term does not exceed the original diameter.

If a particle sticks to the surface, it is necessary to track the mass and energy addition to that surface through a field variable added to the solid object surface set. Mass addition should be tracked on a mass per solid-element surface area basis. Similarly, the energy deposited on a surface should be tracked based on $c_{v,p} (T_p - T_{surface})$ . Erikson and Gill have developed and implemented a model to provide such an interface for a Calore surface [125].