The atomising pulsed jet

Abstract

Direct numerical simulations of the injection of a pulsed round liquid jet in a stagnant gas are performed in a series of runs of geometrically progressing resolution. The Reynolds and Weber numbers and the density ratio are sufficiently large for reaching a complex high-speed atomisation regime but not so large so that the small length scales of the flow are impossible to resolve, except for a very small liquid-sheet thickness. The Weber number based on grid size is then small, an indication that the simulations are very well resolved. Computations are performed using octree adaptive mesh refinement with a finite volume method and height-function computation of curvature, down to a specified minimum grid size $\varDelta$. Qualitative analysis of the flow and its topology reveals a complex structure of ligaments, sheets, droplets and bubbles that evolve and interact through impacts, ligament breakup, sheet rupture and engulfment of air bubbles in the liquid. A rich gallery of images of entangled structures is produced. Most processes occurring in this type of atomisation are reproduced in detail, except at the instant of thin sheet perforation or breakup. We analyse droplet statistics, showing that as the grid resolution is increased, the small-scale part of the distribution does not converge, and contains a large number of droplets close in order of magnitude to the minimum grid size with a significant peak at $d = 3\varDelta$. This non-convergence arises from the numerical sheet breakup effect, in which the interface becomes rough just before it breaks. The rough appearance of the interface is associated with a high-wavenumber oscillation of the curvature. To recover convergence, we apply the controlled ‘manifold death’ numerical procedure, in which thin sheets are detected, and then pierced by fiat before they reach a set critical thickness $h_c$ that is always larger than $6 \varDelta$. This allows convergence of the droplet frequency above a certain critical diameter $d_c$, above and close to $h_c$. A unimodal distribution is observed in the converged range. The number of holes pierced in the sheet is a free parameter in the manifold death procedure; however, we use the Kibble–Zurek theory to predict the number of holes expected on heuristic physical grounds.

1. Introduction

Atomisation simulations have progressed at an amazing rate, however, the topic is still far from mature. As we show in this paper, the prediction of the breakup of liquid masses in a typical large-speed flow is marred by vexing numerical effects and profound physical uncertainties. New developments, such as better codes, rapidly increasing processing power and new numerical methods are however poised to mitigate the difficulties. In this paper we investigate the pulsed jet, a paradigmatic case of atomising flow inspired by diesel engine jets, although we do not aim at solving a particular applied problem or even any experimental configuration, but rather to investigate the potential and limitations of direct numerical simulation (DNS) of an atomising flow, and in particular, the ability of such DNS to reveal physically relevant properties of the flow through statistically converged numerical approaches. The emphasis here is on statistical convergence, rather than ‘trajectory’ convergence, since only the former is conceivable in very complex, irregular flows.

The study of such convergence already has a long history, despite rapid progress. We review specifically the history of round-jet atomising simulations and we refer the reader to reviews on the general topic of atomisation and its simulation (Gorokhovski & Herrmann 2008; Villermaux 2020). The first attempts at testing the convergence of the probability distribution function (PDF) of droplet sizes were those of Herrmann (2011). His simulations of a round jet used a variant of the level-set method and showed that the number of small droplets in the PDF was underestimated by coarse-grid simulations. This can be understood by the fact the the level-set methods tend to eliminate small droplets by ‘evaporating’ them. On the other hand, volume-of-fluid (VOF) methods keep too many droplets in an intermediate range around the grid size $\varDelta$ , and in that range the PDF of droplet sizes is overestimated. This can be seen, for example, in a study of the convergence of the droplet-size distribution in the round jet by Pairetti et al. (2020) where coarse grids overestimate the number of large droplets.

A graphical illustration of the contrasting effects of level-set and VOF methods on the droplet sizes in atomisation is shown on figure 1. Despite the inaccuracy of the distribution for small droplet diameters, both methods may converge as the grid is refined progressively. There is thus some hope that very large simulations will eventually produce converged distributions, a hope we want to explore in this paper.

Figure 1. An illustration of the outcomes for the numerical simulation of a thinning liquid sheet; panel (a) shows the initial configuration, before breakup. The other three images (b–d) show a schematic view of the outcome either in reality or in various types of numerical simulation. It is arbitrarily assumed that there are two breakup locations. Both the VOF and the level-set methods yield topology changes when the sheet thickness reaches the grid size. In (b) fragments larger than the grid size are obtained because of mass conservation in the VOF method. (c) In reality, the sheet thinning continues until much later than in the numerics, unless extremely fine grids are used. The final size of some of the droplets is then much smaller than in the VOF simulation. (d) The level-set or diffuse-interface methods on the other hand evaporate the thin parts of the sheet and loses much more mass.

Figure 2. The increase in two-phase round-jet grid resolution in time. The graph includes two simulations published only on the Gerris and Basilisk websites (and in other channels outside of academic journals) before 2017. The 2024 simulations are those reported in this paper.

Aside from the two investigations cited above, very few studies address the issue of PDF convergence, despite detailed analyses of aspects of the flow. Perhaps the first round, single-jet atomisation simulation in the conditions of a diesel jet was that of Bianchi et al. (2007) using the VOF method, followed by more detailed simulations by Ménard et al. (2007) using the combined level-set VOF method (CLSVOF). This was followed by other simulations using CLSVOF by Lebas et al. (2009), Chesnel et al. (2011) and Anez et al. (2019) and by simulations using the VOF method (Fuster et al. 2009). The latter studied both the coaxial jet cases in the so-called ‘assisted atomisation’ set-up and the conical jet. For the conical jet, Fuster et al. (2009) showed that the PDF changed drastically with grid resolution. The PDF never peaked as the droplet size was decreased towards the grid size. Until 2022, the most detailed (or highest-fidelity) published simulations of a round jet were performed by Shinjo & Umemura (2010) using a combination of level-set and VOF techniques with a ratio of jet diameter to grid size of 286. These simulations, together with those of Herrmann (2011) are the two outliers above the trend line in figure 2. In addition to revealing a wealth of details about ligament and sheet formation, and the perforation of sheets, Shinjo & Umemura have also shown distributions of droplet sizes, but without investigating explicitly their convergence. As in earlier simulations of the round jet, the droplet-size distribution has a single maximum close to the small droplet end of the spectrum, almost at the grid size. Studies by Jarrahbashi & Sirignano (2014) and Jarrahbashi et al. (2016) focused on a round, spatially periodic jet and analysed the effect of vortex dynamics on the development of instabilities along with the jet core. Studies by Zhang et al. (2020) and Pairetti et al. (2020) applied the VOF method with octree adaptive mesh refinement. As already mentioned, the latter paper investigated the distribution of droplet sizes showing some difficulties in reaching statistical convergence. Several authors including Torregrosa et al. (2020) and Salvador et al. (2018) have focused on the effect of injection conditions and turbulence. Khanwale et al. (2022) and Saurabh et al. (2023), using a diffusive interface method with octree adaptive mesh refinement, obtained the most refined simulations published at the time of writing with a special treatment of the refinement in thin regions. It should be noted that all of the round-jet studies cited above involve a variety of physical parameters and grid resolutions, despite the fact that they share similarities, such as moderate density ratios and liquid Reynolds numbers (defined below) in a range $5000 \leqslant Re_l \leqslant 13\,400$ that are attainable by DNS for single-phase flow. For example, the two octree studies of Zhang et al. (2020) and Pairetti et al. (2020) have gas-based Weber numbers (defined below) of respectively 177 and 417. As a general rule, it is better to select relatively small density ratios, Reynolds and Weber numbers to increase the likelihood of reaching convergence in computation. In this work, we thus chose the parameters in the lower range of $Re$ , $\textit{We}$ , as it correspondens to Weber number and density ratio for this type of flow. Similar attempts at obtaining convergence were realised in the simulations of Ling et al. (2017), investigating assisted quasi-planar atomisation after the experiment of Grenoble (Ben Rayana et al. 2006; Fuster et al. 2013). Ling et al. (2017) performed computations on four grids of increasing resolution. Despite the huge computational effort involved in that latter simulation, convergence was clearly not reached, both quantitatively as droplet-size distributions had a very narrow region of overlap, and qualitatively as clearly under-resolved structures in thin sheet-like regions were observed. One of our objectives in this work is to see what happens if ever finer grids are used, until very thin sheets are clearly resolved. In experiments as well as in some simulations it is clear that thin sheets are the site of weak spots (Lohse & Villermaux 2020). The presence of these weak spots as one of the mechanisms leading to atomisation forces the following somewhat sobering conclusion. Numerical simulations of diffuse-interface, level-set or VOF type can never be converged if thin sheets break or perforate only when they reach the grid size. Instead, a physical mechanism for weak spots or perforations must be present to nucleate holes in a manner independent of the numerics. Although such a mechanism may not be known yet, we suggest as a backstop procedure to define a critical sheet thickness $h_c$ beyond which perforations occur with relatively high frequency, that is, at the rate of one or more perforations per connected thin sheet region. Chirco et al. (2022) have called manifold death (MD) this procedure (this choice of words partially echoes that of Debrégeas et al. 1998). It should be noted that bags or membranes are observed at relatively low velocities in many atomisation configurations (droplets, jets etc.) and that at higher velocities, bags are not visible and perhaps non-existent. We note that the literature (Lasheras & Hopfinger 2000; Tolfts et al. 2023) distinguishes two kinds of atomisation, a bag or membrane-type atomisation at low Weber numbers and a fibre-type atomisation at high Weber numbers. The transition between the two behaviours was found by Tolfts et al. (2023) to lie between $\textit{We}_g=70$ and $\textit{We}_g=400$ . From the connection between thin sheets and the lack of statistical convergence one may conjecture that in the fibre-type regime statistical convergence could be achieved. This is however contradicted by the lack of convergence observed by Herrmann (2011) at $\textit{We}_g = 500$ . A definite answer to this argument may have to wait for a study similar to the one in this paper in the fibre-type regime, a study which would probably be even more expensive.

