CONDITIONAL DAMPED RANDOM SURFACE VELOCITY MODEL OF TURBULENT JET BREAKUP

A turbulent jet breakup model is derived using concepts from probability theory. Velocity ﬂuctu-ations at the free surface are hypothesized to be the cause of turbulent jet breakup. This idea is formalized by treating the ﬂuctuations as random variables, subject to damping from the free surface. In contrast to previous theories, this theory uses a conditional ensemble average to determine quantities of interest because not all ﬂuctuations produce droplets. An energy balance and a closure model are used to determine the Sauter mean diameter. Similar approaches are used to determine the breakup onset location, breakup length, and spray angle. To validate the model, data from previous experimental studies using long pipe nozzles was compiled. Data for rough pipes was used to include turbulence intensity in the study.


INTRODUCTION
The main goal of this paper is to develop a model for the breakup of a statistically steady high Weber number turbulent liquid jet injected into a low density quiescent environment.The model is intended to apply when the breakup is caused primarily by the turbulent velocity fluctuations at the free surface, i.e., in the turbulent surface breakup regime, commonly but inaccurately called the "second wind-induced regime" (Trettel, 2020a).I call this model the conditional damped random surface velocity (CDRSV) model.
The breakup of liquid jets under these conditions appears in many applications like fire protection and fuel sprays.Solving this relatively simpler case is necessary to solve more complex cases involving additional physics.I neglect the effects of the transition to turbulence on the jet (as in the case under study the jet is turbulent at the nozzle outlet -see Trettel (2020a) for a discussion of these effects), velocity profile relaxation or boundary layer effects, cavitation, swirl, and any aerodynamic (gas-phase) effects (e.g., the atomization regime, or gas co-or cross-flows).Cavitation and aerodynamic effects are more important in fuel spray problems.
While there have been many previous models of turbulent jet breakup (to be reviewed later in this paper), none of the lower computational cost models are regarded as truly predictive.This is caused primarily by the assumptions and approximations inherent in the models.One important but neglected aspect of the validation of these models is the effect of turbulence, particularly through a measure of the strength of the turbulence like the turbulence intensity.As will be discussed, it is well accepted that the turbulence intensity of the jet is a major factor in the breakup of the jet.Yet relatively few models even consider the turbulence intensity as a factor.And merely considering the turbulence intensity in a model is insufficient -the model must be validated against experimental data with appreciable turbulence intensity variation.This work uses a relatively new experimental database with appreciable turbulence intensity variation (Tu 0 varies from 4.9% to 12.7%).There are many models which consider turbulence intensity, but in ways which are incorrect, as will be discussed.Many models fail to even capture the qualitative trends as the turbulence intensity varies, much less match the data quantitatively.
The paper is organized into two parts.The first part ( § 2) reviews the two main theoretical approaches used to develop turbulent jet breakup models: stability theories and phenomenological theories.Stability theories have received the bulk of previous attention, but unfortunately rarely consider the turbulence intensity as a variable.Various attempts to consider turbulence intensity in stability theories are reviewed.Ultimately these attempts are seen to be problematic, and consequently stability theories are rejected for the time being, though future developments in stability theories may make them viable.Phenomenological theories are more briefly reviewed, with an emphasis on the commonalities between the published models.Phenomenological theories are determined to be more promising for turbulent jet breakup, motivating CDRSV theory.
The second part ( § 3) of the paper develops CDRSV theory, starting with the average velocity of a droplet as it separates from the jet as a motivating example.The reader interested in this aspect of the work can skip the earlier parts if desired as they are merely review and motivation for CDRSV theory.Velocity fluctuations in the jet are treated as random variables and the statistics of the quantities of interest are calculated using a model of the droplet formation process.This is more rigorous than previous phenomenological approaches which typically assume a certain relationship between the turbulent RMS velocity and a quantity of interest.Models are developed for the Sauter mean diameter of a droplet at formation, average droplet radial velocity at formation, average breakup onset location, average breakup length, and spray angle.These models and others are compared against an experimental database including turbulence intensity variation to provide more challenging validation.

Quantities of interest, independent variables, and nomenclature
Figure 1 shows a slice through the center of a statistically steady ensemble averaged circular liquid jet ejected from left to right into a quiescent gas.The nozzle outlet plane is denoted with 0, e.g., the nozzle outlet diameter is d 0 .The x axis starts at the center of the nozzle outlet plane and is oriented with the jet's bulk velocity (U 0 ).The r axis extends radially from the In this work the Reynolds number is Re ij ≡ U j d j /ν i for location j (0 for nozzle outlet) and fluid i (l for liquid, g for gas), e.g., for the liquid phase at the nozzle outlet Re 0 ≡ U 0 d 0 /ν .The Weber number is We ij ≡ ρ i U 2 j d j /σ, e.g., for the liquid phase at the nozzle outlet We 0 ≡ ρ U 2 0 d 0 /σ.The Ohnesorge number can be found given the Reynolds and Weber numbers: 2 )/2 is the turbulent kinetic energy.For simplicity, plane averaged turbulence intensities are used.The spatially averaged equivalent isotropic (u = v = w ) RMS velocity is defined as u j ≡ (2k j /3) 1/2 , not averaging over u j directly.This allows the turbulence intensity to be defined as Tu , which is advantageous in transport equations -like U 0 can be used in mass conservation equations, k can be used in turbulence transport equations -see Trettel (2018) .More simply, Tu 0 is the plane averaged turbulence intensity at the nozzle outlet.The dissipation (ε j ) model also uses k j .The integral length scale is Λ, and does not have a bar for simplicity as it is not a well known quantity in turbulent jet breakup, so whether an average or local value is taken is of little importance.These simplifications can be relaxed in future works, but for the moment they are deemed acceptable given that free turbulent flows tend to become more homogeneous downstream, and that little information is available on the spatial variation of these quantities.It is well accepted that the turbulent kinetic energy peaks near the edge of the jet.A plane homogeneous approximation will be bad in many situations, however, it is a reasonable first approximation.
The average droplet diameter D ij is a diameter which is representative of the spray in some sense.This can be defined in multiple ways by choosing i and j: where f (D) is the droplet size distribution function.While this definition may seem arbitrary at first, various choices of i and j have different physical meanings.For example, the arithmetic mean diameter is D 10 .The volume or mass mean diameter, D 30 , is the diameter of a droplet which has the same volume or mass as the spray as a whole.Similarly, the Sauter mean diameter, D 32 , is the diameter of a droplet with the same surface area to volume ratio as the spray as a whole.
The average droplet velocity is more clear than the average diameter: it is the average velocity of droplets at a particular location.The average droplet velocity at formation from the jet is termed v d in this work.
The average breakup length x b (typically breakup length for short) is defined as the timeaveraged distance from the nozzle where the diameter of the jet core reduces to zero.This typically is found from images or electrical conductivity measurements, which return consistent results -see Trettel (2020a) and Trettel (2019, pp. 3-4).
The average breakup onset location x i (again, typically breakup onset location for short) is defined as the time-averaged distance from the nozzle where breakup first occurs.This is typically measured using imaging techniques.
As has been highlighted by Reitz and Bracco (1986, p. 234-235), at low speeds (the "Rayleigh regime"), the jet breaks up at a single point (i.e., everywhere in a particular cross section), so x i = x b .At higher speeds these two diverge, so that surface breakup can be initiated at a certain distance but the jet core does not end until much farther downstream.Note that some researchers will call x i the "breakup length".I am using the term "onset location" following the work of Wu and Faeth (1995).
The spray angle θ i is roughly defined as the angle produced by the outermost droplets of the spray region at the breakup onset location.Unfortunately there is no standard precise definition of this quantity despite its ubiquity, and consequently there is a large spread in the existing data for the spray angle.See Trettel (2019, p. 3R) for details.
Note that while figure 1 shows the spray cone and jet core as linearly changing over x, this does not necessarily represent the actual shape of the jet.

