[Fluent Inc. Logo] return to home search
next up previous contents index

22.7.2 Droplet Breakup Models

FLUENT offers two droplet breakup models: the Taylor analogy breakup (TAB) model and the wave model. The TAB model is recommended for low-Weber-number injections and is well suited for low-speed sprays into a standard atmosphere. For Weber numbers greater than 100, the wave model is more applicable. The wave model is popular for use in high-speed fuel-injection applications. Details for each model are provided below.

Taylor Analogy Breakup (TAB) Model


The Taylor analogy breakup (TAB) model is a classic method for calculating droplet breakup, which is applicable to many engineering sprays. This method is based upon Taylor's analogy [ 369] between an oscillating and distorting droplet and a spring mass system. Table  22.7.1 illustrates the analogous components.

Table 22.7.1: Comparison of a Spring-Mass System to a Distorting Droplet
Spring-Mass System Distorting and Oscillating Droplet
restoring force of spring surface tension forces
external force droplet drag force
damping force droplet viscosity forces

The resulting TAB model equation set, which governs the oscillating and distorting droplet, can be solved to determine the droplet oscillation and distortion at any given time. As described in detail below, when the droplet oscillations grow to a critical value the "parent'' droplet will break up into a number of smaller "child'' droplets. As a droplet is distorted from a spherical shape, the drag coefficient changes. A drag model that incorporates the distorting droplet effects is available in FLUENT. See Section  22.6 for details.

Use and Limitations

The TAB model is best for low-Weber-number sprays. Extremely high-Weber-number sprays result in shattering of droplets, which is not described well by the spring-mass analogy.

Droplet Distortion

The equation governing a damped, forced oscillator is [ 270]

 F - kx - d \frac{dx}{dt} = m \frac{d^2 x}{dt^2} (22.7-8)

where $x$ is the displacement of the droplet equator from its spherical (undisturbed) position. The coefficients of this equation are taken from Taylor's analogy:

$\displaystyle \frac{F}{m}$ $\textstyle =$ $\displaystyle C_F \frac{\rho_g u^2}{\rho_l r}$ (22.7-9)
$\displaystyle \frac{k}{m}$ $\textstyle =$ $\displaystyle C_k \frac{\sigma}{\rho_l r^3}$ (22.7-10)
$\displaystyle \frac{d}{m}$ $\textstyle =$ $\displaystyle C_d \frac{\mu_l}{\rho_l r^2}$ (22.7-11)

where $\rho_l$ and $\rho_g$ are the discrete phase and continuous phase densities, $u$ is the relative velocity of the droplet, $r$ is the undisturbed droplet radius, $\sigma$ is the droplet surface tension, and $\mu_l$ is the droplet viscosity. The dimensionless constants $C_F$, $C_k$, and $C_d$ will be defined later.

The droplet is assumed to break up if the distortion grows to a critical ratio of the droplet radius. This breakup requirement is given as

 x > C_b r (22.7-12)

where $C_b$ is a constant equal to 0.5 if breakup is assumed to occur when the distortion is equal to the droplet radius, i.e., the north and south poles of the droplet meet at the droplet center. This implicitly assumes that the droplet is undergoing only one (fundamental) oscillation mode. Equation  22.7-8 is nondimensionalized by setting $y=x/(C_b r)$ and substituting the relationships in Equations  22.7-9- 22.7-11:

 \frac{d^2 y}{dt^2} = \frac{C_F}{C_b} \frac{\rho_g}{\rho_l} \... ...ma}{\rho_l r^3} y - \frac{C_d \mu_l}{\rho_l r^2} \frac{dy}{dt} (22.7-13)

where breakup now occurs for $y>1$. For under-damped droplets, the equation governing $y$ can easily be determined from Equation  22.7-13 if the relative velocity is assumed to be constant:

 y(t) = {\rm We}_c + e^{-(t/t_d)}\left[(y_0-{\rm We}_c) \cos ... ...} + \frac{y_0 -{\rm We}_c}{t_d}\right) \sin (\omega t) \right] (22.7-14)