In what follows we first describe the characteristics of our test case. We then describe our numerical method, including the MD procedure. Then we continue with our results that are both qualitative and quantitative. Qualitatively our results reveal new processes of numerical rupture. Quantitatively they allow an analysis of the droplet-size distributions or PDF. A discussion of the phenomenology of hole expansion in a thin sheet is given in Appendix A. We end with a conclusion including perspectives and discussion.

2. Mathematical model and numerical method

Our mathematical model is based on the mass and momentum conservation equations for incompressible and isothermal flow, i.e.

(2.1) \begin{equation} \nabla \cdot \textbf { u} = 0, \end{equation}

(2.2) \begin{equation} \frac {\partial \rho \textbf { u}}{\partial t} + \nabla \cdot (\rho \textbf { u}\textbf { u}) = - \nabla p + \nabla \cdot \left (2\mu \mathbf {D}\right ) + \textbf { f}_\sigma , \end{equation}

where $\textbf { u}(\textbf { x},t)$ is the velocity field and $p(\textbf { x},t)$ is the pressure field. The tensor $\mathbf {D}$ is defined as $ ({1}/{2}) [\nabla \textbf { u} + (\nabla \textbf { u} )^T ]$ . The density and viscosity of the flow are noted as $\rho$ and $\mu$ , respectively. The last term on the right-hand side of the Navier–Stokes equation (2.2) represents the surface tension force

(2.3) \begin{equation} \textbf { f}_\sigma = \sigma \kappa \textbf { n} \delta _s, \end{equation}

which depends on the surface tension $\sigma$ and the interface shape, particularly on its curvature $\kappa$ and unit normal vector $\textbf { n}$ . The Dirac distribution $\delta _s$ indicates that the force only acts at the free surface. We consider the phase distribution function $c(\textbf{x},t)$ that takes value unity in the reference phase and null outside of it. The transport equation for this function is

(2.4) \begin{equation} \frac {\partial c}{\partial t} + \nabla \cdot (c\textbf { u})= c\,\nabla \cdot \textbf { u}. \end{equation}

In the context of incompressible flow, the right-hand side above is equal to zero. The $c$ function implicitly defines the interface at its discontinuity surface, defining also the $\delta _s$ , $\textbf { n}$ and $\kappa$ fields, e.g. $\nabla c = \textbf { n} \delta _s$ .

The VOF method represents the evolution of $c$ using the piecewise linear interface capturing method of Hirt & Nichols (1981) and DeBar (1974). In this context, the mean value of the colour function on a cell is

(2.5) \begin{equation} f = \frac {1}{| \Omega | }\int _\Omega c(\textbf { x},t)\,{\rm d}V, \end{equation}

where $| \Omega |$ is the volume of the cell $\Omega$ . Then $f$ is the volume fraction of the reference phase in the cell. The mixture properties of the cell may then be computed by arithmetic averages

(2.6) \begin{equation} \rho _\Omega = f \rho _l + (1 - f)\rho _g, \qquad \mu _\Omega = f\mu _l + (1 - f)\mu _g. \end{equation}

Spatial discretisation of the above equation is realised on a network of cubic cells obtained by a tree-like subdivision of an initial cubic cell of size $L_0$ . The subdivision is realised using a wavelet-based error estimate and is adapted dynamically as the simulation progresses. At all times the maximum subdivision level is a fixed number $\ell$ and the smallest grid size is

(2.7) \begin{equation} \varDelta _\ell = 2^{-\ell } {L_0}. \end{equation}

From now on, we consider the algebraic equations derived from applying the finite volume method on each cell. In this context, the approximate projection method by Chorin (1968) can be used to solve the coupling between (2.1) and (2.2), considering that the velocity $\textbf { u}$ is staggered in time with respect to the volume fraction $f$ and the pressure $p$ . The discrete set of equations can then be expressed as

(2.8) \begin{equation} \frac {\rho \textbf { u}^{*} - \rho \textbf { u}^{n}}{\Delta t} + \nabla \cdot \left (\rho ^{n + \frac {1}{2}}\,\textbf { u}^{n+1/2}\textbf { u}^{n+1/2}\right ) = \nabla \cdot [\mu ^{n + \frac {1}{2}}(\mathbf {D}^n + \mathbf {D}^*)] + \textbf { f}_{\sigma }^{n+1/2}, \end{equation}

(2.9) \begin{equation} \nabla \cdot \left (\frac {\Delta t}{\rho ^{n+\frac {1}{2}}} \nabla p^{n + \frac {1}{2}}\right ) = \nabla \cdot \textbf { u}^*, \end{equation}

(2.10) \begin{equation} \textbf { u}^{n+1} = \textbf { u}^* - \frac {\Delta t}{\rho ^{n+\frac {1}{2}}} \nabla p^{n + \frac {1}{2}}. \end{equation}