Turbulence effects on breakup
A focus of this paper is on the effect of turbulence properties like the turbulence intensity on the breakup of liquid jets.This is motivated by common fire hose nozzle design guidelines emphasizing the negative effect of turbulence on the performance of water jet systems (Hoyt and Taylor, 1985;McCarthy and Molloy, 1974;Oehler, 1958;and Rouse et al., 1952).While this effect is present in other applications as well, e.g., fuel sprays, it does not appear to have received the same level of attention.A few experiments have clearly shown that many quantities of interest are sensitive to the turbulence intensity, e.g., the effect of the increasing turbulence intensity typically decreases droplet size (Bogdanovich, 1948 andDumouchel et al., 2005), increases spray angle (Ervine andFalvey, 1987 andSkrebkov, 1966), and decreases breakup length (Ervine et al., 1980;Kusui, 1969;and McKeogh and Elsawy, 1980).Kusui (1969) also showed a turbulence intensity effect on the transition to the atomization regime.Most quantities of interest show a turbulence intensity dependence in turbulent jet breakup to my knowledge.
As discussed in Trettel (2019, p. 4), the turbulence intensity is the most natural way to measure the strength of the turbulence, but unfortunately the turbulence intensity is rarely estimated in experiments or considered in models.Frequently, the turbulence intensity is assumed to be roughly constant or a function of only the Reynolds number, making it unnecessary in models.The Reynolds number is often seen as a measure of "how turbulent" a flow is, but this is mistaken.Contrary to what most expect, in fully developed smooth pipe flow, Tu decreases as Re 0 increases -see Trettel (2019).Nozzles do not necessarily have monotonically increasing or decreasing trends in the turbulence intensity (Lebedev, 2019, figs. 2-4).For a particular nozzle, the nozzle outlet turbulence intensity (Tu 0 ) is a function of Re 0 and the nozzle inlet turbulence intensity (Tu in ).Different nozzles have different trends, so both variables are needed.

REVIEW OF PREVIOUS MODELS FOR TURBULENT JET BREAKUP
Broadly, there are two approaches for modeling turbulent jet breakup.These approaches are not necessarily mutually exclusive, and can be viewed as two different ways of describing the same phenomena.The first approach analyzes the stability of the Navier-Stokes equations as applied to a liquid jet.This approach is called "stability theory" in this work.The second approach does not start with the Navier-Stokes equations, and instead assumes that turbulent jet breakup is described through a simplified model of the larger physics.Models taking the Trettel second approach are called "phenomenological" in this work.Some models are hybrids of the two approaches, applying stability theory where it works and phenomenological theory where stability theory fails, e.g., Magnotti (2017) and Som and Aggarwal (2010).
Stability theory has dominated jet breakup modeling from the early work of Rayleigh (1878) to the popular KH-RT breakup model (Beale and Reitz, 1999).Stability theory has proved accurate at low Reynolds and Weber numbers, i.e., in the "Rayleigh" regime (Trettel, 2020a), but has not demonstrated accuracy in other regimes.The alternative phenomenological approach has been rising in popularity due to the shortcomings of stability theory in turbulent regimes.However, phenomenological theories so far have failed to accurately predict all relevant scalings, though they can overcome other shortcomings of existing stability theories (to be discussed) and may be preferred for that reason alone.
A third approach to modeling turbulent jet breakup is detailed computation, whether RANS, LES, DNS, or another approach.RANS models like the ELSA model of Demoulin et al. (2007) tend to treat the liquid and gas phases as an Eulerian continuum with the ability to create Lagrangian droplets when certain criteria are met.LES and DNS models involve less modeling, but have much higher computational costs.However, given the high computational costs associated with these approaches, low-order models will continue to be used for the foreseeable future and consequently will be the focus of this work.A review of previous research on low-order models, much of which is presently obscure, follows.Note that this is not a rejection of computationthe results of DNS studies will be used as appropriate to help construct low-order models in a manner similar to experimental results, but DNS itself is not used as a model in this work.

Have any turbulent jet breakup models been truly validated?
Before describing the models and their qualitative features, it is worth describing in general how successful these models are quantitatively.This would show some of the strengths and weaknesses of each model, which will be elaborated on later.Existing turbulent jet breakup models are not validated in general in my view due to inadequacies in the validation data and tests, despite frequent claims to the contrary.While jet breakup in the laminar Rayleigh regime has had good success in modeling, the history of the validation of turbulent jet breakup has many examples of researchers trying to fit a square peg in a round hole.I am not going to claim success until a model fits the data with little or no special pleading.
Consider table 1, which summarizes how well the most popular turbulent jet breakup models and the models developed in this work (CDRSV) perform when compared against the data compilation for the turbulent surface breakup regime used in this work.Also included are comparisons against power law regressions made in Trettel (2020a).This data compilation is described in § 3.1, and is much more difficult for models to match because it has much more data and also has appreciable variation in the turbulence intensity, unlike most data sets used for validation of turbulent jet breakup models.
Each model has a single coefficient taking the calibrated value on the left.The coefficient of determination (R 2 ), a measure of how well a model fits the data, is on the right, and higher R 2 is better.The R 2 values for the power law regressions (last row) should be higher than the R 2 values for the other models because the power law regressions are more general.However, note that because the power law regressions are linear regressions made in log coordinates (e.g, log We 0 ), the R 2 value in true coordinates (e.g., We 0 ) may have a different ordering.This is why some of the power law regression R 2 values are lower than the theory R 2 values in table 1 -they are actually higher in log coordinates.The most popular model, KH-RT, is frequently misapplied to the turbulent surface breakup regime, as this model is designed for cases where aerodynamic influence is strong (the atomization regime).Consequently, the model performs consistently poorly in this regime, having negative R 2 values for all quantities of interest.Similarly, all theoretical models for the spray angle have negative R 2 values, indicating uniformly bad performance on this quantity of interest.The droplet size model of Huh et al. (1998) has a value of R 2 of zero by definition as it proposes that the droplet size is proportional to the integral scale.Given the lack of existing data, is assumed that the integral scale is proportional to the nozzle diameter, so in effect D 32 /d 0 is a constant in the model of Huh et al.Constants have R 2 values of zero by definition, and in this case a constant is not a good model for other reasons, as will be discussed.

Atomization and Sprays
Models for the breakup length tend to have mixed success at best.The Faeth group model has a slightly better R 2 value than that developed in this work, despite the fact that it has no turbulence intensity variation.This is likely because the Faeth group model has a Weber number exponent closer to that of the empirical regression ( x b /d 0 ∝ ∼ We 1/3 0 in the regression).Unfortunately no turbulent jet breakup model I am aware of has the correct scaling of the breakup length with the liquid Weber number.Some models have better success.The droplet size model of the Faeth group (to be discussed) has the same scaling as the model developed in this work and has a respectable R 2 value.Similarly, calibrated breakup length models from the Faeth group and this work also have respectable R 2 values.The breakup onset location model developed in this work performs significantly better than that of the Faeth group, however.

