# Growth and instability of the liquid rim in the crown splash regime

## Abstract

We study the formation, growth, and disintegration of jets following impact of a drop on a thin film of the same liquid for and using a combination of numerical simulations and linear stability theory (Agbaglah et al., 2013). Our simulations faithfully capture this phenomena and are in good agreement with experimental profiles obtained from high-speed X-ray imaging. We obtain scaling relations from our simulations and use these as inputs to our stability analysis. The resulting prediction for the most unstable wavelength are in excellent agreement with experimental data. Our calculations show that the dominant destabilizing mechanism is a competition between capillarity and inertia but that deceleration of the rim provides an additional boost to growth. We also predict over the entire parameter range of our study the number and timescale for formation of secondary droplets formed during a splash, based on the assumption that the most unstable mode sets the droplet number.

eurm10
\checkfontmsam10
Growth and instability of the liquid rim]Growth and instability of the liquid rim in the crown splash regime
G. Agbaglah and R. D. Deegan]G.\lsAgbaglah^{1}^{2}

reakup/coalescence, drops

## 1 Introduction

The impact of a drop on a film of the same fluid is ubiquitous in nature and arises in many different contexts such as rain and surf interactions with the air-sea interface, fuel injection systems, spray painting, and atomization. Depending on the governing dimensionless parameters – mainly the ratio of the film thickness to the drop diameter, and the Reynolds and Weber numbers – the impact may cause a splash (see e.g. Cossali et al., 1997; Rioboo et al., 2003; Deegan et al., 2008) herein defined as the generation of secondary droplets. Splashing encompasses a broad variety of qualitatively different morphologies distinguished by their regularity and size-distribution of secondary droplets. Whether these differences arise from a single mechanism or multiple mechanisms remains an open question. Indeed, the cause of splashing has been revisited often in century since the pioneering studies of Worthington (1879) and has been answered in many, seemingly contradictory ways: Rieber & Frohn (1999), Bremond & Villermaux (2006), and Zhang et al. (2010) argue for the Rayleigh-Plateau mechanism, Gueyffier & Zaleski (1998) for a Richtmyer-Meshkov mechanism, Krechetnikov & Homsy (2009) for a combination of the Richtmyer-Meshkov and Rayleigh-Taylor mechanism, Krechetnikov (2010) for a combination of the Rayleigh-Plateau and Rayleigh-Taylor mechanism, and Yarin & Weiss (1995) for a nonlinear amplification mechanism governed by the eikonal equation. More broadly, Roisman and collaborators (Roisman et al., 2006, 2007; Roisman, 2010) investigated the linear stability of a receding liquid sheet and concluded that the liquid rim is subject to a Rayleigh-Plateau and Rayleigh-Taylor instability; a followup study by Agbaglah et al. (2013) that included the growth of the liquid rim in the linear stability analysis concluded that the rim is susceptible to both the Rayleigh-Taylor and Rayleigh-Plateau instabilities but at different times during its evolution: at early times when the rim is decelerating sharply the former dominates and at later times the latter dominates.

In the crown splash regime – the focus of this paper – the splashing proceeds through a sequence of clearly distinguishable steps: a sheet-like jet emerges and grows outward from the neck of fluid connecting the drop and the pool; the jet’s leading edge is pulled back into the sheet by surface tension and forms a rim as it moves back into the sheet; lastly, the rim develops transverse corrugations that grow, sharpen into fingers, and ultimately pinch-off into secondary droplets. Because of its high degree of regularity the crown splash offers the simplest scenario in which to examine extant issues on the origin of secondary droplets.

Below we present the results of our numerical study of the initial axisymmetric phase of a crown splash. The validity of these simulations is demonstrated by their faithful reproduction of the finest details resolvable in experiments. Using the characteristics of the jet obtained from simulation as inputs for the linear stability theory of Agbaglah et al. (2013) we are able to reproduce the experimental results of Zhang et al. (2010) on the breakup of the crown. We extend our predictions to regimes for which there are currently no experiments by extracting from simulations scaling relationships for the rim’s characteristics as a function of Weber and Reynolds number and using these to predict the most unstable wavelength throughout the crown splash regime. We show that the Rayleigh-Plateau mechanism is the dominant mechanism for the parameter range of our study though its growth is significantly enhanced by the Raleigh-Taylor mechanism.

