Motion Equation (Spherical Particles)

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}}\).

Motion Equation (Ellipsoidal Particles)

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}

Energy Equation

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}

Spray Modeling

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]].

Kelvin Helmholtz Model

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.

Rayleigh Taylor Model

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.

References

Empirical Drag law correlations