Modeling drop breakage using the full energy spectrum and a specific realization of turbulence anisotropy

The assumption of homogeneous isotropic turbulence when modeling drop breakage in industrially relevant geometries is questionable. We describe the development of an anisotropic breakage model, where the anisotropy is introduced via the inclusion of a perturbed turbulence spectrum. The selection of the perturbed spectrum is itself motivated by our previous large-eddy simulations of high-pressure homogenizers. The model redistributes energy from small to large scales, and assumes that the anisotropic part of the Reynolds stresses is confined to the energy-containing range. The second-order structure function arising from the perturbed spectrum is used in the standard framework of Coulaloglou and Tavlarides to calculate breakage frequency. While the base model exhibits non-monotonic behavior (by predicting a maximum value for a certain drop size), the effect of anisotropy is shown to increase breakage frequency in length scales larger than this peak, thereby reducing non-monotonicity. This effect is more pronounced for small turbulence Reynolds numbers.


| INTRODUCTION
Immiscible liquid-liquid systems are commonly encountered in chemical, petroleum, fast-moving consumer goods and pharmaceutical industries. Process design, scale-up and optimization still rely on timeconsuming and costly experimental programs, 1 motivating in-silico evaluation of emulsification processes via simulation techniques. The drop size distribution (DSD) in such systems is important because it determines critical features of the final product and affects the process dynamics through rheology. 2 The evolution of the DSD is often modeled through a population balance equation (PBE), which is a continuity statement for the population of drops. Several methods have been developed to solve the PBE, such as the classes method, 3 the Monte Carlo method, 4 and the methods of moments 5,6 ; particularly important are the quadrature-based moment methods, 7 such as the quadrature method of moments (QMOM), 8,9 since they are a computationally tractable approach to solve the PBE coupled to CFD. [10][11][12] These methods are suitable when the distribution of the dispersed phase is simple (e.g., monomodal or bimodal), 13 and the main interest is on its average properties rather than on the complete distribution.
For more complex distributions, a large number of moments is required which may impact the stability of the solution when higherorder schemes are employed. 14 A key challenge that applies to all the discussed methods lies in the submodels describing drop breakage and coalescence. 1 Many competing models have been published over many years, and a detailed overview of them is given in the review papers of Lasheras et al, 15 Liao and Lucas, 16 Solsvik et al, 17 Sajjadi et al, 18 Abidin et al, 19 and Falzone et al. 20 In this study, we focus on drop breakage frequency in turbulent flow. A major disagreement between the different models concerns their monotonic or non-monotonic behavior with increasing drop size. Early models and their extensions [21][22][23][24] predict a maximum breakage frequency value for certain drop size, while later models [25][26][27][28] impose monotonicity based on available experimental evidence, 25 arguing that a continuously rising breakage frequency with an increasing drop size is physically more plausible. However, this argument is questionable as recent experimental data 24 have shown non-monotonic behavior in the breakage frequency itself.
Drop breakage in turbulent flows is generally modeled in terms of Kolmogorov's second similarity hypothesis 29 and the second-order longitudinal structure function (hδu 2 i) for the inertial subrange of homogeneous isotropic turbulence (HIT). 30 Following the functional form of Coulaloglou and Tavlarides, 21 a general form of drop breakage frequency may be written as 24 where P b is the breakage probability, ξ is the drop diameter, and the breakage time is given by The underlying assumption is that the motion of the daughter drops is essentially that of turbulent structures in the flow field with the same length scale 21   Another widely used class of models, [35][36][37] which follows a monotonic behavior with increasing drop size is the approach of Narsimhan et al, 27 and its extensions by Alopaeus et al 28 and Laakonen et al. 38 The main assumption in this class is that drop frag- The breakage probability in Equation (4) (where κ = 2π/r) corresponds to a second-order structure function of the form the exponents p and q are 41 In this respect, the implied structure functions of the Coulaloglou and Tavlarides model 21 where and F I G U R E 2 (a) Energy spectrum from Equation (8) account for large and small length scale deviations from the classical inertial (−5/3) behavior. L is the integral length scale, η is the Kolmogorov length scale, β 0 = 5.2, and c L and c η are turbulence Reynolds number dependent parameters. 43 If we calculate the second-order structure function for this spectrum 44 then, for a turbulence Reynolds number Re L = 10 5 , the resulting structure function is shown in Figure 2b. Calculating the breakage time for this structure function using Equation (2), then its slope is shown in Figure 2c; it assumes a value of 2/3 in the inertial subrange, reflecting the scaling of the inertial subrange-based Coulaloglou and Tavlarides model. 21 The decrease of t b toward a zero slope in the dissipation range reflects the scaling of the Alopaeus et al model. 28 In the range ξ/η $ 3 to ξ/η $ 15 it assumes intermediate values between 0 and 2/3.