Stability theory
Linear stability theory examines how a linearized system evolves in time based on a constant background flow with small perturbations added to it.This section is an overview of stability theory based jet breakup models without details of the derivation, highlighting typically ignored weaknesses of stability theory as applied to turbulent jet breakup.The focus is on linear stability theories, however, it is worth noting that most of the mentioned issues apply equally well to non-linear stability theories.
Trettel 2.2.1 Popular stability theory based models Reitz andBracco (1982, 1986) note that many popular linear stability theory based models can be understood as simplifications of a general "dispersion relation" for liquid jets: where ω is the "growth rate", a measure of how quickly a disturbance of wavenumber κ grows, r 0 ≡ d 0 /2 is the nozzle outlet radius, I n are nth order modified Bessel function of the first kind, K n are nth order Bessel functions of the second kind, the primes indicate differentiation, and Equation 2 can be simplified into various limiting forms for different scenarios and solved analytically to obtain the growth rate ω for a particular wavenumber κ, or the full equation could be solved numerically.Traditional linear stability theories then assume that the most unstable mode (i.e., largest ω, denoted ω m ) is responsible for all breakup of the jet.
Rayleigh (1878) developed the first stability analysis of liquid jets.The analysis assumed the jet was inviscid jet and in a vacuum.Rayleigh hypothesized that the most unstable mode would dominate over all others, and numerically calculated the wavelength (λ ≡ 2π/κ) of this mode as λ m = 4.51d 0 . (3) If a cylinder of length λ m with diameter d 0 forms a spherical droplet, it has a diameter of (Lefebvre and McDonell, 2017, p. 25L) Similarly, the breakup length can be computed in Rayleigh's theory.A brief derivation of the breakup length here can be instructive, following McCarthy and Molloy (1974, p. 2R).The amplitude of the perturbation of the radial location of the free surface of the jet traditionally evolves according to so the radial location of the free surface is where κ m is the wavenumber associated with λ m , δ is the amplitude of the free surface perturbation, and δ 0 is "initial disturbance level".The initial disturbance level is simply how disturbed the free surface is at time 0 at x = 0.The breakup length can be found by locating where the jet diameter decreases to zero.This occurs when the disturbance amplitude grows to d 0 /2.Consequently, in Rayleigh's theory Atomization and Sprays which can be used to find an equation for the breakup time t b , The breakup time can be converted into a length through the convective velocity U 0 to find t b ≡ x b /U 0 .Now equation 8 can be rearranged to state Weber (2019, year of translation) in 1931 developed a theory for low-speed viscous jets.In this theory the most unstable wavelength is (Lefebvre andMcDonell, 2017 andWeber, 2019) where the Ohnesorge number is Oh 0 ≡ µ / √ ρ σd 0 .The corresponding droplet diameter is found to be by, again, setting as equivalent the volume of a droplet of diameter D and a cylinder of length λ m with diameter d 0 .
The breakup length in Weber's theory is (McCarthy andMolloy, 1974 andWeber, 2019) Experience has shown that Weber's theory is reasonable for low Weber number jets in general.
Another popular early theory was developed by Taylor (1958) in 1940 and later extended by Ranz (1958) and Ranz and Dreier (1964).Reitz (1978, p. 134) calculated a droplet diameter consistent with Taylor's theory: where B 1 is a model constant and λ * m is a dimensionless wavelength for the most unstable mode, defined as 1, λ * m asymptotically approaches 3/2.Similarly, for the breakup length in the same limit, Taylor's theory states that (Reitz, 1978, p. 135) x where B 4 is a model constant.Equation 15 has proved popular in predicting the breakup length in the atomization regime, though my own regression analysis detailed in Trettel (2020a) suggests Trettel the scaling with ρ /ρ g may be incorrect.Also, the coefficient is known to take wildly varying values, suggesting that the model is not predictive, as discussed in the next section.
In the context of Taylor's theory, Ranz and Dreier (1964, p. 59R) suggested that the spray angle can be modeled as This result applies for linear stability theory in general and can be used with other models.
In 1987, rather than analytically solving limiting cases of the dispersion relation, Reitz (1987) solved the dispersion relation numerically and fitted algebraic equations to it.Changing the notation to match that used in this work, the curve fits for the maximum growth rate ω m and the associated wavelength λ m are These equations fit the numerical solutions well in the ranges 0 ≤ Oh 0 ≤ √ 2/2, ρ /ρ g > 10, and 0 ≤ We g0 < 2000.
Equation 17 is used as a component of the popular KH-RT model (Beale and Reitz, 1999), where the droplet size of an atomizing liquid "blob" is calculated as where B 0 is a constant (presumably universal as the wavelengths should be directly proportional to the droplet size in this theory) which is 0.61.While in principle a breakup length can be calculated with ω m from equation 18, in practice the Ohnesorge number is assumed to be zero and the gas phase Weber number is assumed to be high, leading to simplified expressions for λ m and ω m .Under these approximations Beale (1999, pp. 94-95) finds that which is identical to the result for Taylor's theory, equation 15

Shortcomings of stability theory
Often, the flow exiting a nozzle is turbulent, yet stability theories tend to assume that the flow is initially laminar with only small free surface perturbations (no velocity perturbations).Turbulence influences the breakup of the jet, but how should turbulence properties be included in stability analyses?
Atomization and Sprays

Inclusion of turbulence properties via model coefficients
In stability theory, various model coefficients can vary greatly.For example, in the previously mentioned KH-RT model, one of the coefficients, B 1 , takes calibrated values ranging from 1.73 to 40 (Ning, 2007, p. 14).This has made some researchers conclude that the KH-RT model is not predictive (Magnotti and Genzale, 2017, p. 34L).A widely varying model coefficient suggests that there may be missing variables in the theory.
To be more specific, stability theories offer no clear mechanism by which Tu 0 can be incorporated into the theory beyond the "initial disturbance level" or another similar model coefficient.The physical meaning of this term is unfortunately vague, and how it scales with Tu 0 would at least require additional theoretical work.This fact was unfortunately realized slowly.Phinney (1972) noted that earlier research tended to assume that the initial disturbance level was constant, and suggested instead that the initial disturbance level be considered as a variable.Chen and Davis (1964, p. 191) recommended a study of how the initial disturbance level relates to turbulence quantities at the nozzle outlet, though no such study was ever undertaken to my knowledge.Lin and Reitz (1998, p. 92) note that the initial disturbance level depends on the nozzle geometry, as diesel and water jet cutting systems operate at similar pressures, yet the breakup length for diesel jets is short while that for water jet cutting jets is long.When studying diesel sprays, Reitz and Bracco (1982) discuss how the initial disturbance level changes based on several different nozzle geometries tested.In an earlier work Reitz (1978, p. 133) gives an empirical regression for a model coefficient, A, as a function of the nozzle aspect ratio (L 0 /d 0 ): This model coefficient is higher for more disturbed flow.This equation is only valid for L 0 /d 0 < 10, as the amount of disturbances to the flow should saturate as the flow develops, but the regression suggests that it continues to increase regardless of how developed the flow is.
Assuming that stability analysis still works when the flow is turbulent at the nozzle outlet, an energy flux argument similar to that of Moallemi et al. (2016) can be used to suggest that in Weber's theory the initial disturbance level is a function of the turbulence intensity.The derivation is omitted for brevity but can be found in Trettel (2020b).The result is that the energy-equivalent initial disturbance level is Unfortunately, this equation is a poor fit to the data in the turbulent Rayleigh regime (where it presumably applies as Weber's theory is for the Rayleigh regime), over-predicting the breakup length.
Trettel (2020a) develops a model for an equivalent initial disturbance level which performs better in the turbulent Rayleigh regime than equation 22: where arccsch is the inverse hyperbolic cosecant function.Equation 23 with C v = 0.0615 is an excellent fit to the available data (31 data points, R 2 = 0.961), indicating that modeling the coefficients in stability theory can be viable.However, there are no stability theory coefficient models with similar success in the turbulent surface breakup or atomization regimes.