$\displaystyle {\rm We}$ $\textstyle =$ $\displaystyle \frac{\rho_g u^2 r}{\sigma}$ (22.7-15)
$\displaystyle {\rm We}_c$ $\textstyle =$ $\displaystyle \frac{C_F}{C_k C_b} {\rm We}$ (22.7-16)
$\displaystyle y_0$ $\textstyle =$ $\displaystyle y(0)$ (22.7-17)
$\displaystyle \frac{dy_0}{dt}$ $\textstyle =$ $\displaystyle \frac{dy}{dt}(0)$ (22.7-18)
$\displaystyle \frac{1}{t_d}$ $\textstyle =$ $\displaystyle \frac{C_d}{2} \frac{\mu_l}{\rho_l r^2}$ (22.7-19)
$\displaystyle \omega^2$ $\textstyle =$ $\displaystyle C_k \frac{\sigma}{\rho_l r^3} - \frac{1}{t_d^2}$ (22.7-20)

In Equation  22.7-14, $u$ is the relative velocity between the droplet and the gas phase and We is the droplet Weber number, a dimensionless parameter defined as the ratio of aerodynamic forces to surface tension forces. The droplet oscillation frequency is represented by $\omega$. The default value of $y_0$ is 0, based upon the work of Liu et al. [ 218]. The constants have been chosen to match experiments and theory [ 190]:

\begin{eqnarray*} C_k & = & 8 \\ C_d & = & 5 \\ C_F & = & \frac{1}{3} \\ \end{eqnarray*}

If Equation  22.7-14 is solved for all droplets, those with $y>1$ are assumed to break up. The size and velocity of the new child droplets must be determined.

Size of Child Droplets

The size of the child droplets is determined by equating the energy of the parent droplet to the combined energy of the child droplets. The energy of the parent droplet is [ 270]

 E_{\rm parent} = 4 \pi r^2 \sigma + K \frac{\pi}{5} \rho_l r^5 \left[\left(\frac{dy}{dt}\right)^2 + \omega^2 y^2\right] (22.7-21)

where $K$ is the ratio of the total energy in distortion and oscillation to the energy in the fundamental mode, of the order ( $\frac{10}{3}$). The child droplets are assumed to be nondistorted and nonoscillating. Thus, the energy of the child droplets can be shown to be

 E_{\rm child} = 4 \pi r^2 \sigma \frac{r}{r_{32}} + \frac{\pi}{6}\rho_l r^5 \left(\frac{dy}{dt}\right)^2 (22.7-22)

where $r_{32}$ is the Sauter mean radius of the droplet size distribution. $r_{32}$ can be found by equating the energy of the parent and child droplets (i.e., Equations  22.7-21 and 22.7-22), setting $y=1$, and $\omega^2 = 8\sigma/\rho_l r^3$:

 r_{32} = \frac{r}{1+\frac{8Ky^2}{20} + \frac{\rho_l r^3 (dy/dt)^2}{\sigma} \left(\frac{6K-5}{120}\right)} (22.7-23)

Once the size of the child droplets is determined, the number of child droplets can easily be determined by mass conservation.

Velocity of Child Droplets

The TAB model allows for a velocity component normal to the parent droplet velocity to be imposed upon the child droplets. When breakup occurs, the equator of the parent droplet is traveling at a velocity of $dx/dt = C_b r (dy/dt)$. Therefore, the child droplets will have a velocity normal to the parent droplet velocity given by

 v_{\rm normal} = C_v C_b r \frac{dy}{dt} (22.7-24)

where $C_v$ is a constant of order (1).

Droplet Breakup