The manuscript is divided into a numerical section and a modeling section. In the numerical section we describe our numerical technique; validate the results of the simulations with high speed X-ray images; we compute the thickness of the jet and the radius and position of the rim for a variety of Reynolds and Weber numbers and develop a scaling relationship for the jet thickness at the moment it emerges from the impact. In the modeling section we use these scaling relationship to predict the number of secondary droplets as a function of Weber and Reynolds number.

## 2 Numerical computations

### 2.1 Numerical method

We simulate two incompressible fluids, a liquid and a gas, with constant densities and and constant viscosities and , with the two-fluid Navier-Stokes equations:

In this formulation the density and viscosity are constant within each phase and discontinuous at the interface. The Dirac distribution function expresses the fact that the surface tension term is concentrated at the interface, and n being the curvature and the normal of the interface respectively and the interfacial tension.

Numerical simulations were performed using GERRIS (Popinet, 2003b), an open-source code where the interface is tracked using a Volume-of-Fluid method on an octree structured grid allowing adaptive mesh refinement, and the incompressibility condition is satisfied using a multigrid solver (see Popinet, 2003a, 2009). This numerical code has been validated with numerous examples and used successfully for many different multiphase problems (Fuster et al., 2009; Agbaglah et al., 2011) including splashing (see e.g. Thoraval et al., 2012).

### 2.2 Numerical simulation

We simulated a spherical liquid drop of diameter and velocity impacting normally on a uniform layer of thickness of the same liquid in an initially quiescent gas. The results depend on five parameters: the viscosity and density ratios ( and ), the dimensionless fluid depth , the liquid Reynolds number and Weber number . We used the axisymmetric formulation of the Navier-Stokes equations described above. Typically we stopped our simulation at because experiments (e.g. Zhang et al. (2010)) show that axial symmetry breaks before this point and hence our axisymmetric simulation cannot follow the true dynamics. Our computational domain extends in the radial direction and in the vertical direction. We define as the moment when the drop first touches the pool; this time is slightly later ( s) than the usual experimental choice of as the moment when the bottom of the drop crosses the undisturbed level of the pool.

We assessed the convergence of the code at , , , for various dimensionless pool depths in the range 0.2 to 5 by observing the velocity norm as a function of grid size. We used the , and norms (respectively the average of the absolute values, the root-mean-square norm and the maximum absolute value) of the difference between the velocity at a given grid size () and the velocity at the finest grid size (), at 66 spatial points dispersed over the computational domain and at 20 different times. The average in time of the computed error are shown in Fig. 1. The convergence is at the first order in the spatial resolution.

In the following, a mesh is used in all the simulations. This corresponds to an initial mesh size of , and it is also the mesh size at the maximum refinement. The mesh is adapted based on three criteria: 1) distance to the interface, 2) curvature of the interface, 3) vorticity magnitude.

### 2.3 Validation

Figure 2 compares simulation results with the equivalent experiments obtained by using high speed X-ray imaging as described in Zhang et al. (2011). The simulations at correspond to the two-jets regime near the transition point from the one-jet regime observed experimentally (Zhang et al., 2011). Due to the slenderness of the structures (e.g. the ejecta or the gap between the ejecta and the lamella), at these parameter values insufficient grid resolution is readily apparent and thus our comparison of with experiments is a stringent test of our simulation’s accuracy. We find excellent qualitative agreement between simulations and experiments.

As shown in the merged simulation-experiment panel of Fig. 2, the rim formed at the edge of the crown is accurately captured in the simulation, but there is small shift at the base of the sheet. We attribute the latter to various perturbations (e.g. surface waves generated when the drop detaches from the needle or air resistance as the drop falls) that cause the drop in the experimental system to deviate from the spherical geometry used in simulations (see for instance Thoraval et al. (2013)).

### 2.4 Numerical results

The aim of our numerics was to find a phenomenological law for the evolution of the characteristics of the jet as a prelude to developing a reduced model of the rim’s dynamics. In particular, we measured the size of the jet when it first emerges from neck , the time-dependence of the rim radius , and the vertical and horizontal distance of the rim from the impact center

The time evolution of the rim radius is shown in Fig. 3(a) for six different parameter sets. We varied Re and We (by changing only the liquid properties: and ), keeping and the impact speed constant. These data show that grows linearly in time and the overall scale decreases with increasing Re and We. As shown in 3(b), the radius depends linearly in time; moreover, the slope scales as We:

(1) |

where is a real constant and is the initial diameter of the rim as it emerges from the neck(see Fig. 4). Note that according to Eq. 1 the time at which the jet emerges and are identical; the difference in these times in our numerical calculations is no more than 5 s and thus we assume it is negligible.