Trettel
Additionally, changing the initial disturbance level alone typically can not change the characteristic droplet sizes.The droplet size in linear stability theories is independent of the initial disturbance level -see equation 19.So even if the initial disturbance level increases as turbulence intensity increases, the decrease in droplet size observed experimentally as turbulence intensity increases (Bogdanovich, 1948 andDumouchel et al., 2005) can not be reproduced.In non-linear stability theories the wavelength of the most unstable mode (and consequently, the droplet size) depends on the initial disturbance level, possibly remedying this situation.The only study I am aware of to examine the effect of the initial disturbance level on the wavelength of the most unstable mode is by Wang (1968, p. 312, eqn. 91).Using equation 23, it can be shown that Wang's correction term causes the droplet size to decrease in the turbulent Rayleigh regime.It is difficult to evaluate the accuracy of Wang's theory as the small parameter used in the perturbation analysis is not clearly defined.In the Rayleigh regime in general, the droplet size appears to be insensitive to the initial disturbance level (Trettel, 2020a), though possibly Wang's correction is negligible.In other turbulent regimes the non-linear correction may not be negligible.
As it turns out, modeling the coefficients alone is sufficient for the turbulent Rayleigh regime likely because in this regime the most unstable mode approximation (to be discussed) is valid, the droplet size appears to be roughly independent of the initial disturbance level, and the effect of the turbulence is felt by the jet primarily near the nozzle exit where the initial disturbance level is set (Trettel, 2020a).These characteristics do not extend to other regimes.
In other turbulent regimes, the conclusion is that model parameters like the initial disturbance level are largely empirical at present.If a change is made to a nozzle design, it is typically not possible to estimate the effect on the model coefficients.This contrasts strongly with the phenomenological approach, where a nozzle model can easily be used as an input to a jet breakup model due to the turbulence intensity being explicitly considered.

Inclusion of turbulence properties via uncertainty propagation
Several previous researchers have proposed that existing stability analyses can be used if fluctuations in input quantities (e.g., turbulent velocity fluctuations) are used to generate distributions of output quantities like the breakup length (Lafrance, 1977) or droplet size (Babinsky andSojka, 2002 andSovani et al., 1999).This approach is promising, however, existing implementations have major flaws which must be addressed.
From a model validation perspective, past researchers using this approach have not characterized the input quantities (turbulent velocity fluctuations, again, in our case).This left past researchers unable to evaluate the success of these models more than qualitatively, which has been noted by Babinsky and Sojka (2002, pp. 326-327), proponents of this approach.Fortunately, today we have the information to properly evaluate these models.To use Lafrance (1977) as an example, Lafrance use the stability analysis of Rayleigh and a random nondimensional initial disturbance level to calibrate the RMS nondimensional initial disturbance level to data from Phinney (1973).Lafrance suggests this RMS nondimensional initial disturbance level can be interpreted as a RMS turbulent velocity (presumably a turbulence intensity).The value found, 0.8%, was called "reasonable" by Lafrance, but is actually roughly an order of magnitude too low -see Trettel (2019).Mere assertion is not enough; comparison with unambiguous experimental data is absolutely necessary.The failure of the approach of Lafrance is due to the shortcomings of the model used, not the idea that fluctuations need to be taken into account.Similar criticism also applies to the approach taken by Sovani et al. (1999).
The most significant shortcoming of this approach is that it ignores the fact that velocity fluctuations can appreciably change the mean quantities and do not mainly widen the distributions.Higher turbulence levels typically decrease the breakup length, however, this is not reflected in Lafrance's analysis.Lafrance treated the initial disturbance level and turbulence intensity as independent, when in reality the two are coupled.The mean breakup length in Lafrance's model does not decrease as the turbulence level increases, in contradiction to experimental evidence.The same problem is also demonstrated through characteristic droplet sizes in typical linear stability theories being independent of the initial disturbance level as previously mentioned.In the droplet size case, making the initial disturbance level fluctuate would have no influence at all on the droplet size.The fluctuations are more than just uncertainties to be propagated.Schmid (2007, p. 149) states that stochasticity has two effects on stability in general: 1. the response of a deterministic system disturbed by an external stochastic process, and 2. the modification of the system when it is disturbed by an external stochastic process.
Existing models consider only the first effect, but the second effect is very important in turbulent jet breakup as well.To my knowledge at present no stability-theory-based turbulent jet breakup model (outside of the turbulent Rayleigh regime) considers the second effect.How to incorporate the second effect in stability theory for turbulent jet breakup in general is an open problem.It is possible that combination with an initial disturbance model like equation 23 will partly account for the second effect.However, as the droplet size in linear stability theory is independent of the initial disturbance level, this approach will only be partial.As previously mentioned, a nonlinear theory considering the effect of the initial disturbance level on the breakup process may be necessary to take into account the second effect.Still, at present it is not clear if a non-linear theory is sufficient to consider the second effect.

Inclusion of turbulence properties via using a turbulent viscosity or by analyzing the stability of the RANS equations
Chen and Davis (1964, pp. 190-191) hypothesized that Weber's theory could be used if a turbulent viscosity were used in place of a molecular viscosity.Unfortunately, Chen and Davis found that this assumption did not agree with their experimental data.Use of a turbulent viscosity instead of a molecular viscosity in stability analysis seems unlikely to account for turbulence properties in turbulent jet breakup.Turbulent viscosities are much higher than molecular viscosities.Increased viscosity tends to stabilize liquid jets according to Weber's theory (see the Reynolds number effect in equation 12), but in experiments increased turbulence levels, which would increase turbulent viscosity, tend to destabilize liquid jets.The physical reasons for using a turbulent viscosity instead of the molecular viscosity are not clear.If stability analysis is performed on the RANS equations with a turbulent viscosity model, then neglecting the molecular viscosity, the results will be identical to classical stability theory, albeit with a turbulent viscosity instead of a molecular viscosity.Indeed, Sauerwein (2020) analyzed the stability of the RANS equations in a case similar to that of Rayleigh without using a turbulent viscosity model and obtained similar incorrect results: higher turbulence levels led to more stable jets.
The simplest explanation for the failure is that the breakup stability of the RANS equations is not the same as the breakup stability of a turbulent flow.The perturbations that lead to breakup, even in the Rayleigh case, would seem to be averaged out in a RANS stability approach, leading to a more stable jet.

Assumption that the most unstable mode dominates
The vast majority of previous stability analyses assume that a single wavelength dominates the breakup process.This seems implausible in a broadband phenomena like turbulence.Measurements of surface waves on high Weber number turbulent liquid jets confirm that the breakup process is broadband (Appel andSwenson, 1968 andChen andDavis, 1964).Similarly, droplet size distributions in high Weber number turbulent jet breakup tend to be highly polydisperse, not essentially monodisperse as would be expected if a single wavelength dominated.More recently, Agarwal and Trujillo (2018, pp. 11R-12L) note that in a configuration like typical fuel sprays, the most unstable mode acts primarily to cause breakup of the surface, but the jet core itself is broken up by a different mode downstream.In the laminar and turbulent Rayleigh regimes, the most unstable mode dominance hypothesis appears reasonable.But its success in those regimes does not imply this hypothesis is valid in other regimes.
Consideration of more than one wavelength is necessary to model high Weber number turbulent jet breakup.To my knowledge there have been few attempts to consider more than one disturbance wavelength.Lemberskii and Ferber (1973) propose using the dispersion relationship to calculate a droplet size probability density function.Lemberskii and Ferber assume that the initial disturbance level is independent of the wavenumber.Yi and Reitz (2004) later independently developed a computational model considering multiple different initial disturbances.This is an essentially empirical way to determine the initial disturbance level given a particular turbulent disturbance.Unfortunately, Yi and Reitz (2004, eqn. 22) estimated the nozzle turbulent kinetic energy with a nozzle turbulence model that I showed to be extremely inaccurate (Trettel, 2018).Yi and Reitz also did not consider the possibility that the stochasticity changes the dispersion relationship, like all other precious stochastic models I am aware of.The works of Lemberskii and Ferber and Yi and Reitz are pioneering, but ultimately too flawed to justify using as-is.

