A model for the prediction of the size and velocity distribution of daughter droplets created by the breakup of an unstable parent droplet is proposed. The basis of the model is the maximum entropy formalism, which states that the most probable joint probability density function (JPDF) for the daughter droplet population is the one that maximizes the Bayesian entropy conditional on the enforcement of a set of constraints, which are the conservation laws for the problem. The result is a closed-form expression for the JPDF. Validation against experimental and Direct Numerical Simulations data over the bag, multimode, and sheet-thinning breakup regimes is included. Predictions from the model show that the daughter droplet velocity distribution widens as the droplet size decreases. This result is due to a heightened sensitivity to drag force with lower droplet inertia and coincides with spray behavior. The velocity distribution is found to be near Gaussian. The model does not treat size and velocity as independently distributed, as generally assumed in the literature. In fact, marginal conditional densities derived from JPDF show that the distribution of size and velocity are interrelated.

## I. INTRODUCTION

The conventional picture of spray atomization consists of a primary breakup process, where the liquid core undergoes fragmentation, and a secondary breakup process, where the unstable liquid ligaments or droplets (referred to as parent droplets) are exposed to further breakup.^{1} The scope of the present work is on the development of a new model for predicting the size and velocity distribution of the population of daughter or secondary droplets generated by the breakup of an unstable parent droplet. Existing experimental investigations on droplet breakup (e.g., see Hinze,^{2} Simpkins and Bales,^{3} Krzeczkowski,^{4} Pilch and Erdman,^{5} Wierzba,^{6} Hsiang and Faeth,^{7,8} and Gelfand^{9}) have shown that the droplet breakup process is primarily characterized by the dimensionless gas-based Weber number,

where *ρ _{G}* is the density of the surrounding gas medium, $up,o$ is the initial slip velocity magnitude between the parent droplet and the surrounding gas flows, $Dp,o$ is the initial size of the parent droplet in diameter, and

*σ*is the droplet surface tension coefficient. The

*We*gauges the ratio of the destabilizing drag force to the restoring surface tension force. The fact that this number is a gas-based Weber number rather than liquid-based is a result of dynamic pressure from the gas side, which is responsible for the deformation and breakup of the droplet.

_{G}^{10}If this ratio is above a critical value ($WeG,c$), the drag force outweighs surface tension, inciting the disturbances on the droplet surface, which, in turn, causes the droplet to break up. A value of $WeG,c$ equal to 12 has been experimentally confirmed by Pilch and Erdman

^{5}and Lefebvre and McDonell.

^{10}Depending on the magnitude of the liquid-based Ohnesorge number,

*Oh*, this non-dimensional quantity also affects the breakup characteristics and breakup threshold (Ref. 7, Fig. 1). However, once $OhL\u22720.1$, the influence of

_{L}*Oh*is essentially negligible, and

_{L}*We*controls the breakup process. For all cases considered in the present work,

_{G}*Oh*is at most $O(10\u22122)$, and thus, we limit our considerations to the

_{L}*We*alone.

_{G}Regarding the detailed studies of the droplet breakup dynamics, readers can refer to the review papers by Gelfand,^{9} Guildenbecher *et al.,*^{11} and Theofanous.^{12} A greater value of *We _{G}* is associated with a more violent droplet breakup mode. A classification of the breakup modes includes bag (

*We*$\u2248$ 20), multimode (

_{G}*We*$\u2248$ 60), and sheet-thinning breakups (

_{G}*We*$\u2248$ 100), as documented in a series of experimental studies of droplet breakup by Faeth and his group.

_{G}^{1,7,8,13,14}These works were mainly focused on a single droplet breakup condition. In the recent paper by Wang

*et al.,*

^{15}more complex breakup cases where parent droplets are in tandem formation are considered and are experimentally investigated.

Classical approaches for modeling secondary droplet population in droplet breakup generally treat droplet size and velocity as independently distributed.^{16–19} Besides its use in droplet breakup, the size-velocity independence assumption is also adopted in spray impingement on a flat surface.^{20} To date, numerous models^{17,19,21–40} aimed at describing the droplet size statistics based on various approaches have been proposed. The first set of models is completely empirical, where an optimal fitting function is sought. Prominent examples include Nukiyama-Tanasawa,^{21} Rosin-Rammler,^{22} log-normal,^{23} upper-limit,^{24} chi-squared,^{25} and root-normal distribution functions.^{26} The main problem with this approach is that the model performance can significantly deteriorate if applied outside the experimental regime considered during the curve-fitting process.

The second set of models has been developed within the context of Lagrangian–Eulerian spray simulations. A popular model in this category is the Taylor Analogy Breakup (TAB) model advanced by O'Rourke and Amsden.^{17} Under this model, O'Rourke and Amsden treated the competing mechanisms of destabilizing and stabilizing forces acting on the droplet as a damped/forced harmonic motion. The resulting equation is a linearized ordinary differential equation governing the temporal evolution of droplet deformation. A drawback to the TAB model is that it under-estimates droplet size in diesel sprays.^{27} To remedy this problem, Tanner^{27} proposed the enhanced TAB model (ETAB) by assuming that the ratio of droplet generation rate to the droplet number is constant, and different constants are proposed corresponding to different droplet breakup regimes. Another noticeable drawback of TAB is that it cannot capture the non-linearity in large droplet deformation.^{28} As such, Ibrahim *et al.*^{28} proposed the droplet deformation and breakup (DDB) model, where they assumed that the time rate of change of the droplet internal energy (sum of the kinetic and potential energies) was due purely to work done on the droplet due to pressure and viscous force. Ibrahim *et al.*^{28} showed that for the medium Weber number case, such as Weber number = 50, DDB outperforms TAB. Rimbert *et al.*^{29} later improved DDB by modifying the way of pressure effect treatment. Recent advances in the TAB method can be found in Obenauf and Sojka,^{40} where an improved breakup criterion is introduced along with the inclusion of turbulence effects within the droplet.