| BREAKAGE FREQUENCY USING THE FULL SPECTRUM BASED ON HOMOGENEOUS ISOTROPIC TURBULENCE
To obtain the breakage frequency for the full spectrum, one may use the Coulaloglou and Tavlarides approach by replacing β(ε ξ) 2/3 in the original model (Equation (3)) with hδu 2 i from Equation (11). 45 The breakage frequency is then written as As shown in Figure 2d, the main effect of accounting for the full turbulence spectrum is that the breakage frequency in both the small and the large scales is smaller than that in the inertial subrange, due to the larger breakage times in these regions. The decrease in the small scales is observed for ξ/η ≲ 15, and is driven by the decrease of the breakage time slope from 2/3 to 0 shown in Figure 2c.
The extension of drop breakage modeling to the full spectrum presents a significant development which potentially reduces experimental costs required to refit the parameters in prior breakage frequency models. 20,30 However, the assumptions of homogeneity and isotropy used are still questionable. Solsvik et al 30,46 argued that the effect of anisotropy is negligible compared to the effect of finite Reynolds number on breakage frequency; however, this view was not explored in further detail. As it will be discussed in the last section where we compare breakage frequencies based on isotropic and anisotropic spectra, the relative influence of anisotropy on breakage frequency increases with decreasing Reynolds number.

| DEVIATIONS FROM HOMOGENEOUS ISOTROPIC TURBULENCE
Our previous work on LES of a high-pressure homogenizer 47 showed that the regions in the mixer associated with peak drop fragmentation stresses also exhibit strong turbulence anisotropy; 1-D energy spectra highlight an amplification of energy in the energy-containing range associated with vortex shedding. Kolmogorov's standard 2/3 scaling Perturbed (E 11,ani ) and unperturbed (E 11 ) energy spectra from Equations (13) and (23), respectively, for Re L = 10 5 and the parameters included in Table 1. For E 11,ani : A = 16, κ 1 = 2π=L 11 , L 11 = 0.43L, and σ 1 = κ 1 =10. E 11,ani integrates to the same k as E 11 . The vertical line denotes the low wave number limit of the inertial subrange. (b) Structure functions from Equation (14) corresponding to E 11,ani and E 11 . The vertical line denotes the upper length scale limit of the inertial subrange [Color figure can be viewed at wileyonlinelibrary.com] in the corresponding second-order structure functions is thus masked by vortex shedding and the resulting structure functions exhibit slopes between 2/3 and 2 in the inertial subrange. In a separate study of a rotor-stator mixer, 48 PIV-measured 1-D energy spectra at the stator slot jet reveal deviations from model spectra based on HIT. The deviations manifest themselves as steeper than −5/3 slope in the low wave number end of the inertial subrange and a narrow range of wave numbers with −5/3 slope. A previous PIV study of the same rotor-stator configuration 49 identified periodic vortex shedding inside the stator slot due to the rotor blade passage. The authors in Reference 48 concluded that the deviations between the measured spectra and the model spectra are consistent with finite Reynolds number effects (with Re L $ 1000 − 3000). In the high-pressure homogenizer study, 47 Re L was also moderate 50 ranging from $10 3 − 10 4 .
The amplification of energy in these studies is restricted to the energy-containing range and the low wave number end of the inertial subrange; its effect, however, on the second-order structure function extends to a much broader range of scales, as demonstrated by Portela et al. 51 By assuming a 1-D model energy spectrum 41 of the form (where σ 1 represents the spread of the spectral peak, A is an amplification factor and E 11 is the 1-D model spectrum), the second-order structure function was calculated as 51 where ℱ −1 denotes an inverse Fourier transform. Figure 3a shows a perturbation restricted to the energy-containing range. The low wave number limit of the inertial subrange is indicated by the vertical line. The departure from isotropy can be expressed in terms of the deviatoric part of the Reynolds stress tensor 52 : Assuming that the anisotropy is confined to the large scales and considering only the largest principal stress u 2 1 Ã 2 3 k,2k Â Ã , 52 then we model the 1-D energy spectrum according to Equation (13). κ 1 and σ 1 are suitably chosen such that the perturbation is confined to the energy-containing range. We assume that the contribution to u 2 1 Ã arising from anisotropy is equal to the amount of energy added to the spectrum via the perturbation. u 2 1 Ã is thus written as Since it follows that where f A depends on Re L , κ 1 , and σ 1 (but not on stress anisotropy).
Equation (18) states that the perturbation parameter A is proportional to the departure from turbulence isotropy. To ensure the perturbed spectrum integrates to the same total turbulent kinetic energy as the unperturbed spectrum, Equation (16) is rescaled as follows In the Supporting Information it is shown that and hence, the proposed normalized perturbation to the 1-D energy spectrum is finally written as The energy associated with the perturbation is assumed to be normally distributed around a mean wave number κ 1 = 2π L11 with a SD σ 1 . L 11 is the longitudinal integral length scale which depends on Re L and decreases asymptotically to 0.43L with increasing Re L . 41 For given k, ε, u 2 1 Ã and the above assumptions about κ 1 and σ 1 , Equation (22) returns a spectrum shape which mimics the 1-D spectra observed in our LES of a high-pressure homogenizer. 47 This is achieved by redistributing energy from across the full spectrum to the wave numbers local to the vortex shedding frequency. The model may thus be considered as a specific realization of turbulence anisotropy pertinent to fluid flows dominated by vortex shedding. It is noted that the model is not limited to the configuration studied in our LES: any modeled anisotropy from any known source can in principle be accounted. The second-order structure function is obtained via the Fourier transform of the perturbed spectrum by 41