We measured the size of the jet when it first emerges from the neck connecting the drop and the pool from the distance between the inflection points on the neck as shown in Fig. 4(b). The data for thirty different combinations of We (200, 400, 600, 800, 1000) and Re (500, 1000, 1500, 2000, 3000, 4000) are shown in Fig. 4(b-d). We find empirically that

(2) |

Furthermore, the distance of the neck from the center of the drop when the jet first emerges, usually called contact length, decreases with increasing Re (see Fig. 4(b); the center of the drop is at x/D=0) in agreement with the results of Coppola et al. (2011). Combining the scaling for and yields:

(3) |

where and are constants with values obtained from fits to the data of 13 and 0.35 respectively.

We measured the radius of the crown from the horizontal distance of the leading edge of the rim to the impact center as a function of time and the vertical distance of the crown above the undisturbed pool . These data are plotted in Fig. 5 for various combinations of Re and We; and , largely independent of Re and We. (Liang et al. (2013) observed the same behavior.) Hence the acceleration is primarily along the radial direction.

## 3 Modeling wavelength selection of the rim instability

Zhang et al. (2010) showed good agreement of experiments with a stability analysis based on the Rayleigh-Plateau instability that ignored the connection of the rim to the jet and the deceleration of the rim. Fullana & Zaleski (1999) however argued that the former effect is essential because it saturates the growth of the Rayleigh-Plateau instability and thus prevents the formation of secondary droplets within the relevant timescale. Moreover, Krechetnikov & Homsy (2009) also argued that acceleration is important because their measurements, albeit in a different regime from the crown splash, were consistent with a Ritchmyer-Meshkov instability. Below we examine the effects of the connection of the jet to the rim and the deceleration of the rim with a model that merges the long-wavelength linear stability analysis of Agbaglah et al. (2013) and the results of our numerical calculations.

Agbaglah et al. (2013) analyzed in the inviscid limit the stability of the rim of a flat sheet-like jet of constant thickness and speed. Importantly, it incorporate both acceleration and capillarity. Their model agrees well with full numerical simulations of that scenario. The crown sheet however is curved, its thickness and speed are nonuniform, and its geometry is cylindrical. We gloss over these difficulties in an attempt to construct a minimal model of the splash. The model of Agbaglah et al. (2013) requires as inputs the ratio of the jet thickness to rim radius and the acceleration of the rim. We supply the former from the scaling relationship for under the assumption that this value is representative of and the latter from the rim position. Based on this composite model we predict the most unstable wavelength of a crown splash for a wide range of We and Re numbers. In particular, for and for which there is high quality experimental data available we find excellent agreement.

### 3.1 Wavelength selection and secondary droplet production

Using the rim radius and acceleration obtained in the previous sections as inputs for the linear stability calculation of (Agbaglah et al., 2013), we computed the growth rate as a function of wavenumber , and from this dispersion relationship the power and wavelength of the most unstable mode. The dispersion relationship for We=760 and Re=1060 is plotted in Fig. 6(a). Figure 6(b) shows the equivalent dispersion curves for the classical Rayleigh-Plateau and Rayleigh-Taylor instabilities. The most unstable wavelength of the classical Rayleigh-Plateau most closely tracks the most unstable wavelength of our stability calculation.

From the growth rate we compute the amplitude of the th mode:

(4) |

where is taken as a constant independent of wavelength under the assumption that the initial microscopic corrugations present in the system are distributed like broadband white noise. The peak mode is selected as the mode with the maximum amplitude at any given time. The peak wavelength and its amplitude are compared with the experimental measurements of Zhang et al. (2010) in Fig. 7. In addition, we plot the classical inviscid Rayleigh-Plateau and Rayleigh-Taylor theory (Chandrasekhar, 1981) in Fig. 7 for comparison. In all three calculations is taken as an adjustable parameter chosen to maximize the correspondence between the calculation and the experimental measurements. These data show that our model is in excellent agreement with measurements.

We can predict the number of secondary droplets and the timescale for their onset over the entire parameter range of the crown splash with a few further assumptions. We repeated the above calculation of on a regular grid of parameter values (Re, We) assuming that is the same for all parameters and equal to the value obtained from the comparison in figure 7 (right). We further assume that when the amplitude of the peak corrugation grows to a size comparable to the diameter of the rim (i.e. , the choice of is motivated by the fact that the antisymmetric part of the perturbation is dominant in the nonlinear regime (e.g. see Fig. 9 and 10 of Agbaglah et al. (2013))) nonlinear effects become dominant and lock in the most unstable mode. Lastly, we assume that the number of droplets is given by the peak mode as suggested by the experiments of Zhang et al. (2010). Figure 8 shows the resulting predictions. These data predict that the number of secondary droplets increase with both Re and We. No secondary droplets are predicted for the hashed region of Fig. 8, since for the entire timescale of splashing (i.e. ).