Another route in estimating the secondary droplet size statistics is the Rayleigh–Taylor (RT) model based on Rayleigh–Taylor instability. One early use of the RT model was reported by Bellman and Pennington,^{41} where the surface tension and viscosity effects on the liquid-gas interface were analyzed in detail. Bower *et al.*^{42} later applied the RT instability to study droplet breakup in sprays. A systematic modification of the RT model was subsequently presented by Reitz and his group.^{19,30,31} Among these RT adaptations, the one proposed by Beale and Reitz^{19} is generally regarded as the most polished version.

A modeling strategy that employs the stochastic nature of the droplet breakup process has also been attempted. For instance, Gorokhovski and Saveliev^{32} proposed a stochastic model using Kolmogorov's theory. In this model, the mean number of daughter droplets is assumed to be independent of the parent droplet size, and as such, the droplet size distribution can be shown to take the form of a power function. Apte *et al.*^{33} proposed a hybrid Large-Eddy Simulations (LES)-stochastic model and apply it to spray simulations. The model of Gorokhovski and Saveliev is later revised by the same authors,^{34} with particular attention to the asymptotic behavior of the droplet size evolution.

Models based on the joint use of the probability theories and molecular dynamics are also reported in the literature. These models aim at calculating droplet breakup in turbulent flows. One early model was proposed by Luo and Svendsen.^{35} The model assumes that a breakup event occurs only when turbulent eddies hit a droplet. By analyzing the eddies' collision frequency, a droplet breakup probability function can be obtained. Recent progress in this methodology has extended the models' capabilities to a full spectrum of turbulence, where the energy-containing and the dissipation subranges are also present.^{37–39}

While the models mentioned above provide reasonable estimates of droplet size, they have some fundamental shortcomings. The crucial issue is that they assume a highly idealized physical process in their description of the droplet breakup physics. Specifically, as revealed by Direct Numerical Simulations (DNS),^{43–47} the breakup process begins with severe droplet deformation, leading into a bag/thin sheet formation, interfacial instabilities, ligament destabilization, and fragment dispersion. Moreover, the dynamics are highly dependent on the Weber number. Hence, the actual breakup process is far more complex than what the idealized views suggest. For example, the TAB model^{17} assumes a simplified oscillation treatment, significantly overlooking the underlying complexity.

Another issue is that the models generally treat droplet size and velocity as independently distributed borrowing from ideas used in spray modeling.^{16–20} Again, the reason for this assumption is the reduction of modeling complexity. However, this assumption does not agree with experimental results.^{8,48} For instance, Hsiang and Faeth^{8} show that there is a non-linear correlation between daughter droplet velocity and daughter droplet size. Likewise, the work by Guildenbecher *et al.*^{48} demonstrates that the droplet velocity distribution shrinks as droplet size increases. A detailed analysis of this assumption will be presented in Sec. IV E.

Due to the significant complexity of the breakup process and the need to model it without the overburdening computational expense of DNS, an attractive alternative is to employ the maximum entropy formalism (MEF).^{49} This principle is based on an underlying JPDF that describes the secondary droplet size and velocity distributions. The principle then states that the JPDF that maximizes the statistical entropy while simultaneously obeying a set of physical constraints is the most likely JPDF and the one chosen to characterize a given breakup process. Hence, the application of MEF targets *the results* of hydrodynamic breakup rather than all of the intermediate dynamics. Since the results are of interest in many applications, the use of MEF is well suited for this purpose, and it is the one implemented in the present work.

Thus far, most attempts at employing MEF in breakup problems have targeted the full liquid jet atomization and resulting spray formation problem.^{50–64} In those models, the maximization of statistical entropy is coupled to the constraints imposed by physical conservation laws. For the interested reader, informative reviews are provided by Babinsky and Sojka^{65} and Dumouchel.^{60}

In the case of droplet breakup, the MEF has scarcely been applied. To the authors' best knowledge, only the paper by Bodaghkhani *et al.*^{64} employs this method, in which their treatment inherits the same formulation as in primary atomization of sprays. Hence, the present work represents one of the first studies applying MEF in secondary droplet breakup problems. The new features of our model include the following:

A new statistical description of the breakup process is presented, forming the basis for the ensemble-averaged equations, which are subsequently solved via the MEF method.

The model is derived in a secondary droplet breakup framework, which differs from the existing MEF models derived based on the spray atomization framework.

New sub-models for calculating the source terms included in momentum and energy conservation equations are developed.

In what follows, we begin with a statistical description of a droplet breakup event in Sec. II which forms the foundation for the development of the MEF model, presented in Sec. III. Results consisting of a qualitative assessment of the model performance, experimental and numerical validation, impact of the uncertainty in the sub-model for momentum and energy losses, and an examination of the independence of the size-velocity distribution are all presented in Sec. IV. Finally, a summary of the work together with the concluding thoughts is given in Sec. V.

## II. STATISTICAL DESCRIPTION OF DROPLET BREAKUP

Consider the breakup of an unstable parent droplet into a population of secondary or daughter droplets. Due to the stochastic nature of the breakup process, the discrete size and velocity distribution of the daughter droplet population, and the total number of daughter droplets, are not expected to remain the same from realization to realization (breakup event to breakup event). That is, if we were to consider, for example, three different instances of an unstable parent droplet characterized by some velocity, size, fluid, and ambient conditions, the resulting daughter droplet populations would differ. This result is illustrated in Fig. 1, indicating the randomness of the breakup phenomenon.

Despite the randomness of the breakup event, it must satisfy the following conversation laws for mass, momentum, and energy:

where subscript *q* denotes an arbitrary breakup event, *M _{q}* is the number of secondary droplets created, $Dp,o$ is the initial parent droplet size in diameter,

*D*is the daughter droplet diameter, $up,o$ is the initial velocity of the parent droplet, and