In the above the expression, $\nabla \cdot (\rho ^{n + \frac {1}{2}}\,\textbf { u}^{n+1/2}\textbf { u}^{n+1/2} )$ must be interpreted as the use of a predictor–corrector scheme for advection of the velocity field. Both the predictor and the corrector use the Bell–Collela–Glaz scheme (Popinet 2003, 2009), and involve two projections. A full description of the predictor corrector scheme with projection is best found in the `literate’ source code at http://basilisk.fr/src/navier-stokes/centered.h. In addition to these equations, $f$ is advanced in time using the split-volume-fraction advection scheme of Weymouth & Yue (2010). We use the semi-implicit Crank–Nicholson time stepping to compute the diffusive flux due to the viscous term in (2.8). An important aspect of the method is the numerical approximation of the surface tension force $\textbf { f}_{\sigma }$ . We use the continuous surface force well-balanced formulation $\textbf { f}_{\sigma } = \sigma \kappa \nabla c$ for (2.3), where we discretise $\nabla c$ at cell faces with the same scheme employed to compute the pressure gradient $\nabla p$ . This is useful to reduce spurious currents provided curvature is accurately computed, as explained by Popinet (2009, 2018). The interface curvature is computed using a complex method described in Popinet (2009). When the grid resolution is sufficient, second-order stencils based on height functions are used to compute curvature. In some large-curvature cases the height functions are not defined and fits by quadratic functions are used instead. Equations (2.9) and (2.10) are the projection steps that will ensure continuity for the velocity field at time step $n + 1$ . The mesh is adapted dynamically, by splitting the grid cells to eight smaller cells whenever the estimated local discretisation error exceeds a set threshold (see Appendix B). Details of the procedure can be found on the web site basilisk.fr and, in particular, in the documented free code available at http://basilisk.fr/src/examples/atomisation.c.

In all the data reported in this paper, the velocities are expressed in multiples of the jet velocity $U$ . All the numerical values for time are made dimensionless by $\tau _0 = L_0/(3 U)$ and all the reported lengths are made dimensionless by $L_0/3$ . However in the theoretical discussions the dimensional values are used unless otherwise indicated. To avoid confusion, we note with a $^*$ exponent the dimensionless values, so that

(2.11) \begin{equation} t^* = t/\tau _0 \qquad x^* = 3 x / L_0 \qquad \varDelta ^* = 3 (2^{-\ell }). \end{equation}

In addition to the aforementioned model and procedure, we implement the MD procedure as described by Chirco et al. (2022). This procedure involves the following two primary steps. (i) The identification of thin liquid sheets with a thickness of $h_c$ or less, achieved through a local integration of the characteristic function $c$ . This integration results in the calculation of the signature of a quadratic form, based on the local main inertia moments. (ii) The creation of holes within the so identified thin sheets, using a probabilistic approach. This procedure is applied at a given frequency, defined by the user within this method. From a physical standpoint, this technique aligns with the weak spot model, wherein the likelihood of hole formation is exceedingly low for $h \gt h_c$ and substantially increases, rapidly approaching certainty, for $h \lt h_c$ . Note that this procedure distinguishes thin sheets from other shapes such as ligaments, only the thin sheets are detected and perforated. The specific implementation of the MD method we used is called the signature method. The complete MD code with the signature method can be accessed at http://basilisk.fr/sandbox/lchirco/signature.h. In the present case, hole punching is attempted every time interval of $\tau _m=0.01 \tau _0$ , that is, the raw dimensionless breakup frequency is $\tilde f_b^* = \tau _0/\tau _m = 100$ . At these instants $N_b=200$ holes are attempted. We loosely term the pair of parameters $(\tilde f_b^* ,N_b)$ the breakup frequency. The breakup frequency chosen in this study is probably too high, as will appear in the discussion below. However, it is chosen to maintain the same parameters as in other studies, and possible improvements are discussed in the conclusion section. As we show in Appendix A, the characteristic time for sheet thinning is $\tau _c$ given by

(2.12) \begin{equation} \tau _c = \frac D U \left ( \frac {\rho _l }{\rho _g} \right )^{1/2}, \end{equation}

where $D$ is the jet diameter. From the definition of these quantities, we can obtain the breakup attempt frequency in units of the sheet thinning characteristic time $f_b = \tau _c/\tau _m$ , i.e.

(2.13) \begin{equation} f_b = \frac {3 \tilde f_b^* D}{L_0} \left ( \frac {\rho _l }{\rho _g} \right )^{1/2} , \end{equation}

that is, with the parameters of our simulation

(2.14) \begin{equation} f_b = \frac {50}{3} \left ( \frac {\rho _l }{\rho _g} \right )^{1/2} . \end{equation}

The numerical value is $f_b \simeq 87.9$ , ensuring that punching attempts are sufficiently frequent as not to `miss’ transient thinning sheets. It is large enough to eliminate numerical breakup events but is likely to be too high as discussed below.

Using this numerical method, we analyse the atomisation of a circular jet injected at average velocity $U$ in a gas-filled cubical chamber of edge length $L_0$ . We consider incompressible, isothermal flow. The dimensionless groups most relevant for this problem are the gas-based Weber number, the liquid-based Reynolds number and the ratios

(2.15) \begin{equation} We_g = \frac {\rho _g U^2 D}{\sigma }, \qquad {Re}_l = \frac {\rho _l U D}{\mu _l}, \qquad \rho ^* = \frac {\rho _l}{\rho _g}, \qquad \mu ^* = \frac {\mu _l}{\mu _g}, \qquad \frac {L_0}D, \end{equation}

where $\rho _l$ and $\rho _g$ are the densities of liquid and gas, respectively. The viscosity coefficients are noted by $\mu _l$ and $\mu _g$ , and $\sigma$ is the surface tension coefficient. The boundary condition on the injection plane $x = 0$ imposes the no-slip condition everywhere except on the liquid section ( ${y^2 + z^2} \lt D^2 / 4$ ), where the injection velocity is

(2.16) \begin{equation} u_x(t) = U\,\left [ 1 + A_p\sin \left (\omega _{p} U t/D \right )\right ]. \end{equation}

On the remaining cube sides, we allow free outflow: $\partial _n \textbf { u}_\Gamma = 0$ and $p_\Gamma = 0$ .

The dimensionless parameters of our simulation are given in table 1. These numbers are in the low range of the dimensionless numbers used in the above cited work. They are much smaller than the numbers of the spray G case of the engine combustion network (see Duke et al. 2017) used by Zhang et al. (2020). The Reynolds number and density ratio are identical to those of Ménard et al. (2007) and of the same order as those of Herrmann (2011). The equality of the gas and liquid Reynolds numbers stems from the fact that the viscosity and density ratios are equal, so the ratio of kinematic viscosities is one. However, while the moderate density ratio is realistic in the context of diesel jets, the viscosity ratio is not and would be much larger in a realistic set-up, which would lead to much larger values of ${Re}_g$ that are expensive to resolve by DNS. This motivates the choice of viscosity ratio. We note that the same choice of equal viscosity and density ratios and moderate $\textit{We}_g$ was made by Ling et al. (2017) who also studied atomisation and PDF convergence in a simplified set-up. The Reynolds number is somewhat larger than that of Shinjo & Umemura (2010) while the gas Weber number $\textit{We}_g$ is smaller than that of all the other papers in the literature. The rationale for selecting such a small Weber number is twofold. First the simplest prediction for the droplet size, based on the ratio of Bernoulli pressure $\rho _g U^2$ to surface tension, is $d = D We_g^{-1}$ . Indeed if a purely inviscid flow with a velocity jump is considered, the cutoff wavelength for the Kelvin–Helmholtz instability scales as $D We_g^{-1}$ . The wavelength thus obtained is a natural scale for the subsequent formation of sheets and ligaments. The number of grid points in the diameter of a droplet of this diameter $d$ is then $\textit{We}_\varDelta ^{-1}$ , where