Other assumptions made in stability analyses
Agarwal and Trujillo (2018) criticize some of the assumptions involved in linear stability theory.In particular, the DNS results of Agarwal and Trujillo show that at high Weber numbers, the nonlinear term neglected in the analysis can have significant effects more than 4 nozzle diameters downstream of the nozzle outlet.Agarwal and Trujillo also showed that the Fourier decomposition used in stability analysis quickly becomes a poor approximation because the free-surface is multiple-valued.
Additionally, the vast majority of existing linear stability theories assume that the disturbances are axisymmetric.This assumption is questionable at best -Yang (1992, p. 681) notes that if axisymmetry were correct then rings would break off from the jet, not droplets.Still, Lin and Reitz (1998, p. 101) maintain that axisymmetric disturbances dominate, noting that the growth rate of the axisymmetric disturbances is typically higher than the growth rate for asymmetric disturbances.Lin and Reitz suggest that breakup of the rings into droplets may explain the observation of droplets, but again, droplets are observed to be formed directly from the jet, not rings, so this explanation is implausible.The simplest explanation for why asymmetric disturbances appear in turbulent jet breakup in practice is that the turbulence spectra leads to asymmetric initial disturbances.

Overall assessment of stability theory for turbulent jet breakup
Stability-theory-based models have so far produced results which are implausible with respect to the effect of turbulence intensity.Adding turbulence intensity to stability theories is not trivial.Simple fixes can not solve this problem -substantially new models are needed but do not appear forthcoming.While there are several plausible avenues for research in stability theory, ultimately the alternative phenomenological approach was deemed more promising for developing simple models for turbulent jet breakup as this approach makes avoiding the shortcomings of stability theory easier.

Phenomenological models
In contrast to stability theory models, "phenomenological models" are not based on the Navier-Stokes equations, and instead assume a certain simplified physical process is occurring.I borrow the word "phenomenological" in this context from Wu et al. (1992) and others in the Faeth research group.These models may be developed in mechanistic ways, considering for example an energy or force balance which leads to droplet formation provided certain criteria are met.They also could be developed in heuristic ways like scaling arguments.
It is tempting to think that stability-theory-based models are inherently more credible than phenomenological models, as they are based on the Navier-Stokes equations, which are considered reliable.However, the assumptions and approximations inherent in present stability-theorybased models make these models arguably no more credible than phenomenological models.As has been discussed, these assumptions and approximations are not easily avoided.Further, phenomenological models do not necessarily abandon first principles.Phenomenological models instead choose simplified first principles.The inaccuracy, like in stability theory, comes from the assumptions and approximations.
There are many phenomenological models in the literature.To reduce the scope of this review, only models with turbulence intensity dependence are considered.The turbulence intensity can be more easily incorporated into a phenomenological model than a stability theory model due to the flexibility of construction of a phenomenological model.However, this review should not give the impression that most phenomenological models consider turbulence intensity, as few do.Additionally, in contrast to the stability theory section, some brief derivations of these models will be given as appropriate to give the reader an understanding of how these models work.For detailed derivations, the reader is referred to the original works like before.
It seems obvious that turbulent fluctuations normal to the free surface, i.e., in the radial direction, can perforate the surface and form droplets.An illustration of this process is shown in figure 2. Breakup through this mechanism has been observed in the photographic experiments of Hoyt et al. (1988, pp. 359-360) and Wu et al. (1992), and also the DNS study of Desjardins and Pitsch (2010).Natanzon (2018) developed the earliest quantitative theory of turbulent jet breakup I am aware of in 1938.Natanzon applied the maximum entropy principle with a kinetic energy constraint using k to find the droplet diameter (number) distribution and the mass mean diameter which can be expressed in dimensionless variables: Bogdanovich (1948, pp. 122-123) developed an energy balance argument for an (unspecified) average droplet diameter.This argument proved popular in Eastern Bloc countries during the Cold War (Lebedev, 2019 andSitkei, 1963), and similar arguments (both inspired by and independent of Bogdanovich) were also made by researchers outside of the Eastern Bloc (Fritzsche, 1965 andInoue, 1963).The simplest form of this argument examines the formation of a single spherical droplet of volume -V = π 6 D 3 and surface area SA = πD 2 at the free surface.When inside of the jet, the "droplet" has turbulent kinetic energy k 0 (i.e., the same value as at the nozzle, plane averaged).All of this turbulent energy is consequently converted into surface energy, i.e.: Then, substituting in the definitions of the volume and surface area of a spherical droplet, the droplet diameter obtained is or in dimensionless variables which is the same scaling Natanzon found for D 30 .Lebedev (2019) changes the constant to a model coefficient.Sitkei (1963) modifies the argument to add a term for viscous dissipation.This class of arguments is flawed, however, as the eddies interacting with the free surface do not all have the same kinetic energy as represented by k 0 .Wu et al. (1992) avoided this issue through an energy balance argument combined with inertial range scaling to estimate D 32 for the initial droplets, that is the Sauter mean diameter of the first droplets formed downstream of the nozzle.The difference in this argument is that instead of using the full turbulent kinetic energy k 0 , the turbulence spectrum is introduced as the energy appropriate for an "eddy" of size which is assumed to be in the inertial range.Wu et al. convert the isotropic energy spectrum into an equivalent velocity (see § 3.4 for details).Then, merely assuming proportionality rather than strict equality, Wu et al. state that the surface energy is proportional to the eddy energy: and apply the result below which follows from the inertial range spectrum where Λ 0 is the turbulent integral scale at the nozzle outlet.Wu et al. then assume that D 32 ∝ to obtain the following scaling: where I am using a different notation than Wu et al. which has changed the form of the equation.The use of a spectral energy rather than the full turbulent kinetic energy appreciably changes the form of the droplet size equation.
Broadly, the simple arguments by Bogdanovich and Wu et al. have three problems: 1. the arguments assume that to obtain an average output (e.g., D 32 ), one can replace input quantities with their averages (or a "representative" value).However, this is only true for linear equations and is often false for non-linear equations (e.g., the RANS closure problem), 2. the arguments do not justify which characteristic diameter (D ij ) is appropriate, and 3. the arguments assume that all turbulent surface fluctuations will form droplets.
If strict energy conservation is followed, these arguments assume that the droplets have zero velocity at formation because all input energy (e.g., turbulent kinetic energy) is used to create new free surface.The use of empirical coefficients instead of those implied by strict energy conservation by Lebedev (2019) and Wu et al. (1992) avoids this issue, though I am not sure that these researchers were aware of this.Ultimately, it would be better to follow strict energy conservation and model the distribution of turbulent energy into surface energy and droplet kinetic energy as done later in this work.
Fortunately, Natanzon's theory does not have the first two problems or the droplet velocity problem, as the droplets are assigned velocities through maximum entropy as well, and averages are computed directly from a distribution function rather than assumed to be of a certain form.However, as before, having a model for droplet velocity would be preferably to merely using the maximum entropy principle to assign droplet velocities.
Recently, Schmitz (2011) developed a computational extension of the ideas of Wu et al. (1992) with a more general turbulence spectrum.This allows computation of the full droplet size distribution, avoiding the issue of which characteristic droplet size the model predicts.