_{i}*u*is the daughter droplet velocity. Also,

_{i}*σ*is the surface tension coefficient,

*ρ*is the mass density, and subscripts (

*L*,

*G*) denote liquid and gas. Furthermore, Δ

_{mom}and Δ

_{E}represent the momentum and energy losses during the breakup process due to the aerodynamic drag force exerted by the surroundings.

To arrive at Eqs. (2)–(4), we make the following assumptions:

The lateral velocity of daughter droplet populations is assumed to be negligible. Since the breakup model is targeted for use in spray simulations, the velocity of the liquid in the near-field of the spray, where the secondary droplet process mostly occurs, is heavily weighted toward the initial parent velocity direction with a minor contribution from the radial or lateral directions.

The ambient velocity is treated as zero, i.e., we are considering the general case of spray formation in a quiescent environment (this can be generalized to a non-zero condition without much trouble).

All daughter droplets are assumed to be spherical. As suggested by Helenbrook and Edwards,

^{66}the assumption of sphericity holds for $WeG<O(10)$. In our calculations, we have confirmed that the great majority of daughter droplet populations satisfy this criterion. However, there are secondary droplets with $WeG>12$, and thus further breakup is applied to these unstable droplets.

For the energy conservation [Eq. (4)], both kinetic energy and surface energy are considered as a single constraint. In some MEF papers, such as Sellens and Brzustowski,^{50} Sellens,^{52} and Chin *et al.,*^{53} it is argued that using a single energy constraint cannot provide information on how the total energy is distributed between these two energy modes. Hence, they introduced two separate energy conservation constraints, one for kinetic energy and the other for surface energy. However, employing separate conservation equations for surface energy and kinetic energy suggests a violation of the conservation laws established in continuum mechanics; thus, this approach is *not* employed in the present work.

We begin the statistical description by introducing the JPDF, which is denoted by $fDu(D,u)$. Following convention, it obeys the normalization requirement,

The expected value or ensemble average of any quantity $\u27e8Drus\u27e9$ (*r* and $s\u2208\mathbb{Z}$) can be obtained from

where the sample space, $\Omega \u2282\mathbb{R}2$, is spanned by $D\u2208[Dmin,Dmax]$ and $u\u2208[umin,umax]$. Here, *D _{min}*,

*D*,

_{max}*u*, and

_{min}*u*are, respectively, the lower and upper bounds of

_{max}*D*and

*u*, both of which are discussed below in Sec. III A. For secondary droplet statistics, an equivalent manner of expressing an ensemble average is

where *M*_{1}, *M*_{2}, and *M _{q}* are the total number of secondary droplets in realization 1, 2, or generally in realization

*q*. The expected value of

*M*can be similarly obtained from

To obtain an ensemble-averaged mass conservation equation, we begin with Eq. (7) evaluated at *r *=* *3 and *s *=* *0, namely, the expected value of *D*^{3},

Since the mass conservation expression in Eq. (2) holds for any breakup event, i.e.,

this allows to rephrase Eq. (9) as

Noting that the parent droplet size $Dp,o$ does not change from realization to realization yields

Introducing Eq. (8) in this expression gives the ensemble-averaged form of mass conservation as

For momentum, we start with the expected value of $D3u$ given by

From Eq. (3), the conservation of momentum holds for all realizations, i.e.,

where Δ_{mom} is a deterministic quantity that is only dependent on parent droplet conditions, and thus it remains unchanged in all breakup events. Substituting Eq. (15) into Eq. (14) gives

Using the expression for $\u27e8M\u27e9$ in Eq. (8) yields the ensemble-averaged form for momentum conservation,

Last, for energy conservation, the goal is to obtain an expected value expression for the kinetic and surface energy term. Beginning with

we note that energy is conserved for each breakup event, namely,

Again, since the parent droplet conditions do not change from realization to realization, this yields the ensemble-averaged form of energy conservation,

where Eq. (8) has been used.

In the MEF model, explicit use of $fDu(D,u)$ is made. Hence, the expected values appearing in Eqs. (13), (17), and (21) are expressed in terms of this JPDF, namely,

where $\Delta \xafmom=\Delta mom/[(\pi /6)\rho LDp,o3up,o]$ and $\Delta \xafE=\Delta E/[(\pi /12)\rho LDp,o3up,o2]$ are momentum loss and energy loss over the breakup process in dimensionless form. These terms are discussed in Sec. II A, Also, $WeL=(\rho Lup,o2Dp,o)/\sigma $ is the liquid-based Weber number. Equations (25)–(27), along with the normalization constraint of $fDu(D,u)$ defined in Eq. (5), constitute a set of constraints that $fDu(D,u)$ must satisfy.

### A. Model for estimating momentum loss ($\Delta \xafmom$) and energy loss ($\Delta \xafE$)

We begin with the equation of motion for a droplet,^{18}

where only the drag and gravitational terms are accounted for since $\rho L\u226b\rho G$. In this expression, $up(t)$ denotes the parent droplet velocity, $Dp(t)$ the parent droplet size, $CD(t)$ the drag coefficient, $uEul(xd(t))$ the gas velocity at the droplet position $xd(t)$, and **g** the gravity. In the current model, droplet size and drag coefficient are treated as functions of time. Considering the breakup of a droplet into a quiescent environment, such that $uEul(xd(t))$ is set to zero (non-zero values for gas velocity can also be considered), Eq. (28) can be simplified to

where $Dp(t),\u2009up(t)$, and $CD(t)$ are unknowns. In this expression, it is being acknowledged that drag forces dominate, and gravitational effects are almost negligible. This material is illustrated in greater detail in Appendix A.

Address first $Dp(t)$. This parameter is obtained by matching the experimental findings of a droplet breakup reported by Hsiang and Faeth,^{7} where this quantity is described by