(2.17) \begin{equation} \ {We}_l = {\rho _g U^2 l }/{\sigma } \end{equation}

is the Weber number based on length scale $l$ and $\varDelta$ is the grid size. Similarly,

(2.18) \begin{equation} \ {Oh}_l = \mu _l (\sigma \rho _l l)^{-1/2} \end{equation}

is the Ohnesorge number based on scale $l$ . The values of $\textit{We}_\varDelta$ for two of the grids we use are given in table 2. We also give the value of the Ohnesorge numbers ${Oh}_{\varDelta }$ and ${Oh}_{h_c}$ , which characterise the interplay of viscous and surface tension effects. In particular, they typify the regime of the Taylor–Culick rims (Song & Tryggvason 1999; Savva & Bush 2009) occurring at the edge of a sheet of thickness $\varDelta$ or $h_c$ . The $Oh$ number of the thin sheets at the moment of hole formation is relatively large, implying a viscous Taylor–Culick flow and the absence of unsteady fragmentation as in Wang & Bourouiba (2018) and Kant et al. (2023).

Table 1. The dimensionless numbers characterising our simulation. The Reynolds number ${Re}_l$ based on the liquid is rather moderate. The density and viscosity ratios are identical, which implies that ${Re}_g = \ {Re}_l$ .

Table 2. Characteristic numbers related to the grid size. The level $\ell$ is defined in (2.7).

Table 3. Characteristic numbers based on the critical sheet thickness in our simulations. The MD level $m$ is defined in the text.

The second important consequence of selecting the relatively small value $\textit{We}_g =200$ is that it makes the simulation sit on the boundary between membrane type and fibre type defined in Lasheras & Hopfinger (2000). Membrane type in that later paper refers to the formation of bags and sheets. Thus, the prevalence of thin sheets discussed in this paper may be less marked at higher Weber numbers, the sheets being replaced by fibres, and also less important at smaller Weber numbers where the flow may be less unstable.

3. Results

3.1. Uncontrolled perforation: numerical sheet breakup