Trettel
The computational KH-RT model (Beale and Reitz, 1999) mentioned in the stability theory section has seen many improvements over the years.One notable recent phenomenological improvement by Magnotti (2017) to the model has been the incorporation of a droplet diameter model from Wu et al. (1992, eqn. 14).The model is a hybrid of the stability and phenomenological approaches, picking the approach that works best in different regimes.Ignoring cavitation and droplet/secondary breakup, Magnotti's model essentially switches between the stability theory "KH" model for the atomization regime to the phenomenological model by Wu et al. outside of the atomization regime.One major criticism of Magnotti's model is that it does not consider the effect of turbulence intensity as currently implemented.
Another computational approach is one-dimensional turbulence (Movaghar et al., 2017), abbreviated ODT.ODT is based on random mixing events in a simplified one-dimensional (transverse) domain that moves at the local (constant) convection velocity U 0 .When the mixing event causes liquid to no longer be attached to the core, breakup occurs.The accuracy of ODT for turbulent jet breakup has not yet been fully demonstrated to my knowledge; the model appears capable of predicting jet breakup lengths reasonably well, though droplet size appears to have some issues downstream, possibly due to a change in the mechanism of droplet formation from surface breakup to column (Rayleigh) breakup.I also am concerned that the spatial grid resolution would impact the estimation of the smallest droplets.ODT is also considerably more complex and computationally expensive than the model proposed by Schmitz (2011).
The CDRSV model developed in this work is based on energy balance ideas, like in many previous models, however, averages are computed explicitly to make clear which characteristic droplet size is implied by the argument.This model also considers the fact that not all turbulent surface fluctuations will lead to droplet formation, which is essential for the prediction of quantities like the spray angle and droplet radial velocity.
Separate from energy balance arguments, Huh et al. (1998) and Skrebkov (1966) assume that a characteristic droplet diameter is proportional to the integral scale Λ: This does not avoid the previously mentioned problems.And, as the integral scale is the largest turbulent scale, this selection seems plausible for only the largest droplet diameter.Still, the model of Huh et al. has proved popular, being developed further into in part of the hybrid model of Som and Aggarwal (2010).
The previous discussion in this section focused mainly on droplet size, which is the primary concern of most previous researchers, however, other quantities of interest may also be computed with phenomenological approaches.Natanzon (2018, p. 6) assumed all breakup occurs at the nozzle outlet, so x b = 0. Wu and Faeth (1995, p. 2916R) assumed the jet core ends where the local Sauter mean droplet diameter increases to the local jet diameter, returning which is identical to Weber's theory at high Reynolds number (equation 12).As discussed in Trettel (2020a), the We 1/2 0 scaling is not appropriate for liquid jets in the turbulent surface breakup regime or the atomization regime.This model also has no dependence on the turbulence intensity, making it implausible as-is.
The approach of Wu and Faeth does not directly predict where the jet core ends.In contrast, I take a more direct approach in this work by estimating the surface mass flux to find where the jet core ends on average.
The experiments of Wu et al. (1992, p. 305) suggest that v d may scale with the radial turbulent RMS velocity v .Natanzon (2018) and Huh et al. (1998) developed spray angle θ i models with this assumption.In the limit of high ρ /ρ g , the model of Huh et al. is where C θ i is a model constant.Rather than assuming this scaling, Skrebkov (1966) used a energy balance including v to determine v d and θ i .I use a force balance in this work.Kerstein et al. (2017) and Wu et al. (1992) estimate the breakup onset location x i by finding the time required for breakup to occur, but they estimate this time differently.Wu et al. (1992, eqn. 10) find that Kerstein et al. (2017, eqn.2.5) instead find that Kerstein et al. is more consistent with the data.My model in this work is similar to that of Kerstein et al., however, it has a more detailed justification and replaces the Re 0 dependence with a Tu 0 dependence.

Experimental data compilation
An experimental database was compiled for all quantities of interest for model validation.The database is limited to non-cavitating liquid jets ejected from long pipes ("pipe jets") at low Mach numbers (< 0.4).Excluding non-pipe jets reduces the impact of factors which are typically unknown but roughly constant for fully developed pipes, e.g., the velocity profile and the integral scale.Using solely pipe nozzles also allows us to estimate Tu 0 .For that purpose, I developed a regression between the friction factor and Tu 0 (≡ (2k 0 /3) 1/2 /U 0 ) for fully developed pipe flows: Tu 0 = 0.366f 0.459 (9 smooth and 8 rough points, R 2 = 0.975) (Trettel, 2018, p. 6).As Tu 0 in smooth fully developed pipe flows is a function of only Re 0 , rough pipes are needed to avoid confounding between Re 0 and Tu 0 .Unfortunately I am aware of only two rough pipe jet breakup studies.Skrebkov (1966) has 3 measurements of θ i and Kusui (1969) has over 150 measurements of x b .Kusui had a 8.75d 0 smooth section after their rough pipe, complicating estimating Tu 0 .Presumably Tu decays in the smooth section.However, a power law regression of pipe data for x b fits non-pipe data (Ervine et al., 1980 andMcKeogh andElsawy, 1980) best with no decay -see Trettel (2020a).While I assume there is no decay in Tu 0 in this work, ultimately, Kusui's data is imprecise about how x b varies with Tu 0 .
Only data in the "turbulent surface breakup" regime is used in this work.In this regime the breakup is caused by turbulence with negligible aerodynamic effects (Trettel, 2020a).
FIG. 3: A Gaussian velocity probability density function, with the fluctuations that can lead to breakup shaded.The vertical dashed line is the average velocity fluctuation which causes breakup, which in the simplified model is the same as the average radial droplet velocity, v d .

Turbulence evolution in liquid jets
Understanding how turbulence quantities (k and Λ) evolve spatially in the jet is required to develop models of the breakup of the entire jet.Kim (1983, p. 23) and Huh et al. (1998, p. 458) used turbulence models to estimate the decay of turbulence in the jet.Experiments show that turbulence in a liquid jet decays at the centerline initially (Mansour and Chigier, 1994, p. 3390).However, shear at the jet surface causes production of turbulence, such that k can increase.As droplets are formed at the free surface, using solely decay is not necessarily correct if production is significant.The measurements of Mansour and Chigier (1994, p. 3389) suggest that k at the jet boundary grows slowly downstream.This is inconsistent with the measurements of Wolf et al. (1995, p. 402L), which suggest that k only decays at the boundary.Given the complexity of turbulence modeling, I will use the approximation of Wu et al. (1992, p. 308): k and Λ do not vary downstream.The turbulence will also be approximated as homogeneous in the radial and angular directions and isotropic.Spatial averaged k will approximate the k profile.In reality, k peaks near the free surface, becoming more homogeneous downstream.

Simplified example of CDRSV theory
Jet breakup models often assume that the average radial droplet velocity at formation is proportional to the turbulent RMS velocity, v d ∝ v (Huh et al., 1998;Natanzon, 2018;and Wu et al., 1992).CDRSV theory can easily mathematically derive this result and predict the constant of proportionality.
Consider an eddy with radial velocity v approaching the free surface as seen in figure 2. Additionally assume that the free surface returns to a straight line once all interactions with eddies are complete and that no more than one eddy interacts with the free surface at any time.The velocity fluctuations at the free surface are assumed to be accurately described by a Gaussian probability density function (Gaussian PDF for short).
In this example, the free surface presents no obstacle to the eddy, i.e., the surface tension is zero and there is no "damping".So v d = v.Or, in words, the droplet velocity equals the fluctuation velocity.A droplet is formed if the droplet velocity is greater than zero, i.e., v d > 0, Atomization and Sprays so if v > v min = 0 in this case -v min will not be zero in the more complex cases considered later in this paper.I abbreviate this condition as DF (droplet formation).The term DF will be dropped if redundant.
Without conditioning, v d = 0.With the v > 0 condition, the average radial droplet velocity has to be greater than zero.Graphically, the effect of the conditioning can be see in figure 3. The shaded portion of the Gaussian PDF is the part conditioned on.The average of the conditioned portion is given by the dashed vertical line.Obviously, by excluding all fluctuations in the negative direction, the droplet velocity must move to greater than zero, as the math shows.This basic procedure can be applied to calculate various quantities of interest in turbulent jet breakup, not just the droplet radial velocity.