| Model assumptions
To describe the evolution of energy spectra and second-order structure functions, a model able to capture spectral dynamics is We have thus restricted ourselves to a specific shape of the spectrum and a simple algebraic expression to link anisotropy to the perturbation in the spectrum which is tractable for use in CFD-PBM simulations.

| Comparison of model spectra with LES spectra
shown in Figure 4b, reduces rapidly with increasing streamwise distance; in the region of peak dissipation rate the magnitude of a 33 Ã remains considerable. Re L at the same locations (shown in Figure 4c) is O 10 3 −10 4 .

| Comparison of structure functions based on 1-D and 3-D energy spectra
Before implementing the proposed model (Equations (22)- (24)) into a breakage frequency model it is instructive to examine its behavior in the limit of isotropic turbulence. In the isotropic limit, 2 3 , and Equation (22) simplifies to Equation (23). We can then compare this with the 3-D energy spectrum-structure function model in HIT 30,42 (Equation (8) and (11)). Figure 5a shows a 3-D spectrum from Equation (8) and a 1-D spectrum from Equation (22) for 3 . The corresponding structure functions, obtained from Equations (11) and (24), respectively, are shown in Figure 5b; the structure functions coincide for the entire dissipation range and inertial subrange. There is a minor difference in the energy-containing range shown in the inset of Figure 5b which, however, is relatively less important as the DSDs typically lie within the dissipation range and the inertial subrange. We thus conclude that the isotropic limit of our 1-D energy spectrum-structure function approach using Equations (23) and (24) is equivalent to the HIT 3-D energy spectrum-structure function approach. 30,42 F I G U R E 5 (a) 3-D and 1-D model spectra in HIT from Equations (8) and (23), respectively, and; (b) their corresponding second-order structure functions from Equations (11) and (24) [Color figure can be viewed at wileyonlinelibrary.com]

| Comparison between isotropic and anisotropic models
We now examine the effect of turbulence anisotropy on breakage time and frequency by implementing the 1-D energy spectrumstructure function in a breakage frequency model. The breakage frequency is calculated in the same manner as with other full spectrum approaches (Equation (12)). Figure 6 shows energy spectra, structure functions, breakage times and frequencies for different values of In the small length scales of the inertial subrange, the breakage time of the full spectrum isotropic approach (Equations (11) and (12)) and that of the Alopaeus et al model 28 may be viewed as two extremes of the length scale exponent ( 2 3 and 0, respectively) since from Equations (2) and (6) F I G U R E 6 (a) Energy spectra from Equation (22) where q ≤ 2 (Equation (7)) is the slope of the structure function. Our

