Air-blast atomization of liquid films, when exposed to highspeed co-flowing gas streams, leads to the production of fine droplets. This phenomenon has a wide range of applications across various industries, including agricultural irrigation, spray painting, gas turbines and spray coating 1,2 . Given the importance of accurately predicting spray characteristics and optimizing industrial devices, considerable research has been devoted to understanding the breakup dynamics and the mechanisms underlying this process [3][4][5][6][7] . It is shown that the liquid film undergoes oscillations in both the streamwise and spanwise directions, forming longitudinal and transverse waves, driven by Kelvin-Helmholtz (KH) and Rayleigh-Taylor (RT) type instabilities, respectively. As these instabilities grow, the gas flow stretches the liquid film, leading to the formation of hollow bags between transverse waves. These bags eventually rupture into multiple holes, which collide with one another, fragmenting the thin film forming the bag structures into fine droplets and ligaments. This atomization process spans multiple timescales and length scales, including the atomizer lip thickness, film thickness, instability wavelengths, boundary layer thickness, and turbulence scale. As a result, gaining a detailed understanding of these mechanisms is complex, yet essential.
The stability of a thin, sinusoidally fluctuating liquid sheet moving through quiescent air was analyzed theoretically by Squire 8 , showing that the longitudinal wavelength of maximum instability, λ lon , is governed by a balance between surface tension and inertial forces. Hagerty and Shea 9 demonstrated that planar liquid sheets exhibit two distinct wave a) Electronic mail: gilou.agbaglah@wayne.edu modes: sinuous (in-phase surface oscillations) and dilational (out-of-phase) modes. However, under identical conditions, the sinuous mode consistently grows faster, making it the dominant breakup dynamics. Models for droplet diameter, based on theoretical analysis, were developed by examining the disintegration of the liquid film into ligaments and, subsequently, into droplets due to Rayleigh instability 10,11 . A review of liquid film rupture was proposed 12 , exploring various factors such as laser pulses and surface tension heterogeneity. The fragmentation phenomena of films and ligaments are summarized in 13 , which also discusses the droplet diameter distribution resulting from ligament breakup. Additionally, a compound gamma distribution with two parameters is used to model the droplet diameter distribution generated by the breakup of the liquid film produced by swirl and fan spray nozzles 14 . Empirical correlations for droplet diameter in airblast atomizers are proposed by several researchers 15,16 . However, existing models and empirical correlations are often limited in their applicability, as they fail to capture the full complexity of atomization processes, such as gas flow-induced transverse fluctuations, hole collisions, and the formation of bag-like structures. As a result, these correlations require parameter tuning for each specific application and set of conditions. A recent approach to modeling the mean droplet size resulting from the fragmentation of liquid films involves considering two distinct types of film acceleration: first, due to the KH instability, and second, due to the drag force exerted by the co-current gas flow 7 . Large droplets are primarily attributed to ligament breakup, while small droplets result from bag breakup. The mean errors obtained when comparing the predictions to the experiment, using a liquid film thickness of 200 µm, are 47% for the ligament size, 31% for the mean diameter of large droplets, and 60% for the mean diameter of small droplets. This highlights the need to improve the prediction of the droplet size distribution in the atomization of air-blasted liquid films, which is crucial for industrial applications.
In a previous study 17 , we investigated the aerodynamicallydriven rupture of a liquid film through three-dimensional numerical simulations. A thin liquid film with a thickness of h l = 25µm was used to accelerate the fragmentation process. It is worth noting that this thickness is approximately ten times smaller than the values typically used in experimental studies, e.g. 7,18 . A fixed gas Reynolds number, Re g = ρ g U g h l /µ g = 5000, was considered, where ρ g and µ g represent the gas density and viscosity, respectively, and U g is the gas velocity. We demonstrated that the liquid film is entrained by the fastmoving gas streams, causing it to oscillate and stretch in the longitudinal direction. Small transverse deformations at the liquid interface were observed and attributed to the influence of rib-like vortical structures in the gas flow. Due to the small film thickness considered, bag structures do not form between the transverse waves, as seen in experimental observations. Instead, the liquid film between the transverse waves fragments through hole nucleation and subsequent collisions.
In the present study, we use the same liquid film thickness of 25µm but vary the gas Reynolds numbers (Re g = 300, 500, 1000) to examine the inertial and viscous effects on the gas flow dynamics and liquid film fragmentation. For the gas flow dynamics, distinct three-dimensional vortex shedding modes are observed depending on Re g . These modes significantly influence the liquid film fragmentation process, controlling the transverse wavelength and the localization of hole formation, which ultimately determines the number of ligaments and the droplet size distribution. This study clearly highlights the pivotal role of instability modes in the gas flow, establishing a direct connection between the gas vortical structures and the breakup process of the entrained liquid film. This novel insight opens the door to enhanced predictive models for air-blast atomization of liquid films.
We consider a two-phase (air/water) system composed of a liquid and a gas, each with constant density, ρ g for the gas phase and ρ l for the liquid phase, and dynamic viscosities µ g and µ l , respectively. The dynamics of the system are governed by the incompressible Navier-Stokes equations, solved using the one-fluid approach, as outlined below:
where u is the velocity vector, p the pressure, D = [∇u + (∇u) T ]/2 the deformation tensor. µ ≡ µ(x,t), ρ ≡ ρ(x,t) and σ represent the viscosity, density and the surface tension respectively.
The surface tension force is applied exclusively at the interface and is modeled using the Dirac delta function δ s . The curvature of the interface is denoted by κ, and the unit normal vector perpendicular to the interface is given by n.
For simulations, we use the open-source solver Basilisk 19 which has been widely used for a variety of multi-phase flow problems [20][21][22][23][24] . In Basilisk, the liquid-gas interface is tracked using the Volume-of-Fluid (VOF) method on an octree grid, which enables adaptive mesh refinement. The physical properties, such as µ and ρ, are determined based on the volume fraction of the liquid phase f using linear interpolation:
The transport equation of f is given by:
Surface tension is computed using a balanced-force approach, while the interface curvature is accurately determined using the height-function method. Adaptive mesh refinement is applied using a wavelet-based estimation of discretization error, as described in 25 . Further details on the numerical methods implemented in Basilisk can be found in 26 .
We first perform simulations of a double-shear layer gas flow, in which a gas layer of thickness h gm = 25 µm is sandwiched between two high-speed gas streams, each with a thickness of h g = 10h gm . The inflow velocity field (u,0,0) is described by the following equations:
U g tanh
where the domain subdivisions is described below:
(a) bottom co-flow: -L y /2 ≤ y < -0.5h gmh g (b) bottom fast gas: -0.5h gmh g ≤ y < -0.5h gm (c) middle gas layer: -0.5h gm ≤ y < 0.5h gm
The middle gas layer has a velocity U gm = U g /26, where U g represents the velocity of the two co-flowing high-speed gas streams. This velocity ratio is chosen for a reason that will be explained later. U g is defined based on the corresponding Reynolds number Re g = ρ g U g h gm /µ g = 300, 500, 1000. Additionally, we account for the entrainment of the gas by adding This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. FIG. 2. Mean kinetic energy across the entire computational domain as a function of dimensionless time for ( Re g = 500). E k0 represents the mean kinetic energy at tU g /δ g = 0.
5% of the inflow velocity between the gas streams and the upper and lower walls of the simulation domain, following a similar approach to previous studies [27][28][29] . The gas boundary layer is set to δ g = h gm /2 and L y = 50h gm represents the vertical distance of the computational domain. It is worth noting that the velocity profile described above is similar to the one used in a previous study 30 to induce the flapping of the liquid sheet. Except for the finite thickness of the two co-flowing high-speed gas streams, the double-shear layer flow defined by this velocity profile resembles the velocity distribution observed in the wake of bluff bodies, such as cylinders, at comparable Reynolds numbers 31 . The simulation domain extends over a region of 50h gm × 50h gm × 50h gm in the streamwise (x), normal (y), and spanwise (z) directions, respectively, see Fig. 1(a). The domain is discretized using adaptive mesh refinement, with up to 2 11 grid points along each axis. This corresponds to a minimum grid size of approximately ∆ = h gm /41 = δ g /20.5 ≈ 0.6 µm, which would result in approximately 8.6 billion cells if a uniform grid were used. It is worth noting that lower mesh resolutions, specifically ∆ = h gm /6 and h gm /22, have previously been used and shown to accurately capture holes expansion and collision in thin liquid sheets 17,32 . The mesh is automatically adapted using the VOF field f , the velocity components, and the vorticity magnitude, with refinement tolerances set to 0.01, 1, 1, 1 and 4000 respectively. A Dirichlet boundary condition is applied at the domain inlet, with the velocity field (u, 0, 0). The inlet velocity profile u, defined by the equations above, is illustrated in Fig. 1(b). The downstream boundary uses an outflow condition, while symmetry conditions are applied at the upper and lower boundaries of the domain. Additionally, periodic boundary conditions are imposed in the spanwise direction.
The gas flow simulations are run until the system reaches a statistically steady state, which is determined by the temporal evolution of the average kinetic energy across the computational domain, e.g. see Fig. 2. Once steady-state conditions are achieved, the results are used to initialize the two-phase flow simulation described in section II C. Different vortex shedding modes are observed depending on Re g . For low Reynolds numbers, Re g = 300, 500, oblique cellular vortex shedding modes are observed, while for Re g = 1000, parralel vortices are seen, with three-dimensional streamwise vortical structures resembling the rib-like patterns observed in the wake of cylinders [33][34][35][36] , see Fig. 3 (c1). As shown in Fig. 3 (a1) and (b1), two cells are observed along the span for Re g = 300 while only one cell forms at Re g = 500. Cellular oblique vortex shedding has been previously observed in studies of wake flow past cylinders [37][38][39] . The formation of multiple cells is often linked to vortex dislocations at the junctions between cells. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
We g (= ρ g U 2 g δ g /σ ) 300 0.8 0.5 500 0.9 0.8 1000 1 3 TABLE I. Simulation parameters.
Two-phase flow simulations are subsequently performed, in which a thin, low-velocity liquid sheet is subjected to shear forces from faster-moving air streams, a scenario commonly referred to as an air-blasted liquid sheet. The gas flows simulated earlier in Section II B are used to initialize the twophase flow simulations. Dimensionless parameters used to control the simulations are listed in Table I, where M represents the dynamic pressure ratio and We g the Weber number of the gas stream. The liquid sheet is positioned midway between the fast-moving gas streams, with a thickness h l = h gm = 25 µm and a velocity U l = U gm = U g /26. This velocity ratio is kept constant across the three different Reynolds numbers and is chosen to approximate the convective velocity of liquid waves 40 . The initial streamwise length of the liquid sheet is set to 7.5h l . For the two-phase flow simulations, the origin of time is defined at the end of the gas flow simulation, t fg ; that is, the dimensionless time is given by t * =(tt fg )U g /δ g . Fig. 3(a2)-(a3), (b2)-(b3) and (c2)-(c3) illustrate the dynamics of the liquid sheet for Re g = 300, 500 and 1000 respectively. The liquid sheet undergoes longitudinal oscillations and deforms in the transverse direction. These deformations closely follow those observed in the gas flow as highlighted by the vortical structures in Fig. 3(a1), (b1) and (c1). The three-dimensional oblique waves and rib-like structures of the gas flow are clearly imprinted on the liquid sheet, providing strong evidence that its dynamics are governed by flow structures already present in the gas phase. The angle of the oblique waves, as well as the longitudinal and spanwise wavelengths, match those observed in the gas flow. The velocity distribution within the liquid sheet is shown in Fig. 4, where the instantaneous streamwise velocity U, Fig. 4 (a1) (a2) (a3) (a4) (a5) (a6) FIG. 5. Zoomed-in view of hole collisions and the subsequent formation of ligaments for Re = 500. (a1)-(a6) present a time sequence of snapshots at t * = 29, 36, 42, 49, 63, and 77. The color indicates velocity magnitude ranging from 5U l to 15U l .
is accelerated downstream by the gas flow, with the streamwise velocity at the free end of the sheet being approximately 10 times larger than the liquid inlet velocity. The normal velocity alternates between positive and negative values in the downstream direction, indicating the longitudinal oscillations of the liquid sheet. While the spanwise velocity appears noisy and lacks clear patterns for Re g = 300 and 500, alternating positive and negative velocities are observed in Fig. 4(c3) for Re g = 1000. This spanwise velocity pattern is clearly a result of the rib-like vortex structures present in the gas flow.
The fragmentation process of the liquid sheet varies depending on Re g and the vortical structures present in the gas flow. For Re g = 300, the liquid sheet fragments at its free end due to longitudinal rim instability 41 , forming liquid fingers, as shown in Fig. 3(a3). The liquid fingers subsequently break up through the end-pinching 42 or Rayleigh-Plateau-type mechanism. At Re g = 500, fragmentation occurs through hole formation and subsequent collision. Specifically, the sheet thins near the troughs and crests of the longitudinal waves and even-tually ruptures, forming multiple holes that expand over time and merge, in a manner similar to that described in 32,43 . A time sequence of hole growth and collision is shown in Fig. 5. Following the collision of holes, ligaments are formed, which then stretch and break up into droplets. For Re g = 1000, holes form between the spanwise waves, as shown in Fig. 3(c3), in a manner similar to that reported in 17 . Note that for this particular Re g , complete fragmentation of the liquid sheet occurs more rapidly than for the other two cases, due to the formation of stronger transverse deformations. This enables more detailed analysis and examination of droplet statistics, as presented in the following section.
In this section, we analyze the wavelength of the rib-like structures formed in the prior gas flow simulation for Re g = 1000, the transverse waves that develop on the liquid sheet, This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Transverse Wavelengths λtra 𝐴 = 𝜋𝑟 2 Ligaments λlon λ FIG. 7. Top view of the air-blasted liquid sheet fragmentation for Re g = 1000 at t * = 102, showing transverse waves and ligaments (left). The right image presents a side view, showing the longitudinal oscillation.
and the droplet size distribution resulting from its complete fragmentation. The droplet size distributions at low Reynolds numbers, Re g = 300, 500 are also analyzed.
Fig. 6(a) shows the magnitude of the streamwise vorticity in the gas flow along a vertical zy plane at a fixed downstream location x/h l = 10 at t * = 0. The obtained spatial undulation, exhibiting alternating positive and negative vorticity in the z-direction, clearly reflects the rib-like structures shown in Fig. 3(c1). A Fast Fourier Transform (FFT) algorithm is applied to the spatial undulation to identify the most representative wavelength. The process begins with the application of a Hamming window to minimize spectral leakage, followed by signal filtering using fourth-order Butterworth filters to separate the low-and high-frequency components. The dominant wavelength is then determined by applying the FFT to the high-frequency component and taking the inverse of the resulting wavenumbers. As shown in Fig. 6(b), the wavelength corresponding to the maximum amplitude is λ tra /h l = 2.1. It is worth noting that a transverse wavelength of 1.5 was obtained at a higher Reynolds number, Re g = 5000, in 17 . This suggests that λ tra /h l decreases with increasing Reynolds number. A similar trend is reported in the experimental study by Oshima & Sou 44 and aligns with the rib-like structures observed in the wake of a circular cylinder. For example, Williamson 45 demonstrates that smaller-scale structures be-come more prominent at higher Reynolds numbers (see their Figure 5). Given that the computational domain spans a distance of 50h l , this suggests that approximately 24 transverse waves are expected to form on the liquid sheet. This number closely matches the number of ligaments (≃ 26) formed following hole collisions, as illustrated in the top-view image in Fig. 7(left), which shows transverse waves near the inlet and ligaments downstream. The longitudinal oscillation of the liquid sheet is also visible in the side-view image in Fig. 7(right). The agreement between the wavelength of the streamwise vortical structures in the gas flow and those observed on the liquid sheet and downstream ligaments underscores that the breakup process of the thin liquid sheet can be predicted by analyzing the characteristics of the vortical structures in the gas phase.
The droplet size distribution resulting from the fragmentation of an air-blasted liquid sheet is critical for various industrial applications. We compare the droplet size probability density functions (PDFs) obtained from simulations at the specified Reynolds numbers: Re g = 300, 500 and 1000. A built-in tool of the Gerris flow solver 46 (the predecessor of Basilisk) is used to detect and compute the size distribution of droplets. This tool identifies each connected liquid region using the tracer field f , integrates f dV over that region to obtain its volume (dV is the volume of a computational cell), and then assigns a size-ranked index. An equivalent spher-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0278367
FIG. 8. Probability Density Function (PDF) of the droplet size distribution for Re g = 300 at t * = 298 (a1), Re g = 500 at t * = 183 (b1) and Re g = 1000 at t * = 115 (c1). The log-normal distribution (red line) and the two-parameter gamma distribution (blue line) are plotted based on Equations ( 6)- (8). 〈d〉 denotes the average droplet diameter, and the vertical black dashed line indicates the minimum mesh size, ∆. (a2), (b2), and (c2) show the same data as in (a1), (b1), and (c1), but plotted on a logarithmic scale.
ical diameter is subsequently computed from the calculated volume. The droplet size distribution is commonly analyzed using two primary distribution functions found in the literature: the log-normal distribution, Equation ( 6), and the twoparameter gamma distribution, Equation (7), 13,[47][48][49] .
where d ′ = d/ d is the normalized droplet diameter, with d the mean droplet diameter, µ and σ are the distribution's mean and standard deviation,
where
are the shape and rate parameters, respectively, estimated using the method of moments. Γ denotes the gamma function. It should be noted that there is ongoing debate regarding the most appropriate distribution function for representing sprays. The log-normal distribution is commonly used when sprays result from a combination of various breakup mechanisms, whereas the gamma distribution is typically assumed when the spray is produced by a single or dominant breakup mechanism 49,50 . Figure 8 presents the obtained PDFs along with their comparison to log-normal and gamma distributions. The bottom row displays the same data on a logarithmic y-scale to emphasize the distribution of larger droplets. The mean and standard deviation are computed from the simulation data, while the shape and rate parameters for the gamma distribution are determined using Equation (8). It is important to note that no curve fitting is applied. The results indicate that the log-normal distribution provides a better fit to the droplet size data, particularly around the peak, which represents the most probable droplet size, see Figures 8 (a1), (b1), and (c1) corresponding to Re g = 300, 500 and 1000, respectively. However, the tail of the distribution, shown in Figures 8 (a2), (b2), and (c2), is not well captured by either distribution function at the low Reynolds number Re g = 300, although the log-normal distribution provides a closer approximation. At higher Reynolds numbers, Re g = 500 and 1000, the lognormal distribution shows better agreement with the distribution of larger droplets. This trend is consistent with the involvement of multiple breakup mechanisms, such as hole collision, rim and ligament breakup, during the fragmentation of the liquid sheet.
To examine the temporal evolution of the droplet size distribution, a time-resolved diagram is constructed, as shown in Figure 9. Note that the x-axis range differs in Figures 9 (a FIG. 10. Histogram of the droplet size distribution for Re g = 1000 at t * = 115. Results without the Manifold Death (MD) algorithm are shown in blue, while simulations with the MD algorithm using hole creation intervals of δt * = 7 and δt * = 27 are shown in green and red, respectively. The vertical black dashed line indicates the minimum mesh size, ∆.
mains consistent, aligning with recent experimental observations of aerodynamic droplet breakup 51 .
It is well known that numerical simulations of thin liquid sheet fragmentation face significant challenges due to meshinduced spurious breakup. This artifact leads to the premature rupture of liquid films when their thickness approaches the grid resolution limit, making it difficult to achieve mesh convergence in the droplet size distribution. To address this issue, a Death (MD) algorithm is recently developed 52 , in which thin sheets are artificially perforated upon reaching a predefined critical thickness, independent of the grid resolution. Using the MD algorithm to define a preset cut-off film thickness, a grid convergence of the droplet size distribution was achieved and validated against experimental data 23,53 . The critical film thickness, h c , is defined as h c = 3∆ sig , where ∆ sig denotes the grid resolution used to calculate the quadratic moments of the VOF color function, which serves to identify thin liquid structures. h c represents the smallest length scale at which liquid films can be accurately represented, as below this threshold, they would be artificially perforated by the MD algorithm. To verify that the most representative droplet size observed in the distributions of Figure 8 is not a numerical artifact, simulations are performed both with and without the MD algorithm. These simulations are conducted at Re g = 1000, as this case requires less computational effort compared to the other Reynolds numbers considered in this study. When the MD algorithm is applied, ∆ sig is set to 2∆, and the formation of up to 10 holes is permitted. To investigate the effect of perforation frequency on the droplet size distribution, two different time intervals for hole creation are examined: δt * = 7 and δt * = 27. Figure 10 presents the droplet count as a function of spherical diameter, with blue representing the results without MD, green for simulations using MD with δt * = 7, and red for δt * = 27. Interestingly, the dominant droplet diameter in the histogram remains unchanged with the use of the MD algorithm, suggesting that hole formation and collisions are well resolved by the simulation grid. However, a reduction in the number of droplets, particularly those of smaller sizes, is observed when the MD algorithm is applied. This effect becomes more pronounced with increased perforation frequency, indicating that the MD algorithm helps suppress the formation of small fragments that may arise from numerical artifacts. Figure 11 compares the droplet size PDFs obtained from simulations with the MD algorithm (with δt * = 27) to the log-normal and gamma distributions. While the agreement of the log-normal distribution with the small droplet sizes is similar to that obtained without the MD algorithm (see Fig. 8(c1)), the use of the MD algorithm yields better agreement for larger droplets, particularly for normalized droplet diameters between 3 and 8.
The aerodynamically driven fragmentation of a thin liquid sheet with a thickness of h l = 25 µm is investigated using three-dimensional numerical simulations. The liquid sheet is positioned between fast-moving gas flows on both its top and bottom surfaces, which are first simulated until the sys-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. tem reaches a statistically steady state. Three-different gas Reynolds numbers, namely Re g = 300, 500 and 1000 are considered to investigate the effect of gas flow on liquid sheet breakup. The dynamics of the preceding gas flow simulations reveal distinct downstream vortical structures for each Re g value. Specifically, an oblique cellular vortex shedding mode is observed for Re g = 300 and 500, featuring two and one cells along the span, respectively. In contrast, parallel rib-like vortical structures are observed for Re g = 1000. It is worth noting that similar vortical structures have previously been observed in the wake of bluff bodies [33][34][35][36][37][38][39] . The liquid sheet exhibits longitudinal and transverse oscillations that closely mirror those in the surrounding gas flow. For Re g = 300 and 500, oblique flapping motions are observed, while at Re g = 1000, transverse undulations emerge with wavelengths comparable to those of the rib-like vortical structures in the gas flow. At Re g = 300, the liquid sheet fragments at its free end due to longitudinal rim instability 41 , forming liquid fingers that subsequently break up via the end-pinching 42 or Rayleigh-Plateau-type mechanism. For Re g = 500, the sheet thins near the troughs and crests of the longitudinal waves and eventually ruptures, leading to the formation of multiple holes. These holes expand and collide over time, generating ligaments that break up into droplets, in a manner consistent with observations reported in 32,43 . For Re g = 1000, holes develop between the transverse waves, and the sheet fragments through hole collision and ligament breakup. The number of ligaments formed following hole collisions closely matches the number of transverse waves observed in the preceding gas flow structures.
The droplet size distribution is analyzed by comparing simulations with and without the Manifold Death (MD) algorithm 52 , which imposes a preset cutoff film thickness that is six times larger than the finest grid size. The dominant droplet diameter in the histogram remains unchanged with the use of the MD algorithm, while a reduction in the number of small droplets is observed. This suggests that the MD algorithm effectively suppresses the formation of small fragments likely caused by numerical artifacts. Finally, a comparison with two commonly used distribution functions in the literature, the log-normal and two-parameter gamma distributions, shows that the log-normal distribution provides a better fit to the droplet size distribution.
This study provides evidence that the air-blast atomization of thin liquid sheets is governed by the underlying gas flow structures, which can ultimately inform and enhance predictive models for liquid sheet fragmentation.
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0278367