Figure 3 shows a global view of the jet evolution (see also supplementary movies 1–6 available at https://doi.org/10.1017/jfm.2025.218). The dynamics is initialised with a slightly penetrated jet at $t=0$ . At $t=0.1$ , a mushroom-like structure starts rolling up. The annular structure stretches and extends in the form of a `corona flap’ that develops behind the mushroom head. A rim forms at the edge of this flap then detaches. It is visible at $t=0.2$ as a corrugated and partially fragmented annular ligament. As the jet evolves, pulsation results in a periodic series of corona flaps seen at $t=0.8$ . These flaps interact with the ligaments coming from the mushroom head. At $t=2$ , the jet has evolved further and we can identify an interval in the core structure where the effect of pulsation is lost. Several mechanisms of droplet formation like sheet rupture, ligament breakup and their interactions with the jet core and corona flaps are seen. At $t=3$ , the jet is fully developed and we see long axial ligaments with a marked velocity gradient made visible by the colour change along the axis, circular ligaments, corona fingers merging and sheet rupture. A large number of droplets and ligaments are produced in the mushroom tip region. Most of the phenomena illustrated above are also observed with controlled perforation, except the initiation of sheet rupture and annular ligament detachment.

Figure 3. The advancing pulsed jet at various time instants $t$ and level $\ell =14$ . The fluid interface is coloured by the axial velocity and the background is coloured by the vorticity. The background also shows the mesh refinement. (a) The pulsed jet develops a mushroom head and a rim. In (b) the rim detaches. (c) Development of flaps coming from the sinusoidal pulsation. (d) A jet entering in a regime where the effect of pulsation is lost at the mushroom head. (e) A fully developed jet and a rich spectrum of droplets and ligaments. Since $Re_g$ is rather low at $5800$ , there is relatively little vorticity away from the interface unlike in the case of Kant et al. (2023).

Figure 4. View from the inlet at time $t=3.04$ showing the inner region of the central core liquid jet. Plot (a) is coloured by the curvature showing the encapsulation of gas bubbles identified by the negative curvature (blue) in the liquid core encircled in the black circle. The droplets have positive curvature (red). The entrained bubbles travel with the core jet velocity and could also result in the formation of a few compound droplets during atomisation or provide a physical breakup mechanism for thin sheets. Plot (b) is the same as (a) but coloured by the axial velocity. The simulation corresponds to level $\ell =14$ with the MD method applied at level $m=13$ .

Figure 4 shows a contrasting view of the jet with the camera direction aligned with the jet axis and the view positioned at the inlet. A full evolution of the jet, as seen from the inlet, is available in supplementary movie 6. In figure 4(a) the objects are coloured by curvature, so small droplets are red and small bubbles, with negative curvature, are blue. The bubbles have been trapped in the bulk or core of the jet or inside other liquid masses. The likely mechanism for the formation of a bubble is liquid mass collision, for example, droplets and wavelets of fluid mass impacting on the jet core. Many such droplet or ligament impacts are seen in supplementary movie 1. A surprising aspect of the display is the relatively large number of droplets and bubbles seen. In figure 4(b) the objects are coloured by axial velocity, as with supplementary movies 1 and 3. A further interesting aspect of this view is to show the large dispersion in the axial velocities of the droplets.

At this point in the exposition it may be useful to define the number distribution and the PDF and to recall the connection between the two. We consider a fixed instant of time $t_0$ . We partition the interval of droplet sizes $(0,d_M)$ into a discrete set of $M$ intervals $[d_i, d_{i+1} = d_i + \delta _i]$ . The index $i$ varies from 0 to $i_m$ . The intervals are geometrically distributed, with $d_{i+1} = 1.047 d_i$ , so $\log(d_i)$ is evenly distributed in log space. The number of droplets with diameter in the interval $[d_i, d_i + \delta _i]$ is $n_{d,i} \delta _i$ . The PDF $p(d)$ is proportional to the expected value $\langle n_{d,i} \rangle$ as

(3.1) \begin{equation} p(d_i) = \frac {1}{\mathcal N} \langle n_{d,i} \rangle , \end{equation}

where the normalising factor is

(3.2) \begin{equation} {\mathcal N} = \sum _{i=0}^{M} \langle n_{d,i} \rangle \delta _i. \end{equation}

We show two kinds of plots, one with the number $N(d_i)$ of droplets in each bin $i$ , that is,

(3.3) \begin{equation} N(d_i) = n(d_i) \delta _i , \end{equation}

and the other with the PDF $p(d)$ . Because of the difference between the definitions of $N$ (which is bin-size dependent) and $n$ and $p$ (which are both bin-size independent), there is a scaling relation

(3.4) \begin{equation} N(d) \sim {\rm d} p(d) \sim {\rm d} n(d) \end{equation}

between the two distributions. (We note that in the literature both $n$ and $N$ are reported with the name ‘number distribution’.) The pulsed jet process is both deterministic (except for the random character of sheet perforations) and transient, so the value of $n_{d,i}$ obtained at $t=t_0$ in a simulation is well defined. The expected value $\langle n_{d,i} \rangle$ differs very little from the numerically obtained value $n_{d,i}$ . Moreover, in such a deterministic process no difference is expected between two realisations, that is, two identical simulations should yield identical droplet counts. On the other hand, in a jet with a random upstream injection condition instead of a pulsed injection condition, the simulation is not deterministic and the value of $n_{d,i}$ obtained in a simulation at some time $t$ is fluctuating. The transient nature of the jet is also important. In a steady-state simulation (which should be obtained some time after the jet starts flowing out of the simulation domain), initial perturbations are amplified chaotically and the system is again probabilistic. Thus, our pulsating initial condition simplifies somewhat the analysis by yielding a deterministic instead of a probabilistic number density.

Figure 5. (a) The bubble size distribution for the image shown in figure 4. The number $N$ is defined in (3.3). (b) A 2-D histogram for the same image showing bubble size distribution on the transverse coordinates. The jet is advancing in the x direction, which is the axial direction. The solid white circle is the inlet region. The curly dotted white lines are drawn to highlight pale patches to indicate that some bubbles also exist outside the jet core, indicating possible compound drops.

Figure 5 shows the number distribution of bubbles. The ordinate $N$ is the number $n_i$ of droplets defined in the text. We see that above the critical thickness, a distribution $N \sim r^{-\gamma }$ is seen. This distribution may have consequences in the sheet rupture and, hence, the droplet-size distribution. We have discussed the physical implications of bubble size distribution in Appendix A. In figure 5(b) we plot a two-dimensional (2-D) distribution of the bubble size on a y–z plane. The white circle represents the inlet. The plot shows that most of the bubbles are near the outer boundary of the jet core where we see wavelet interactions and drop impacts the most (figure 3(c) and supplementary movie 1). The distribution also shows eight characteristic high bubble-density zones. The occurrence of these zones and what decides their number is beyond the scope of the current paper. We also see some bubbles present outside the jet core encircled in dotted white lines. This indicates that we do have a relatively small number of compound droplets.

Figure 6. (a) The droplet-size number distribution ( $N$ is defined in (3.3)) and (b) the probability density function of the droplet diameter (defined in (3.1)) at various time instants. The plot shows that the PDF has converged in time $t \sim 2.5$ . This time convergence determines the choice of the end time of the simulation at t = 3.5. Both (a) and (b) have 200 bins. The simulation corresponds to level 14 and the vertical dashed line represents the grid size.

We now discuss the statistical properties of the jet, in particular, the droplet-size distribution in the case of uncontrolled sheet perforation or ‘no-MD’ case. As the jet advances, the droplet-size distribution evolves. This is shown in figure 6. One can see that the droplet-size distribution has started to converge at $t=1.8$ and is approximately converged at $t=2.8$ . The histogram thus only shifts in the vertical direction (the number of droplets), while maintaining the same qualitative nature. This histogram is bimodal with two peaks or local maxima. There is a large amplitude peak corresponding to a smaller droplet size (referred to as peak 1) and a smaller amplitude one corresponding to a larger droplet size (referred to as peak 2).

Figure 7. (a) The droplet-size distribution and (b) the probability density function of the droplet diameter $d$ at the final time $t=3.5$ for various grid resolutions. The grid sizes are shown with vertical dashed lines and $\varDelta _{\ell }$ is defined as in (2.7). The inset in (a) shows the converged tail region. All the plots are done with 200 bins. A Pareto distribution $N(d) \sim d^{-2}$ is added to compare with the converged region of the distributions.

We investigate the statistical convergence of the droplet-size distribution by plotting in figure 7 the histogram at $t=3.5$ when the jet is fully developed. We plot various grid resolutions or levels on that same figure. The highest resolution in the no-MD case corresponds to the $\varDelta _{14}$ line. This implies $D/\varDelta = 920$ grid points per initial diameter. Figure 7 has two important characteristics: (i) the histograms at various resolutions are qualitatively similar, that is, bimodal with a major peak 1 at small scales and a minor peak 2 at large scales; and (ii) despite such high resolution, it is clear that there is no convergence for the droplet-size distribution. However, the tail region of the distribution, to the right of peak 2, shows some degree of convergence. In that tail region a Pareto distribution $n(d) \sim d^{-2}$ is seen reminiscent of such distributions found in other contexts (Balachandar et al. 2020; Pairetti et al. 2021). The dependence of the two peaks with the grid size is shown explicitly in figure 8(b) where we plot the ratios $d_i/\varDelta _\ell$ corresponding to the peaks as a function of the maximum level $\ell$ . We observe that the ratio for peak 1, $d_1/\varDelta _\ell$ , remains constant around a value of 3, while for peak 2, $d_2/\varDelta _\ell$ weakly increases with $\ell$ . The grid dependence of these two peaks is synonymous with the absence of statistical convergence. We discuss possible explanations for the behaviour of these two peaks below and, in particular, the relation between peak 1 and the curvature ripples appearing just before sheet rupture. We also discuss a mechanism and a possible scaling for peak 2 in Appendix A. Since the grid dependence behaviour is maintained at all resolutions, the constant scaling of the diameter $d_1$ of peak 1 with the grid size implies that no matter how much we refine the grid, the distribution of droplet sizes is grid dependent.

Figure 8. (a) Schematic of the bimodal droplet-size distribution. The arrows indicate the logarithmic distance from the grid size $\log d_i/\varDelta = \log d_i - \log \varDelta$ . (b) The values of $d_i/\varDelta$ as a function of the maximum grid refinement level $\ell$ for the distribution of figure 7. The dashed line is at $3\varDelta$ , implying the droplets at peak 1 have a grid-dependent diameter of $d_1 \sim 3\varDelta$ .

Figure 9. Droplet impacts on the frontal liquid sheet. The points of interest are indicated by the arrows. Droplet impacts can result in holes with characteristic ligaments as seen at $t=3.10$ . In some cases droplets coalesce into the sheet seen at $t=3.12$ and $t=3.14$ . At $t=3.16$ , the sheet is ruptured but the droplet is still identifiable. All simulations are for the maximum level $\ell = 13$ . The MD method is not applied. The interface is coloured by axial velocity. The artefacts or corrugated surfaces particularly visible at t = 3.02 and t = 3.10 at the top right are not aliasing artefacts (due to an approximate isosurface interpolation) but are representative of the curvature oscillations seen in the next figures.

Figure 10. Appearance and evolution of the curvature ripples on the interface. The interface is coloured by the interface curvature. These ligaments are coloured a darker red and the interfaces on the thin sheets are closer to white. Ligament rupture is encircled in blue. The rupture or ligament pinch off appears in a darker red colour. Sheet rupture, encircled in green, displays curvature ripples in the form of red–blue oscillations. At time $t=2.82$ , the green-circled rupture region displays a grid-dependent ligament network that has evolved from these ripples. These images correspond to a level $\ell = 14$ simulation. (See also supplementary movie 5.)

The mechanism leading to peak 1 is related to the perforation of thin sheets. This perforation is observed to occur in two contrasted ways. A simple mechanism is the impact of small droplets on the liquid sheet. The droplets are formed by breakup previously occurring elsewhere in the simulation. The impacts yield in some cases the formation of an expanding hole, as shown on figure 9. We note that such impacts have already been observed in previous studies such as those of Ménard et al. (2007) and Shinjo & Umemura (2010). Not all droplet impacts create perforations. In figure 9, at time $t=2.94$ , we track the droplets indicated by arrows. We see that droplet impacts often create holes with a characteristic finger-like ligament inside the hole. Impacts are shown by arrows at time $t=3.10$ in figure 9. In the lower panel at $t=3.12$ and $t=3.14$ some droplet impacts do not result in holes but instead droplets merge with the sheet. This hypothesis was confirmed by repeated observations of the sheets at various angles and by visioning supplementary movie 5. The red arrow from $t=3.14$ and $t=3.16$ captures an interesting moment, where the droplet impact and sheet rupture happen simultaneously and it is unclear if the hole was created by the numerical rupture process discussed below or by droplet impact.

Another mode of thin sheet perforation, much more complex and purely numerical, is illustrated on figure 10 as well as supplementary movie 5. The curvature colouring of the interface helps us identify the perforation spots. We see high-frequency oscillations of the curvature field in the form of alternating red and blue colours. These indicate that the ripples originate at the fluid interface before the sheet rupture happens. Finally, at $t=2.82$ , we see a dense network of small diameter ligaments, scaling with grid size. Another mechanism that may contribute to peaks 1 and 2 is the rupture of ligaments, also shown on figure 10. Small diameter ligaments have large curvature and are thus easily identified by their colour. They break physically and not numerically by the Rayleigh–Plateau instability. The origin of these ligaments themselves is often an expanding hole whose rim collides into other rims.

Figure 11. The fully developed jet at $t= 2.8$ . The inset shows a zoom in at the sheet rupture spot showing the curvature ripples in the weak spot about to be punctured and the resulting ligament network. This ligament network eventually produces grid-dependent droplets. The simulation shown here is at level $\ell =14$ .

Figure 12. The fully developed jet at $t= 3.03$ . The inset shows a zoom in at the sheet rupture spot showing the curvature ripples in the weak spot about to be punctured. The simulation shown here is at level $\ell =13$ .

Figure 13. The fully developed jet at $t= 3.01$ . The inset shows a zoom in at the sheet rupture spot showing the curvature ripples in the weak spot about to be punctured. The simulation shown here is at level $\ell =14$ .

Figure 14. The $(14,13)$ jet at $t= 3.1$ with MD applied. The basilisk level $\ell =14$ and MD level is $13$ . The insets and the boxes indicate droplets that are near the maximum size of the distribution in the converged region.

We note that the high-frequency pressure oscillations are due in part to the use of the height-function method to compute curvature. The height-function method ceases to work when the sheet is too thin. However, other methods will also fail as represented in figure 1. To investigate the small-scale ligament networks, we zoom around such a breakup event in figure 11. We see that a weak spot develops and curvature ripples appear at $t=2.78$ as soon as the sheet reaches a thickness of a few times the cell size, as shown in the inset. These curvature ripples then give rise to ligament networks inside the expanded hole eventually leading to the grid-dependent droplets. To summarise, as the sheet reaches a thickness of order $ 3 \varDelta$ , the rupture mechanism is as follows: thin sheet regions $\gt$ curvature ripples $\gt$ ligament networks and tiny droplets $(d \sim 3\varDelta )$ . Since this sequence of events is inherent to the numerical breakup of thin sheets, it is impossible to escape its occurrence by increasing the grid refinement, although refinement can delay it. The numerical sheet breakup identified here is the direct cause of the lack of convergence of peak 1. Indirectly, the formation of droplets and ligaments of a given size in a non-converged manner leads to all the droplet size counts being unconverged as these other size droplets are formed through coalescence and breakup events from the peak 1 droplets and similarly scaled ligaments.

We also note that the point of view and region of figure 9 were selected to show a large number of drop impact situations, which are however scarcer than the numerical breakup of sheets.

To give an overall idea of the evolution of the jet, we show the fully evolved jet at $\ell = 13$ in figure 12 and the same fully evolved jet but at $\ell = 14$ in figure 13. The mist of tiny droplets near the mushroom head and near the inlet flaps is more pronounced in the $\ell = 14$ case than in the $\ell = 13$ case. The $\ell =14$ case retains more large drops and deformed ligaments at mid-length compared with the $\ell =13$ case. The inset zoom shows the numerical sheet rupture. The ligament network structure is more clearly seen at $\ell = 14$ than at $\ell =13$ .

3.2. Controlled perforation by the MD method

We now present the results of atomisation when we apply the MD method. Here we adopt the following conventions.

1. (i) The octree maximum refinement level $\ell$ defines the finest grid size as given in (2.7).

2. (ii) The MD level $m$ implies that the critical thickness for punching holes is given by $h_c = 3 \varDelta _{m}$ .

As an example, the notation $(\ell ,m)=(13,12)$ would mean that the finest cell size in the simulation is $\varDelta _{13}$ and the critical thickness at which holes are being punched is $h_c = 3\varDelta _{12} = 3\times (2\varDelta _{13}) = 6 \varDelta _{13}$ . This convention for the definition of the ‘MD level $m$ ’ maintains consistency with the comments in the code available at http://basilisk.fr/sandbox/lchirco/signature.h.

Figure 15. The probability density function at various times $t$ for the droplet-size distribution when MD is applied. The PDF has converged in time at $t=2.5$ . (a) The PDF at maximum level $\ell =13$ and MD level $m=12$ and (b) the PDF at maximum level $\ell =14$ and MD level $m=12$ . The critical hole punching thickness $h_c = 3 \varDelta _{m=12}$ is shown as the black dashed line and is the same for both plots.

The extent of the integrated region in the quadratic form computation leads to the correspondence $h_c = 3 \varDelta _m$ . The overall evolution of the jet shows important variations with the no-MD case. Comparing figure 14 MD $(14,13)$ to figures 12 ( $\ell =13$ ) and 13 ( $\ell =14$ ), we note that (i) the mist of tiny droplets near the mushroom head has disappeared in figure 14, (ii) the mushroom-head corona-flap sheet is more affected by holes and breakups and is closer to figure 13 $\ell =14$ than to figure 12 $\ell =13$ , (iii) the cloud of droplets at mid-length is similar to both no-MD cases but is more uniform and closer to the figure 12 $\ell =13$ case.

Figure 15 shows the evolution of the PDF of the droplet-size distribution with time. Similar to the no-MD plots, we see that here also the statistics are converged in time as soon as $t=2.5$ . When compared with figure 7, we see that the peak corresponding to the smaller droplet size has a smaller amplitude.

Figure 16. The droplet-size distribution from simulations using the MD method. Each curve is labelled $\varDelta _{\ell }-\varDelta _{MD=m}$ , where $\ell$ is the grid level and $m$ is the MD level. The dashed black vertical line represents $h_c = 3 \varDelta _{m}$ , that is, the critical thickness of punching holes. The solid vertical lines indicate the values of the converged diameter $d_c$ , corresponding to the smallest diameter above which the difference between subsequent distributions is not significant. All curves in (a) have fixed $h_c = 3\varDelta _{12}$ and all curves in (b) have fixed $h_c = 3 \varDelta _{13}$ . Note that the $\varDelta _{11}$ curve in (a) and the $\varDelta _{12}$ curve in (b) indicate the no-MD method and roughly correspond to respectively $h_c = 3\varDelta _{MD=12}$ and $h_c = 3 \varDelta _{MD=13}$ , given the breakup due to a lack of sheet resolution in VOF simulations. The shading highlights the departure from the converged region. We see that $d_c$ approaches $h_c$ as the grid size is reduced. The inset shows the important change in the typical sizes of droplets in the first peak above $h_c$ as MD is applied. Peak A corresponds to peak 2 in the no-MD distributions of figure 7 while peak B, with much smaller droplets, corresponds to droplets seen in the MD figure 14.

We show the number density $n(d)$ in figure 16. Figure 16(a) shows the histograms at various grid resolutions for a fixed MD threshold $h_c = 3 \varDelta _{12}$ . The distribution is again bimodal, but the peak corresponding to the smaller droplets is weaker. We also observe a clear converged region in the interval $[d_c, d_{i_m+1})$ , where $d_c$ is the minimum converged equivalent droplet diameter or ‘converged diameter’ for simplicity. The extent of the converged region increases, that is, $d_c$ decreases as the grid size $\varDelta _l$ decreases. In figure 16(a) the vertical line marks $d_c$ for two consecutive grid sizes. We have shaded the onset of the non-converged region between two consecutive grid refinements for clarity.

In practice, the ‘converged diameter’ $d_c$ is determined from the comparison of a pair of distributions from simulations of two different grid sizes $\varDelta _l$ . We note this pair  $\varDelta _{l_1} - \varDelta _{l_2}$ . In figure 16(a) the orange vertical line marks $d_c$ for the pair { $\varDelta _{14} - \varDelta _{13}$ }. Similarly, the purple vertical line marks $d_c$ for the pair { $\varDelta _{15} - \varDelta _{14}$ }. Some additional information may be obtained by investigating a no-MD case. In such a case, there is a typical sheet thickness at breakup $h_e$ that is comparable to the controlled breakup for $h_c \simeq h_e$ . Indeed, the critical breakup film thickness $h_e$ for no-MD is $h_e \simeq 1.5 \varDelta$ , Then the above case $h_c = 3\varDelta _{12}$ with MD has a critical sheet thickness equal to, per expression (2.7), $h_e = 1.5 \varDelta _{11}$ in the no-MD case at level 11. Thus, the no-MD simulation at $\ell =11$ is a rough equivalent to the three simulations with MD level $m=12$ . The vertical black line is the converged thickness $d_c$ for the pair { $\varDelta _{13}, \varDelta _{11}$ } where the first simulation is MD and the second is no-MD. One sees on figure 16(a) that the converged diameters, indicated by the solid vertical lines, decrease and tend to approach the critical sheet thickness $h_c$ for MD breakup, represented by the dashed vertical line. We now show the droplet-size distribution and converged region for a smaller critical sheet thickness $h_c = 3 \varDelta _{13}$ . In figure 16(b) we see similar results but with only two instead of three levels, $l=14$ and 15. Since a simulation at $\varDelta _{15}$ is already quite expensive (it implies an equivalent of $35$ trillion cells for a uniform grid with $\varDelta _{15}$ ), we did not perform a third simulation at $l=16$ in that case. A zoomed-in view for the converged region $d_c$ for each consecutive pair is shown in Appendix E.

The closely packed peaks around $\log(d)=-3$ in figure 16(a), mostly visible for $(\ell ,m)=(15,12)$ , are related to a process that produces a narrow spectrum of small droplets of size $\bar d$ with $\varDelta _{15} \lt \bar d \lt h_c$ . These small droplets coalesce at least two or three times to produce additional small droplets of size $2\bar d, 3\bar d$ , etc. We are unsure of the mechanism producing these initial droplets of well-defined diameter $\bar d$ and why it appears most clearly for the pair ( $15,12$ ).

Figure 17. (a) Schematic of the comparison of two simulations and how it allows us to define the `converged diameter’ $d_c$ . It is seen that $d_c$ approaches $h_c$ as the simulations are refined. (b) Plot of $d_c/h_c$ versus $h_{c}$ corresponding to the $d_c$ values of figure 16. Along a fixed $h_c$ , the $d_c/h_c$ points move closer to unity as the maximum level is increased. The black squares corresponds to the no-MD case, while all circles are for the cases with MD.

The variation of the converged diameter $d_c$ as a function of $h_c$ and $\varDelta _\ell$ is shown explicitly in figure 17. Figure 17(a) illustrates schematically how the converged diameter approaches $h_c$ as $\varDelta _\ell$ is reduced. Figure 17(b) shows the data in the $d_c/h_c, h_c$ plane. Each data point in that plane corresponds to a simulation pair discussed above such as { $\varDelta _{14} - \varDelta _{13}$ }, and to a vertical line in figure 16. Following a vertical line in figure 17(b), one can see $d_c$ tending to $h_c$ as the grid is refined for a fixed $h_c$ .

Figure 18. Sheet perforation when the MD method is applied. Unlike figure 11, the holes are not preceded by large regions containing curvature ripples. Holes are punched at $t=3.01$ and expand with time as shown. The interface in the left part is coloured by curvature and on the right is coloured by the axial velocity. The simulation shown here corresponds to $(\ell ,m)=(14,13)$ .

Figure 18 shows a zoom in at a hole expansion that was initiated by the MD method. The phase modification region is a cube but in less than 5–10 time steps the hole takes a rounded shape due to surface tension and starts expanding. When these holes expand, we do not see the noisy curvature ripples that resulted in grid-dependent droplets in the no-MD simulation of figure 11. This qualitative improvement in sheet collapse helps to obtain grid convergence.

Figure 19. Further evolution of the holes of figure 18. The image is coloured by axial velocity. We see that near the rims of the main expanding holes, additional holes are punched, encircled in black. (Manifold death is applied.)

The MD procedure we use is tuned to a relatively large frequency of hole formation. This results in features such as those shown in figure 19. We see that once we detect a thin sheet and the MD process generates holes, these expand. However, it often happens that before the holes have expanded over the entire sheet more holes are punched just outside the rim. These kind of holes are identified as black circles in figure 19. This effect is also seen in supplementary movie 4 that shows the evolution of thin sheets with MD. Experimental photographs and the reasoning in Appendix A point to a smaller number of holes. This smaller number could be achieved by an improved MD algorithm that would involve a more realistic model of the physical process of hole formation.

Figure 20. Two distributions of droplet sizes at two MD levels, $m=12$ and $m=13$ , with the same ratio of the grid size to the MD threshold $h_c/\varDelta _\ell$ . We only plot the converged region, that is, $d\gt d_c(h_c)$ . The most refined simulation $(l,m)=(15,13)$ yields a distribution somewhat but not exactly similar to the less refined. Both distributions show a sudden increase in droplet number (a `ridge’) as the droplet size is decreased at the edge of the converged region. The vertical line shows the limit of convergence $d_c$ . Plot (a) shows the un-rescaled number distribution, while plot (b) shows the PDFs. Plot (a) is clipped to show the PDFs $p(d)$ for $d\gt d_c(m=12)$ , where $d_c$ is the edge of the converged region for $m=12$ . The Pareto distribution $n(d) \sim d^{-2}$ is also plotted.