| Comparison between anisotropic and intermittency models
In the previous section it was shown that the introduction of anisotropy via a perturbation to the energy spectrum alters the slope of the structure function in the inertial subrange, thereby affecting the shape of the breakage frequency curve. A different phenomenon which also has an impact on that slope is intermittency. In this section, we compare the proposed anisotropic model with an intermittency model in order to assess their relative importance.
Intermittency of turbulence refers to spatial fluctuations of quantities based on velocity differences or velocity gradients, for example vorticity or energy dissipation rate. Based on the scale of the fluctuations, intermittency is distinguished as external or internal. 44 External intermittency is related to the large scales and is not considered to be universal. Internal intermittency is related to the small scales and is considered to affect the universality of small-F I G U R E 7 (a) Energy spectra from Equation (22); (b) structure functions from Equation (24); (c) breakage times and; (d) breakage frequency functions from Equation (12) for the parameters included in Tables 1 and 2 where ζ n is the revised exponent in δu n h i/r ζ n and μ = 0.25 is a model parameter fitted with experimental data. 41 The correction is considerable for higher-order structure functions. The revised exponent of the second-order structure function, relevant to droplet breakage modeling, is $4% larger than the value predicted by the original theory. This is considered a small change taking account of the other uncertainties involved in droplet breakage. 63,64 Meyers and Menevau 65 introduced the effects of intermittency by modifying Pope's 3-D model spectrum. 41 Their model is written as where F I G U R E 8 (a) Isotropic and anisotropic spectra from Equation (22), and intermittency-corrected spectrum from Equations (23) and (28); (b) structure functions from Equation (24); (c) breakage times and; (d) breakage frequency functions from Equation (12) for the parameters included in Tables 1 and 2 and Re L = 10 5 [Color figure can be viewed at wileyonlinelibrary.com] and are reshape functions for the energy containing range and the dissipation range, respectively, and the factor is a bottleneck correction used to shape the spectral bump on the onset of the dissipation range; β = μ/9 is the intermittency correction for the inertial subrange and μ = 0.25, defined in Equation (27); α 1 − α 5 are model parameters, fitted with experimental data, and we use their respective values reported by Solsvik and Jakobsen. 64 An intermittency-corrected 1-D spectrum is obtained using Equation (23) by numerically integrating the 3-D spectrum, Equation (28), similar to the approach of Meyers and Menevau. 65 The result is shown in Figure 8a along with isotropic and anisotropic 1-D spectra.
Comparing to the baseline isotropic spectrum, the intermittency model redistributes energy to the high wavenumbers of the inertial subrange, to account for the bottleneck effect, while the anisotropy model redistributes energy to the energy containing range. The amount of redistributed energy with the anisotropy model is much larger than that with the intermittency model, even though we have used a moderate value for the stress anisotropy,

| CONCLUSION AND FUTURE WORK
The majority of turbulent drop breakage in emulsification devices is expected to occur in locations where turbulence production and dissipation peak. These regions are characterized by turbulence anisotropy and deviations of the energy distribution from model spectra assuming homogeneous isotropic turbulence.
We have presented a novel approach aiming to capture the effect of a specific realization of turbulence anisotropy on the second-order structure function used in drop breakage modeling. The specific realization is related to spectral peaks and slopes in the inertial subrange steeper than Kolmogorov's −5/3 scaling, which have been observed in industrially relevant geometries. Our approach incorporates the recent extension 30 of drop breakage modeling to the full spectrum of isotropic turbulence, as a limiting case of zero anisotropy. By using the new structure function in the model of Coulaloglou and Tavlarides, 21 we have shown that the main effect of anisotropy is to alter the non-monotonic behavior of an HIT-based breakage frequency toward a monotonic one; this effect is more pronounced for smaller turbulence Reynolds numbers. Apart from drop breakage processes, the proposed model may also be used in rate functions of other particulate processes related to the turbulent stresses, such as drop coalescence.
In future work, this model is to be implemented in a CFD solver and used in CFD-PBM simulations of turbulent emulsification in highpressure homogenizers. Instead of two-equation RANS models, typically used in such studies, a Reynolds-stress transport model is employed to account for turbulence anisotropy, which is then introduced in drop breakage modeling through Equations (22) and (24).