over the range of ($10<WeG<105$). Note that our interest is for $WeG>12$, and thus the above relationship applies. Here, *D _{max}* is the maximum droplet size at the onset of the breakup, which can be correlated with

*We*

_{G}^{7}as

The characteristic time for breakup has been reported by Nicholls and Ranger^{67} as

The change in drag coefficient due to a deforming droplet, $CD(t)$, is given empirically by Hsiang and Faeth^{7} as

The relation in Eq. (34) holds for *Re _{G}* in the range of 1000–2500, where $ReG=(\rho GDp,oup,o)/\mu G$, and

*μ*is the dynamic viscosity of the gas. In the limits of no deformation, i.e., $Dp(t)=Dp,o$, Eq. (34) converts back to the drag coefficient for a solid sphere.

_{G}^{68}In our calculations, the test

*Re*range is 773–2446, corresponding to a range for

_{G}*We*from 12 to 120. Hence, Eq. (34) is mostly applicable.

_{G}Returning to the parent droplet's equation of motion, substitutions for $Dp(t)$ [Eq. (33)] and $CD(t)$ [Eq. (34)] are incorporated into Eq. (29), giving

This expression can be manipulated to yield

Integrating from *t *=* *0 to an arbitrary value of *t* gives

Solving this equation for $up(t)$ gives

Let *t _{break}* denote the time required to reach the droplet breaking point. Then the corresponding velocity at this time is

Pilch and Erdman's experimental correlation of *t _{break}*

^{5}is

where *t _{char}* was previously defined in Eq. (32). Substituting this empirical correlation into the equation for $up(tbreak)$ [Eq. (39)] results in

where $\Phi =(WeG\u221212)\u22120.25$.

With expressions for the final droplet velocity at the breakup point, the momentum loss (Δ_{mom}) is given by

or in dimensionless form as

Since Eq. (44) shows that $\Delta \xafmom$ is a function of *We _{G}* and the liquid-gas density ratio ($\rho L/\rho G$), an instructive exercise is to plot a 2D iso-surface contour of $\Delta \xafmom$ in terms of

*We*and $\rho L/\rho G$. The result is shown in Fig. 2. The parameter range considered for

_{G}*We*is 20–100, covering the bag, multimode, and sheet-thinning breakup regimes. For $\rho L/\rho G$, the range is between 20 and 1000, which is the expected range employed in engineering sprays. The results from Fig. 2 demonstrate that $\Delta \xafmom$ is not sensitive to changes in

_{G}*We*, remaining nearly unchanged as this parameter is varied. However, it is highly sensitive to $\rho L/\rho G$, particularly as this quantity decreases from approximately 265 to 1000. This behavior is attributed to the functional dependence of momentum loss to

_{G}*We*and $\rho L/\rho G$, which from Eq. (44) can be estimated to be $\Delta \xafmom\u223c(\rho L/\rho G)1/2$ and $\Delta \xafmom\u223cWeG1/4$.

_{G}For Δ_{E}, we consider a kinetic energy form, which is related to Δ_{mom} via

or in dimensionless form as

Thus,

Based on the same parameter ranges considered in Fig. 2, the global maximum of $\Delta \xafE$ is only around 0.04.

An unlikely, but potential issue, exists with the proposed equations for $\Delta \xafmom$ and $\Delta \xafE$ at the precise stability limit, i.e., at *We _{G}* = 12. At this boundary, $\Delta \xafmom=1$ and $\Delta \xafE=1$. This result indicates that 100% of the parent droplet momentum is consumed by drag force during the breakup process. In this circumstance, the momentum equation defined in Eq. (26) indicates that the expected value of daughter droplet momentum is zero, which is unphysical. This problem is caused by the singularity of the breakup time, $tbreak=1.9(WeG\u221212)\u22120.25tchar$, which diverges to infinity at

*We*= 12. For

_{G}*We*sufficiently close to 12, it will take an impractically long time for the breakup to occur, therefore depleting all the parent droplet momentum. To avoid this unphysical situation, another empirical time suggested by Guildenbecher

_{G}*et al.*

^{11}in initiating a bag breakup event, 2

*t*, is also considered. The revised expression for

_{char}*t*is

_{break}## III. ESTABLISHING THE MEF MODEL

The maximum entropy principle states that potentially infinite numbers of distributions can satisfy a set of physical constraints. Among these, the most likely probability distribution is the one that maximizes the measure of entropy as explained by Kapur.^{49} The most classical measure is the Shannon entropy,^{56,69,70}

Unfortunately, the Shannon entropy does not take any underlying physics-based or empirical information into account. For hydrodynamic breakup problems, exercising the maximum entropy principle does not generally yield a globally proper JPDF. A noticeable issue has been reported by Babinsky and Sojka^{65} in the limit as $D\u21920$, where $fDu(D,u)$ is expected to be zero, but an implementation based on the maximization of the Shannon entropy produces an $fDu(D,u)$ that does not approach zero in this limit. The expectation of this zero limit comes from the fact that below some minuscule droplet size, there are no more droplets. For circumventing this issue, the Bayesian entropy is used instead of the Shannon entropy, which is defined^{57} as

where $fo(D)$ is *a priori* droplet size distribution known empirically for small droplet populations characterized by $D\u226aD30$. Here, *D*_{30} is the mass mean diameter and is defined as

Comparing Shannon entropy with Bayesian entropy shows that Shannon entropy is a particular case of Bayesian entropy with $fo(D)=1$, implying that all *small* droplets have the same size, which is generally not valid in breakup problems.^{65} The details of $fo(D)$, *D*_{30}, and their relation to the present MEF work will be discussed later in Sec. III A.

To obtain an $fDu(D,u)$ that maximizes the Bayesian entropy while satisfying the physical conservation laws, we employ the Lagrange multiplier method described in Kapur.^{49} This JPDF is further refined by enforcing that the small size droplet distribution ($D\u226aD30$) must follow a prescribed probability density function (PDF) given by $fo(D)$. First, the constraints imposed on $fDu(D,u)$[see Eqs. (5) and (25)–(27)] are re-cast in the form of $gi(fDu(D,u))=0,i=$ 0, 1, 2, and 3, namely,