Figure 21. Same as figure 20 but with a linear abcissa and a fit to gamma and log-normal distributions.

The final result of our manifold-death-enhanced study of droplet sizes is summarised on figure 20. In this figure we select two distributions of droplet sizes at two MD levels, $m=12$ and $m=13$ , and we keep the ratio of the grid size to the MD threshold $h_c/\varDelta _\ell$ constant so the level difference is always $\ell - m= 2$ . The two obtained distributions are similar to each other, with a sudden increase in droplet number at the edge of the converged region. Moreover, these plots show that a sufficiently thin grid produces a bump-shaped distribution with a clear maximum (this is the case at $m=13$ but not at $m=12$ ). This maximum is located at $d_{\textrm {max}}/D = 4.3 \, 10^{-2}$ in our case. It is interesting to try to understand the origin of this bump. By dimensional analysis the maximum should be expressed as

(3.5) \begin{equation} d_{max} = D f(h_c/D,N_1,\ldots, N_n), \end{equation}

where the $N_i$ are the dimensionless numbers of the problem other than $h_c/D$ . It is quite tempting to assume that in some range of $h_c/D$ and of the dimensionless numbers a proportionality relation holds. From the existing data it would be

(3.6) \begin{equation} d_{max} = 4.3 \, h_c. \end{equation}