Droplet radial velocity v d for a particular eddy and the Hinze scales -the damping
A model of the droplet formation process considering surface tension is needed.Consider a random turbulent velocity fluctuation v (mean zero) at the free surface at time 0 (so ṽ(t = 0) = v).A droplet forms if the radial velocity ṽ(t) > 0 when a droplet detachment condition is met.Surface tension opposes/damps the turbulent fluctuations.This force F σ = A • p σ where A is the cross-sectional area of the surface perturbation and p σ = 2σ/R is the capillary pressure, where R is the radius of curvature.I assume that the surface perturbations are spherical, with a radius of curvature R equal to the distance δ the eddy penetrates outside the free surface (see figure 4).
Multiplying by an arbitrary constant, I find that F σ = 2πC F σδ.I assume that the eddy has a diameter proportional to ≡ 2π/κ, where κ is the wavenumber of the turbulence associated with the velocity fluctuation v. (Note that despite the eddy's nominal diameter being , I select the radius of curvature as δ for simplicity.)The eddy's mass then is C - V ρ π 3 /6, with another arbitrary constant.The equations of motion of the eddy as it penetrates the surface are dδ dt = ṽ and which have the solutions If it is assumed that the droplet detaches after traveling a distance δ = C lig (C lig 2, so that detachment occurs when the lower end of the ligament is beyond the original free surface location), then the breakup time t b can be found and the droplet velocity at detachment (v d = ṽ(t = t b )) is: This model is oversimplified, but it has the desired features.The last term is an inverse eddy Weber number, As such, minimum scales for droplet formation exist.An arbitrary eddy velocity v can be related to a corresponding eddy wavenumber with v = κE(κ) (Hinze, 1975, p. 222).If I assume that the minimum scales are in the inertial range and apply v = κE(κ) to the inertial range spectrum E(κ) = C K ε 2/3 κ −5/3 , I find that = 2πv 3 /(C 1/2 K ε).Note that C K = 1.5 is recommended by Pope (2000, p. 231).From there I can calculate the Hinze scales (Hinze, 1955 andKolmogorov, 1991), the smallest for which droplet formation can occur (v d = 0): The velocity v σ is the minimum for droplets to form if surface tension dominates.At high Weber numbers, v σ may decrease below v K , the Kolmogorov velocity scale (Pope, 2000, p. 185), and in that case v K will be the minimum.I use the term v min for whichever minimum applies.Because not all fluctuations produce droplets, the ensemble averages I calculate will be conditioned on droplet formation, abbreviated DF.The condition notation will be dropped for terms which imply breakup occurs, e.g., v d | DF would be redundant.Additionally, I'll use v min = v σ for simplicity in this paper.Analogous expressions for v min = v K are easily found.I am unaware of data capable of validating the minimum droplet velocity and diameter estimates.The smallest droplet observed by Wu (1983, p. 36) was 3 µm in diameter (< 0.5 µm uncertainty, ρ /ρ g < 40, likely in the atomization regime), but insufficient detail was provided to estimate σ or K , the Kolmogorov length scale (Pope, 2000, p. 185), for this case.The smallest droplets measured by Wu et al. (1992, p. 307) (ρ /ρ g > 500) were said to be much larger than K in the turbulent surface breakup regime.The DNS study of McCaslin and Desjardins (2015, p. 5, fig. 2b) suggests that surface perturbations are suppressed for scales smaller than σ if σ > K .Also, note that while the DNS study of Ling et al. (2019) suggests that the Hinze scale is larger than the smallest observed droplets, this study is not in the turbulent surface breakup regime studied here, as that work has very strong shear that is absent in the problem studied in this work.
3.5 Sauter mean diameter, D 32 , and average droplet radial velocity, v d Wu et al. (1992, p. 312) assume that D 32 scales with a representative length (in the terminology used here, D 32 ∝ | DF ), however, I will not assume this.D 32 is controlled by the surface energy, not the size of the eddies directly.Energy conservation suggests (assuming the process Atomization and Sprays is adiabatic and neglecting rotational and other energies): or simplified: where in the first equation the left side is before breakup and right side is after breakup.I assumed that only one droplet is formed per eddy event.The eddy/droplet has volume -V , and the formed droplet has surface area SA.The model also implicitly assumes that v d and droplet diameter D are perfectly correlated.For simplicity I assume that the r direction is always normal to the liquid surface, accurate for large x b /d 0 .Like v (= ṽ(t = 0) as before), u and w are turbulent velocity fluctuations with mean zero defined in the streamwise and angular directions, respectively.The mean velocities in the radial and angular directions are zero.(If the jet is decelerating, there is a mean V as well, however, I neglect this as I assume aerodynamic drag is negligible.)I assume that the free surface does not affect streamwise or angular velocities such that u d ≡ U 0 + u and w d ≡ w.These cancel with the input energy, leaving the surface area to volume ratio to be determined by the energy left over from the damping.Now, I apply the conditional average and the model for v d (equation 40), and note that by hypothesis SA/-V ≈ SA / -V = 6/D 32 : which returns Contrary to what one might expect, D 32 is proportional to the harmonic mean −1 DF −1 , not the arithmetic mean | DF .The two terms are the same to first-order, but not identical.This term is unclosed, so it requires a model.The concept of an "eddy" in this work will be clarified.The length associated with a particular velocity fluctuation v is ambiguous.The energy spectrum as used by Wu andFaeth (1995, p. 2916) can relate v and , but this is only a heuristic.More than one "eddy" can contribute to velocity fluctuations at a particular location.Smaller lengths likely have only one eddy contribution, making the idea behind the Hinze scales reasonable.Larger velocity fluctuations may involve more than one eddy, making the spectrum heuristic incorrect.I'll use the functional form of the average to inform the choice of the model.For D 32 specifically, I'll use the inertial range spectrum, as the average is more strongly influenced by the smallest scales.Averages controlled by larger scales require a different length scale specification.The inertial range spectrum with the dissipation model ε 0 = C ε k 3/2 0 /Λ 0 suggests D 32 ∝ v −3 DF , which can be computed with a prescribed PDF.The value C ε = 0.43 is recommended by Pope (2000, p. 244).
To maintain analytical tractability, a power law velocity PDF (f v (v) = Cv −α ) will be used.A Gaussian PDF would be more realistic, but will be used in future work to keep this work simple.Both are compared in figure 5 implying v d ∝ v as hypothesized by Wu et al. (1992, p. 305), but for power law PDFs min with no v dependence.(Again, DF means v > v min here.)Using a power law PDF, I find that v −3 v > v min = (α − 1)v −3 min /(α − 2).To compute D 32 , I start with equation 45, then use the inertial range spectrum to eliminate , substitute in the dissipation and v −3 v > v min models, and choose v min = v σ (equation 41) to find which has a similar scaling to Wu et al. (1992, p. 308) for the initial value of D 32 , despite the difference in the definition (see equation 32).This is a consequence of the power law PDF.
Alternative choices could make how | DF and −1 DF −1 scale differ.To find the average droplet velocity v d I start with equation 40 and apply an approach similar to that for D 32 , noting that v −5 v > v min = (α − 1)v −5 min /(α + 4) for a power law PDF.I find that The theory will now be calibrated against experimental data.Only initial droplet diameter and velocity measurements are compared because the constant k and Λ approximations may be 10 inaccurate downstream.For initial D 32 , three data sources are available (Wu and Faeth, 1993;Wu et al., 1995Wu et al., , 1992)).For initial droplet radial velocity, I use data from Wu et al. (1992, p. 305).
None of these sources have rough tubes, so the data has almost no variation in Tu 0 .Fitting the models to the data, the coefficient C D32 = 0.522 (29 points, R 2 = 0.730) and the coefficient C vd = 0.254 (17 points, R 2 = −0.625).Figure 6 shows a comparison of the D 32 theory against the data.The measurement error in v d is large -60% according to Wu (1992, p. 129)making a close fit impossible for any reasonable model.Given the small variation in Tu 0 for the data, for the moment the most that can be said is that the theory is not inconsistent with the data for v d .Due to the high measurement error and consequential poor fit, no plot comparing the v d theory and the data is presented.
Similar procedures can find other diameters.The mass mean diameter has the same average mass as the ensemble averaged spray at that location, so m d = ρ πD 3 30 /6.Consistent with the model used to find v d , This term is more strongly influenced by the larger scales than D 32 . Trettel

