Equations
Spray Equations
This section outlines the modeling and mathematics for the spray routines in PeleMP. Firstly, spray modeling in PeleMP relies on the following assumptions:
Dilute spray: droplet volume inside an Eulerian cell is much smaller than the volume of the gas phase; the droplets can be modeled as Lagrangian point source terms relative to the Eulerian gas phase
Infinite conductivity model: temperature within a droplet is temporally varying but spatially uniform
One-third rule: the thermophysical properties in the film of an evaporating droplet can be approximated as a weighted average of the state at the droplet surface (weighted as 2/3) and the state of the surrounding gas (weighted as 1/3)
Ideal equilibrium: the liquid and vapor state at the surface of the droplet are in equilibrium
The radiation, Soret, and Dufour effects are neglected
The evaporation models follow the work by Abramzon and Sirignano [1] and the multicomponent evaporation is based on work by Tonini. [2] Details regarding the energy balance are provided in Ge et al. [3]
The subscript notation for this section is: \(d\) relates to the liquid droplet, \(v\) relates to the vapor state that is in equilibrium with the liquid and gas phase, \(L\) relates to the liquid phase, and \(g\) relates to the gas phase. The subscript \(r\) relates to the reference state with which to approximate the thermophysical and transport properties. This reference state is assumed to be in the evaporating film that surrounds the droplet state and is approximated as
where \(A = 1/3\) according the the one-third rule. Additional nomenclature: \(M_n\) is the molar mass of species \(n\), \(\overline{M}\) is the average molar mass of a mixture, \(\mathcal{R}\) is the universal gas constant, \(N_L\) is the number of liquid species, and \(N_s\) is the number of gas phase species. \(Y_n\) and \(\chi_n\) are the mass and molar fractions of species \(n\), respectively. The user is required to provide a reference temperature for the liquid properties, \(T^*\), the critical temperature for each liquid species, \(T_{c,n}\), the boiling temperature for each liquid species at atmospheric pressure, \(T^*_{b,n}\), the latent heat and liquid specific heat at the reference temperature, \(h_{L,n}(T^*)\) and \(c_{p,L,n}(T^*)\), respectively. Note: this reference temperature is a constant value for all species and is not related to the reference state denoted by the subscript \(r\).
The equations of motion, mass, momentum, and energy for the Lagrangian spray droplet are:
where \(\mathbf{X}_d\) is the spatial vector, \(\mathbf{u}_d\) is the velocity vector, \(T_d\) is the droplet temperature, \(m_d\) is the mass of the droplet, \(\mathbf{g}\) is an external body force (like gravity), \(\dot{m}\) is evaporated mass, \(\mathcal{Q}_d\) is the heat transfer between the droplet and the surrounding gas, and \(\mathbf{F}_d\) is the momentum source term. The density of the liquid mixture, \(\rho_d\), depends on the liquid mass fractions of the dropet, \(Y_{d,n}\),
The droplets are assumed to be spherical with diameter \(d_d\). Therefore, the mass is computed as
The procedure is as follows for updating the spray droplet:
Interpolate the gas phase state to the droplet location using a trilinear interpolation scheme.
Compute the boiling temperature for species \(n\) at the current gas phase pressure using the Clasius-Clapeyron relation
\[T_{b,n} = \left(\log\left(\frac{p_{\rm{atm}}}{p_g}\right) \frac{\mathcal{R}}{M_n h_{L,n}(T^*_{b,n})} + \frac{1}{T^*_{b,n}}\right)\]The boiling temperature of the droplet is computed as
\[T_{d,b} = \sum^{N_L}_{n=0} Y_{d,n} T_{b,n}\]Since we only have the latent heat at the reference condition temperature, we estimate the enthalpy at the boiling condition using Watson’s law
\[h_{L,n}(T^*_{b,n}) = h_{L,n}(T^*) \left(\frac{T_{c,n} - T^*}{T_{c,n} - T^*_{b,n}} \right)^{-0.38}\]Compute the latent heat of the droplet using
\[h_{L,n}(T_d) = h_{g,n}(T_d) - h_{g,n}(T^*) + h_{L,n}(T^*) - c_{p,L,n}(T^*) (T_d - T^*) \,.\]and the saturation pressure using either the Clasius-Clapeyron relation
\[p_{{\rm{sat}}, n} = p_{\rm{atm}} \exp\left(\frac{h_{L,n}(T_d) M_n}{\mathcal{R}} \left(\frac{1}{T^*_{b,n}} - \frac{1}{T_d}\right)\right)\]or the Antoine curve fit
\[p_{{\rm{sat}},n} = d 10^{a - b / (T_d + c)}\]Estimate the mass fractions in the vapor state using Raoult’s law
\[ \begin{align}\begin{aligned}Y_{v,n} &= \frac{\chi_{v,n} M_n}{\overline{M}_v + \overline{M}_g (1 - \chi_{v,{\rm{sum}}})} \; \forall n \in N_L\\\chi_{v,{\rm{sum}}} &= \sum^{N_L}_{n=0} \chi_{v,n}\\\chi_{v,n} &= \frac{\chi_{d,n} p_{{\rm{sat}},n}}{p_g}\\\chi_{d,n} &= \frac{Y_{d,n}}{M_n}\left(\sum^{N_L}_{k=0} \frac{Y_{d,k}}{M_k}\right)^{-1}\\\overline{M}_v &= \sum^{N_L}_{n=0} \chi_{v,n} M_n\end{aligned}\end{align} \]If \(\chi_{g,n} p_g > p_{{\rm{sat}},n}\), then \(\chi_{v,n} = Y_{v,n} = 0\) for that particular species in the equations above, since that means the gas phase is saturated. The mass fractions in the reference state for the fuel are computed using the one-third rule and the remaining reference mass fractions are normalized gas phase mass fractions to ensure they sum to 1
\[\begin{split}Y_{r,n} = \left\{\begin{array}{c l} \displaystyle Y_{v,n} + A (Y_{g,n} - Y_{v,n}) & {\text{If $Y_{v,n} > 0$}}, \\ \displaystyle\frac{1 - \sum^{N_L}_{k=0} Y_{v,k}}{1 - \sum^{N_L}_{k=0} Y_{g,k}} Y_{g,n} & {\text{Otherwise}}. \end{array}\right. \; \forall n \in N_s.\end{split}\]The average molar mass, specific heat, and density of the reference state in the gas film are computed as
\[ \begin{align}\begin{aligned}\overline{M}_r &= \left(\sum^{N_s}_{n=0} \frac{Y_{r,n}}{M_n}\right)^{-1},\\c_{p,r} &= \sum^{N_s}_{n=0} Y_{s,n} c_{p,g,n}(T_r),\\\rho_r &= \frac{\overline{M}_r p_g}{\mathcal{R} T_r}.\end{aligned}\end{align} \]Transport properties are computed using the reference state: dynamic viscosity, \(\mu_r\), thermal conductivity, \(\lambda_r\), and mass diffusion coefficient for species \(n\), \(D_{r,n}\).
It is important to note that PelePhysics provides mixture averaged mass diffusion coefficient \(\overline{(\rho D)}_{r,n}\), which is converted into the binary mass diffusion coefficient using
\[(\rho D)_{r,n} = \overline{(\rho D)}_{r,n} \overline{M}_r / M_n.\]Mass diffusion coefficient is then normalized by the total fuel vapor molar fraction
\[(\rho D)^*_{r,n} = \frac{\chi_{v,n} (\rho D)_{r,n}}{\chi_{v,{\rm{sum}}}} \; \forall n \in N_L\]and the total is
\[(\rho D)_r = \sum_{n=0}^{N_L} (\rho D)_{r,n}^*\]The momentum source is a function of the drag force
\[\mathbf{F}_d = \frac{1}{2} \rho_r C_D A_d \left\|\Delta \mathbf{u}\right\| \Delta \mathbf{u}\]where \(\Delta \mathbf{u} = \mathbf{u}_g - \mathbf{u}_d\), \(A_d = \pi d_d^2/4\) is the frontal area of the droplet, and \(C_D\) is the drag coefficient for a sphere, which is estimated using the standard drag curve for an immersed sphere
\[\begin{split}C_D = \frac{24}{{\rm{Re}}_d}\left\{\begin{array}{c l} 1 & {\text{If Re$_d$ < 1}}, \\ \displaystyle 1 + \frac{{\rm{Re}}^{2/3}_d}{6} & {\text{Otherwise}}. \end{array}\right.\end{split}\]The droplet Reynolds number is defined as
\[{\rm{Re}}_d = \frac{\rho_r d_d \left\|\Delta \mathbf{u}\right\|}{\mu_r}\]The mass source term is modeled according to Abramzon and Sirignano (1989). The following non-dimensional numbers and factors are used:
\[ \begin{align}\begin{aligned}F(B) &= (1 + B)^{0.7}\frac{\log(1 + B)}{B}\\F_2 &= \max(1, \min(400, {\rm{Re}}_d)^{0.077})\\{\rm{Pr}}_r &= \frac{\mu_r c_{p,r}}{\lambda_r}\\{\rm{Sc}}_r &= \frac{\mu_r}{(\rho D)_r}\\{\rm{Sh}}_0 &= 1 + (1 + {\rm{Re}}_d {\rm{Sc}}_r)^{1/3} F_2\\{\rm{Nu}}_0 &= 1 + (1 + {\rm{Re}}_d {\rm{Pr}}_r)^{1/3} F_2\\{\rm{Sh}}^* &= 2 + \frac{{\rm{Sh}}_0 - 2}{F(B_M)}\\{\rm{Nu}}^* &= 2 + \frac{{\rm{Nu}}_0 - 2}{F(B_T)}\end{aligned}\end{align} \]The Spalding numbers for mass transfer, \(B_M\), and heat transfer, \(B_T\), are computed using
\[ \begin{align}\begin{aligned}B_M &= \displaystyle\frac{\sum^{N_L}_{n=0} Y_{v,n} - \sum^{N_L}_{n=0} Y_{g,n}}{1 - \sum^{N_L}_{n=0} Y_{v,n}}\\B_T &= \left(1 + B_M\right)^{\phi} - 1\end{aligned}\end{align} \]where
\[\phi = \frac{c_{p,r} (\rho D)_r {\rm{Sh}}^*}{\lambda_r {\rm{Nu}}^*}\]Note the dependence of \({\rm{Nu}}^*\) on \(B_T\) means an iterative scheme is required to solve for both. The droplet vaporization rate and heat transfer become
\[ \begin{align}\begin{aligned}\dot{m}_n &= -\pi (\rho D)_{r,n}^* d_d {\rm{Sh}}^* \log(1 + B_M). \; \forall n \in N_L\\\mathcal{Q}_d &= \pi \lambda_r d_d (T_g - T_d) {\rm{Nu}}^* \frac{\log(1 + B_T)}{B_T}\end{aligned}\end{align} \]If the gas phase is saturated for all liquid species, the equations for heat and mass transfer become
\[ \begin{align}\begin{aligned}\dot{m}_n &= 0\\\mathcal{Q}_d &= \pi \lambda_r d_d (T_g - T_d) {\rm{Nu}}_0\end{aligned}\end{align} \]
To alleviate conservation issues at AMR interfaces, each parcel only contributes to the gas phase source term of the cell containing it. The gas phase source terms for a single parcel to the cell are
\[ \begin{align}\begin{aligned}S_{\rho} &= \mathcal{C} \sum^{N_L}_{n=0} \dot{m}_n,\\S_{\rho Y_n} &= \mathcal{C} \dot{m}_n,\\\mathbf{S}_{\rho \mathbf{u}} &= \mathcal{C} \mathbf{F}_d,\\S_{\rho h} &= \mathcal{C}\left(\mathcal{Q}_d + \sum_{n=0}^{N_L} \dot{m}_n h_{g,n}(T_d)\right),\\S_{\rho E} &= S_{\rho h} + \frac{1}{2}\left\|\mathbf{u}_d\right\| S_{\rho} + \mathcal{C} \mathbf{F}_d \cdot \mathbf{u}_d\end{aligned}\end{align} \]where
\[\mathcal{C} = -\frac{N_{d}}{V_{\rm{cell}}},\]\(N_{d}\) is the number of droplets per computational parcel, and \(V_{\rm{cell}}\) is the volume for the cell of interest. Note that the cell volume can vary depending on AMR level and if an EB is present.
Soot Equations
In PeleMP, soot formation, growth, and oxidation is modeled using the hybrid-method of moments (HMOM) model developed by Mueller et al. [4]. This approach combines the numerical ease of the method of moments with interpolative closure (MOMIC) with the ability to capture the bimodal nature of the soot number density function (NDF) provided by the direct quadrature method of moments (DQMOM). \(M_{x,y}\) is the moment of the soot NDF, where \(x\) is the order for volume and \(y\) for surface area; these terms are modeled according to
where \(R\) is the order of the polynomial interpolation, \(a_{r,k}\) are the interpolation coefficients, \(N_0\) is the weight of the delta function, and \(V_0\) and \(S_0\) are the volume and surface area of the nucleated spherical soot particles, respectively. The location of the delta function is fixed at coordinates \(V_0\) and \(S_0\) and assumed to be equal to the nucleated particle size. The nucleated particle volume and surface area are fixed according to \(V_0 = 2 W_C C_{\rm{dimer}} / \rho_{\rm{soot}}\) and \(S_0 = (36 \pi)^{1/3} V_0^{2/3}\), where \(W_C\) is the molar mass of carbon, \(C_{\rm{dimer}}\) is the average number of carbon atoms per dimer, and \(\rho_{\rm{soot}}\) is the density of soot (\(\rho_{\rm{soot}} = 1800 {\text{ kg/m}}^3\)). If the first-order polynomial interpolation of the moments (\(R=1\)) is used, the above equation reduces to
where \(V_L\) and \(S_L\) are the mean volume and surface area of the second mode (large particles). For first-order polynomial interpolation, four transport equations are solved: \(M_{0,0}\) (number density), \(M_{1,0}\) (volume fraction, also denoted as \(f_v\)), \(M_{0,1}\), and \(N_0\). The governing equations for the soot moments are [5]
where \(\boldsymbol{J}_{M}\) is the soot mass flux and \(\dot{M}_{x,y}\) is the soot source term. The current formulation ignores the molecular diffusion and thermophoretic effects of the soot, making the first term of the right-hand side zero.
For more details regarding the HMOM model, users are encouraged to consult the references cited.