corresponding to the normalization requirement, mass, momentum, and energy conservation.

Next, four Lagrange multipliers, *λ _{i}*, $i=$ 0, 1, 2, and 3, are introduced into the following function:

Under the Lagrange multiplier method, the optimum $fDu(D,u)$ is obtained from

After some manipulations the above equation becomes

Letting $\lambda 0*=1+\lambda 0$, the optimum JPDF that maximizes Bayesian entropy in the presence of the constraints is found to be

Overall, the MEF system consists of four equations, Eqs. (52)–(55) with six unknowns, namely, [$\lambda 0*,\lambda 1,\lambda 2,\lambda 3,fo(D)$, and $\u27e8D3\u27e9$]. To have a unique solution to this system requires introducing two more independent equations that $fo(D)$ and $\u27e8D3\u27e9$ have to satisfy.

### A. Development for $fo(D)$

A mathematical form for $fo(D)$ that matches the empirical size distributions for $D\u226aD30$ is sought. Dumouchel^{59} indicates that $fo(D)$ should be related to *D*, due partially to the role of surface tension forces in characterizing liquid-gas interfaces. Continuing this idea, Dumouchel^{59} suggests that $fo(D)$ should follow a power-law distribution with respect to *D* for small droplets. Following this suggestion, we assign $fo(D)$ the well-established Nukiyama–Tanasawa distribution,^{10} given by

where *B* is a constant, *C* is the size parameter, and *p* and *q* are distribution parameters.

We will treat $fo(D)$ as a PDF, which automatically implies obeying the normalization requirement. González-Tello *et al.*^{71} suggest that in spray problems, the size of small droplets characterized by $D\u226aD30$ is much smaller than the size parameter *C*, i.e., $D\u226aC$. Expanding out the exponential term in Eq. (61) yields

Also, Li and Tankin^{51} point out that the constant *q* is typically greater than one. Hence, the entire denominator in Eq. (62) can be approximated with good accuracy to be equal to one. Furthermore, the parameter *p* is determined by matching the empirical observations of Li and Tankin^{51} and Lefebvre and McDonell^{10} in liquid jet breakup problems, which gives $p=$ 2. With all of this information, $fo(D)=BD2$, and from the normalization requirement

gives

Following the procedure documented by Movahednejad *et al.*^{61} for spray problems, we assign the following respective bound values:

Hence, the final form of $fo(D)$ is

### B. Development for $\u27e8D3\u27e9$

Regarding $\u27e8D3\u27e9$, it is constructed by performing a non-linear regression to the DNS data from Jain *et al..*^{46} The correlation is then calibrated in line with spray measurements of *D*_{30}, i.e., $\u27e8D3\u27e91/3$, from Kim *et al.,*^{57} yielding

Equation (67) shows that $\u27e8D3\u27e91/3$ is correlated with $up,o$ via $\u27e8D3\u27e91/3\u221dup,o\u22121.48$. The droplet breakup experiments from Wolfe and Andersen^{72} also report a similar empirical correlation of $\u27e8D3\u27e91/3\u221dup,o\u22121.33$, which provides some degree of reassurance. Combining $fo(D)$ together with the general JPDF in Eq. (60) gives the final expression for $fDu(D,u)$ as

## IV. RESULTS AND DISCUSSION

The results are divided into five sections. In Sec. IV A, a qualitative assessment of $fDu(D,u)$ is completed to confirm whether the key features of this distribution resemble those predicted in previous MEF spray breakup studies, namely, in the works of Sellens,^{52} Li *et al.,*^{73} Ayres *et al.,*^{56} Movahednejad *et al.,*^{61} Yan *et al.,*^{62} and Fu *et al.*^{63} In Sec. IV B, the MEF model is compared to experiments from Hsiang and Faeth.^{7,8} A similar comparison is reported in Sec. IV C, but this time the source is from a DNS by Jalaal and Mehravaran.^{43} Since the MEF model includes source terms to address the effects of drag force on droplet momentum and energy, the sensitivities of uncertainties in these source terms to $fDu(D,u)$ are studied and presented in Sec. IV D. Finally, the independence of size and velocity distributions is examined.

### A. Qualitative assessment of $fDu(D,u)$

The characteristics of $fDu(D,u)$ are examined for a case consisting of the breakup of a droplet within the multimode breakup regime. This regime is chosen because of its relatively complex breakup mechanism compared to the bag and sheet-thinning regimes.^{74} As such, it can serve as a good benchmark for assessing the model performance. The droplet, in this case, has a diameter ($Dp,o$) of 200 *μ*m and a velocity ($up,o$) of 135.07 m/s, with a respective Weber number (*We _{G}*) of 60. The physical properties are listed in Table I. The predicted $fDu(D,u)$ using the MEF model is shown as a 2D surface plot in Fig. 3.

Liquid density (water) | ρ = 997 kg/m_{L}^{3} | ||

Gas density (air) | ρ = 1.1839 kg/m_{G}^{3} | ||

Coefficient of surface tension | σ = 0.072 kg/s^{2} |

Liquid density (water) | ρ = 997 kg/m_{L}^{3} | ||

Gas density (air) | ρ = 1.1839 kg/m_{G}^{3} | ||

Coefficient of surface tension | σ = 0.072 kg/s^{2} |