Average breakup onset location, x i
I define the breakup onset location as the average distance eddies travel from the nozzle outlet in the time it takes for breakup to occur: x i ≡ (U 0 + u 0 )t b,0 DF ≈ U 0 t b,0 assuming that uv is small (because u and v are correlated) and that t b,0 is small (or else an integral with a random integrand would need to be computed).To second-order t b = C lig /v (see equation 40) so t b,0 ∝ /v DF , which is difficult to model.The term is not influenced by the smallest scales as much as D 32 .As such, I assume that the conditioning has little effect.By hypothesis, the parameters influencing the breakup time are σ (N/m), ρ (kg/m 3 ), and v 0 (m/s), from which a unique time scale can be formed: t b,0 ∝ σ/(ρ v 0 3 ), leading to This result is equivalent that of Kerstein et al. (2017) if one replaces their u τ,0 with v 0 .Their model would have no Re 0 dependence with this modification.The v 0 ∝ u τ,0 scaling implies Tu 0 ∝ √ f , similar to the regression for fully developed pipe flows (Tu 0 ∝ f 0.4587 ).This scaling is consistent with multiple physical pictures, not just the boundary layer scaling described by Kerstein et al.The available pipe jet data is unfortunately not able to distinguish between the two theories.Fitting the theory to the data returns C x i = 20.4(52 data points, R 2 = 0.905).A comparison of the theory and experimental data is in figure 7.

Average breakup length, x b
To determine the breakup length, I first calculate the average surface mass flux of droplets from the jet, ṁ .I decompose the surface into waves of wavenumbers κ ∝ 1/ in the streamwise and angular directions.I assume droplets are formed with frequency v/ and mass proportional to ρ 3 .I ensemble average to determine ṁ : which is constant downstream because I take k and Λ as constant.Similarly, the dimensionless quantity a constant.But the experiments of Sallam et al. (2002, fig. 10) show that this quantity increases with x from O(10 −2 ) to O(1).Consequently the model is not correct.The inaccuracy could be due to the ṁ model, v d model, or both.For simplicity, I assume that x i = 0 for the derivation of x b .Otherwise, a delay differential equation would be required to account for the delay between an eddy impacting the surface and droplet formation.After applying mass conservation for a particular realization of the jet to a differential element, I find that where U 0 is the (constant) jet convection velocity, the jet is assumed to have a circular cross section, d j (x) is the diameter of the jet at x, A(x) = πd 2 j /4 is the cross sectional area, and P (x) = πd j is the perimeter.Consistent with how x b is measured, I define x b with d j (x b ) ≡ 0, so to first-order d j ( x b ) = 0. Solving equation 52 for x b with the d j ( x b ) = 0 approximation using the ṁ model (equation 49), I obtain In the case of the breakup length, comparison against an empirical regression is worthwhile.The power law regression for the turbulent surface breakup regime developed in a companion paper (Trettel, 2020a) is (R 2 = 0.958): The signs of the Tu 0 and We 0 exponents are correct, but the magnitudes are in error.The most likely cause of the error may be the model for ṁ , as use of the ṁ correlation from Sallam et al. (2002, p. 446) for the turbulent surface breakup regime, ṁ /(ρ v d ) ∝ x/[Λ 0 (We 0 Λ 0 /d 0 ) 1/2 ], returns x b /d 0 ∝ Tu −3/10 0 We 3/10 0 .This suggests that the breakup length in the turbulent surface breakup regime is controlled mainly by turbulent primary breakup at the free surface rather than the Rayleigh mechanism, which breaks up the entire jet core.

Spray angle, θ i
Similar to previous works (Huh et al., 1998;Natanzon, 2018;and Skrebkov, 1966), I define the spray angle through tan θ i /2 ∝ v d /u d (at x = x i ).In other words, the spray angle is determined through simple geometry via the ratio of the radial to streamwise droplet velocities.The spray angle is a maximum angle rather than an average angle, as the observed boundary of the spray is the maximum extent of the spray, so the spray angle is not written as an average.However, it is assumed that the maximum is proportional to the average, similar to Markov's inequality.
As u d = U 0 +u, then u d = U 0 because there is an additional term with the correlation uv .I assume this effect is negligible as I did for x i , so u d = U 0 .Then tan θ i /2 = C θ i v d /U 0 , so Atomization and Sprays

Trettel
Fitting the theory against the experimental data returns C θi = 0.584 (5 data points, R 2 = −0.889).As shown in table 1, all models tested in this work fit the spray angle data very poorly.
It is worth examining a more general regression of the data to see what causes the poor fit.
The available spray angle data in the turbulent surface breakup regime was very noisy and lacked appreciable turbulence intensity variation.Consequently, a regression was made from two studies from the Faeth group (Ruff, 1990 andSallam, 2002) which were less noisy, likely due to using a more consistent definition of the spray angle (R 2 = 0.983): We 0.621 0 . (56) See Trettel (2020a) for more detail on this regression, including how the turbulence intensity term was developed.One major source of the poor fit is the variation with the Weber number.In CDRSV theory, tan θ i /2 decreases with We 0 , contrary to the regression.The only model I am aware of where tan θ i /2 increases with We 0 is that of Skrebkov (1966, p. 145), who suggests that (tan θ i /2) 2 = Tu 2 0 + 12Cρ g /ρ − 12/(DWe 0 ) for high Re 0 .The model of Huh et al. (1998) has no We 0 variation at all.Why CDRSV theory obtains the wrong scaling with the Weber number for the spray angle is unclear at present.

CONCLUSIONS
Conventional stability theory has, so far, failed to work in the turbulent surface breakup regime.Phenomenological theories like that developed in this work appear to have more promise due to their flexibility.However, future stability theories may completely supplant phenomenological theories if the issues identified in this work are solved.
Jet breakup does not occur for all surface fluctuations, so it is inappropriate for a model to imply that it does.Conditional averages must be computed to account for this feature of turbulent jet breakup.
While excellent agreement between current CDRSV theory and measurements was found for D 32 and x i , the theory has only modest success for x b , and none for θ i -see table 1. Ultimately no theory for turbulent jet breakup has been fully validated, in part due to the failures of these theories with existing data, and also because Tu 0 varies little in existing data.Alternative modeling choices might improve accuracy.
One possible avenue for improvement is using a more general turbulence spectrum like Schmitz (2011) rather than the Kolmogorov inertial range spectrum.This would likely require a computational model, which is why this approach was not used in this work.This work focused on the issues of the definitions of quantities of interest and analytical modeling rather than detailed modeling of each quantity of interest.
Another way to improve CDRSV theory would be to use a more accurate velocity probability density function function.A Gaussian probability density function would be more accurate.A Gaussian PDF was used in an example in § 3.3, and the theory can easily be extended to use Gaussian PDFs.
of Texas at Austin's Interlibrary Services are also thanked for help obtaining many obscure papers used in this work.
FIG.1: Jet breakup variables labeled on a schematic liquid jet.d 0 is the nozzle outlet diameter, x i is the average breakup onset location, θ i is the spray angle, and x b is the breakup length.

FIG. 2 :
FIG. 2:There are two possibilities when a turbulent event occurs at the free surface: breakup (top right) and no breakup (bottom right), depending on whether the velocity fluctuation exceeds the critical velocity v min which could be a function of the length scale associated with the fluctuation.
FIG. 5:A comparison of Gaussian and power law probability density functions for the tails of the velocity fluctuations.Both are conditioned on v > v min = 0.3v .
FIG. 7: Comparison of the x i theory (equation 48) against experimental data.
the result for a power law PDF,v | DF = v | v > v min = (α−1)v min /(α−2),and also chose v min = v σ (equation 41).The theory was fitted to the data, returning C x b = 5.62 (193 data points, R 2 = 0.719).The theory and experimental data is compared in figure8.
Trettelcenter.The quantities of interest are the average droplet diameter at formation (D ij , e.g., D 32 for the Sauter mean diameter), average droplet radial velocity at formation ( v d ), average breakup length ( x b ), average breakup onset location ( x i ), and average (full) spray angle (θ i ).I will typically drop the phrase "average" for the quantity of interest.Bars denote spatial averages, and angle brackets denote ensemble averages.