We note that this is to our knowledge the only example of a converged bump-shaped (unimodal) distribution in simulations of jet atomisation at a large Weber number. We attempt to fit this bump shape to standard distributions in the literature, namely the gamma and the log-normal. Both distributions are acceptable in the large $d$ range, above 0.005, as seen in figure 21. The fitting distributions are defined as follows:

(3.7) \begin{align} \text {Log-normal: } \, P\left ( x ; \mu , \sigma \right ) &= \frac {1}{x \sigma \sqrt {2\pi }} \,\textrm {exp}\left [-\frac {1}{2}\left (\frac {\text {ln}( x) - \mu }{\sigma }\right )^2\right ] \,,\\ \nonumber \text {Gamma: } \, P\left ( x ; \alpha , \beta \right ) &= \frac {\beta ^{\alpha }}{\Gamma (\alpha )} x^{\alpha - 1} \textrm {exp}\left (-\beta x\right ) \,. \end{align}

4. Conclusions

We have studied an atomising, pulsed round jet at a succession of grid sizes with unprecedented resolution. Visual inspection of the resulting interface topology reveals a new numerical phenomenon preceding hole formation in the VOF method: high-frequency curvature oscillations followed by the appearance of a small-scale network of ligaments. While the experimental observation and theoretical analysis of droplet formation by expansion of holes is not new, this is the first time that these high-frequency oscillations of curvature are observed in VOF simulations. The inspection of droplet-size frequency distributions reveals the presence of two large peaks tied to the grid size. It is clear that the small- $d$ peak and probably the second one are related to the unphysical breakup.