## 4 Discussion and Conclusion

We performed axisymmetric simulations of drop impact on a thin film in the crown splash regime using a volume-of-fluid implementation of the Navier-Stokes equations and validated these by comparing with experimental profiles obtained using high speed X-ray imaging. By combining the result of our axisymmetric simulations with the theory of Agbaglah et al. (2013) for symmetry breaking of a flat sheet, we predict the most unstable mode of the axial symmetry breaking instability in a crown splash. With the further assumption, supported by experiments, that the crests of the most unstable mode are the origin of secondary droplets, we predict the secondary droplet production throughout the crown splash regime (, ).

From simulations we obtain scaling relations for the position ( and ) and the size () of the rim, and the initial thickness () of the jet. These are then used in the linear stability theory of Agbaglah et al. (2013) to obtain the most unstable wavelength and hence the number of secondary droplets. No adjustable parameters are used to predict the most unstable wavelength; a single adjustable parameter, the initial amplitude of the corrugation that initiates the instability, is used to predict the magnitude of the instability. Our predictions are in excellent agreement with the experimental measurements of Zhang et al. (2010).

The agreement between experiments and theory could be perhaps improved if we allowed for a dependence of on Re and We. However, in absence of experimental evidence that would warrant this more complicated scenario, we chose to be independent of wavelength, Re and We.

The debate about the cause of splashing has largely revolved around the relative importance of capillarity versus acceleration, or – as it is frequently phrased in the literature – whether the instability is primarily Rayleigh-Plateau-like or Rayleigh-Taylor-like. With our simulation results we can quantitatively examine this question in the crown splash regime. As shown by Fig. 6, our results favor a Rayleigh-Plateau-like instability as the primary mechanism for wavelength selection.

Our calculation however, differs from the classical Rayleigh-Plateau in that it includes the effects of deceleration of the rim and the connection of the rim to the jet. The combination of these two effects yields a better fit to the experimental data. The rim-jet connection weakens the destabilizing role of capillarity. For example, consider the extreme case that prevails when the jet first forms: the jet ends in a semicircle with a diameter equal to the sheet thickness. Unlike the case of a full cylinder, there is no surface energy gain for transverse corrugation for the free end of the jet. Thus, we expect that the classical Rayleigh-Plateau theory which is for a full cylinder will produce a greater growth rate than our theory; this is indeed the case as shown by Fig. 6. The deceleration of the rim however, counteracts the dampening effect of the rim-jet connection. As shown in Fig. 9(a), the amplitude of the most unstable mode is significantly more amplified when acceleration is included. Roisman (2010) reached a similar conclusion in the context of flat sheets.

How can we reconcile the obvious effect of acceleration on the magnitude of the mode with its absence in wavelength selection? The deceleration phase is too short-lived to effect wavelength selection. As the rim approaches the Taylor-Culick limiting speed (relative to the jet speed), acceleration and its destabilizing effect vanish. We find in our simulations that this happens rapidly. The legacy of this initial acceleration phase is a greater initial amplitude for the subsequent evolution of the rim due to the capillary instability. Moreover, since the Rayleigh-Plateau and Rayleigh-Taylor instabilities have similar peak wavelengths because they both originate from the same driving force, the principal contribution of the acceleration instability is centered around the modes susceptible to the capillary instability.

## 5 Acknowledgements

The authors thank the James S. McDonnell Foundation for support through a 21st Century Science Initiative in Studying Complex Systems Research Award; L. V. Zhang, P. Ray and C. Josserand for valuable discussions; and S. Veerapaneni, S. Zaleski and the Institut Jean le Rond d’Alembert for sharing their computational resources.

### Footnotes

- thanks: Email address for correspondence: agbagla@umich.edu
- thanks: Email address for correspondence: rddeegan@umich.edu

### References