To model droplet breakup, the TAB model first determines the amplitude for an undamped oscillation $(t_d \approx \infty$) for each droplet at time step $n$ using the following:

 A = \sqrt{(y^n - {\rm We}_c)^2 + \left(\frac{(dy/dt)^n}{\omega}\right)^2} (22.7-25)

According to Equation  22.7-25, breakup is possible only if the following condition is satisfied:

 {\rm We}_c + A > 1 (22.7-26)

This is the limiting case, as damping will only reduce the chance of breakup. If a droplet fails the above criterion, breakup does not occur. The only additional calculations required then, are to update $y$ using a discretized form of Equation  22.7-14 and its derivative, which are both based on work done by O'Rourke and Amsden [ 270]:

 y^{n+1} = {\rm We}_c + e^{-(\Delta t/t_d)} \left\{(y^n - {\r... ... + \frac{y^n-{\rm We}_c}{t_d} \right] \sin (\omega t) \right\} (22.7-27)

\left(\frac{dy}{dt}\right)^{n+1} = \frac{{\rm We}_c - y^{n+1}... ...[\left(\frac{dy}{dt}\right)^n + \frac{y^n-{We}_c}{t_d} \right]}

 \omega e^{-(\Delta t/t_d)} \left\{\frac{1}{\omega} \left[\le... ...\Delta t) - (y^n - {\rm We}_c) \sin (\omega \Delta t) \right\} (22.7-28)

All of the constants in these expressions are assumed to be constant throughout the time step.

If the criterion of Equation  22.7-26 is met, then breakup is possible. The breakup time, $t_{\rm bu}$, must be determined to see if breakup occurs within the time step $\Delta t$. The value of $t_{\rm bu}$ is set to the time required for oscillations to grow sufficiently large that the magnitude of the droplet distortion, $y$, is equal to unity. The breakup time is determined under the assumption that the droplet oscillation is undamped for its first period. The breakup time is therefore the smallest root greater than $t^n$ of an undamped version of Equation  22.7-14:

 {\rm We}_c + A\cos [\omega(t-t^n) + \phi] = 1 (22.7-29)


 \cos \phi = \frac{y^n - {\rm We}_c}{A} (22.7-30)


 \sin \phi = -\frac{(dy/dt)^n}{A \omega} (22.7-31)

If $t_{\rm bu} > t^{n+1}$ , then breakup will not occur during the current time step, and $y$ and $(dy/dt)$ are updated by Equations  22.7-27 and 22.7-28. The breakup calculation then continues with the next droplet. Conversely, if $t^n < t_{\rm bu} < t^{n+1}$, then breakup will occur and the child droplet radii are determined by Equation  22.7-23. The number of child droplets, $N$, is determined by mass conservation:

 N^{n+1} = N^n \left(\frac{r^n}{r^{n+1}}\right)^3 (22.7-32)

It is assumed that the child droplets are neither distorted nor oscillating; i.e., $y = (dy/dt) = 0$. The child droplets are represented by a number of child parcels which are created from the original parcel. These child parcels are distributed equally along the equator of the parent droplet in a plane normal to the parent relative velocity vector. The diameter of each of the child parcels is sampled from a Rosin Rammler distribution based on the Sauter mean radius (Equation  22.7-23) and a spread parameter of 3.5.

A velocity component normal to the relative velocity vector, with magnitude computed by Equation  22.7-24, is imposed upon the child droplets. It is decomposed at the equator into components pointing radially outward.


A large number of child parcels ensures a smooth distribution of particle diameters and source terms which is needed when simulating, for example, evaporating sprays.

Wave Breakup Model


An alternative to the TAB model that is appropriate for high-Weber-number flows is the wave breakup model of Reitz [ 299], which considers the breakup of the droplets to be induced by the relative velocity between the gas and liquid phases. The model assumes that the time of breakup and the resulting droplet size are related to the fastest-growing Kelvin-Helmholtz instability, derived from the jet stability analysis described below. The wavelength and growth rate of this instability are used to predict details of the newly-formed droplets.

Use and Limitations

The wave model is appropriate for high-speed injections, where the Kelvin-Helmholtz instability is believed to dominate droplet breakup ( ${\rm We} > 100$). Because this breakup model can increase the number of computational parcels, you may wish to inject a modest number of droplets initially.

Jet Stability Analysis

The jet stability analysis described in detail by Reitz and Bracco [ 301] is presented briefly here. The analysis considers the stability of a cylindrical, viscous, liquid jet of radius $a$ issuing from a circular orifice at a velocity $v$ into a stagnant, incompressible, inviscid gas of density $\rho_2$. The liquid has a density, $\rho_1$, and viscosity, $\mu_1$, and a cylindrical polar coordinate system is used which moves with the jet. An arbitrary infinitesimal axisymmetric surface displacement of the form

 \eta = \eta_0 e^{ikz + \omega t} (22.7-33)

is imposed on the initially steady motion and it is thus desired to find the dispersion relation $\omega = \omega(k)$ which relates the real part of the growth rate, $\omega$, to its wave number, $k=2\pi/\lambda$.

In order to determine the dispersion relation, the linearized equations for the hydrodynamics of the liquid are solved assuming wave solutions of the form

$\displaystyle \phi_1$ $\textstyle =$ $\displaystyle C_1 I_0 (kr) e^{ikz + \omega t}$ (22.7-34)
$\displaystyle \psi_1$ $\textstyle =$ $\displaystyle C_2 I_1 (Lr) e^{ikz + \omega t}$ (22.7-35)

where $\phi_1$ and $\psi_1$ are the velocity potential and stream function, respectively, $C_1$ and $C_2$ are integration constants, $I_0$ and $I_1$ are modified Bessel functions of the first kind, $L^2 = k^2 + \omega/\nu_1$, and $\nu_1$ is the liquid kinematic viscosity [ 299]. The liquid pressure is obtained from the inviscid part of the liquid equations. In addition, the inviscid gas equations can be solved to obtain the fluctuating gas pressure at $r = a$:

 -p_{21} = -\rho_2(U - i \omega k)^2 k \eta \frac{K_0 (ka)}{K_1 (ka)} (22.7-36)

where $K_0$ and $K_1$ are modified Bessel functions of the second kind and $u$ is the relative velocity between the liquid and the gas. The linearized boundary conditions are

$\displaystyle v_1$ $\textstyle =$ $\displaystyle \frac{\partial \eta}{\partial t}$ (22.7-37)
$\displaystyle \frac{\partial u_1}{\partial r}$ $\textstyle =$ $\displaystyle -\frac{\partial v_1}{\partial z}$ (22.7-38)


 -p_1 + 2\mu_1 - \frac{\sigma}{a^2}\left(\eta + a^2 \frac{\partial^2 \eta}{\partial z^2}\right) + p_2 = 0 (22.7-39)

which are mathematical statements of the liquid kinematic free surface condition, continuity of shear stress, and continuity of normal stress, respectively. Note that $u_1$ is the axial perturbation liquid velocity, $v_1$ is the radial perturbation liquid velocity, and $\sigma$ is the surface tension. Also note that Equation  22.7-38 was obtained under the assumption that $v_2=0$.

As described by Reitz [ 299], Equations  22.7-37 and 22.7-38 can be used to eliminate the integration constants $C_1$ and $C_2$ in Equations  22.7-34 and 22.7-35. Thus, when the pressure and velocity solutions are substituted into Equation  22.7-39, the desired dispersion relation is obtained:

\omega^2 + 2 \nu_1 k^2 \omega \left[\frac{I'_1 (ka)}{I_0 (ka)... ...2}\frac{I_1 (ka)}{I_0 (ka)}\frac{I'_1 (La)}{I_1 (La)}\right] =

 \frac{\sigma k}{\rho_1 a^2}(1 - k^2 a^2)\left(\frac{L^2 -a^2... ...a^2}\right) \frac{I_1 (ka)}{I_0 (ka)}\frac{K_0 (ka)}{K_1 (ka)} (22.7-40)

As shown by Reitz [ 299], Equation  22.7-40 predicts that a maximum growth rate (or most unstable wave) exists for a given set of flow conditions. Curve fits of numerical solutions to Equation  22.7-40 were generated for the maximum growth rate, $\Omega$, and the corresponding wavelength, $\Lambda$, and are given by Reitz [ 299]:

$\displaystyle \frac{\Lambda}{a}$ $\textstyle =$ $\displaystyle 9.02 \frac{(1 + 0.45{\rm Oh}^{0.5})(1 + 0.4{\rm Ta}^{0.7})}{(1 + 0.87{\rm We}_2^{1.67})^{0.6}}$ (22.7-41)
$\displaystyle \Omega \left(\frac{\rho_1 a^3}{\sigma}\right)$ $\textstyle =$ $\displaystyle \frac{(0.34 + 0.38{\rm We}_2^{1.5})}{(1 + {\rm Oh})(1 + 1.4{\rm Ta}^{0.6})}$ (22.7-42)

where ${\rm Oh} = \sqrt{{\rm We}_1}/{\rm Re}_1$ is the Ohnesorge number and ${\rm Ta} = {\rm Oh} \sqrt{{\rm We}_2}$ is the Taylor number. Furthermore, ${\rm We}_1 = \rho_1 U^2a/\sigma$ and ${\rm We}_2 = \rho_2 U^2a/\sigma$ are the liquid and gas Weber numbers, respectively, and ${\rm Re}_1 = U a/\nu_1$ is the Reynolds number.

Droplet Breakup

In the wave model, breakup of droplet parcels is calculated by assuming that the radius of the newly-formed droplets is proportional to the wavelength of the fastest-growing unstable surface wave on the parent droplet. In other words,

 r = B_0 \Lambda (22.7-43)

where $B_0$ is a model constant set equal to 0.61 based on the work of Reitz [ 299]. Furthermore, the rate of change of droplet radius in the parent parcel is given by

 \frac{da}{dt} = -\frac{(a - r)}{\tau}, \; \; r \leq a (22.7-44)

where the breakup time, $\tau$, is given by

 \tau = \frac{3.726 B_1 a}{\Lambda \Omega} (22.7-45)

and $\Lambda$ and $\Omega$ are obtained from Equations  22.7-41 and 22.7-42, respectively. The breakup time constant, $B_1$, is set to a value of 1.73 as recommended by Liu et al. [ 218]. Values of $B_1$ can range between 1 and 60, depending on the injector characterization.

In the wave model, mass is accumulated from the parent drop at a rate given by Equation  22.7-45 until the shed mass is equal to 5% of the initial parcel mass. At this time, a new parcel is created with a radius given by Equation  22.7-43. The new parcel is given the same properties as the parent parcel (i.e., temperature, material, position, etc.) with the exception of radius and velocity. The new parcel is given a component of velocity randomly selected in the plane orthogonal to the direction vector of the parent parcel, and the momentum of the parent parcel is adjusted so that momentum is conserved. The velocity magnitude of the new parcel is the same as the parent parcel.

You must also specify the model constants which determine how the gas phase interacts with the liquid droplets. For example, the breakup time constant B1 is the constant multiplying the time scale which determines how quickly the parcel will loose mass. Therefore, a larger number means that it takes longer for the particle to loose a given amount. A larger number for B1 in the context of interaction with the gas phase would mean that the interaction with the subgrid is less intense. B0 is the constant for the drop size and is generally taken to be 0.61.

next up previous contents index Previous: 22.7.1 Droplet Collision Model
Up: 22.7 Spray Model Theory
Next: 22.8 Atomizer Model Theory
© Fluent Inc. 2006-09-20