The first noticeable characteristic from Fig. 3 is that the velocity distribution becomes wider with decreasing drop size. This behavior is attributed to an increased sensitivity to drag force with small droplet inertia. Similar findings have also been reported for sprays by Sellens.^{52} A second characteristic is that $fDu(D,u)$ decays significantly for droplets below approximately 2 *μ*m, reaching zero at D $\u226a$ 2 *μ*m. This result is attributed to the $fo(D)$, i.e., small droplet distribution. Another feature is that $fDu(D,u)$ exhibits the near-symmetric velocity distribution about $up,o$. With regard to spray problems, similar behavior was also documented in Sellens and Brzustowski,^{50} where the authors state that the symmetric behavior is very similar to a Gaussian profile. Sellens and Brzustowski^{50} explain this observation by arguing that the liquid momentum is a linear function in the velocity space, and the liquid kinetic energy is a square function of velocity. These two constraints taken in combination should allow for the estimation of the mean and variance of velocity. With the information of mean and variance, a Gaussian distribution is the most likely distribution that can be deduced. As a means of confirming this near-Gaussian behavior for the velocity distribution, the respective marginal velocity density function,

is compared against a Gaussian distribution in Fig. 4. For this Gaussian distribution, the mean value is equal to the parent droplet velocity $up,o$, and the standard deviation is 0.26 $up,o$. Overall, good qualitative matching is observed, confirming expectations from Sellens and Brzustowski.^{50}

Since the droplet breakup severity is heavily influenced by *We _{G}*, it is expected that $fDu(D,u)$ would reflect this sensitivity. To examine this aspect, we consider a parent water droplet having a size of 200

*μ*m and impose different initial velocities of 77.98, 135.07, and 174.38 m/s, corresponding to

*We*of 20, 60, and 100, respectively. The choice of these Weber numbers ensures that the three regimes of the breakup, namely, bag, multimode, and sheet-thinning breakup, are considered. The surrounding medium is air, and physical properties are the same as those listed in Table I. Surface plots of $fDu(D,u)$ pertaining to the three different cases are shown Fig. 5. Similar to the previous case, the density function is symmetric about the parent droplet velocity. The maximum value of the JPDF, $fDu(D,u)|max$, is located at (

_{G}*D*,

*u*) = (12.86

*μ*m, 74.71 m/s), (5.75

*μ*m, 131.01 m/s), and (3.95

*μ*m, 169.74 m/s) corresponding, respectively, to

*We*= 20, 60, and 100. Hence, there is a natural shift to higher velocities and smaller drop sizes with increasing

_{G}*We*.

_{G}A related trend displayed in Fig. 5 is the concentration of progressive smaller daughter droplets with increasing breakup intensity, i.e., larger values for *We _{G}*. This pattern is quantified by calculating the size variance of $fDu(D,u)$, i.e., $\sigma D2$ = $\u222bDminDmax\u222buminumax(D\u2212\u27e8D\u27e9)2fDu(D,u)dudD$. This variance, $\sigma D2$, is found to be equal to 22.43, 4.56, and 2.14 $\mu m2$, corresponding again to the respective

*We*= 20, 60, and 100 cases. This characteristic is one of the most distinct features of the breakup process with increasing

_{G}*We*.

_{G}### B. Experimental validation

Breakup predictions from the MEF model in the bag, multimode, and sheet-thinning breakup regimes are compared to the experiments from Hsiang and Faeth.^{7,8} Since in Hsiang and Faeth's work,^{7} a variable drop size is used for different test fluids, which include water, n-heptane, ethyl alcohol, mercury, and glycerol solutions, in the model predictions, we chose to let the initial parent droplet size match the one for water drop conditions, which is 1000 *μ*m. The results are quantified using a normalized cumulative volume distribution denoted by CDF(*D*^{3}) as a function of droplet size for different *We _{G}*. The parent droplet is formed using water and its initial velocity is systematically varied to let the parent droplet have a respective

*We*equal to 20, 60, and 100. The surrounding medium is air, and physical properties are the same as those listed in Table I.

_{G}The comparisons between model solutions and the experimental measurements are shown in Fig. 6. In the figure, the *x*-axis is the droplet diameter normalized by the mass medium diameter *MMD*, which is implicitly defined by

The *y*-axis of Fig. 6 shows the values of CDF(*D*^{3}) calculated by

where $fD(D)$ is the marginal size density function given by

Considering the results shown in Fig. 6, at least two observations can be made. The first is that for values *D*/*MMD* below one, the agreement with experimental data is appreciably better than for the range *D*/*MMD* greater than one. The exception to this trend is for the sheet-thinning regime employing *We _{G}* = 100, where the level of agreement between larger and smaller droplets relative to

*MMD*is approximately the same. A second observation concerns the noticeable deviation between the MEF predicted curves and the experiments for the larger secondary droplet populations. MEF does not predict as many larger droplets as experiments would indicate.

#### 1. MMD/$D32$

Another metric that can be used in validating the droplet size characteristics is the ratio of *MMD* to *D*_{32}, where *D*_{32} is the Sauter mean diameter defined as

This metric was originally advanced by Simmons,^{26} where the author showed that $MMD/D32$ was universally equal to 1.2 under various spray operating conditions. Hsiang and Faeth^{7,8} and Guildenbecher *et al.*^{48} later show that this conclusion still holds in droplet breakup experiments in the bag, the multimode, and the sheet-thinning breakup regimes. Our calculations show that the value of $MMD/D32$ is 1.1 in each *We _{G}* case reported, agreeing with the experimental values within 10%.

### C. Numerical validation

We compare our predicted marginal droplet size distributions, $fD(D)$, to results from DNS reported by Jalaal and Mehravaran.^{43} In these simulations, the authors^{43} employed the open-source solver Gerris, a geometric Volume-of-Fluid (gVoF) algorithm with adaptive mesh refinement (AMR) capabilities developed by Popinet,^{75} to capture the liquid–gas interface. A total of three cases are considered having Weber numbers and parent droplet sizes given, respectively, by (*We _{G}*,

*D*) = (37.6, $3.2\xd7104$

*μ*m), (

*We*,

_{G}*D*) = (47.8, $3.9\xd7104$

*μ*m), and (

*We*,

_{G}*D*) = (59, $4.5\xd7104$

*μ*m). The physical properties are the same as listed in Table I. In the work of Jalaal and Mehravaran,