- Agbaglah, G., Delaux, S., Fuster, D., Hoepffner, J., Josserand, C., Popinet, S. Ray, P., Scardovelli, R. & Zaleski, S. 2011 Parallel simulation of multiphase flows using octree adaptivity and the volume-of-fluid method. C. R. Mecanique. 339, 194–207.
- Agbaglah, G., Josserand, C. & Zaleski, S. 2013 Longitudinal instability of a liquid rim. Phys. Fluids 25, 022103.
- Bremond, N. & Villermaux, E. 2006 Atomization by jet impact. Journal of Fluid Mechanics 549, 273–306.
- Chandrasekhar, S. 1981 Hydrodynamic and hydromagnetic stability. Dover, New York .
- Coppola, G., Rocco, G. & de Luca, L. 2011 Insights on impact of a plane drop on a thin liqui film. Phys. Fluids. 23, 022105.
- Cossali, G. E., Coghe, A. & Marengo, M. 1997 Impact of a single drop on a wetted solid surface. Experiments in Fluids 22 (6), 463–472.
- Deegan, R. D., Brunet, P. & Eggers, J. 2008 Complexities of splashing. Nonlinearity 21 (1), C1–C11, 0951-7715.
- Fullana, J. M. & Zaleski, S. 1999 Stability of a growing end rim in a liquid sheet of uniform thickness. Physics of Fluids 11 (5), 952–954.
- Fuster, Daniel, Agbaglah, Gilou, Josserand, Christophe, Popinet, Stephane & Zaleski, Stephane 2009 Numerical simulation of droplets, bubbles and waves: state of the art. Fluid Dynamics Research 41 (6), 065001.
- Gueyffier, D. & Zaleski, S. 1998 Finger formation during droplet impact on a liquid film. C.R. Acad. Sci., Ser. IIb: Mec., Phys., Chim., Astron 326 (12), 839–844.
- Krechetnikov, R. 2010 Stability of liquid sheet edges. Physics of Fluids 22, 092101.
- Krechetnikov, Rouslan & Homsy, George M. 2009 Crown-forming instability phenomena in the drop splash problem. J Colloid Interface Sci 331 (2), 555–9.
- Liang, G., Guo, Y., Shen, S. & Yang, Y. 2013 Crown behavior and bubble entrainment during a drop impact on a liquid film. Theor. Comput. Fluid Dyn. DOI: 10.1007/s00162-013-0308-z.
- Popinet, S. 2003a Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries. J. Comp. Phys. 190, 572–600.
- Popinet, S. 2003b Gerris flow solver. http://gfs.sourceforge.net/.
- Popinet, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. J. Comput. Phys. 228, 5838–5866.
- Rieber, M. & Frohn, A. 1999 A numerical study on the mechanism of splashing. Int. J. Heat Fluid Flow 20 (5), 455–461.
- Rioboo, R., Bauthier, C., Conti, J., Voue, M. & De Coninck, J. 2003 Experimental investigation of splash and crown formation during single drop impact on wetted surfaces. Experiments in Fluids 35 (6), 648–652.
- Roisman, I. V. 2010 On the instability of a free viscous rim. Journal of Fluid Mechanics 661, 206–228.
- Roisman, I. V., Gambaryan-Roisman, T., Kyriopoulos, O., Stephan, P. & Tropea, C. 2007 Breakup and atomization of a stretching crown. Physical Review E 76 (2).
- Roisman, I. V., Horvat, K. & Tropea, C. 2006 Spray impact: Rim transverse instability initiating fingering and splash, and description of a secondary spray. Physics of Fluids 18 (10).
- Thoraval, Marie-Jean, Takehara, Kohsei, Etoh, Takeharu Goji, Popinet, Stephane, Ray, Pascal, Josserand, Christophe, Zaleski, Stephane & Thoroddsen, Sigurdur T. 2012 von karman vortex street within an impacting drop. Physical Review Letters 108 (26).
- Thoraval, M. J., Takehara, K., Etoh, T. G. & Thoroddsen, S. T. 2013 Drop impact entrapment of bubble rings. Journal of Fluid Mechanics 724, 234–258.
- Worthington, A. M. 1879 On the spontaneous segmentation of a liquid annulus. Proc. Phys. Soc. London 30, 49–60.
- Yarin, A. L. & Weiss, D. A. 1995 Impact of drops on solid-surfaces - self-similar capillary waves, and splashing as a new-type of kinematic discontinuity. Journal of Fluid Mechanics 283, 141–173.
- Zhang, L. V., Brunet, P., Eggers, J. & Deegan, R. D. 2010 Wavelength selection in the crown splash. Physics of Fluids 22 (12).
- Zhang, L. V., Toole, J., Fezzaa, K. & Deegan, R. D. 2011 Evolution of the ejecta sheet from the impact of a drop with a deep pool. J. Fluid Mech. 396, 1–11.