The equation of motion for a discrete particle p with velocity \mv{u}_p in non-dimensional form reads
\begin{align} \frac{d\mv{u}_p}{dt} &= (1 - \frac{\rho}{\rho_p} ) \frac{1}{Fr_0^2} \mv{e_g} + \frac{C_{D}}{24} \frac{Re_p}{\tau_p} (\mv{u} - \mv{u}_p), \\ \frac{d\mv{x}_p}{dt} &= \mv{u}_p \end{align}
with the Froude-number Fr_0 = \frac{\mdim{a_{0}}}{\sqrt{\mdim{g}\, \mdim{L}}}, unity gravity vector \mv{e_g}, particle Re number Re_p, and particle relaxation time \tau_p given as
\begin{align} Re_p = \frac{\rho||\mv{u} - \mv{u}_p|| d_p}{\mu}, \quad \tau_p = \frac{\rho_p d^2_p}{18\mu}. \end{align}
The drag coefficient C_D is given by a empirical relations, such as
\begin{align} C_D = \begin{cases} \frac{24}{Re'} & \text{for } Re'\leq 0.1\\ \frac{24}{Re'}(1+\frac{1}{6} (Re') ^{\frac{2}{3}}) & \text{for } Re' \leq 1000\\ 0.424 & \text{for } Re' > 1000\\ \end{cases} \end{align}
from[[putnam1961integratable]] with Re number Re' = Re_p Re_0. The additional Reynolds number Re_0 appears due to the non-dimensionalization of the fluid properties by the stagnation point reference value, i.e., \mu = \frac{\mdim{\mu}}{\mdim{\mu_0}}.
The motion of a discrete rotational symmetric ellipsoidal particle p with principal semi-minor axes a=b=c/\beta where \beta refers to the aspect ratio is characterized by its velocity \mv{u}_p and angular velocity \hat{\mv{\omega}}. Quantities in the particle fixed coordinate system are denoted by \hat{\circ}.
In the Stokes limit, the equation of motion for the translational movement of an ellipsoidal particle [[Oberbeck1876]] reads
\begin{align} \frac{d\mv{u}_p}{dt} &= (1 - \frac{\rho}{\rho_p} ) \frac{1}{Fr_0^2} \mv{e_g} + \frac{24}{9} \frac{\beta^{2/3}}{\tau_p} \mm{R^{-1}} \begin{pmatrix} \frac{1}{\chi_0/a^2 + \alpha_0} & 0 & 0 \\ 0 & \frac{1}{\chi_0/a^2 + \alpha_0} & 0 \\ 0 & 0 & \frac{1}{\chi_0/a^2 + \beta^2 \gamma_0} \end{pmatrix} \mm{R} (\mv{u} - \mv{u}_p), \\ \frac{d\mv{x}_p}{dt} &= \mv{u}_p \end{align}
with the Froude-number Fr_0 = \frac{\mdim{a_{0}}}{\sqrt{\mdim{g}\, \mdim{L}}}, unity gravity vector \mv{e_g}, rotation matrix \mm{R}, and the particle relaxation time \tau_p given as
\begin{align} \tau_p = \frac{\rho_p d^2_p}{18 \rho \mu}. \end{align}
\alpha_0, \gamma_0, and \chi_0 are shape factors that depend on the aspect ratio \beta. Analytical expressions can be found in [[Happel&Brenner1965]].
For particles with \beta=1 the equation reduces to the simplified Maxey-Riley equation for spherical particles.
The rotational motion equation depends on the hydrodynamic torque [[Jeffrey1992]] and reads
\begin{align} \begin{pmatrix} \frac{d\omega_{\hat{x}}}{dt} \\ \frac{d\omega_{\hat{y}}}{dt} \\ \frac{d\omega_{\hat{z}}}{dt} \end{pmatrix} = \begin{pmatrix} \omega_{\hat{y}}\omega_{\hat{z}}\frac{\beta^2-1}{1+\beta^2} \\ \omega_{\hat{z}}\omega_{\hat{x}}\frac{\beta^2-1}{1+\beta^2} \\ 0 \end{pmatrix} + \frac{40}{9} \frac{\beta^{2/3}}{\tau_p} \begin{pmatrix} \frac{1}{\alpha_0 + \beta^2\gamma_0} \\ \frac{1}{\alpha_0 + \beta^2\gamma_0} \\ \frac{1}{2\alpha_0} \end{pmatrix} \begin{pmatrix} \frac{1-\beta^2}{1+\beta^2} \tau_{\hat{z}\hat{y}} + (\zeta_{\hat{z}\hat{y}}-\omega_{\hat{x}}) \\ \frac{\beta^2-1}{1+\beta^2} \tau_{\hat{x}\hat{z}} + (\zeta_{\hat{x}\hat{z}}-\omega_{\hat{y}}) \\ (\zeta_{\hat{y}\hat{x}}-\omega_{\hat{z}}) \\ \end{pmatrix} \end{align}
where the fluid shear stress and vorticity are defined as
\begin{align} \tau_{ij} = \frac{1}{2} \left( \fracpart{u_i}{x_j} + \fracpart{u_j}{x_i} \right), \quad \zeta_{ij} = \frac{1}{2} \left( \fracpart{u_i}{x_j} - \fracpart{u_j}{x_i} \right) . \end{align}
For the total particle mass and temperature temporal rate of change ( \dot{m}_p , \dot{T}_p ) two differential equations following the evaporation model of Miller and Bellan~([[miller1999]]-[[Miller1998]]) are included in non-dimensional form
\begin{align} \label{parteqm} \frac{d m_p}{d t} &= -\frac{Sh'}{3\,Sc' \, Re_0} (\frac{m_p}{\tau_p})\ln(1+B_{M,neq})\\ \label{parteqt} \frac{dT_p}{dt} &=\frac{Nu'_{neq}}{3\,Pr'\,Re_0}\theta_1 \frac{T-T_p}{\tau_p} +L_{ev} (\gamma - 1) \theta_1 \frac{\dot{m_p}}{m_p} \quad , \end{align}
where \theta_1 is the ratio of fluid-to-liquid heat capacity \theta_1=\frac{\mdim{C_p}}{\mdim{C_{p,p}}} and L_{ev}=\frac{\mdim{L_{ev}}}{\mdim{a_{0}}} the non-dimensional latent heat of vaporization. The Schmidt and Prandtl number of the fluid phase are given by
\begin{align} \label{eqSchmidt} Sc' &= \frac{\mu}{\rho \Gamma} Sc_0 &\quad \text{with} \quad Sc_0 =\frac{\mdim{\mu_0}}{\mdim{\rho_0} \mdim{\Gamma_0}} \\ Pr' &= \frac{\mu C_p}{\lambda} Pr_0 &\quad \text{with} \quad Pr_0 = \frac{\mdim{\mu_0} \mdim{C_p}}{\mdim{\lambda_0}} \\ \end{align}
with the dimensionless binary diffusion coefficient \Gamma, the thermal conductivity \lambda, and the respective dimensional reference state values. For the convective heat and mass transfer, the empirical Ranz-Marshall correlations ([[RanzMarshall1]]-[[RanzMarshall2]]) for the Sherwood and Nusselt number
\begin{align} \label{eqRanzMarshall} Sh' &= 2+0.552 Re'^{\frac{1}{2}} Sc'^{\frac{1}{3}} \\ Nu' &= 2+0.552 Re'^{\frac{1}{2}} Pr'^{\frac{1}{3}} \end{align}
are used. The equilibrium vapor mole fraction is obtained from the Clausius-Clapeyron equation
\begin{align} \chi_{s,eq} = \frac{p_{B}}{p}\exp\left(\frac{\gamma L_{ev}}{\theta_2}\left(\frac{1}{T_B}-\frac{1}{T_p}\right) \right) \label{eq::clausius} \end{align}
with the liquid boiling point ( T_B at p_B) and the mole-fraction ratio \theta_2 = \frac{M}{M_p}. Non-equilibrium effects are considered by a correction of the Nusselt number with the non-dimensional evaporation parameter \beta and by a correction of the vapor mole fraction at the surface by the Langmuir-Knudsen law
\begin{align} \label{eqNeq} Nu'_{neq} &= Nu' \frac{\beta}{e^{\beta}-1} &\quad{\text{with} } \quad &\beta = 0.5 Pr' Re_b'\\ \chi_{s,neq} &= \chi_s - 2 \beta L_k, &\quad{ \text{with} } \quad &L_k = \frac{\mu \sqrt{2 \pi T_p \frac{\theta_2}{\gamma}}}{p d_p Sc' Re_0} \end{align}
with Re_b' = \frac{\dot{m}_p}{\mu \pi d_p} Re_0 as
the Reynolds number based on the blowing velocity u_b, which is obtained from mass conservation at the droplet surface with \dot{m}_p=-\pi d_p^2 u_b \rho as in[[miller1999]].
The non-equilibrium species vapor mass fraction at the particle surface Y_{s,neq} and the Spalding transfer number for mass B_{M,neq} are defined as
\begin{align} \label{eqNeq2} Y_{s,neq} = \frac{\chi_{s,neq}}{\chi_{s,neq} + (1-\chi_{s,neq})\theta_2} \, \text{and} \quad B_{M,neq} = \frac{Y_{s,neq} - Y}{1 - Y_{s,neq}}. \end{align}
The "1/3 rule" following [[Hubbard]] is used as a reference state in the boundary layer around the particle to obtain fluid properties. The Wilke mixing rule [[Wilke]] is applied for the fuel vapor and fluid mixing in the particle boundary layer.
The following source terms are used for the mass, momentum and energy exchange between the particle and the fluid phase
\begin{align} \begin{bmatrix} S_\rho \\ \mv{S_u} \\[6pt] S_E \\[6pt] S_Y \end{bmatrix} = \begin{bmatrix} \dot{m}_p \\ m_p \frac{C_{D}}{24} \frac{Re_p}{\tau_p Re_0} (\mv{u} - \mv{u}_p) + \dot{m}_p \mv{u}_p \\[6pt] \mv{S_i} \cdot \mv{u}_p - 0.5 \dot{m}_p \mv{u}_p \cdot \mv{u}_p + \frac{C_{p,p}}{\gamma -1} \left(m_p \dot{T}_p + \dot{m}_p T_p \right) \\[6pt] \dot{m}_p \end{bmatrix}. \end{align}
To reduce the computational effort generated by the spray injection, a number of fuel droplets N_d are grouped into a parcel. All droplets in a parcel are assumed to have the same physical values and suffer from the same change in their properties, i.e., they are tracked as a single Lagrangian particle.
The spray atomization and secondary breakup process, leading to decreasing particle diameter, is modeled by a combination of the Kelvin-Helmholtz (KH) and Rayleigh-Taylor (RT) model. The present implementation is based on the formulation of Reitz [[reitz1987]] and Beale and Reitz [[Reitz1999]].
The Kelvin-Helmholtz (KH) wave model assumes droplet break up due to wave growth excited by aerodynamic forces based on the relative velocity between liquid and gas phase. Continuous mass stripping from the parent droplet into smaller child droplets occurs. A correlation for the most unstable wave growth rate \Omega_{KH} and wavelength \Lambda_{KH} can be obtained from first order linear analysis such that
\begin{align} \Omega_{KH} &= \frac{0.34 + 0.38 We'^{1.5}}{(1 + Z')(1 + 1.4 T'^{0.6})} \sqrt{\frac{8 \sigma_p}{\rho_p d^3_p We_0}} \\ \Lambda_{KH} &= \frac{4.51 d_p (1 + 0.45 \sqrt{Z'})(1 + 0.4 T'^{0.7})} {(1 + 0.865 We'^{1.67})^{0.6}}, \end{align}
with the Ohnesorge number Z'={\sqrt{We_l'}}/{Re_{l}'} and the Taylor number T' = Z'\sqrt{We'}. The Weber numbers are defined as We'=\frac{0.5 \rho ||\mv{u}-\mv{u}_p||^2 d_p}{\sigma_p} We_0, the liquid Weber number is We_l'=\frac{0.5 \rho_p||\mv{u}-\mv{u}_p||_2^2d_p}{\sigma_p} We_0 and the liquid Reynolds number is Re_{l}' = 0.5 \frac{\rho_p||\mv{u} - \mv{u}_p||_2 d_p}{\mu_p} Re_0. The breakup time \tau_{KH} and droplet size d_{KH} can be determined as
\begin{align} \tau_{KH} &= \frac{1.863 B_1 d_p}{\Omega_{KH} \Lambda_{KH}} \quad \text{and} \quad d_{KH} = 2 B_0 \Lambda_{KH}. \end{align}
The model adjustable control parameters are B_0 and B_1. If a KH breakup occurs, the shedded mass of child droplets is accumulated until a mass limit m_{KH,lim} or the parent droplet mass has been reached. If the mass limit is reached, a new parcel with the mass of the shedded droplets and d_{KH} is added. The number of droplets N_d in the parent parcel remains the same, and a new parent diameter is computed based on mass conservation. If the parent mass is reached, the parent diameter is set to d_{KH} and the number of droplets is adapted such that mass conservation is achieved. See Aguerre et. al.[[Aguerre]] for further details of this implementation.
The Rayleigh-Taylor (RT) model describes the catastrophic breakup mechanism caused by deceleration effects on the liquid droplet due to aerodynamic drag. The growth rate \Omega_{RT} and wavelength \Lambda_{RT} of the RT wave are defined as
\begin{align} \Omega_{RT} &= \sqrt{\frac{2}{3 \sqrt{3\sigma_p}} \frac{(||\mv{a}_p|| (\rho_p - \rho))^{\frac{3}{2}}}{\rho_p + \rho} } We_0^{0.25}\\ \Lambda_{RT} &= 2 \pi \sqrt{\frac{3 \sigma_p} {||\mv{a}_p|| (\rho_p - \rho)}}, \end{align}
with \mv{a}_p being the droplet acceleration and \sigma_p the surface tension. The droplet size after breakup d_{RT} and the RT breakup time \tau_{RT} are given by
\begin{align} \tau_{RT} =\frac{1}{\Omega_{RT}} \quad \text{and} \quad d_{RT} = C_{RT} \Lambda_{RT} \label{eq_RT} \end{align}
with C_{RT} as the primary model control parameter. If RT breakup occurs the parcel diameter is reduced to d_{RT} and the number of droplets is increased based on mass conservation.
The interaction of the two breakup models is based on Beale and Reitz[[Reitz1999]], where it is assumed that within a certain distance L_{bu} the spray consists of a liquid core
\begin{align} L_{bu} = \frac{C_{bu}}{2} d_i \sqrt{\frac{\rho_p}{\rho_i}} \quad , \end{align}
where C_{bu} is a model parameter, d_i the injector orifice diameter, and \rho_i the initial fluid density. Within L_{bu} only KH breakup is considered, while beyond both mechanisms are applied with the RT model as the predominant mechanism.
1.9.5