^{43}the breakup is initiated by the falling action of the droplet due to gravity.

Results and comparisons are shown in Fig. 7 according to each of the three conditions considered. The overall agreement is reasonable for the larger secondary droplet, i.e., $D\u226b103\mu m$. For the smaller droplet distribution, the deviation between the model and the DNS data is substantial. Besides recognizing the possibility of modeling error in this part of the droplet size spectrum, other sources of discrepancy should be considered. For instance, in all Weber numbers provided, the $fD(D)$ from the DNS goes to zero at roughly half $\Delta xmin$, where $\Delta xmin$ is the minimum mesh size. Hence, there is not enough spatial resolution to capture these small-scale features. Previous work in our group performing similar DNS^{76} indicates that at least $O(10)$ cells should be employed across the diameter of a droplet to capture its large scale dynamics. Thus, it can be argued that insufficient grid resolution is being employed by the authors^{43} for the smallest part of the daughter droplet population.

### D. Analysis of the source term modeling uncertainties

Consider the momentum and energy conservation equations. Losses in momentum and energy are accounted for through $\Delta \xafmom$ and $\Delta \xafE$, respectively. Since these terms are obtained from models, they are guaranteed not to be perfect. For the sake of argument, let us assume that we have a ±30% uncertainty in $\Delta \xafmom$, which will transcend to ±69% uncertainty in $\Delta \xafE$ since both terms are related by $\Delta \xafE=\Delta \xafmom2$.

A plot of $\Delta \xafmom$ shown in Fig. 2 (see Sec. II A) indicates that its greatest sensitivity is due to changes in liquid-gas density ratio instead of the Weber number. Hence, in the following analysis, we consider three different scenarios corresponding to $\rho L/\rho G$ = 1000, 100, and 50. All of which correspond to *We _{G}* = 60, i.e., the same multimode regime studied previously in Sec. IV A. In each case, the parent droplet has an initial size of 200

*μ*m, a density of 1000 kg/m

^{3}, and a surface tension coefficient of 0.072 kg/s

^{2}.

At each density ratio, the nominal value for $\Delta \xafmom$ is varied ± 30%, and its effects on the resulting secondary droplet population are examined. Here, we are interested in the effects of the uncertainty of $\Delta \xafmom$ on $fDu(D,u)$ as a function of $\rho L/\rho G$, not on the changes of $fDu(D,u)$ as a result of variations of density ratio. This point is subtle but important. The results, expressed in terms of mean values with respect to size and velocity and the respective standard deviations, are included in Table II. The mean values of these statistics are normalized by the respective nominal quantities. For instance, at $\rho L/\rho G=1000$, the mean diameter corresponding to $\Delta \xafmom\xb1(0.3)\Delta \xafmom$ is divided by the mean diameter evaluated at $\Delta \xafmom$. For standard deviations, they are normalized by the mean values. This operation gives a direct indication of the sensitivity of each statistic to the uncertainty associated with $\Delta \xafmom$.

. | $\rho L/\rho G$ = 1000 . | $\rho L/\rho G$ = 100 . | $\rho L/\rho G$ = 50 . | ||||||
---|---|---|---|---|---|---|---|---|---|

$\Delta \xafmom$ . | 30% . | Nominal . | −30% . | 30% . | Nominal . | −30% . | 30% . | Nominal . | −30% . |

$\u27e8D\u27e9/\u27e8D\u27e9nom$ | 1.006 | 1 | 0.991 | 1.014 | 1 | 0.967 | 1.032 | 1 | 0.881 |

$\sigma D/\u27e8D\u27e9$ | 0.458 | 0.467 | 0.473 | 0.439 | 0.463 | 0.529 | 0.452 | 0.504 | 0.790 |

$\u27e8u\u27e9/\u27e8u\u27e9nom$ | 1.003 | 1 | 0.996 | 1.003 | 1 | 0.997 | 0.999 | 1 | 1.018 |

$\sigma u/\u27e8u\u27e9$ | 0.452 | 0.423 | 0.386 | 0.550 | 0.505 | 0.446 | 0.565 | 0.513 | 0.479 |

. | $\rho L/\rho G$ = 1000 . | $\rho L/\rho G$ = 100 . | $\rho L/\rho G$ = 50 . | ||||||
---|---|---|---|---|---|---|---|---|---|

$\Delta \xafmom$ . | 30% . | Nominal . | −30% . | 30% . | Nominal . | −30% . | 30% . | Nominal . | −30% . |

$\u27e8D\u27e9/\u27e8D\u27e9nom$ | 1.006 | 1 | 0.991 | 1.014 | 1 | 0.967 | 1.032 | 1 | 0.881 |

$\sigma D/\u27e8D\u27e9$ | 0.458 | 0.467 | 0.473 | 0.439 | 0.463 | 0.529 | 0.452 | 0.504 | 0.790 |

$\u27e8u\u27e9/\u27e8u\u27e9nom$ | 1.003 | 1 | 0.996 | 1.003 | 1 | 0.997 | 0.999 | 1 | 1.018 |

$\sigma u/\u27e8u\u27e9$ | 0.452 | 0.423 | 0.386 | 0.550 | 0.505 | 0.446 | 0.565 | 0.513 | 0.479 |

The general trend observed from the results in Table II is one of increasing uncertainty impacts with decreasing density ratio. Only the mean velocity shows negligible variation with different values of $\rho L/\rho G$, but, even in this case, there is a slighter larger variation at $\rho L/\rho G=50$ in comparison to the other cases. This heightened sensitivity to uncertainties in $\Delta \xafmom$ at the low end of the density ratio spectrum is a direct result of increased momentum exchange between the droplet and its surroundings due to large gas density. Specifically, the low-density ratio is governed by a higher gas density since *ρ _{L}* remains the same. Hence, any variations in $\Delta \xafmom$ at higher