To mitigate this effect, we apply the MD procedure of Chirco et al. (2022). The procedure forces hole formation when sheets reach a set thickness $h_c$ . Number frequency plots show convergence in a range of droplet sizes starting at a critical diameter $d_c$ . The extent of this converged range increases when either the grid size $\varDelta _\ell$ or the critical thickness $h_c$ are decreased. For a fixed $h_c$ , when we reduce $\varDelta _\ell$ , the converged diameter $d_c$ approaches the critical thickness $h_c$ . We thus recover the convergence properties observed by Chirco et al. (2022) and by Tang et al. (2023). We characterise them in the $h_c$ – $d_c$ plot of figure 17(b). We note that the statistical accuracy is much stronger in the present case since our pulsed jet produces three orders of magnitude more droplets than the phase inversion case of Chirco et al. (2022). We do not have a definite explanation for the lack of convergence below the diameter $h_c$ , that is, we are not sure why $d_c$ does not decrease below $h_c$ as $\varDelta _\ell$ is decreased. We can however offer several avenues for improvement of this feature and for future exploration. One is the smoothness of the MD procedure itself and the choice of its parameters such as the number of holes punched or the punching frequency $1/\tau _m$ . For the first time in VOF simulations of atomisation, a bump-shaped, converged droplet-size distribution is obtained. Even without convergence, bump-shaped distributions are harder to obtain in VOF simulations than in level-set or diffuse-interface simulations where they are common (Herrmann 2010; Khanwale et al. 2022; Saurabh et al. 2023), because these two other methods tend to ‘evaporate’ small droplets, as shown on figure 1. A scaling relation for the position of the bump is suggested. The bump only appears at the highest level. Whether the proposed scaling (3.6) for the bump would persist if yet more refined simulations are performed is an interesting open question.

Beyond the physics of hole formation and expansion, the simulations reveal a rich spectrum of elementary phenomena that are the building blocks of atomisation. These include networks of well-resolved ligaments in and around holes, ligaments detached from the main liquid core and either aligned or at an angle with it, bubbles trapped in liquid regions, droplets impacting on sheets, etc. This rich physics is seen in increasing detail as grid resolution is increased. For example, the sharpness and apparent realism of the various expanding holes, with attached ligaments or networks of ligaments inside, was not seen in previous simulations, for example, in Ling et al. (2017). Moreover, a large number of each type of object or ‘building block’ are visible at each instant of time, leading to more significant and converged statistics.

A tantalising question is that of the physical, as opposed to numerical, origin of the perforations. We attempt some discussion of this issue in Appendix A. Two types of rate processes for heterogeneous nucleation are postulated, one based on in-sheet germs or nuclei, such as small bubbles, the other on external perturbations such as floating droplets or particles. Our simulations do not allow us to decide at present between these two processes for hole formation. However, they may in the future helpus to understand how to discriminate between them.

Unlike most of the other round-jet simulations in the literature, ours is bimodal and even trimodal after the application of MD. It seems however that a single mode emerges in the converged range of the best resolved MD simulations, indicating that the appearance of a maximum of droplet size well in the converged region is possible only if the physics is correctly captured down to scales much smaller than the integral scales such as the jet diameter $D$ .

Another important characteristic of simulations and experiments are the state of the fluids at the inlet. The primary mixing-layer instability in our simulations is initiated not by upstream turbulence, as in most round-jet simulations in the literature, but by a pulsating flow. In future investigations it would be interesting to compare the visual aspects of the fragmenting interface as well as number distributions and PDFs, between simulations with or without noisy/turbulent inflow conditions.

The other perspectives opened by this work are numerous. Clearly, simulations of comparable quality at a higher Weber number may allow investigation of the ‘fibre-type’ regime observed in many experiments. Currently, our simulations sit on the boundary between membrane type and fibre type. Simulations of comparable resolution should also be performed in the planar mixing-layer case investigated experimentally by the Grenoble group (Ben Rayana et al. 2006; Fuster et al. 2013) and numerically by Jiang & Ling (2021). The latter authors performed very high resolution simulations with hole formation and interestingly obtained a bimodal distribution, but they did not apply a controlled perforation scheme as in this paper. Simulations of secondary atomisation, that is, the atomisation of a droplet suddenly plunged in a high-speed flow, has already been performed together with controlled perforation by Tang et al. (2023) achieving convergence of the PDF. The secondary atomisation case is probably the most promising for an eventual comparison of the distribution of droplet sizes obtained by simulation with the rich experimental data reviewed and modelled, for example, by Jackiw & Ashgriz (2022).

More distant perspectives or more difficult problems remaining to solve are in our opinion twofold. One is clearly the improvement of the MD procedure, to have a better control of the number of holes punched in each thin sheet region, and a smoother manner to ‘cut out’ the initial hole. The breakup frequency and number of attempted breakups were chosen in a rather arbitrary way. These values are too high as shown above. An adequate breakup frequency would correspond to smaller values of $f_b$ and $N_b$ . A caveat is however that if the breakup frequency ( $f_b$ , $N_b$ ) is too low, some bags may escape MD and breakup numerically. An expensive trial-and-error procedure could be attempted to tune the breakup frequency. Another avenue of improvement would be to modify the MD method to perform an analysis of the bag region to identify the connected thin sheet regions. Then a controlled number of perforations could be performed in any of these regions. This controlled number of perforations per thin region could be inspired by experimental observations or by a theory such as the one in Appendix A. A final, intriguing possibility is that if a level-set or diffuse-interface method is used in conjunction with MD, and the frequency is on the low side, numerical breakup will be avoided in the regions that would not be sufficiently perforated because of the low frequency. There would still be mass loss because of the evaporation of thin sheets, but no curvature oscillations nor an excessive number of small droplets scaling with grid size.

Another challenging future task is to use this type of simulation to infer, by careful matching of simulation results and experiments, the physical mechanisms for weak spots and sheet perforations. Given the wide variety of such mechanisms that have been advanced, this is a formidable enterprise for which we hope this paper could provide partial guidance.