*ρ*will have greater impact on the characteristics of $fDu(D,u)$.

_{G}### E. Examining the droplet size-velocity independence assumption

As discussed, previous spray models often assume that the velocity and size distributions of daughter droplets are independent of each other^{16–20,53} without actually examining the validity of this assumption. To inspect this treatment in the present work, we simply consider the respective conditional densities. For reference, if *D* and *u* are independently distributed, then the corresponding JPDF is simply a product of the marginal densities, i.e., $hDu(D,u)=hD(D)hu(u)$. Consequently, the distribution of velocity conditional on *D* = *D _{a}* is given by

which simply gives the same $hu(u)$ for any value of *D _{a}* due to independence. Similar dependencies can also be shown for the marginal size distribution conditional on a specified velocity of daughter droplet, i.e., $hD|u=ua(D)$.

Regarding the present JPDF corresponding to the MEF, its conditional marginal density is

where $fD(D)$ has a closed-form expression given by (see Appendix B for the derivation)

Here, $\xi (u)=(\lambda 3D3)/\u27e8D3\u27e9[u/up,o+\lambda 2/(2\lambda 3)],\u2009\xi max=\xi (umax),\xi min=\xi (umin)$, and $erf(\xi (u))=(2/\pi )\u222b0\xi (u)\u2009exp\u2009(\u2212x2)dx$ is the error function of $\xi (u)$.

which after some manipulations can be simplified as (see Appendix C)

This expression clearly shows that the marginal density conditional on *D* = *D _{a}* is explicitly dependent on

*D*, i.e., the independently distributed condition does not apply. In particular, $fu|D=Da(u)$ is exponentially correlated with $\u2212Da3$. A numerical illustration is given in Fig. 8. Similar conclusions also hold for the marginal density of size conditional on a prescribed velocity.

_{a}## V. CONCLUSIONS

Motivated by the complexity of the droplet breakup phenomena and the associated burdensome costs of performing DNS, the present effort provides an alternative treatment offered by the maximum entropy formalism. Hence, the model is not targeted at all of the complicated intermediate dynamics of the breakup process; instead, it is focused on the final outcome in terms of secondary droplet distribution with respect to velocity and size. In various practical applications, it is precisely the results of the breakup process that are of interest, and thus the present model suits this purpose. This model would also be useful in situations where many droplets exist under hydrodynamically unstable conditions and are too expensive to be resolved by a DNS. The model accounts for all the conservation laws and includes submodels for momentum and energy losses during the breakup. It provides in *closed form* an expression for the most probable joint probability density function in size and velocity space based on the maximization of the Bayesian entropy. While we would prefer an empiricism-free solution, we acknowledge that due to the complexity of the phenomena, which is breakup regime dependent and which is characterized by pronounced changes in droplet shape and fragmentation, some empirical relations have been introduced. This was primarily done in the closure of $fo(D)$ and $\u27e8D3\u27e9$.

Predictions from the model show that the droplet velocity distribution widens as the droplet size decreases. This result is due to a heightened sensitivity to drag force with lower droplet inertia and coincides with spray behavior documented by Sellens.^{52} In addition, the velocity distribution is found to be near the Gaussian distribution. For the size distribution, it is similar to a root-normal distribution, agreeing with Guildenbecher *et al.*'s experimental observations on droplet breakup.^{48} As the Weber number increases, the model shows that the size distribution skews toward more concentrated smaller droplet sizes in agreement with experiments. Further analysis shows that the size and velocity are dependently distributed based on an inspection of the conditional marginal densities. This last feature contrasts with common assumptions made in the literature claiming that size and velocity are independently distributed.

## ACKNOWLEDGMENTS

Support from Caterpillar Inc. and the Army Research Laboratory is gratefully acknowledged. A portion of this research was sponsored by the Army Research Laboratory and was accomplished under Cooperative Agreement No. W911NF-20-2-0181. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Laboratory or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein. CEI is acknowledged for granting use of their post-processing software, EnSight.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: ANALYSIS OF GRAVITATIONAL EFFECTS

Let us return to the droplet momentum equation defined in Eq. (28). The drag force term is

while the gravitational term is

Recalling $uEul(xd(t))$ = 0 in the present work, the difference in the magnitude between these two terms is

Since

and

it then follows from Eq. (A3) that the gravitational term is nearly negligible in the droplet momentum equation.

The other metric for evaluating the influence of the gravitational effects is the Bond number (*Bo*),

which quantifies the ratio of gravitational force to surface tension force. A 2D surface plot of *Bo* is shown in Fig. 9 as a function of *ρ _{L}* and $Dp,o$. In this calculation, the surface tension coefficient is given as 0.072 kg/m

^{3}; the maximum values of $Dp,o$ and

*ρ*are given as 1000

_{L}*μ*m and 1000 kg/m

^{3}. The results show that Bond number is $O(10\u22122)$–$O(10\u22121)$, which confirms the almost negligible effect of gravity.

### APPENDIX B: DERIVATION OF THE ANALYTICAL EXPRESSION FOR $fD(D)$

The procedure begins with substituting the final form of $fDu(D,u)$ [Eq. (68)] into the definition of $fD(D)$ [Eq. (72)], yielding

The integrand of Eq. (B2) after some manipulations can be re-cast as

Perform change of variable by defining

whose derivative is

where $\xi max=\xi (umax),\u2009\xi min=\xi (umin)$, and $erf(\xi (u))=(2/\pi )\u222b0\xi (u)\u2009exp\u2009(\u2212x2)dx$ is the error function of $\xi (u)$.

### APPENDIX C: SIMPLIFIED EXPRESSION FOR $fu|D=Da(u)$

The expression for $fu|D=Da(u)$ [Eq. (77)] can be simplified as follows: